Deploying to gh-pages from @ 6cee125 🚀

Author: florence-bockting

dev/articles/loo2-elpd.html

@@ -18,6 +18,7 @@ matching*](https://mc-stan.org/loo/articles/loo2-moment-matching.html).
 Here is the Stan code for fitting a Poisson regression model:
 
 ``` r
+
 # Note: some syntax used in this Stan program requires RStan >= 2.26 (or CmdStanR)
 # To use an older version of RStan change the line declaring `y` to: int y[N];
 stancode <- "
@@ -60,6 +61,7 @@ In addition to **loo**, we load the **rstan** package for fitting the
 model. We will also need the **rstanarm** package for the data.
 
 ``` r
+
 library("rstan")
 library("loo")
 seed <- 9547
@@ -76,6 +78,7 @@ is evaluated on the held-out “test” data.
 The data is divided between train (80% of the data) and test (20%):
 
 ``` r
+
 # Prepare data
 data(roaches, package = "rstanarm")
 roaches$roach1 <- sqrt(roaches$roach1)
@@ -111,15 +114,16 @@ Next we fit the model to the “test” data in Stan using the **rstan**
 package:
 
 ``` r
+
 # Compile
 stanmodel <- stan_model(model_code = stancode)
 ```
 
     Trying to compile a simple C file
 
-    Running /opt/R/4.5.3/lib/R/bin/R CMD SHLIB foo.c
+    Running /opt/R/4.6.0/lib/R/bin/R CMD SHLIB foo.c
     using C compiler: ‘gcc (Ubuntu 13.3.0-6ubuntu2~24.04.1) 13.3.0’
-    gcc -std=gnu2x -I"/opt/R/4.5.3/lib/R/include" -DNDEBUG   -I"/home/runner/work/_temp/Library/Rcpp/include/"  -I"/home/runner/work/_temp/Library/RcppEigen/include/"  -I"/home/runner/work/_temp/Library/RcppEigen/include/unsupported"  -I"/home/runner/work/_temp/Library/BH/include" -I"/home/runner/work/_temp/Library/StanHeaders/include/src/"  -I"/home/runner/work/_temp/Library/StanHeaders/include/"  -I"/home/runner/work/_temp/Library/RcppParallel/include/"  -I"/home/runner/work/_temp/Library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/home/runner/work/_temp/Library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include    -fpic  -g -O2  -c foo.c -o foo.o
+    gcc -std=gnu2x -I"/opt/R/4.6.0/lib/R/include" -DNDEBUG   -I"/home/runner/work/_temp/Library/Rcpp/include/"  -I"/home/runner/work/_temp/Library/RcppEigen/include/"  -I"/home/runner/work/_temp/Library/RcppEigen/include/unsupported"  -I"/home/runner/work/_temp/Library/BH/include" -I"/home/runner/work/_temp/Library/StanHeaders/include/src/"  -I"/home/runner/work/_temp/Library/StanHeaders/include/"  -I"/home/runner/work/_temp/Library/RcppParallel/include/"  -I"/home/runner/work/_temp/Library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/home/runner/work/_temp/Library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include    -fpic  -g -O2  -c foo.c -o foo.o
     In file included from /home/runner/work/_temp/Library/RcppEigen/include/Eigen/Core:19,
                      from /home/runner/work/_temp/Library/RcppEigen/include/Eigen/Dense:1,
                      from /home/runner/work/_temp/Library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
@@ -128,9 +132,10 @@ stanmodel <- stan_model(model_code = stancode)
       679 | #include <cmath>
           |          ^~~~~~~
     compilation terminated.
-    make: *** [/opt/R/4.5.3/lib/R/etc/Makeconf:202: foo.o] Error 1
+    make: *** [/opt/R/4.6.0/lib/R/etc/Makeconf:190: foo.o] Error 1
 
 ``` r
+
 # Fit model
 fit <- sampling(stanmodel, data = data_train, seed = seed, refresh = 0)
 ```
@@ -142,6 +147,7 @@ using independent data, the log predictive density coincides with the
 log likelihood of the test data.
 
 ``` r
+
 gen_test <- gqs(stanmodel, draws = as.matrix(fit), data= data_test)
 log_pd <- extract_log_lik(gen_test)
 ```
@@ -152,9 +158,11 @@ Now we evaluate the predictive performance of the model on the test data
 using [`elpd()`](https://mc-stan.org/loo/dev/reference/elpd.md).
 
 ``` r
+
 (elpd_holdout <- elpd(log_pd))
 ```
 
+
     Computed from 4000 by 52 log-likelihood matrix using the generic elpd function
 
          Estimate    SE
@@ -177,6 +185,7 @@ We use the data that is already pre-processed and we divide it in 10
 random folds using `kfold_split_random`
 
 ``` r
+
 # Prepare data
 roaches$fold <- kfold_split_random(K = 10, N = nrow(roaches))
 ```
@@ -192,6 +201,7 @@ this loop is a matrix of the log pointwise predictive densities of all
 the observations.
 
 ``` r
+
 # Prepare a matrix with the number of post-warmup iterations by number of observations:
 log_pd_kfold <- matrix(nrow = 4000, ncol = nrow(roaches))
 # Loop over the folds
@@ -226,9 +236,11 @@ Now we evaluate the predictive performance of the model on the 10 folds
 using [`elpd()`](https://mc-stan.org/loo/dev/reference/elpd.md).
 
 ``` r
+
 (elpd_kfold <- elpd(log_pd_kfold))
 ```
 
+
     Computed from 4000 by 262 log-likelihood matrix using the generic elpd function
 
          Estimate     SE

dev/articles/loo2-elpd.md

@@ -28,6 +28,7 @@ In addition to the **loo** package, we’ll also be using **rstanarm** and
 **bayesplot**:
 
 ``` r
+
 library("rstanarm")
 library("bayesplot")
 library("loo")
@@ -67,6 +68,7 @@ were used is not the same for all apartments in the sample, we use the
 the linear predictor.
 
 ``` r
+
 # the 'roaches' data frame is included with the rstanarm package
 data(roaches)
 str(roaches)
@@ -80,6 +82,7 @@ str(roaches)
      $ exposure2: num  0.8 0.6 1 1 1.14 ...
 
 ``` r
+
 # rescale to units of hundreds of roaches
 roaches$roach1 <- roaches$roach1 / 100
 ```
@@ -90,6 +93,7 @@ We’ll fit a simple Poisson regression model using the `stan_glm`
 function from the **rstanarm** package.
 
 ``` r
+
 fit1 <-
   stan_glm(
     formula = y ~ roach1 + treatment + senior,
@@ -127,6 +131,7 @@ We’ll also use the argument `save_psis = TRUE` to save some intermediate
 results to be re-used later.
 
 ``` r
+
 loo1 <- loo(fit1, save_psis = TRUE)
 ```
 
@@ -139,9 +144,11 @@ correct the difference. We can see more details by printing the `loo`
 object.
 
 ``` r
+
 print(loo1)
 ```
 
+
     Computed from 4000 by 262 log-likelihood matrix.
 
              Estimate     SE
@@ -186,6 +193,7 @@ to fit the model) with horizontal lines corresponding to the same
 categories as in the printed output above.
 
 ``` r
+
 plot(loo1)
 ```
 
@@ -209,6 +217,7 @@ QQ-plot the LOO-PIT values for our model (y-axi) is compared to standard
 uniform distribution (x-axis).
 
 ``` r
+
 yrep <- posterior_predict(fit1)
 
 ppc_loo_pit_qq(
@@ -238,19 +247,23 @@ Unlike the Poisson distribution, the negative binomial distribution
 allows the conditional mean and variance of \\y\\ to differ.
 
 ``` r
+
 fit2 <- update(fit1, family = neg_binomial_2)
 ```
 
 ``` r
+
 loo2 <- loo(fit2, save_psis = TRUE, cores = 2)
 ```
 
     Warning: Found 1 observation(s) with a pareto_k > 0.7. We recommend calling 'loo' again with argument 'k_threshold = 0.7' in order to calculate the ELPD without the assumption that these observations are negligible. This will refit the model 1 times to compute the ELPDs for the problematic observations directly.
 
 ``` r
+
 print(loo2)
 ```
 
+
     Computed from 4000 by 262 log-likelihood matrix.
 
              Estimate   SE
@@ -269,6 +282,7 @@ print(loo2)
     See help('pareto-k-diagnostic') for details.
 
 ``` r
+
 plot(loo2, label_points = TRUE)
 ```
 
@@ -292,6 +306,7 @@ recombined with the approximate LOO calculations already carried out for
 the observations without problematic \\k\\ values:
 
 ``` r
+
 if (any(pareto_k_values(loo2) > 0.7)) {
   loo2 <- loo(fit2, save_psis = TRUE, k_threshold = 0.7)
 }
@@ -300,12 +315,15 @@ if (any(pareto_k_values(loo2) > 0.7)) {
     1 problematic observation(s) found.
     Model will be refit 1 times.
 
+
     Fitting model 1 out of 1 (leaving out observation 93)
 
 ``` r
+
 print(loo2)
 ```
 
+
     Computed from 4000 by 262 log-likelihood matrix.
 
              Estimate   SE
@@ -330,6 +348,7 @@ is much better than the `p_loo` estimate for the Poisson model.
 For further model checking we again examine the LOO-PIT values.
 
 ``` r
+
 yrep <- posterior_predict(fit2)
 ppc_loo_pit_qq(roaches$y, yrep, lw = weights(loo2$psis_object))
 ```
@@ -346,12 +365,18 @@ We can use the `loo_compare` function to compare our two models on
 expected log predictive density (ELPD) for new data:
 
 ``` r
+
 loo_compare(loo1, loo2)
 ```
 
-         elpd_diff se_diff
-    fit2     0.0       0.0
-    fit1 -5352.0     709.2
+     model elpd_diff se_diff p_worse diag_diff       diag_elpd
+      fit2       0.0     0.0      NA                      <NA>
+      fit1   -5352.0   709.2    1.00           17 k_psis > 0.7
+
+
+    Diagnostic flags present.
+    See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
+    or https://mc-stan.org/loo/reference/loo-glossary.html.
 
 The difference in ELPD is much larger than several times the estimated
 standard error of the difference again indicating that the