large-data.Rmd: specicy r_eff argument

Author: jgabry

vignettes/large-data.Rmd

@@ -91,10 +91,12 @@ log-likelihood values in memory.
 
 The log-likelihood in R can be coded as follows:
 ```{r llfun_logistic}
-llfun_logistic <- function(data_i, draws) {
+# we'll add an argument log to toggle whether this is a log-likelihood or 
+# likelihood function. this will be useful later in the vignette.
+llfun_logistic <- function(data_i, draws, log = TRUE) {
   x_i <- as.matrix(data_i[, which(grepl(colnames(data_i), pattern = "X")), drop=FALSE])
   logit_pred <- draws %*% t(x_i)
-  dbinom(x = data_i$y, size = 1, prob = 1/(1 + exp(-logit_pred)), log = TRUE)
+  dbinom(x = data_i$y, size = 1, prob = 1/(1 + exp(-logit_pred)), log = log)
 }
 ```
 
@@ -133,15 +135,28 @@ function we wrote is working as it should. The `loo_i()` function is a helper
 function that can be used to test a log-likelihood function on a single observation.
 
 ```{r, eval=FALSE}
+# used for draws argument to loo_i
 parameter_draws_1 <- extract(fit_1)$beta
+
+# used for data argument to loo_i
 stan_df_1 <- as.data.frame(standata)
-loo_i(1, llfun_logistic, data = stan_df_1, draws = parameter_draws_1)
+
+# compute relative efficiency (this is slow and optional but is recommended to allow 
+# for adjusting PSIS effective sample size based on MCMC effective sample size)
+r_eff <- relative_eff(llfun_logistic, 
+                      log = FALSE, # relative_eff wants likelihood not log-likelihood values
+                      chain_id = rep(1:4, each = 1000), 
+                      data = stan_df_1, 
+                      draws = parameter_draws_1, 
+                      cores = 2)
+
+loo_i(i = 1, llfun_logistic, r_eff = r_eff, data = stan_df_1, draws = parameter_draws_1)
 ```
 
 ```
 $pointwise
-    elpd_loo mcse_elpd_loo        p_loo     looic
-1 -0.3310342  0.0002908997 0.0003487243 0.6620683
+    elpd_loo mcse_elpd_loo        p_loo     looic influence_pareto_k
+1 -0.3314552  0.0002887608 0.0003361772 0.6629103        -0.05679886
 ...
 ```
 
@@ -151,16 +166,16 @@ We can then use the `loo_subsample()` function to compute the efficient PSIS-LOO
 approximation to exact LOO-CV using subsampling:
 
 ```{r, eval=FALSE}
-parameter_draws_1 <- extract(fit_1)$beta
-stan_df_1 <- as.data.frame(standata)
-
 set.seed(4711)
 loo_ss_1 <-
   loo_subsample(
     llfun_logistic,
+    observations = 100, # take a subsample of size 100
+    cores = 2,
+    # these next objects were computed above
+    r_eff = r_eff, 
     draws = parameter_draws_1,
-    data = stan_df_1,
-    observations = 100 # take a subsample of size 100
+    data = stan_df_1
   )
 print(loo_ss_1)
 ```
@@ -205,8 +220,14 @@ simply add more samples until we are satisfied using the `update()` method.
 
 ```{r, eval=FALSE}
 set.seed(4711)
-loo_ss_1b <- update(loo_ss_1, draws = parameter_draws_1, data = stan_df_1, 
-                    observations = 200) # subsample 200 instead of 100
+loo_ss_1b <-
+  update(
+    loo_ss_1,
+    observations = 200, # subsample 200 instead of 100
+    r_eff = r_eff,
+    draws = parameter_draws_1,
+    data = stan_df_1
+  ) 
 print(loo_ss_1b)
 ```
 
@@ -241,12 +262,14 @@ set.seed(4711)
 loo_ss_1c <-
   loo_subsample(
     x = llfun_logistic,
+    r_eff = r_eff,
     draws = parameter_draws_1,
     data = stan_df_1,
     observations = 100,
     estimator = "hh_pps", # use Hansen-Hurwitz
     loo_approximation = "lpd", # use lpd instead of plpd
-    loo_approximation_draws = 100
+    loo_approximation_draws = 100,
+    cores = 2
   )
 print(loo_ss_1c)
 ```
@@ -297,7 +320,8 @@ loo_ap_1 <-
     draws = parameter_draws_laplace,
     data = stan_df_1,
     log_p = log_p,
-    log_g = log_g
+    log_g = log_g,
+    cores = 2
   )
 print(loo_ap_1)
 ```
@@ -340,7 +364,8 @@ loo_ap_ss_1 <-
     data = stan_df_1,
     log_p = log_p,
     log_g = log_g,
-    observations = 100
+    observations = 100,
+    cores = 2
   )
 print(loo_ap_ss_1)
 ```
@@ -387,23 +412,46 @@ parameter_draws_2 <- extract(fit_2)$beta
 stan_df_2 <- as.data.frame(standata)
 
 # recompute subsampling loo for first model for demonstration purposes
+
+# compute relative efficiency (this is slow and optional but is recommended to allow 
+# for adjusting PSIS effective sample size based on MCMC effective sample size)
+r_eff_1 <- relative_eff(
+  llfun_logistic,
+  log = FALSE, # relative_eff wants likelihood not log-likelihood values
+  chain_id = rep(1:4, each = 1000),
+  data = stan_df_1,
+  draws = parameter_draws_1,
+  cores = 2
+)
+
 set.seed(4711)
-loo_ss_1 <-
-  loo_subsample(
-    x = llfun_logistic,
-    draws = parameter_draws_1,
-    data = stan_df_1,
-    observations = 200
-  )
+loo_ss_1 <- loo_subsample(
+  x = llfun_logistic,
+  r_eff = r_eff_1,
+  draws = parameter_draws_1,
+  data = stan_df_1,
+  observations = 200,
+  cores = 2
+)
 
 # compute subsampling loo for a second model (with log-arsenic)
-loo_ss_2 <-
-  loo_subsample(
-    x = llfun_logistic,
-    draws = parameter_draws_2,
-    data = stan_df_2,
-    observations = 200
-  )
+
+r_eff_2 <- relative_eff(
+  llfun_logistic,
+  log = FALSE, # relative_eff wants likelihood not log-likelihood values
+  chain_id = rep(1:4, each = 1000),
+  data = stan_df_2,
+  draws = parameter_draws_2,
+  cores = 2
+)
+loo_ss_2 <- loo_subsample(
+  x = llfun_logistic,
+  r_eff = r_eff_2, 
+  draws = parameter_draws_2,
+  data = stan_df_2,
+  observations = 200,
+  cores = 2
+)
 
 print(loo_ss_2)
 ```
@@ -455,9 +503,11 @@ we can simply extract the observations used in `loo_ss_1` and use them in
 loo_ss_2 <-
   loo_subsample(
     x = llfun_logistic,
+    r_eff = r_eff_2,
     draws = parameter_draws_2,
     data = stan_df_2,
-    observations = loo_ss_1
+    observations = loo_ss_1,
+    cores = 2
   )
 ```
 
@@ -466,13 +516,14 @@ helper function:
 
 ```{r, eval=FALSE}
 idx <- obs_idx(loo_ss_1)
-loo_ss_2 <-
-  loo_subsample(
-    x = llfun_logistic,
-    draws = parameter_draws_2,
-    data = stan_df_2,
-    observations = idx
-  )
+loo_ss_2 <- loo_subsample(
+  x = llfun_logistic,
+  r_eff = r_eff_2, 
+  draws = parameter_draws_2,
+  data = stan_df_2,
+  observations = idx,
+  cores = 2
+)
 ```
 
 ```
@@ -510,7 +561,13 @@ It is also possible to compare a subsampled loo computation with a full loo obje
 
 ```{r, eval=FALSE}
 # use loo() instead of loo_subsample() to compute full PSIS-LOO for model 2
-loo_full_2 <- loo(x = llfun_logistic, draws = parameter_draws_2, data = stan_df_2)
+loo_full_2 <- loo(
+  x = llfun_logistic,
+  r_eff = r_eff_2,
+  draws = parameter_draws_2,
+  data = stan_df_2,
+  cores = 2
+)
 loo_compare(loo_ss_1, loo_full_2)
 ```