> 
> # entering data (p. 434)
> id<-(1:15)
> libido<-c(3,2,1,1,4,5,2,4,2,3,7,4,5,3,6)
> dose<-gl(3,5, labels = c("Placebo", "Low Dose", "High Dose"))
> viagraData<-data.frame(dose, libido)
> 
> 
> ###################### regular, non Bayesian analyses of the data ######################
> 
> # conventional ANOVA (p. 439)
> viagraModel <- aov(libido ~ dose, data = viagraData)
> summary.lm(viagraModel)

Call:
aov(formula = libido ~ dose, data = viagraData)

Residuals:
   Min     1Q Median     3Q    Max 
  -2.0   -1.2   -0.2    0.9    2.0 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)     2.2000     0.6272   3.508  0.00432 **
doseLow Dose    1.0000     0.8869   1.127  0.28158   
doseHigh Dose   2.8000     0.8869   3.157  0.00827 **
---
Signif. codes:  0 Ô***Õ 0.001 Ô**Õ 0.01 Ô*Õ 0.05 Ô.Õ 0.1 Ô Õ 1

Residual standard error: 1.402 on 12 degrees of freedom
Multiple R-squared:  0.4604,	Adjusted R-squared:  0.3704 
F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469

> 
> 
> 
> ###################### Bayesian ANOVA using MCMCglmm ######################
> 
> model1 = MCMCglmm(libido ~ dose,  data=viagraData, verbose=F,
+                   nitt=50000, burnin=3000, thin=10)
> summary(model1)

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

 DIC: 57.58071 

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units     2.329   0.8284    4.453     4700

 Location effects: libido ~ dose 

              post.mean l-95% CI u-95% CI eff.samp   pMCMC   
(Intercept)      2.2000   0.8032   3.4958     4700 0.00553 **
doseLow Dose     0.9949  -0.9133   2.8483     4700 0.28553   
doseHigh Dose    2.7779   0.9522   4.7755     4700 0.01064 * 
---
Signif. codes:  0 Ô***Õ 0.001 Ô**Õ 0.01 Ô*Õ 0.05 Ô.Õ 0.1 Ô Õ 1
> autocorr(model1$VCV)
, , units

                units
Lag 0    1.0000000000
Lag 10   0.0127178114
Lag 50  -0.0047286002
Lag 100 -0.0001693583
Lag 500 -0.0277471522

> plot(model1)
Hit <Return> to see next plot: 
> 
> # running the same analyses using different priors, & then plotting the various posteriors for "dose"
> 
> 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(viagraModel)
> 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
 3.04  1.82 4.40 1.49  4.20
 1.61 -0.24 2.70 3.29  2.54
 2.83  2.21 4.86 1.32 -0.47
> 
> # (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:

6.29333
12.58667
12.58667
> 
> 
> for (lupeSets in 1:Ndatasets) {
+ 	prior.list <- list(B = list(mu = MUs[,lupeSets], V = Vlmx2))
+ 	model2 = MCMCglmm(libido ~ dose,  data=viagraData, 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 <- viagraData$libido   # 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(libido ~ dose, data = viagraData,
+                     warmup = 1000, iter = 2000, sparse = FALSE, seed = 123)

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

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 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.034771 seconds (Warm-up)
               0.032755 seconds (Sampling)
               0.067526 seconds (Total)


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

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


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.031963 seconds (Warm-up)
               0.045802 seconds (Sampling)
               0.077765 seconds (Total)


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

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


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.033683 seconds (Warm-up)
               0.047593 seconds (Sampling)
               0.081276 seconds (Total)


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

Gradient evaluation took 9e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.035613 seconds (Warm-up)
               0.053674 seconds (Sampling)
               0.089287 seconds (Total)

> 
> summary(model10, digits = 3)

Model Info:

 function:     stan_glm
 family:       gaussian [identity]
 formula:      libido ~ dose
 algorithm:    sampling
 priors:       see help('prior_summary')
 sample:       4000 (posterior sample size)
 observations: 15
 predictors:   3

Estimates:
                mean    sd      2.5%    25%     50%     75%     97.5%
(Intercept)     2.273   0.681   0.933   1.842   2.271   2.695   3.610
doseLow Dose    0.909   0.954  -0.991   0.302   0.906   1.513   2.827
doseHigh Dose   2.666   0.926   0.778   2.089   2.681   3.264   4.424
sigma           1.496   0.329   1.004   1.263   1.442   1.662   2.296
mean_PPD        3.456   0.568   2.306   3.112   3.460   3.814   4.572
log-posterior -33.898   1.624 -37.995 -34.698 -33.534 -32.699 -31.893

Diagnostics:
              mcse  Rhat  n_eff
(Intercept)   0.013 1.001 2624 
doseLow Dose  0.018 1.001 2728 
doseHigh Dose 0.018 1.000 2525 
sigma         0.007 1.003 2319 
mean_PPD      0.009 1.000 3754 
log-posterior 0.047 1.003 1218 

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)
> 
> 
>