@@ -56,7 +56,7 @@ weights or the similar pseudo-BMA weights.
5656In addition to the __ loo__ package we will also load the __ rstanarm__ package
5757for fitting the models.
5858
59- ``` {r, setup, message=FALSE}
59+ ``` {r setup, message=FALSE}
6060library(rstanarm)
6161library(loo)
6262```
@@ -68,7 +68,7 @@ example as follows:
6868
6969> A popular hypothesis has it that primates with larger brains produce more energetic milk, so that brains can grow quickly. ... The question here is to what extent energy content of milk, measured here by kilocalories, is related to the percent of the brain mass that is neocortex. ... We'll end up needing female body mass as well, to see the masking that hides the relationships among the variables.
7070
71- ``` {r}
71+ ``` {r data }
7272data(milk)
7373d <- milk[complete.cases(milk),]
7474d$neocortex <- d$neocortex.perc /100
7878We repeat the analysis in Chapter 6 of _ Statistical Rethinking_ using the
7979following four models (here we use the default weakly informative priors in __ rstanarm__ , while flat priors were used in _ Statistical Rethinking_ ).
8080
81- ``` {r, results="hide"}
81+ ``` {r fits , results="hide"}
8282fit1 <- stan_glm(kcal.per.g ~ 1, data = d, seed = 2030)
8383fit2 <- update(fit1, formula = kcal.per.g ~ neocortex)
8484fit3 <- update(fit1, formula = kcal.per.g ~ log(mass))
@@ -94,7 +94,7 @@ by the __rstanarm__ package (a wrapper around `waic` from the __loo__ package),
9494which allows us to just pass in our fitted model objects instead of first
9595extracting the log-likelihood values.
9696
97- ``` {r}
97+ ``` {r waic }
9898waic1 <- waic(fit1)
9999waic2 <- waic(fit2)
100100waic3 <- waic(fit3)
@@ -117,7 +117,7 @@ needed values to pass to the __loo__ package. (Like __rstanarm__, some other R
117117packages for fitting Stan models, e.g. __ brms__ , also provide similar methods
118118for interfacing with the __ loo__ package.)
119119
120- ``` {r}
120+ ``` {r loo }
121121# note: the loo function accepts a 'cores' argument that we recommend specifying
122122# when working with bigger datasets
123123
@@ -136,7 +136,7 @@ lpd_point <- cbind(
136136With ` loo ` we don't get any warnings for models 3 and 4, but for illustration of
137137good results, we display the diagnostic details for these models anyway.
138138
139- ``` {r}
139+ ``` {r print-loo }
140140print(loo3)
141141print(loo4)
142142```
@@ -149,7 +149,7 @@ Next we compute and compare 1) WAIC weights, 2) Pseudo-BMA weights without
149149Bayesian bootstrap, 3) Pseudo-BMA+ weights with Bayesian bootstrap, and 4)
150150Bayesian stacking weights.
151151
152- ``` {r}
152+ ``` {r weights }
153153waic_wts <- exp(waics) / sum(exp(waics))
154154pbma_wts <- pseudobma_weights(lpd_point, BB=FALSE)
155155pbma_BB_wts <- pseudobma_weights(lpd_point) # default is BB=TRUE
@@ -173,7 +173,7 @@ bunch of times, but you can imagine that instead we would have ten alternative
173173models with about the same predictive performance. WAIC weights for such a
174174scenario would be close to the following:
175175
176- ``` {r}
176+ ``` {r waic_wts_demo }
177177waic_wts_demo <-
178178 exp(waics[c(1,1,1,1,1,1,1,1,1,1,2,3,4)]) /
179179 sum(exp(waics[c(1,1,1,1,1,1,1,1,1,1,2,3,4)]))
@@ -192,7 +192,7 @@ very similar models (in this toy example repeated models) to share their weight
192192while more unique models keep their original weights. In our example
193193we can see this difference clearly:
194194
195- ``` {r}
195+ ``` {r stacking_weights }
196196stacking_weights(lpd_point[,c(1,1,1,1,1,1,1,1,1,1,2,3,4)])
197197```
198198Using stacking, the weight for the best model stays essentially unchanged.
@@ -208,7 +208,7 @@ McElreath describes as follows:
208208We build models predicting the total number of tools given the log
209209population size and the contact rate (high vs. low).
210210
211- ``` {r}
211+ ``` {r Kline }
212212data(Kline)
213213d <- Kline
214214d$log_pop <- log(d$population)
@@ -219,7 +219,7 @@ str(d)
219219We start with a Poisson regression model with the log population size,
220220the contact rate, and an interaction term between them (priors are
221221informative priors as in _ Statistical Rethinking_ ).
222- ``` {r, results="hide"}
222+ ``` {r fit10 , results="hide"}
223223fit10 <-
224224 stan_glm(
225225 total_tools ~ log_pop + contact_high + log_pop * contact_high,
@@ -234,7 +234,7 @@ fit10 <-
234234Before running other models, we check whether Poisson is good choice
235235as the conditional observation model.
236236
237- ``` {r}
237+ ``` {r loo10 }
238238loo10 <- loo(fit10)
239239print(loo10)
240240```
@@ -249,7 +249,7 @@ We can compute LOO more accurately by running Stan again for the
249249leave-one-out folds with high $k$ estimates. When using __ rstanarm__
250250this can be done by specifying the ` k_threshold ` argument:
251251
252- ``` {r}
252+ ``` {r loo10-threshold }
253253loo10 <- loo(fit10, k_threshold=0.7)
254254print(loo10)
255255```
@@ -258,7 +258,7 @@ In this case we see that there is not much difference, and thus it is relatively
258258safe to continue.
259259
260260As a comparison we also compute WAIC:
261- ``` {r}
261+ ``` {r waic10 }
262262waic10 <- waic(fit10)
263263print(waic10)
264264```
@@ -269,7 +269,7 @@ optimistic. We recommend using the PSIS-LOO results instead.
269269To assess whether the contact rate and interaction term are useful, we can make
270270a comparison to models without these terms.
271271
272- ``` {r, results="hide"}
272+ ``` {r contact_high , results="hide"}
273273fit11 <- update(fit10, formula = total_tools ~ log_pop + contact_high)
274274fit12 <- update(fit10, formula = total_tools ~ log_pop)
275275```
@@ -288,7 +288,7 @@ lpd_point <- cbind(
288288```
289289
290290For comparison we'll also compute WAIC values for these additional models:
291- ``` {r}
291+ ``` {r waic-contact_high }
292292waic11 <- waic(fit11)
293293waic12 <- waic(fit12)
294294waics <- c(
@@ -304,7 +304,7 @@ Finally, we compute 1) WAIC weights, 2) Pseudo-BMA weights without
304304Bayesian bootstrap, 3) Pseudo-BMA+ weights with Bayesian bootstrap, and
3053054 ) Bayesian stacking weights.
306306
307- ``` {r}
307+ ``` {r weights-contact_high }
308308waic_wts <- exp(waics) / sum(exp(waics))
309309pbma_wts <- pseudobma_weights(lpd_point, BB=FALSE)
310310pbma_BB_wts <- pseudobma_weights(lpd_point) # default is BB=TRUE
@@ -336,7 +336,7 @@ model objects from __rstanarm__) as well as fitted model objects from
336336other packages (e.g. __ brms__ ) that do the preparation work for the user
337337(see, e.g., the examples at ` help("loo_model_weights", package = "rstanarm") ` ).
338338
339- ``` {r}
339+ ``` {r loo_model_weights }
340340# using list of loo objects
341341loo_list <- list(loo10, loo11, loo12)
342342loo_model_weights(loo_list)
0 commit comments