@@ -15,6 +15,9 @@ params:
1515``` {r, child="children/SETTINGS-knitr.txt"}
1616```
1717
18+ ``` {r, child="children/SEE-ONLINE.txt", eval = if (isTRUE(exists("params"))) !params$EVAL else TRUE}
19+ ```
20+
1821# Introduction
1922
2023This vignette demonstrates how to use the __ loo__ package to carry out
@@ -35,7 +38,7 @@ encourage readers to refer to the following papers for more details:
3538In addition to the __ loo__ package, we'll also be using __ rstanarm__ and
3639__ bayesplot__ :
3740
38- ``` {r, setup, message=FALSE}
41+ ``` {r setup, message=FALSE}
3942library("rstanarm")
4043library("bayesplot")
4144library("loo")
@@ -72,7 +75,7 @@ Because the number of days for which the roach traps were used is not the same
7275for all apartments in the sample, we use the ` offset ` argument to specify that
7376` log(exposure2) ` should be added to the linear predictor.
7477
75- ``` {r}
78+ ``` {r data }
7679# the 'roaches' data frame is included with the rstanarm package
7780data(roaches)
7881str(roaches)
@@ -86,7 +89,7 @@ roaches$roach1 <- roaches$roach1 / 100
8689We'll fit a simple Poisson regression model using the ` stan_glm ` function from
8790the __ rstanarm__ package.
8891
89- ``` {r, count-roaches-mcmc, results="hide"}
92+ ``` {r count-roaches-mcmc, results="hide"}
9093fit1 <-
9194 stan_glm(
9295 formula = y ~ roach1 + treatment + senior,
@@ -121,7 +124,7 @@ to the `loo` function (see, e.g. `help("loo.array", package = "loo")`). We'll
121124also use the argument ` save_psis = TRUE ` to save some intermediate results to be
122125re-used later.
123126
124- ``` {r}
127+ ``` {r loo1 }
125128loo1 <- loo(fit1, save_psis = TRUE)
126129```
127130
@@ -130,7 +133,7 @@ some observations the leave-one-out posteriors are different enough from the
130133full posterior that importance-sampling is not able to correct the difference.
131134We can see more details by printing the ` loo ` object.
132135
133- ``` {r}
136+ ``` {r print-loo1 }
134137print(loo1)
135138```
136139
@@ -159,7 +162,7 @@ Using the `plot` method on our `loo1` object produces a plot of the $k$ values
159162with horizontal lines corresponding to the same categories as in the
160163printed output above.
161164
162- ``` {r, out.width = "70%"}
165+ ``` {r plot-loo1 , out.width = "70%"}
163166plot(loo1)
164167```
165168
@@ -181,7 +184,7 @@ the LOO-PIT values for our model (thick curve) is compared to many
181184independently generated samples (each the same size as our dataset) from the
182185standard uniform distribution (thin curves).
183186
184- ``` {r}
187+ ``` {r ppc_loo_pit_overlay }
185188yrep <- posterior_predict(fit1)
186189
187190ppc_loo_pit_overlay(
@@ -202,17 +205,17 @@ regression, which is commonly used for overdispersed count data.
202205Unlike the Poisson distribution, the negative binomial distribution
203206allows the conditional mean and variance of $y$ to differ.
204207
205- ``` {r, count-roaches-negbin, results="hide"}
208+ ``` {r count-roaches-negbin, results="hide"}
206209fit2 <- update(fit1, family = neg_binomial_2)
207210```
208211
209- ``` {r}
212+ ``` {r loo2 }
210213loo2 <- loo(fit2, save_psis = TRUE, cores = 2)
211214print(loo2)
212215```
213216
214217
215- ``` {r}
218+ ``` {r plot-loo2 }
216219plot(loo2, label_points = TRUE)
217220```
218221
@@ -232,7 +235,7 @@ calculations are performed exactly for that observation. The results are then
232235recombined with the approximate LOO calculations already carried out for the
233236observations without problematic $k$ values:
234237
235- ``` {r}
238+ ``` {r reloo }
236239if (any(pareto_k_values(loo2) > 0.7)) {
237240 loo2 <- loo(fit2, save_psis = TRUE, k_threshold = 0.7)
238241}
@@ -248,7 +251,7 @@ still some degree of model misspecification, but this is much better than the
248251` p_loo ` estimate for the Poisson model.
249252
250253For further model checking we again examine the LOO-PIT values.
251- ``` {r}
254+ ``` {r ppc_loo_pit_overlay-negbin }
252255yrep <- posterior_predict(fit2)
253256ppc_loo_pit_overlay(roaches$y, yrep, lw = weights(loo2$psis_object))
254257```
@@ -263,7 +266,7 @@ the data.
263266We can use the ` loo_compare ` function to compare our two models on
264267expected log predictive density (ELPD) for new data:
265268
266- ``` {r, count-roaches-loo }
269+ ``` {r loo_compare }
267270loo_compare(loo1, loo2)
268271```
269272
0 commit comments