> 
> # enter data (p. 377)
> Group<-gl(2, 12, labels = c("Picture", "Real Spider"))
> Anxiety<-c(30, 35, 45, 40, 50, 35, 55, 25, 30, 45, 40, 50, 40, 
+            35, 50, 55, 65, 55, 50, 35, 30, 50, 60, 39)
> spiderLong<-data.frame(Group, Anxiety)
> 
> 
> 
> ###################### regular, non Bayesian analyses of the data ######################
> 
> # Conventional t-test as GLM (Field online script)
> 
> t.test.GLM <- lm(Anxiety ~ Group, data = spiderLong)
> summary(t.test.GLM)

Call:
lm(formula = Anxiety ~ Group, data = spiderLong)

Residuals:
   Min     1Q Median     3Q    Max 
 -17.0   -8.5    1.5    8.0   18.0 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        40.000      2.944  13.587 3.53e-12 ***
GroupReal Spider    7.000      4.163   1.681    0.107    
---
Signif. codes:  0 Ô***Õ 0.001 Ô**Õ 0.01 Ô*Õ 0.05 Ô.Õ 0.1 Ô Õ 1

Residual standard error: 10.2 on 22 degrees of freedom
Multiple R-squared:  0.1139,	Adjusted R-squared:  0.07359 
F-statistic: 2.827 on 1 and 22 DF,  p-value: 0.1068

> 
> 
> 
> 
> ###################### Bayesian t-test using MCMCglmm ######################
> 
> model1 = MCMCglmm(Anxiety ~ Group, data = spiderLong, verbose = F,
+                   nitt=50000, burnin=3000, thin=10)
> 
> summary(model1)

 Iterations = 3001:49991
 Thinning interval  = 10
 Sample size  = 4700 

 DIC: 183.4878 

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units     114.7    54.55    189.7     4700

 Location effects: Anxiety ~ Group 

                 post.mean l-95% CI u-95% CI eff.samp  pMCMC    
(Intercept)         40.001   34.217   46.303     4700 <2e-04 ***
GroupReal Spider     6.910   -1.508   15.434     4233  0.097 .  
---
Signif. codes:  0 Ô***Õ 0.001 Ô**Õ 0.01 Ô*Õ 0.05 Ô.Õ 0.1 Ô Õ 1
> autocorr(model1$VCV)
, , units

               units
Lag 0    1.000000000
Lag 10  -0.009766803
Lag 50   0.011863046
Lag 100 -0.012402356
Lag 500 -0.005309054

> plot(model1)
Hit <Return> to see next plot: 
> 
> # running the same analyses using different priors, & then plotting the various posteriors for "Group"
> 
> 
> estimateNum = 2 # the number of regression estimate to be used for the analyses e.g., the 1st (intercept), or 2nd, etc
> 
> Ndatasets = 5   # the number of additional analyses using different priors
> 
> # plot data for the first Bayes model
> densdat  <- density(model1$Sol[,estimateNum])
> plotdatx <- densdat$x
> plotdaty <- densdat$y
> 
> # for MCMCglmm priors, B is for fixed effects. B is a list containing the expected value (mu) 
> # and a (co)varcv    bgiance matrix (V) representing the strength of belief: the defaults are B$mu=0
> # and B$V=I*1e+10, where where I is an identity matrix of appropriate dimension
> 
> # the different priors will be derived from the lm, non-Bayesian regression coefficients
> # the mu values will be random values from a 4 SE range around the lm estimates
> # the V values will be set at (4 * the SEs)**2
> 
> fit.list <- summary.lm(t.test.GLM)
> coeffs <- fit.list$coefficients
> coefnames <- rownames(coeffs)
> 
> # ranges around the lm estimates
> MUranges <- cbind( (coeffs[,1] - (4 * coeffs[,2])), (coeffs[,1] + (4 * coeffs[,2]))  )
> 
> # random values from the MUranges for each coefficient
> MUs <- matrix(-9999,nrow(coeffs),Ndatasets)
> for (lupe in 1:nrow(coeffs)) { MUs[lupe,] <- runif(Ndatasets, min=MUranges[lupe,1], max=MUranges[lupe,2]) }
> dimnames(MUs) <-list(rep("", dim(MUs)[1]))
> colnames(MUs) <- seq(1:Ndatasets)
> cat('\n\nPrior regression coefficient estimates for each additional analysis:')


Prior regression coefficient estimates for each additional analysis:
> print(round(MUs,2))
     1     2     3     4     5
 42.63 40.35 34.86 43.46 34.46
 -4.47 -0.50 21.79  2.04 16.95
> 
> # (co)variance matrix (V) representing the strength of belief - which will be (2 * the lm SEs)**2
> Vlmx2 <- diag((4 * coeffs[,2])**2)
> colnames(Vlmx2) <- seq(1:nrow(MUs))
> cat('\n\nPrior variance (strength of belief) value used for each regression estimate in the additional analyses:\n', 
+     round(diag(Vlmx2),5), sep="\n")


Prior variance (strength of belief) value used for each regression estimate in the additional analyses:

138.6667
277.3333
> 
> 
> for (lupeSets in 1:Ndatasets) {
+ 	prior.list <- list(B = list(mu = MUs[,lupeSets], V = Vlmx2))
+ 	model2 = MCMCglmm(Anxiety ~ Group, data = spiderLong, prior = prior.list,
+                   nitt=50000, burnin=3000, thin=10, verbose=F)
+ 
+ 	densdat  <- density(model2$Sol[,estimateNum])
+ 	plotdatx <- cbind( plotdatx, densdat$x)
+ 	plotdaty <- cbind( plotdaty, densdat$y)
+ }
> 
> matplot( plotdatx, plotdaty, main="", xlab=paste('Estimate for', coefnames[estimateNum]), ylab='Density',  
+          font.lab=1.8, type='l', cex.lab=1.2, col=c(2,rep(1,ncol(plotdatx))), lty=1, lwd=1 )
> title(main=paste('Density Plots of the Posteriors for', coefnames[estimateNum], '\nProduced Using Differing Priors'))
> legend("topright", c("broad priors","random priors"), bty="n", lty=c(1,1),
+        lwd=2, cex=1.2, text.col=c(2,1), col=c(2,1) )
> 
> 
> 
> ################ posterior predictive checks using the bayesplot package ######################
> 
> # from the bayesplot package documentation:
> # The bayesplot package provides various plotting functions for graphical posterior predictive checking, 
> # that is, creating graphical displays comparing observed data to simulated data from the posterior 
> # predictive distribution. The idea behind posterior predictive checking is simple: if a model is a 
> # good fit then we should be able to use it to generate data that looks a lot like the data we observed. 
> 
> y <- spiderLong$Anxiety   # the DV that was used for the MCMCglmm model
> 
> # generating the simulated data using the simulate.MCMCglmm function from the MCMCglmm package
> yrep <- simulate.MCMCglmm(model1, 25) 
> 
> yrep <- t( yrep ) # transposing yrep for the bayesplot functions
> 
> 
> # ppc_stat: A histogram of the distribution of a test statistic computed by applying stat to each 
> #           dataset (row) in yrep. The value of the statistic in the observed data, stat(y), is 
> #           overlaid as a vertical line.
> ppc_stat(y, yrep, binwidth = 1)
> 
> 
> # ppc_dens_overlay: Kernel density or empirical CDF estimates of each dataset (row) in yrep are 
> #                   overlaid, with the distribution of y itself on top (and in a darker shade).
> ppc_dens_overlay(y, yrep[1:25, ])
> 
> 
> # ppc_scatter_avg: A scatterplot of y against the average values of yrep, i.e., the 
> #                  points (mean(yrep[, n]), y[n]), where each yrep[, n] is a vector of length equal 
> #                  to the number of posterior draws.
> ppc_scatter_avg(y, yrep)
> 
> 
> # ppc_hist: A separate histogram, shaded frequency polygon, smoothed kernel density estimate, 
> #           or box and whiskers plot is displayed for y and each dataset (row) in yrep. 
> #           For these plots yrep should therefore contain only a small number of rows.
> ppc_hist(y, yrep[1:8, ], binwidth = .3)
> 
> 
> 
> ###################### now conduct the Bayesian analyses using the rstanarm package ######################
> 
> 
> model10 <- stan_glm(Anxiety ~ Group, data = spiderLong,
+ 	warmup = 3000, iter = 50000, sparse = FALSE, seed = 123)

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).

Gradient evaluation took 3.9e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
Adjust your expectations accordingly!


Iteration:     1 / 50000 [  0%]  (Warmup)
Iteration:  3001 / 50000 [  6%]  (Sampling)
Iteration:  8000 / 50000 [ 16%]  (Sampling)
Iteration: 13000 / 50000 [ 26%]  (Sampling)
Iteration: 18000 / 50000 [ 36%]  (Sampling)
Iteration: 23000 / 50000 [ 46%]  (Sampling)
Iteration: 28000 / 50000 [ 56%]  (Sampling)
Iteration: 33000 / 50000 [ 66%]  (Sampling)
Iteration: 38000 / 50000 [ 76%]  (Sampling)
Iteration: 43000 / 50000 [ 86%]  (Sampling)
Iteration: 48000 / 50000 [ 96%]  (Sampling)
Iteration: 50000 / 50000 [100%]  (Sampling)

 Elapsed Time: 0.083947 seconds (Warm-up)
               1.65063 seconds (Sampling)
               1.73458 seconds (Total)


SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).

Gradient evaluation took 1.4e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
Adjust your expectations accordingly!


Iteration:     1 / 50000 [  0%]  (Warmup)
Iteration:  3001 / 50000 [  6%]  (Sampling)
Iteration:  8000 / 50000 [ 16%]  (Sampling)
Iteration: 13000 / 50000 [ 26%]  (Sampling)
Iteration: 18000 / 50000 [ 36%]  (Sampling)
Iteration: 23000 / 50000 [ 46%]  (Sampling)
Iteration: 28000 / 50000 [ 56%]  (Sampling)
Iteration: 33000 / 50000 [ 66%]  (Sampling)
Iteration: 38000 / 50000 [ 76%]  (Sampling)
Iteration: 43000 / 50000 [ 86%]  (Sampling)
Iteration: 48000 / 50000 [ 96%]  (Sampling)
Iteration: 50000 / 50000 [100%]  (Sampling)

 Elapsed Time: 0.081407 seconds (Warm-up)
               1.54933 seconds (Sampling)
               1.63074 seconds (Total)


SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).

Gradient evaluation took 1.4e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
Adjust your expectations accordingly!


Iteration:     1 / 50000 [  0%]  (Warmup)
Iteration:  3001 / 50000 [  6%]  (Sampling)
Iteration:  8000 / 50000 [ 16%]  (Sampling)
Iteration: 13000 / 50000 [ 26%]  (Sampling)
Iteration: 18000 / 50000 [ 36%]  (Sampling)
Iteration: 23000 / 50000 [ 46%]  (Sampling)
Iteration: 28000 / 50000 [ 56%]  (Sampling)
Iteration: 33000 / 50000 [ 66%]  (Sampling)
Iteration: 38000 / 50000 [ 76%]  (Sampling)
Iteration: 43000 / 50000 [ 86%]  (Sampling)
Iteration: 48000 / 50000 [ 96%]  (Sampling)
Iteration: 50000 / 50000 [100%]  (Sampling)

 Elapsed Time: 0.084535 seconds (Warm-up)
               1.66287 seconds (Sampling)
               1.7474 seconds (Total)


SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).

Gradient evaluation took 1.6e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
Adjust your expectations accordingly!


Iteration:     1 / 50000 [  0%]  (Warmup)
Iteration:  3001 / 50000 [  6%]  (Sampling)
Iteration:  8000 / 50000 [ 16%]  (Sampling)
Iteration: 13000 / 50000 [ 26%]  (Sampling)
Iteration: 18000 / 50000 [ 36%]  (Sampling)
Iteration: 23000 / 50000 [ 46%]  (Sampling)
Iteration: 28000 / 50000 [ 56%]  (Sampling)
Iteration: 33000 / 50000 [ 66%]  (Sampling)
Iteration: 38000 / 50000 [ 76%]  (Sampling)
Iteration: 43000 / 50000 [ 86%]  (Sampling)
Iteration: 48000 / 50000 [ 96%]  (Sampling)
Iteration: 50000 / 50000 [100%]  (Sampling)

 Elapsed Time: 0.096829 seconds (Warm-up)
               1.53848 seconds (Sampling)
               1.63531 seconds (Total)

> 
> summary(model10, digits = 3)

Model Info:

 function:     stan_glm
 family:       gaussian [identity]
 formula:      Anxiety ~ Group
 algorithm:    sampling
 priors:       see help('prior_summary')
 sample:       188000 (posterior sample size)
 observations: 24
 predictors:   2

Estimates:
                   mean     sd       2.5%     25%      50%      75%      97.5% 
(Intercept)        40.067    3.053   34.012   38.068   40.067   42.073   46.079
GroupReal Spider    6.814    4.263   -1.615    4.018    6.808    9.610   15.219
sigma              10.547    1.655    7.901    9.373   10.346   11.499   14.357
mean_PPD           43.470    3.075   37.401   41.464   43.465   45.481   49.540
log-posterior     -97.952    1.294 -101.314  -98.538  -97.617  -97.008  -96.481

Diagnostics:
                 mcse  Rhat  n_eff 
(Intercept)      0.007 1.000 167803
GroupReal Spider 0.011 1.000 164773
sigma            0.004 1.000 145428
mean_PPD         0.007 1.000 175162
log-posterior    0.005 1.000  73524

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
> plot(model10)
> 
>