# t-TEST in chapter 9 of Field, Miles, & Field, 2012 #Install packages install.packages ("MCMCglmm") #Initiate packages library(MCMCglmm) # 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) ###################### Bayesian t-test using MCMCglmm ###################### model1 = MCMCglmm(Anxiety ~ Group, data = spiderLong, verbose = F, nitt=50000, burnin=3000, thin=10) summary(model1) autocorr(model1$VCV) plot(model1) # 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:') print(round(MUs,2)) # (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") 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 ###################### install.packages("bayesplot"); library(bayesplot) # 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 ###################### install.packages("httpuv") install.packages("xts") install.packages("rstanarm") library(httpuv) library(xts) library(rstanarm) model10 <- stan_glm(Anxiety ~ Group, data = spiderLong, warmup = 3000, iter = 50000, sparse = FALSE, seed = 123) summary(model10, digits = 3) plot(model10)