Fix demos to use data frames

Closes #188

Author: jgabry

demo/ARM_Ch03.R

@@ -2,13 +2,14 @@
 demo("SETUP", package = "rstanarm", verbose = FALSE, echo = FALSE, ask = FALSE)
 
 source(paste0(ROOT, "ARM/Ch.3/kidiq.data.R"), local = DATA_ENV, verbose = FALSE)
+dat <- with(DATA_ENV, data.frame(kid_score, mom_hs, mom_iq))
 
 # Estimate four contending models
-post1 <- stan_glm(kid_score ~ mom_hs, data = DATA_ENV, 
+post1 <- stan_glm(kid_score ~ mom_hs, data = dat, 
                   family = gaussian(link = "identity"), 
                   seed = SEED, refresh = REFRESH)
 post2 <- update(post1, formula = kid_score ~ mom_iq)
-post3 <- stan_lm(kid_score ~ mom_hs + mom_iq, data = DATA_ENV,
+post3 <- stan_lm(kid_score ~ mom_hs + mom_iq, data = dat,
                  prior = R2(location = 0.25, what = "mean"), 
                  seed = SEED, refresh = REFRESH)
 post4 <- update(post3, formula = kid_score ~ mom_hs * mom_iq,
@@ -45,15 +46,22 @@ source(paste0(ROOT, "ARM/Ch.3/kids_before1987.data.R"),
        local = DATA_ENV, verbose = FALSE)
 source(paste0(ROOT, "ARM/Ch.3/kids_after1987.data.R"), 
        local = DATA_ENV, verbose = FALSE)
-post5 <- stan_lm(ppvt ~ hs + afqt, data = DATA_ENV,
+
+fit_data <- with(DATA_ENV, data.frame(ppvt, hs, afqt))
+pred_data <- with(DATA_ENV, data.frame(ppvt_ev, hs_ev, afqt_ev))
+
+post5 <- stan_lm(ppvt ~ hs + afqt, data = fit_data,
                  prior = R2(location = 0.25, what = "mean"), 
                  seed = SEED, refresh = REFRESH)
-y_ev <- posterior_predict(post5, newdata = 
-                          data.frame(hs = DATA_ENV$hs_ev, afqt = DATA_ENV$afqt_ev))
+y_ev <- posterior_predict(
+  post5, 
+  newdata = with(pred_data, data.frame(hs = hs_ev, afqt = afqt_ev))
+)
 par(mfrow = c(1,1))
-hist(-sweep(y_ev, 2, STATS = DATA_ENV$ppvt_ev, FUN = "-"), prob = TRUE,
+hist(-sweep(y_ev, 2, STATS = pred_data$ppvt_ev, FUN = "-"), prob = TRUE,
      xlab = "Predictive Errors in ppvt", main = "", las = 2)
 
+
 ANSWER <- tolower(readline("Do you want to remove the objects this demo created? (y/n) "))
 if (ANSWER != "n") {
   rm(IQ_SEQ, y_nohs, y_hs, y_ev, ANSWER)

demo/ARM_Ch04.R

@@ -2,6 +2,7 @@
 demo("SETUP", package = "rstanarm", verbose = FALSE, echo = FALSE, ask = FALSE)
 
 source(paste0(ROOT, "ARM/Ch.4/earnings.data.R"), local = DATA_ENV, verbose = FALSE)
+earnings_dat <- with(DATA_ENV, data.frame(earn, height, male))
 
 # The stuff in sections 4.0 -- 4.3 is not very relevant 
 # Moreover, centering predictors is NOT recommended in the rstanarm package
@@ -10,16 +11,16 @@ source(paste0(ROOT, "ARM/Ch.4/earnings.data.R"), local = DATA_ENV, verbose = FAL
 
 # These two models are essentially equivalent in the likelihood
 # But the "same" priors affect the posterior differently
-post1 <- stan_glm(log(earn) ~ height, data = DATA_ENV, 
+post1 <- stan_glm(log(earn) ~ height, data = earnings_dat, 
                   family = gaussian(link = "identity"), 
                   seed = SEED, refresh = REFRESH)
-# post2 <- stan_glm(earn ~ height, data = DATA_ENV, 
+# post2 <- stan_glm(earn ~ height, data = earnings_dat, 
 #                   family = gaussian(link = "log"), 
 #                   seed = SEED, refresh = REFRESH)
 # and this does not even converge
 
 # These models add terms to the right-hand side
-post3 <- stan_lm(log(earn) ~ height + male, data = DATA_ENV,
+post3 <- stan_lm(log(earn) ~ height + male, data = earnings_dat,
                  prior = R2(location = 0.3, what = "mean"), 
                  seed = SEED, refresh = REFRESH)
 post4 <- update(post3, formula = log(earn) ~ height * male)
@@ -49,15 +50,18 @@ boxplot(y_men, outline = FALSE, col = "red", axes = FALSE, log = "y", ylim = YLI
 axis(1, at = 1:ncol(y_men), labels = MEN_SEQ, las = 3)
 
 # Prediction of the weight of mesquite trees 
+DATA_ENV <- new.env()
 source(paste0(ROOT, "ARM/Ch.4/mesquite.data.R"), local = DATA_ENV, verbose = FALSE)
+tree_dat <- as.data.frame(do.call(cbind, as.list(DATA_ENV)))
+
 CONTINUE1 <- tolower(readline(
   paste("A heads up: the next part of the demo (Predicting weight of mesquite trees )",
         "prints many lines \nto the console as it runs many models and compares the results", 
         "Proceed? (y/n)")
 ))
 if (CONTINUE1 != "n") {
   post5 <- stan_lm(weight ~ diam1 + diam2 + canopy_height + total_height +
-                     density + group, data = DATA_ENV,
+                     density + group, data = tree_dat,
                    prior = R2(0.9), seed = SEED, refresh = REFRESH)
   post6 <- update(post5, formula = log(weight) ~ log(diam1) + log(diam2) + 
                     log(canopy_height) + log(total_height) + log(density) + group)
@@ -80,11 +84,13 @@ CONTINUE2 <- tolower(readline(
         "Proceed? (y/n)")
 ))
 if (CONTINUE2 != "n") {
-  YEARS <- as.character(seq(from = 1972, to = 2000, by = 4))
+  YEARS <- as.character(seq(from = 1972, to = 1980, by = 4))
   round(digits = 2, x = sapply(YEARS, FUN = function(YEAR) {
+    DATA_ENV <- new.env()
     source(paste0(ROOT, "ARM/Ch.4/nes", YEAR, ".data.R"), local = DATA_ENV, verbose = FALSE)
+    pid_dat <- as.data.frame(do.call(cbind, as.list(DATA_ENV)))
     coef(stan_lm(partyid7 ~ real_ideo + I(race_adj == 1) + as.factor(age_discrete) + 
-                   educ1 + gender + income, data = DATA_ENV, prior = R2(0.5),
+                   educ1 + gender + income, data = pid_dat, prior = R2(0.5),
                  seed = SEED, refresh = 0))
   }))
 }

demo/ARM_Ch07.R

@@ -2,28 +2,32 @@
 demo("SETUP", package = "rstanarm", verbose = FALSE, echo = FALSE, ask = FALSE)
 
 source(paste0(ROOT, "ARM/Ch.7/congress.data.R"), local = DATA_ENV, verbose = FALSE)
+cong_dat <- with(DATA_ENV, data.frame(incumbency_88, vote_88, vote_86))
 
 # The stuff in sections 7.0 -- 7.2 is not very relevant 
 
-post1 <- stan_lm(vote_88 ~ vote_86 + incumbency_88, data = DATA_ENV, 
+post1 <- stan_lm(vote_88 ~ vote_86 + incumbency_88, data = cong_dat, 
                  prior = R2(0.9, what = "mean"), 
                  seed = SEED, refresh = REFRESH)
 post1 # badly underfitting
 y_tilde <- posterior_predict(post1) # incumbency_90 is not available
 summary(rowSums(y_tilde > 0.5))
 
-source(paste0(ROOT, "ARM/Ch.6/wells.data.R"), local = DATA_ENV, verbose = FALSE)
-post2 <- stan_glm(switch ~ I(dist / 100), data = DATA_ENV, family = "binomial", 
+
+data(wells, package = "rstanarm")
+wells$dist100 <- with(wells, dist / 100)
+post2 <- stan_glm(switch ~ dist100, data = wells, family = "binomial", iter = 100, chains = 1,
                   seed = SEED, refresh = REFRESH)
-prop.table(table(c(ppd)))
+prop.table(table(c(posterior_predict(post2))))
+
 
 # the compound model is not good because it assumes the two errors are 
 # independent. rstanarm will eventually support Heckman models, which
 # would be a better choice here.
 
 ANSWER <- tolower(readline("Do you want to remove the objects this demo created? (y/n) "))
 if (ANSWER != "n") {
-  rm(y_tilde, ppd, ANSWER)
+  rm(y_tilde, wells, ANSWER)
   # removes stanreg and loo objects, plus what was created by STARTUP
   demo("CLEANUP", package = "rstanarm", verbose = FALSE, echo = FALSE, ask = FALSE)
 }

demo/ARM_Ch08.R

@@ -2,10 +2,10 @@
 demo("SETUP", package = "rstanarm", verbose = FALSE, echo = FALSE, ask = FALSE)
 
 source(paste0(ROOT, "ARM/Ch.8/lightspeed.data.R"), local = DATA_ENV, verbose = FALSE)
-
+light_dat <- with(DATA_ENV, data.frame(y))
 # The stuff in sections 8.0 -- 8.2 is not very relevant 
 
-(post1 <- stan_glm(y ~ 1, data = DATA_ENV, seed = SEED, refresh = REFRESH))
+(post1 <- stan_glm(y ~ 1, data = light_dat, seed = SEED, refresh = REFRESH))
 y_rep <- posterior_predict(post1)
 
 pp_check(post1, plotfun = "stat", stat = "min") + 
@@ -25,20 +25,21 @@ pp_check(post1, plotfun = "hist") + ggtitle(ttl)
 # Make similar plot manually but combine all y_rep
 op <- par('mfrow')
 par(mfrow = 1:2, mar = c(5,4,1,1) + .1)
-hist(DATA_ENV$y, prob = TRUE, main = "", las = 1,
+hist(light_dat$y, prob = TRUE, main = "", las = 1,
      xlab = "Measurement Error for the Speed of Light")
 hist(y_rep, prob = TRUE, main = "", las = 1,
      xlab = "Predicted Measurement Error")
 par(mfrow = op)
 
-source(paste0(ROOT, "ARM/Ch.8/roaches.data.R"), local = DATA_ENV, verbose = FALSE)
-post2 <- stan_glm(y ~ roach1 + treatment + senior, data = DATA_ENV, 
+# Roaches example
+data(roaches, package = "rstanarm")
+post2 <- stan_glm(y ~ roach1 + treatment + senior, data = roaches, 
                   family = poisson(link = "log"), seed = SEED, refresh = REFRESH)
 y_rep <- posterior_predict(post2)
 
 # Compare observed proportion of zeros to predicted proportion of zeros
 mean(y_rep == 0)
-mean(DATA_ENV$y == 0)
+mean(roaches$y == 0)
 summary(apply(y_rep == 0, 1, mean))
 prop0 <- function(x) mean(x == 0)
 pp_check(post2, plotfun = "stat", stat = "prop0") # model doesn't predict enough zeros

demo/ARM_Ch09.R

@@ -3,30 +3,32 @@ demo("SETUP", package = "rstanarm", verbose = FALSE, echo = FALSE, ask = FALSE)
 # read data into DATA_ENV environment
 source(paste0(ROOT, "ARM/Ch.9/electric_grade4.data.R"), local = DATA_ENV, 
        verbose = FALSE)
+dat <- with(DATA_ENV, data.frame(post_test, grade, pre_test, treatment))
 
-post1 <- stan_lm(post_test ~ treatment * pre_test, data = DATA_ENV, 
+post1 <- stan_lm(post_test ~ treatment * pre_test, data = dat, 
                  prior = R2(0.75), seed = SEED, refresh = REFRESH)
 post1 # underfitting but ok because it is an experiment
 plot(post1)
 
-y_0 <- posterior_predict(post1, data.frame(treatment = 0, pre_test = DATA_ENV$pre_test))
-y_1 <- posterior_predict(post1, data.frame(treatment = 1, pre_test = DATA_ENV$pre_test))
+y_0 <- posterior_predict(post1, data.frame(treatment = 0, pre_test = dat$pre_test))
+y_1 <- posterior_predict(post1, data.frame(treatment = 1, pre_test = dat$pre_test))
 diff <- y_1 - y_0
 mean(diff)
 sd(diff) # much larger than in ARM
 hist(diff, prob = TRUE, main = "", xlab = "Estimated Average Treatment Effect", las = 1)
 
 
-stopifnot(require(gridExtra))
+stopifnot(require(bayesplot))
 plots <- sapply(1:4, simplify = FALSE, FUN = function(k) {
-  source(paste0(ROOT, "ARM/Ch.9/electric_grade", k, "_supp.data.R"), 
-         local = DATA_ENV, verbose = FALSE)
-  out <- plot(stan_lm(post_test ~ supp + pre_test, data = DATA_ENV, 
+  dat$supp <-
+    source(paste0(ROOT, "ARM/Ch.9/electric_grade", k, "_supp.data.R"),
+           verbose = FALSE)$value
+  out <- plot(stan_lm(post_test ~ supp + pre_test, data = dat, 
                     seed = SEED, refresh = REFRESH,
                     prior = R2(0.75, what = "mean")))
   out + ggtitle(paste("Grade =", k))
 })
-marrangeGrob(plots, nrow = 2, ncol = 2)
+bayesplot_grid(plots = plots, grid_args = list(nrow = 2, ncol = 2))
 
 ANSWER <- tolower(readline("Do you want to remove the objects this demo created? (y/n) "))
 if (ANSWER != "n") {