1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
params <-
list(EVAL = TRUE)
## ----SETTINGS-knitr, include=FALSE--------------------------------------------
stopifnot(require(knitr))
opts_chunk$set(
comment=NA,
eval = if (isTRUE(exists("params"))) params$EVAL else FALSE,
dev = "png",
dpi = 150,
fig.asp = 0.618,
fig.width = 5,
out.width = "60%",
fig.align = "center"
)
## ----setup, message=FALSE-----------------------------------------------------
library(rstanarm)
library(loo)
## ----data---------------------------------------------------------------------
data(milk)
d <- milk[complete.cases(milk),]
d$neocortex <- d$neocortex.perc /100
str(d)
## ----fits, results="hide"-----------------------------------------------------
fit1 <- stan_glm(kcal.per.g ~ 1, data = d, seed = 2030)
fit2 <- update(fit1, formula = kcal.per.g ~ neocortex)
fit3 <- update(fit1, formula = kcal.per.g ~ log(mass))
fit4 <- update(fit1, formula = kcal.per.g ~ neocortex + log(mass))
## ----waic---------------------------------------------------------------------
waic1 <- waic(fit1)
waic2 <- waic(fit2)
waic3 <- waic(fit3)
waic4 <- waic(fit4)
waics <- c(
waic1$estimates["elpd_waic", 1],
waic2$estimates["elpd_waic", 1],
waic3$estimates["elpd_waic", 1],
waic4$estimates["elpd_waic", 1]
)
## ----loo----------------------------------------------------------------------
# note: the loo function accepts a 'cores' argument that we recommend specifying
# when working with bigger datasets
loo1 <- loo(fit1)
loo2 <- loo(fit2)
loo3 <- loo(fit3)
loo4 <- loo(fit4)
lpd_point <- cbind(
loo1$pointwise[,"elpd_loo"],
loo2$pointwise[,"elpd_loo"],
loo3$pointwise[,"elpd_loo"],
loo4$pointwise[,"elpd_loo"]
)
## ----print-loo----------------------------------------------------------------
print(loo3)
print(loo4)
## ----weights------------------------------------------------------------------
waic_wts <- exp(waics) / sum(exp(waics))
pbma_wts <- pseudobma_weights(lpd_point, BB=FALSE)
pbma_BB_wts <- pseudobma_weights(lpd_point) # default is BB=TRUE
stacking_wts <- stacking_weights(lpd_point)
round(cbind(waic_wts, pbma_wts, pbma_BB_wts, stacking_wts), 2)
## ----waic_wts_demo------------------------------------------------------------
waic_wts_demo <-
exp(waics[c(1,1,1,1,1,1,1,1,1,1,2,3,4)]) /
sum(exp(waics[c(1,1,1,1,1,1,1,1,1,1,2,3,4)]))
round(waic_wts_demo, 3)
## ----stacking_weights---------------------------------------------------------
stacking_weights(lpd_point[,c(1,1,1,1,1,1,1,1,1,1,2,3,4)])
## ----Kline--------------------------------------------------------------------
data(Kline)
d <- Kline
d$log_pop <- log(d$population)
d$contact_high <- ifelse(d$contact=="high", 1, 0)
str(d)
## ----fit10, results="hide"----------------------------------------------------
fit10 <-
stan_glm(
total_tools ~ log_pop + contact_high + log_pop * contact_high,
family = poisson(link = "log"),
data = d,
prior = normal(0, 1, autoscale = FALSE),
prior_intercept = normal(0, 100, autoscale = FALSE),
seed = 2030
)
## ----loo10--------------------------------------------------------------------
loo10 <- loo(fit10)
print(loo10)
## ----loo10-threshold----------------------------------------------------------
loo10 <- loo(fit10, k_threshold=0.7)
print(loo10)
## ----waic10-------------------------------------------------------------------
waic10 <- waic(fit10)
print(waic10)
## ----contact_high, results="hide"---------------------------------------------
fit11 <- update(fit10, formula = total_tools ~ log_pop + contact_high)
fit12 <- update(fit10, formula = total_tools ~ log_pop)
## ----loo-contact_high---------------------------------------------------------
(loo11 <- loo(fit11))
(loo12 <- loo(fit12))
## ----relo-contact_high--------------------------------------------------------
loo11 <- loo(fit11, k_threshold=0.7)
loo12 <- loo(fit12, k_threshold=0.7)
lpd_point <- cbind(
loo10$pointwise[, "elpd_loo"],
loo11$pointwise[, "elpd_loo"],
loo12$pointwise[, "elpd_loo"]
)
## ----waic-contact_high--------------------------------------------------------
waic11 <- waic(fit11)
waic12 <- waic(fit12)
waics <- c(
waic10$estimates["elpd_waic", 1],
waic11$estimates["elpd_waic", 1],
waic12$estimates["elpd_waic", 1]
)
## ----weights-contact_high-----------------------------------------------------
waic_wts <- exp(waics) / sum(exp(waics))
pbma_wts <- pseudobma_weights(lpd_point, BB=FALSE)
pbma_BB_wts <- pseudobma_weights(lpd_point) # default is BB=TRUE
stacking_wts <- stacking_weights(lpd_point)
round(cbind(waic_wts, pbma_wts, pbma_BB_wts, stacking_wts), 2)
## ----loo_model_weights--------------------------------------------------------
# using list of loo objects
loo_list <- list(loo10, loo11, loo12)
loo_model_weights(loo_list)
loo_model_weights(loo_list, method = "pseudobma")
loo_model_weights(loo_list, method = "pseudobma", BB = FALSE)
|
