Jump to main navigation


Workshop 9.6b - Factorial ANOVA (Bayesian)

23 April 2011

Factorial ANOVA references

  • McCarthy (2007) - Chpt 6
  • Kery (2010) - Chpt 10
  • Gelman & Hill (2007) - Chpt 4
  • Logan (2010) - Chpt 12
  • Quinn & Keough (2002) - Chpt 9
View R code for preliminaries.
> library(R2jags)
> library(ggplot2)
> library(grid)
> #define my common ggplot options
> murray_opts <- opts(panel.grid.major=theme_blank(),
+                            panel.grid.minor=theme_blank(),
+                            panel.border = theme_blank(),
+                            panel.background = theme_blank(),
+                            axis.title.y=theme_text(size=15, vjust=0,angle=90),
+                            axis.text.y=theme_text(size=12),
+                            axis.title.x=theme_text(size=15, vjust=-1),
+                            axis.text.x=theme_text(size=12),
+                            axis.line = theme_segment(),
+                            legend.position=c(1,0.1),legend.justification=c(1,0),
+                            legend.key=theme_blank(),
+                            plot.margin=unit(c(0.5,0.5,1,2),"lines")
+ )
> 
> #create a convenience function for predicting from mcmc objects
> predict.mcmc <- function(mcmc, terms=NULL, mmat, trans="identity") {
+   library(plyr)
+   library(coda)
+   if(is.null(terms)) mcmc.coefs <- mcmc
+   else  {
+      iCol<-match(terms, colnames(mcmc))
+          mcmc.coefs <- mcmc[,iCol]
+   }
+   t2=")"
+   if (trans=="identity") {
+         t1 <- ""
+         t2=""
+   }else if (trans=="log") {
+         t1="exp("
+   }else if (trans=="log10") {
+         t1="10^("
+   }else if (trans=="sqrt") {
+         t1="("
+         t2=")^2"
+   }
+   eval(parse(text=paste("mcmc.sum <-adply(mcmc.coefs %*% t(mmat),2,function(x) {
	data.frame(mean=",t1,"mean(x)",t2,",
    median=",t1,"median(x)",t2,",									  
	",t1,"coda:::HPDinterval(as.mcmc(x), p=0.68)",t2,",
	",t1,"coda:::HPDinterval(as.mcmc(x))",t2,"
    )})", sep="")))
+ 
+   mcmc.sum
+ }
> 
> 
> #create a convenience function for summarizing mcmc objects
> summary.mcmc <- function(mcmc, terms=NULL, trans="identity") {
+   library(plyr)
+   library(coda)
+   if(is.null(terms)) mcmc.coefs <- mcmc
+   else  {
+      iCol<-grep(terms, colnames(mcmc))
+          mcmc.coefs <- mcmc[,iCol]
+   }
+   t2=")"
+   if (trans=="identity") {
+         t1 <- ""
+         t2=""
+   }else if (trans=="log") {
+         t1="exp("
+   }else if (trans=="log10") {
+         t1="10^("
+   }else if (trans=="sqrt") {
+         t1="("
+         t2=")^2"
+   }
+   eval(parse(text=paste("mcmc.sum <-adply(mcmc.coefs,2,function(x) {
	data.frame(mean=",t1,"mean(x)",t2,",
    median=",t1,"median(x)",t2,",									  
	",t1,"coda:::HPDinterval(as.mcmc(x), p=0.68)",t2,",
	",t1,"coda:::HPDinterval(as.mcmc(x))",t2,"
    )})", sep="")))
+ 
+   mcmc.sum$Perc <-100*mcmc.sum$median/sum(mcmc.sum$median)
+   mcmc.sum
+ }

Two-factor ANOVA

A biologist studying starlings wanted to know whether the mean mass of starlings differed according to different roosting situations. She was also interested in whether the mean mass of starlings altered over winter (Northern hemisphere) and whether the patterns amongst roosting situations were consistent throughout winter, therefore starlings were captured at the start (November) and end of winter (January). Ten starlings were captured from each roosting situation in each season, so in total, 80 birds were captured and weighed.

Download Starling data set
Format of starling.csv data files
SITUATIONMONTHMASSGROUP
S1November78S1Nov
........
S2November78S2Nov
........
S3November79S3Nov
........
S4November77S4Nov
........
S1January85S1Jan
........
SITUATIONCategorical listing of roosting situations
MONTHCategorical listing of the month of sampling.
MASSMass (g) of starlings.
GROUPCategorical listing of situation/month combinations - used for checking ANOVA assumptions
Starlings

Open the starling data file.

Show code
> starling <- read.table("../downloads/data/starling.csv",
+     header = T, sep = ",", strip.white = T)
> head(starling)
  SITUATION    MONTH MASS GROUP
1        S1 November   78 S1Nov
2        S1 November   88 S1Nov
3        S1 November   87 S1Nov
4        S1 November   88 S1Nov
5        S1 November   83 S1Nov
6        S1 November   82 S1Nov

Exploratory data analysis did not reveal any issues with normality or homogeneity of variance.

  1. Fit the JAGS model including all pairwise comparisons AND (just because we can), compare the bare surfaces to the surfaces with an algal biofilm. Use with either:
    • Effects parameterization
      massi = β0 + β1j(i)*monthi + β2j(i)*situationi + β3j(i)*monthi*situationi + εi where ε ∼ N(0,σ2)
      View full code
      > modelString="
      model {
         #Likelihood
         for (i in 1:n) {
            mass[i]~dnorm(mean[i],tau.res)
            mean[i] <- beta0+beta.month[month[i]]+beta.situation[situation[i]]+
      	             beta.int[month[i],situation[i]]
            y.err[i] <- mass[i] - mean[i]
         }
      
         #Priors and derivatives
         beta0 ~ dnorm(0,1.0E-6)
         beta.month[1] <- 0
         beta.int[1,1] <- 0
         for (i in 2:nMonth) {
           beta.month[i] ~ dnorm(0.001, 1.0E-6) #prior
           beta.int[i,1] <- 0
         }
         beta.situation[1] <- 0
         for (i in 2:nSituation){
           beta.situation[i] ~ dnorm(0.001, 1.0E-6) #prior cannot be 0, or else will have an 
           ##init of zero and then will be dividing by zero when calculating cohen's D
           beta.int[1,i] <- 0
         }
         for (i in 2:nMonth) {
           for (j in 2:nSituation){
             beta.int[i,j] ~ dnorm(0, 1.0E-6)
           }
         }
      
         #Means
         for (i in 1:nMonth) {
           for (j in 1:nSituation) {
             Mean[i,j] <- beta0 + beta.month[i] + beta.situation[j] + beta.int[i,j]
           }
         } 
      
         sigma.res ~ dgamma(0.001, 0.001) #prior on sd residuals
         tau.res <- 1 / (sigma.res * sigma.res)
         sd.Month <- sd(beta.month[])
         sd.Situation <- sd(beta.situation[])
         sd.Int <- sd(beta.int[,])
         sd.Resid <- sd(y.err[])
      
         #Pairwise comparisons of situations
         for (i in 1:nSituation) {
           for (j in 1:nSituation) {
      	   eff[i,j] <- (beta.situation[i] - beta.situation[j])
      	 }
         }
         #Effects of month within each situation
         for (i in 1:nSituation) {
            eff2[i] <- Mean[1,i] - Mean[2,i]
         }
       }
      "
      > ## write the model to a text file (I suggest you alter the path to somewhere more relevant 
      > ## to your system!)  
      > writeLines(modelString,con="../downloads/BUGSscripts/ws9.6bQ1.1a.txt")
      > starling.list <- with(starling,
      +                          list(mass=MASS,
      +                          month=as.numeric(MONTH), situation=as.numeric(SITUATION),
      +                                  n=nrow(starling), nMonth=length(levels(MONTH)),
      +                                  nSituation=length(levels(SITUATION))
      +                         )
      +            )
      > 
      > params <- c("beta0","beta.month","beta.situation","beta.int","sigma.res","sd.Month","sd.Situation","sd.Int","sd.Resid", "Mean","eff","eff2")
      > adaptSteps = 1000
      > burnInSteps = 200
      > nChains = 3
      > numSavedSteps = 5000
      > thinSteps = 10
      > nIter = ceiling((numSavedSteps * thinSteps)/nChains)
      > 
      > library(R2jags)
      > ## 
      > starling.r2jags <- jags(data=starling.list,
      +           inits=NULL, #or inits=list(inits,inits,inits) # since there are three chains
      +           parameters.to.save=params,
      +           model.file="../downloads/BUGSscripts/ws9.6bQ1.1a.txt",
      +           n.chains=3,
      +           n.iter=nIter,
      +           n.burnin=burnInSteps,
      +       n.thin=thinSteps
      +           )
      
      Compiling model graph
         Resolving undeclared variables
         Allocating nodes
         Graph Size: 367
      
      Initializing model
      
      
      > print(starling.r2jags)
      
      Inference for Bugs model at "../downloads/BUGSscripts/ws9.6bQ1.1a.txt", fit using jags,
       3 chains, each with 16667 iterations (first 200 discarded), n.thin = 10
       n.sims = 4941 iterations saved
                        mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
      Mean[1,1]          90.770   1.368  88.088  89.864  90.750  91.695  93.470 1.001  4900
      Mean[2,1]          83.630   1.342  80.970  82.743  83.623  84.520  86.276 1.001  4900
      Mean[1,2]          90.159   1.348  87.482  89.264  90.161  91.045  92.819 1.001  4900
      Mean[2,2]          79.401   1.348  76.754  78.510  79.401  80.284  82.053 1.001  4900
      Mean[1,3]          88.211   1.348  85.561  87.329  88.206  89.106  90.850 1.001  4200
      Mean[2,3]          78.593   1.341  75.931  77.692  78.588  79.492  81.244 1.001  4900
      Mean[1,4]          84.206   1.354  81.575  83.272  84.202  85.124  86.843 1.002  2300
      Mean[2,4]          75.438   1.353  72.744  74.560  75.443  76.336  78.079 1.001  4900
      beta.int[1,1]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,1]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[1,2]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,2]      -3.618   2.741  -9.114  -5.424  -3.608  -1.854   1.970 1.001  4900
      beta.int[1,3]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,3]      -2.478   2.745  -7.906  -4.302  -2.465  -0.639   2.888 1.001  4800
      beta.int[1,4]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,4]      -1.628   2.720  -6.953  -3.497  -1.603   0.259   3.625 1.001  4900
      beta.month[1]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.month[2]      -7.140   1.943 -10.964  -8.441  -7.137  -5.837  -3.421 1.001  4900
      beta.situation[1]   0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.situation[2]  -0.611   1.917  -4.434  -1.894  -0.582   0.672   3.168 1.001  4900
      beta.situation[3]  -2.559   1.900  -6.357  -3.805  -2.557  -1.335   1.222 1.001  2800
      beta.situation[4]  -6.564   1.938 -10.448  -7.838  -6.562  -5.261  -2.807 1.001  4900
      beta0              90.770   1.368  88.088  89.864  90.750  91.695  93.470 1.001  4900
      eff[1,1]            0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[2,1]           -0.611   1.917  -4.434  -1.894  -0.582   0.672   3.168 1.001  4900
      eff[3,1]           -2.559   1.900  -6.357  -3.805  -2.557  -1.335   1.222 1.001  2800
      eff[4,1]           -6.564   1.938 -10.448  -7.838  -6.562  -5.261  -2.807 1.001  4900
      eff[1,2]            0.611   1.917  -3.168  -0.672   0.582   1.894   4.434 1.001  4900
      eff[2,2]            0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[3,2]           -1.948   1.906  -5.640  -3.226  -1.960  -0.673   1.807 1.001  3200
      eff[4,2]           -5.953   1.897  -9.722  -7.190  -5.958  -4.696  -2.243 1.001  3800
      eff[1,3]            2.559   1.900  -1.222   1.335   2.557   3.805   6.357 1.001  2800
      eff[2,3]            1.948   1.906  -1.807   0.673   1.960   3.226   5.640 1.001  3200
      eff[3,3]            0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[4,3]           -4.005   1.896  -7.681  -5.283  -3.979  -2.766  -0.270 1.002  2000
      eff[1,4]            6.564   1.938   2.807   5.261   6.562   7.838  10.448 1.001  4900
      eff[2,4]            5.953   1.897   2.243   4.696   5.958   7.190   9.722 1.001  3800
      eff[3,4]            4.005   1.896   0.270   2.766   3.979   5.283   7.681 1.002  2000
      eff[4,4]            0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff2[1]             7.140   1.943   3.421   5.837   7.137   8.441  10.964 1.004  4900
      eff2[2]            10.758   1.931   6.988   9.493  10.754  12.027  14.569 1.001  4900
      eff2[3]             9.618   1.918   5.815   8.290   9.647  10.909  13.365 1.001  4900
      eff2[4]             8.768   1.878   5.189   7.496   8.742   9.981  12.466 1.002  4900
      sd.Int              1.803   0.803   0.516   1.206   1.702   2.306   3.631 1.001  4900
      sd.Month            3.570   0.972   1.710   2.918   3.568   4.221   5.482 1.004  4900
      sd.Resid            4.185   0.108   4.035   4.106   4.166   4.242   4.446 1.002  1700
      sd.Situation        2.741   0.657   1.475   2.296   2.735   3.175   4.089 1.001  3300
      sigma.res           4.257   0.358   3.649   4.003   4.235   4.475   5.033 1.001  4900
      deviance          457.997   4.558 451.319 454.640 457.201 460.513 468.982 1.002  1700
      
      For each parameter, n.eff is a crude measure of effective sample size,
      and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
      
      DIC info (using the rule, pD = var(deviance)/2)
      pD = 10.4 and DIC = 468.4
      DIC is an estimate of expected predictive error (lower deviance is better).
      
    • Means parameterization
      massi = αj(i)*monthi*situationi + εi where ε ∼ N(0,σ2)
      View full code
      > modelString="
      model {
         #Likelihood
         for (i in 1:n) {
            mass[i]~dnorm(mean[i],tau.res)
            mean[i] <- alpha[month[i],situation[i]]
            y.err[i] <- mass[i] - mean[i]
         }
      
         #Priors and derivatives
         for (i in 1:nMonth) {
           for (j in 1:nSituation){
             alpha[i,j] ~ dnorm(0, 1.0E-6)
           }
         }
      
         sigma.res ~ dgamma(0.001, 0.001) #prior on sd residuals
         tau.res <- 1 / (sigma.res * sigma.res)
      
         #Effects
         beta0 <- alpha[1,1]
         for (i in 1:nMonth) {
           beta.month[i] <- alpha[i,1]-beta0
         }
         for (i in 1:nSituation) {
           beta.situation[i] <- alpha[1,i]-beta0
         }
      
         for (i in 1:nMonth) {
           for (j in 1:nSituation) {
      	   beta.int[i,j] <- alpha[i,j]-beta.month[i]-beta.situation[j]-beta0
      	 }
         }
         sd.Month <- sd(beta.month[])
         sd.Situation <- sd(beta.situation[])
         sd.Int <- sd(beta.int[,])
         sd.Resid <- sd(y.err[])
      
         #Pairwise comparisons of situations
         for (i in 1:nSituation) {
           for (j in 1:nSituation) {
      	   eff[i,j] <- beta.situation[i] - beta.situation[j]
      	 }
         }
         #Effects of month within each situation
         for (i in 1:nSituation) {
            eff2[i] <- alpha[1,i] - alpha[2,i]
         }
       }
      "
      > ## write the model to a text file (I suggest you alter the path to somewhere more relevant 
      > ## to your system!)  
      > writeLines(modelString,con="../downloads/BUGSscripts/ws9.6bQ1.1b.txt")
      > starling.list <- with(starling,
      +                          list(mass=MASS,
      +                          month=as.numeric(MONTH), situation=as.numeric(SITUATION),
      +                                  n=nrow(starling), nMonth=length(levels(MONTH)),
      +                                  nSituation=length(levels(SITUATION))
      +                         )
      +            )
      > 
      > params <- c("beta.month","beta.situation","beta.int","sigma.res","sd.Month","sd.Situation","sd.Int","sd.Resid", "eff","eff2")
      > adaptSteps = 1000
      > burnInSteps = 200
      > nChains = 3
      > numSavedSteps = 5000
      > thinSteps = 10
      > nIter = ceiling((numSavedSteps * thinSteps)/nChains)
      > 
      > library(R2jags)
      > ## 
      > starling.r2jags2 <- jags(data=starling.list,
      +           inits=NULL, #or inits=list(inits,inits,inits) # since there are three chains
      +           parameters.to.save=params,
      +           model.file="../downloads/BUGSscripts/ws9.6bQ1.1b.txt",
      +           n.chains=3,
      +           n.iter=nIter,
      +           n.burnin=burnInSteps,
      +       n.thin=thinSteps
      +           )
      
      Compiling model graph
         Resolving undeclared variables
         Allocating nodes
         Graph Size: 381
      
      Initializing model
      
      
      > print(starling.r2jags2)
      
      Inference for Bugs model at "../downloads/BUGSscripts/ws9.6bQ1.1b.txt", fit using jags,
       3 chains, each with 16667 iterations (first 200 discarded), n.thin = 10
       n.sims = 4941 iterations saved
                        mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
      beta.int[1,1]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,1]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[1,2]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,2]      -3.566   2.721  -8.826  -5.398  -3.568  -1.710   1.735 1.001  4900
      beta.int[1,3]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,3]      -2.402   2.666  -7.569  -4.189  -2.387  -0.634   2.796 1.002  1100
      beta.int[1,4]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,4]      -1.568   2.725  -6.805  -3.450  -1.564   0.321   3.712 1.002  1500
      beta.month[1]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.month[2]      -7.203   1.919 -10.918  -8.490  -7.208  -5.906  -3.403 1.002  2200
      beta.situation[1]   0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.situation[2]  -0.648   1.921  -4.359  -1.946  -0.637   0.616   3.119 1.001  4900
      beta.situation[3]  -2.604   1.892  -6.204  -3.883  -2.642  -1.322   1.198 1.001  4900
      beta.situation[4]  -6.639   1.913 -10.306  -7.943  -6.645  -5.333  -2.831 1.002  2300
      eff[1,1]            0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[2,1]           -0.648   1.921  -4.359  -1.946  -0.637   0.616   3.119 1.001  4900
      eff[3,1]           -2.604   1.892  -6.204  -3.883  -2.642  -1.322   1.198 1.001  4900
      eff[4,1]           -6.639   1.913 -10.306  -7.943  -6.645  -5.333  -2.831 1.002  2300
      eff[1,2]            0.648   1.921  -3.119  -0.616   0.637   1.946   4.359 1.001  4900
      eff[2,2]            0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[3,2]           -1.956   1.914  -5.707  -3.223  -1.957  -0.648   1.726 1.001  4900
      eff[4,2]           -5.991   1.894  -9.677  -7.263  -5.999  -4.747  -2.186 1.002  1800
      eff[1,3]            2.604   1.892  -1.198   1.322   2.642   3.883   6.204 1.001  4900
      eff[2,3]            1.956   1.914  -1.726   0.648   1.957   3.223   5.707 1.001  4900
      eff[3,3]            0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[4,3]           -4.035   1.904  -7.732  -5.325  -4.031  -2.746  -0.327 1.001  3900
      eff[1,4]            6.639   1.913   2.831   5.333   6.645   7.943  10.306 1.002  2300
      eff[2,4]            5.991   1.894   2.186   4.747   5.999   7.263   9.677 1.002  1800
      eff[3,4]            4.035   1.904   0.327   2.746   4.031   5.325   7.732 1.001  3900
      eff[4,4]            0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff2[1]             7.203   1.919   3.403   5.906   7.208   8.490  10.918 1.001  3300
      eff2[2]            10.769   1.898   7.117   9.512  10.798  12.029  14.455 1.001  4900
      eff2[3]             9.605   1.904   5.859   8.338   9.627  10.853  13.318 1.004  1200
      eff2[4]             8.772   1.916   5.019   7.474   8.779  10.084  12.475 1.001  4400
      sd.Int              1.779   0.787   0.508   1.195   1.698   2.287   3.499 1.001  4900
      sd.Month            3.602   0.960   1.701   2.953   3.604   4.245   5.459 1.001  3300
      sd.Resid            4.184   0.104   4.037   4.108   4.167   4.237   4.431 1.001  4900
      sd.Situation        2.766   0.644   1.503   2.339   2.768   3.199   4.013 1.002  1800
      sigma.res           4.253   0.358   3.626   4.001   4.229   4.476   5.011 1.002  1400
      deviance          458.035   4.463 451.313 454.807 457.362 460.524 468.601 1.001  4900
      
      For each parameter, n.eff is a crude measure of effective sample size,
      and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
      
      DIC info (using the rule, pD = var(deviance)/2)
      pD = 10.0 and DIC = 468.0
      DIC is an estimate of expected predictive error (lower deviance is better).
      
    • Before fully exploring the parameters, it is prudent to examine the convergence and mixing diagnostics. Chose either the effects or means parameterization model (they should yield much the same).
      • Trace plots
        View code
        > plot(as.mcmc(starling.r2jags))
        
      • Raftery diagnostic
        View code
        > raftery.diag(as.mcmc(starling.r2jags))
        
        [[1]]
        
        Quantile (q) = 0.025
        Accuracy (r) = +/- 0.005
        Probability (s) = 0.95 
        
        You need a sample size of at least 3746 with these values of q, r and s
        
        [[2]]
        
        Quantile (q) = 0.025
        Accuracy (r) = +/- 0.005
        Probability (s) = 0.95 
        
        You need a sample size of at least 3746 with these values of q, r and s
        
        [[3]]
        
        Quantile (q) = 0.025
        Accuracy (r) = +/- 0.005
        Probability (s) = 0.95 
        
        You need a sample size of at least 3746 with these values of q, r and s
        
        
      • Autocorrelation
        View code
        > autocorr.diag(as.mcmc(starling.r2jags))
        
                Mean[1,1]  Mean[2,1] Mean[1,2] Mean[2,2] Mean[1,3] Mean[2,3]  Mean[1,4] Mean[2,4]
        Lag 0    1.000000  1.0000000  1.000000  1.000000 1.0000000  1.000000  1.0000000  1.000000
        Lag 10  -0.001908  0.0039394  0.007642 -0.003132 0.0008062 -0.009733  0.0015500  0.004218
        Lag 50  -0.010543  0.0007551 -0.023796  0.007433 0.0175925 -0.017918 -0.0002732 -0.018483
        Lag 100 -0.017432  0.0037170  0.011733 -0.020159 0.0175249 -0.019843  0.0045985 -0.000722
        Lag 500 -0.027325 -0.0155207  0.002354 -0.007411 0.0041462  0.002988 -0.0086342  0.011634
                beta.int[1,1] beta.int[2,1] beta.int[1,2] beta.int[2,2] beta.int[1,3] beta.int[2,3]
        Lag 0             NaN           NaN           NaN      1.000000           NaN      1.000000
        Lag 10            NaN           NaN           NaN      0.004752           NaN      0.003448
        Lag 50            NaN           NaN           NaN     -0.002842           NaN     -0.024150
        Lag 100           NaN           NaN           NaN     -0.017860           NaN     -0.004386
        Lag 500           NaN           NaN           NaN     -0.009378           NaN     -0.015808
                beta.int[1,4] beta.int[2,4] beta.month[1] beta.month[2] beta.situation[1] beta.situation[2]
        Lag 0             NaN       1.00000           NaN     1.0000000               NaN          1.000000
        Lag 10            NaN       0.01550           NaN    -0.0045391               NaN          0.007261
        Lag 50            NaN       0.01528           NaN     0.0006787               NaN         -0.007635
        Lag 100           NaN      -0.02872           NaN    -0.0114574               NaN          0.005118
        Lag 500           NaN      -0.02084           NaN    -0.0224774               NaN         -0.003564
                beta.situation[3] beta.situation[4]     beta0   deviance eff[1,1]  eff[2,1]  eff[3,1]
        Lag 0            1.000000           1.00000  1.000000  1.0000000      NaN  1.000000  1.000000
        Lag 10           0.001950           0.00691 -0.001908 -0.0098722      NaN  0.007261  0.001950
        Lag 50          -0.003427           0.01735 -0.010543  0.0200729      NaN -0.007635 -0.003427
        Lag 100          0.004746          -0.01890 -0.017432  0.0006352      NaN  0.005118  0.004746
        Lag 500         -0.019668          -0.03489 -0.027325  0.0051727      NaN -0.003564 -0.019668
                eff[4,1]  eff[1,2] eff[2,2]  eff[3,2]  eff[4,2]  eff[1,3]  eff[2,3] eff[3,3] eff[4,3]
        Lag 0    1.00000  1.000000      NaN  1.000000  1.000000  1.000000  1.000000      NaN  1.00000
        Lag 10   0.00691  0.007261      NaN -0.008402 -0.004394  0.001950 -0.008402      NaN  0.01358
        Lag 50   0.01735 -0.007635      NaN  0.022240 -0.027287 -0.003427  0.022240      NaN  0.01526
        Lag 100 -0.01890  0.005118      NaN  0.008424  0.020021  0.004746  0.008424      NaN  0.02826
        Lag 500 -0.03489 -0.003564      NaN -0.019342  0.009797 -0.019668 -0.019342      NaN -0.02784
                eff[1,4]  eff[2,4] eff[3,4] eff[4,4]    eff2[1]   eff2[2]   eff2[3]   eff2[4]    sd.Int
        Lag 0    1.00000  1.000000  1.00000      NaN  1.0000000  1.000000  1.000000  1.000000  1.000000
        Lag 10   0.00691 -0.004394  0.01358      NaN -0.0045391  0.007980  0.008606  0.007452  0.007769
        Lag 50   0.01735 -0.027287  0.01526      NaN  0.0006787 -0.001802 -0.006377 -0.016473 -0.006266
        Lag 100 -0.01890  0.020021  0.02826      NaN -0.0114574 -0.006551  0.008250 -0.003164 -0.025175
        Lag 500 -0.03489  0.009797 -0.02784      NaN -0.0224774 -0.015014  0.015986  0.002236  0.005423
                  sd.Month  sd.Resid sd.Situation sigma.res
        Lag 0    1.0000000  1.000000    1.0000000  1.000000
        Lag 10  -0.0045391 -0.001852   -0.0066975 -0.002443
        Lag 50   0.0006787 -0.002548    0.0002704 -0.004719
        Lag 100 -0.0114574  0.002006   -0.0066467  0.004449
        Lag 500 -0.0224774  0.017954   -0.0203869 -0.006040
        
  2. A summary figure would be nice
    View code
    > starling.mcmc <- as.mcmc(starling.r2jags$BUGSoutput$sims.matrix)
    > #get the column indexes for the parameters we are interested in
    > starling.sum <- summary.mcmc(starling.mcmc, terms=c("Mean")) #or "alpha"
    > starling.sum <- cbind(starling.sum, expand.grid(Month=c("November","January"), Situation=paste("S",1:4,sep="")))
    > 
    > 
    > #starling.sum$name <- factor(rownames(starling.sum), levels=rownames(starling.sum))
    > 
    > p1<-ggplot(starling.sum,aes(y=mean, x=Situation, group=Month))+
    + geom_errorbar(aes(ymin=lower, ymax=upper),width=0,size=1.5)+
    + geom_errorbar(aes(ymin=lower.1, ymax=upper.1),width=0)+
    + geom_point(aes(shape=Month, fill=Month),size=3)+
    + scale_y_continuous("Starling mass (g)")+
    + scale_x_discrete("Roosting situation")+
    + scale_shape_manual(name="Month",values=c(21,16), breaks=c("November","January"), labels=c("November","January"))+
    + scale_fill_manual(name="Month",values=c("white","black"), breaks=c("November","January"), labels=c("November","January"))+
    + murray_opts+opts(legend.position=c(1,0.99),legend.justification=c(1,1))
    > 
    > 
    > #Finite-population standard deviations
    > starling.sd <- summary.mcmc(starling.mcmc, terms=c("sd"))
    > rownames(starling.sd) <- c("Month:Situation","Month","Residuals", "Situation")
    > starling.sd$name <- factor(rownames(starling.sd), levels=c("Residuals", "Month:Situation", "Situation","Month"))
    > 
    > p2<-ggplot(starling.sd,aes(y=name, x=median))+
    + geom_vline(xintercept=0,linetype="dashed")+
    + geom_hline(xintercept=0)+
    + scale_x_continuous("Finite population \nvariance components (sd)")+
    + geom_errorbarh(aes(xmin=lower.1, xmax=upper.1), height=0, size=1)+
    + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0, size=1.5)+
    + geom_point(size=3)+
    + geom_text(aes(label=sprintf("(%4.1f%%)",Perc),vjust=-1))+
    + murray_opts+
    + opts(axis.line=theme_blank(),axis.title.y=theme_blank(),
    + axis.text.y=theme_text(size=12,hjust=1))
    > 
    > grid.newpage()
    > pushViewport(viewport(layout=grid.layout(1,2,5,5,respect=TRUE)))
    > pushViewport(viewport(layout.pos.col=1, layout.pos.row=1))
    > print(p1,newpage=FALSE)
    > popViewport(1)
    > pushViewport(viewport(layout.pos.col=2, layout.pos.row=1))
    > print(p2,newpage=FALSE)
    > popViewport(0)
    
    > 
    + 
    + #eff and eff2
    > starling.eff <- summary.mcmc(starling.mcmc, terms=c("eff"))
    > rows <- c(matrix(upper.tri(diag(paste("S",1:4,sep=""))),ncol=1)[,1],rep(TRUE,4))
    > nms <- paste("Situation",1:4)
    > nms <- with(expand.grid(S1=nms, S2=nms),paste(S1,S2,sep=" vs "))
    > nms <- c(nms,"Nov - Jan (S1)","Nov - Jan (S2)","Nov - Jan (S3)","Nov - Jan (S4)")
    > starling.eff$name <- nms
    > starling.eff <- starling.eff[rows,]
    > 
    > p3 <- ggplot(starling.eff, aes(y=name, x=mean))+
    + geom_vline(xintercept=0,linetype="dashed")+
    + geom_hline(xintercept=0)+
    + geom_errorbarh(aes(xmin=lower.1, xmax=upper.1), height=0, size=1)+
    + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0, size=1.5)+
    + geom_point(size=3)+murray_opts+
    + scale_x_continuous(name="Effect size")+
    + opts(axis.line=theme_blank(),axis.title.y=theme_blank(),
    + axis.text.y=theme_text(size=12,hjust=1))
    > p3
    
    > 
    + 
    + 
    + #cohen's D
    > eff<-NULL
    > eff.sd <- NULL
    > for (i in seq(1,7,b=2)) {
    +   eff <-cbind(eff,starling.mcmc[,i] - starling.mcmc[,i+1])
    +   eff.sd <- cbind(eff.sd,mean(apply(starling.mcmc[,i:(i+1)],2,sd)))
    + }
    > cohen <- sweep(eff,2,eff.sd,"/")
    > adply(cohen,2,function(x) {
    +    data.frame(mean=mean(x), median=median(x),
    +                  HPDinterval(as.mcmc(x), p=0.68), HPDinterval(as.mcmc(x)))
    + })
    
    Error: unable to find an inherited method for function "HPDinterval", for signature "mcmc"

Unbalanced Two-factor ANOVA

Here is a modified example from Quinn and Keough (2002). Stehman and Meredith (1995) present data from an experiment that was set up to test the hypothesis that healthy spruce seedlings break bud sooner than diseased spruce seedlings. There were 2 factors: pH (3 levels: 3, 5.5, 7) and HEALTH (2 levels: healthy, diseased). The dependent variable was the average (from 5 buds) bud emergence rating (BRATING) on each seedling. The sample size varied for each combination of pH and health, ranging from 7 to 23 seedlings. With two factors, this experiment should be analyzed with a 2 factor (2 x 3) ANOVA.

Download Stehman data set
Format of stehman.csv data files
PHHEALTHGROUPBRATING
3DD30.0
........
3HH30.8
........
5.5DD5.50.0
........
5.5HH5.50.0
........
7DD70.2
........

PHCategorical listing of pH (not however that the levels are numbers and thus by default the variable is treated as a numeric variable rather than a factor - we need to correct for this)
HEALTHCategorical listing of the health status of the seedlings, D = diseased, H = healthy
GROUPCategorical listing of pH/health combinations - used for checking ANOVA assumptions
BRATINGAverage bud emergence rating per seedling
Starlings

Open the stehman data file.

Show code
> stehman <- read.table("../downloads/data/stehman.csv", header = T,
+     sep = ",", strip.white = T)
> head(stehman)
  PH HEALTH GROUP BRATING
1  3      D    D3     0.0
2  3      D    D3     0.8
3  3      D    D3     0.8
4  3      D    D3     0.8
5  3      D    D3     0.8
6  3      D    D3     0.8

The variable PH contains a list of pH values and is supposed to represent a factorial variable. However, because the contents of this variable are numbers, R initially treats them as numbers, and therefore considers the variable to be numeric rather than categorical. In order to force R/JAGS to treat this variable as a factor (categorical) it must be in categorical form (yet as numbers). Confused? What this means, is that to be treated as a factor, its levels must be indices and therefore the PH variable needs to be converted into a factor before being converted to a numeric. That way the levels 3, 5.5 and 7 will be coded as 1, 2 and 3.

Exploratory data analysis did not reveal any issues with normality or homogeneity of variance.

  1. Fit the JAGS model including all pairwise comparisons AND (just because we can), compare the bare surfaces to the surfaces with an algal biofilm. Use with either:
    • Effects parameterization
      bratingi = β0 + β1j(i)*phi + β2j(i)*healthi + β3j(i)*phi*healthi + εi where ε ∼ N(0,σ2)
      View full code
      > modelString="
      model {
         #Likelihood
         for (i in 1:n) {
            brating[i]~dnorm(mean[i],tau.res)
            mean[i] <- beta0+beta.ph[ph[i]]+beta.health[health[i]]+
      	             beta.int[ph[i],health[i]]
            y.err[i] <- brating[i] - mean[i]
         }
      
         #Priors and derivatives
         beta0 ~ dnorm(0,1.0E-6)
         beta.ph[1] <- 0
         beta.int[1,1] <- 0
         for (i in 2:nPh) {
           beta.ph[i] ~ dnorm(0.001, 1.0E-6) #prior
           beta.int[i,1] <- 0
         }
         beta.health[1] <- 0
         for (i in 2:nHealth){
           beta.health[i] ~ dnorm(0.001, 1.0E-6) #prior cannot be 0, or else will have an 
           ##init of zero and then will be dividing by zero when calculating cohen's D
           beta.int[1,i] <- 0
         }
         for (i in 2:nPh) {
           for (j in 2:nHealth){
             beta.int[i,j] ~ dnorm(0, 1.0E-6)
           }
         }
      
         #Means
         for (i in 1:nPh) {
           for (j in 1:nHealth) {
             Mean[i,j] <- beta0 + beta.ph[i] + beta.health[j] + beta.int[i,j]
           }
         } 
      
         sigma.res ~ dgamma(0.001, 0.001) #prior on sd residuals
         tau.res <- 1 / (sigma.res * sigma.res)
         sd.Ph <- sd(beta.ph[])
         sd.Health <- sd(beta.health[])
         sd.Int <- sd(beta.int[,])
         sd.Resid <- sd(y.err[])
      
         #Pairwise comparisons of PH
         for (i in 1:nPh) {
           for (j in 1:nPh) {
      	   eff[i,j] <- (beta.ph[i] - beta.ph[j])
             eff2[i,j] <- ((Mean[i,1]+Mean[i,2])/2) - ((Mean[j,1]+Mean[j,2])/2)
      	 }
         }
      #   #Effects of ph within each health
      #   for (i in 1:nHealth) {
      #      eff2[i] <- Mean[1,i] - Mean[2,i]
      #   }
       }
      "
      > ## write the model to a text file (I suggest you alter the path to somewhere more relevant 
      > ## to your system!)  
      > writeLines(modelString,con="../downloads/BUGSscripts/ws9.6bQ2.1a.txt")
      > stehman.list <- with(stehman,
      +                          list(brating=BRATING,
      +                          ph=as.numeric(as.factor(PH)), health=as.numeric(HEALTH),
      +                                  n=nrow(stehman), nPh=length(levels(as.factor(PH))),
      +                                  nHealth=length(levels(HEALTH))
      +                         )
      +            )
      > 
      > params <- c("beta0","beta.ph","beta.health","beta.int","sigma.res","sd.Ph","sd.Health","sd.Int","sd.Resid", "Mean","eff","eff2")
      > #params <- c("beta0","beta.ph","beta.health","beta.int","sigma.res","Mean", "sd.Ph","sd.Health","sd.Int","sd.Resid", "eff","eff2")
      > 
      > adaptSteps = 1000
      > burnInSteps = 200
      > nChains = 3
      > numSavedSteps = 5000
      > thinSteps = 10
      > nIter = ceiling((numSavedSteps * thinSteps)/nChains)
      > 
      > library(R2jags)
      > ## 
      > stehman.r2jags <- jags(data=stehman.list,
      +           inits=NULL, #or inits=list(inits,inits,inits) # since there are three chains
      +           parameters.to.save=params,
      +           model.file="../downloads/BUGSscripts/ws9.6bQ2.1a.txt",
      +           n.chains=3,
      +           n.iter=nIter,
      +           n.burnin=burnInSteps,
      +       n.thin=thinSteps
      +           )
      
      Compiling model graph
         Resolving undeclared variables
         Allocating nodes
         Graph Size: 385
      
      Initializing model
      
      
      > print(stehman.r2jags)
      
      Inference for Bugs model at "../downloads/BUGSscripts/ws9.6bQ2.1a.txt", fit using jags,
       3 chains, each with 16667 iterations (first 200 discarded), n.thin = 10
       n.sims = 4941 iterations saved
                     mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
      Mean[1,1]        1.193   0.107   0.979   1.120   1.194   1.266   1.400 1.001  2700
      Mean[2,1]        0.807   0.107   0.592   0.735   0.807   0.879   1.016 1.001  3400
      Mean[3,1]        1.125   0.112   0.905   1.050   1.124   1.199   1.348 1.001  4900
      Mean[1,2]        1.618   0.156   1.311   1.516   1.616   1.725   1.920 1.001  4900
      Mean[2,2]        1.229   0.192   0.849   1.099   1.228   1.359   1.597 1.002  2100
      Mean[3,2]        1.337   0.163   1.012   1.228   1.337   1.447   1.656 1.002  2200
      beta.health[1]   0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.health[2]   0.425   0.189   0.059   0.299   0.425   0.556   0.796 1.002  2400
      beta.int[1,1]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,1]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[3,1]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[1,2]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,2]   -0.003   0.292  -0.561  -0.201  -0.001   0.195   0.573 1.001  2600
      beta.int[3,2]   -0.214   0.273  -0.738  -0.399  -0.211  -0.028   0.333 1.002  1600
      beta.ph[1]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.ph[2]      -0.386   0.152  -0.683  -0.488  -0.387  -0.288  -0.081 1.001  4900
      beta.ph[3]      -0.068   0.156  -0.374  -0.172  -0.069   0.034   0.246 1.001  4700
      beta0            1.193   0.107   0.979   1.120   1.194   1.266   1.400 1.001  2700
      eff[1,1]         0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[2,1]        -0.386   0.152  -0.683  -0.488  -0.387  -0.288  -0.081 1.001  4900
      eff[3,1]        -0.068   0.156  -0.374  -0.172  -0.069   0.034   0.246 1.001  4700
      eff[1,2]         0.386   0.152   0.081   0.288   0.387   0.488   0.683 1.001  4900
      eff[2,2]         0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[3,2]         0.318   0.153   0.013   0.215   0.320   0.423   0.621 1.001  4900
      eff[1,3]         0.068   0.156  -0.246  -0.034   0.069   0.172   0.374 1.001  4700
      eff[2,3]        -0.318   0.153  -0.621  -0.423  -0.320  -0.215  -0.013 1.001  4900
      eff[3,3]         0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff2[1,1]        0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff2[2,1]       -0.388   0.147  -0.676  -0.486  -0.389  -0.289  -0.099 1.001  4900
      eff2[3,1]       -0.175   0.135  -0.442  -0.265  -0.175  -0.083   0.091 1.001  4900
      eff2[1,2]        0.388   0.147   0.099   0.289   0.389   0.486   0.676 1.001  4900
      eff2[2,2]        0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff2[3,2]        0.213   0.149  -0.072   0.114   0.213   0.313   0.504 1.001  4900
      eff2[1,3]        0.175   0.135  -0.091   0.083   0.175   0.265   0.442 1.001  4900
      eff2[2,3]       -0.213   0.149  -0.504  -0.313  -0.213  -0.114   0.072 1.001  4900
      eff2[3,3]        0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      sd.Health        0.214   0.093   0.036   0.150   0.212   0.278   0.398 1.004  1500
      sd.Int           0.145   0.074   0.027   0.090   0.137   0.193   0.310 1.001  4900
      sd.Ph            0.181   0.059   0.064   0.141   0.181   0.220   0.297 1.001  4900
      sd.Resid         0.504   0.009   0.492   0.498   0.502   0.509   0.525 1.001  4900
      sigma.res        0.512   0.039   0.442   0.484   0.509   0.537   0.595 1.001  4900
      deviance       141.478   3.856 135.933 138.589 140.846 143.614 150.620 1.001  4900
      
      For each parameter, n.eff is a crude measure of effective sample size,
      and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
      
      DIC info (using the rule, pD = var(deviance)/2)
      pD = 7.4 and DIC = 148.9
      DIC is an estimate of expected predictive error (lower deviance is better).
      
    • Means parameterization
      bratingi = αj(i)*phi*healthi + εi where ε ∼ N(0,σ2)
      View full code
      > modelString="
      model {
         #Likelihood
         for (i in 1:n) {
            brating[i]~dnorm(mean[i],tau.res)
            mean[i] <- alpha[ph[i],health[i]]
            y.err[i] <- brating[i] - mean[i]
         }
      
         #Priors and derivatives
         for (i in 1:nPh) {
           for (j in 1:nHealth){
             alpha[i,j] ~ dnorm(0, 1.0E-6)
           }
         }
      
         sigma.res ~ dgamma(0.001, 0.001) #prior on sd residuals
         tau.res <- 1 / (sigma.res * sigma.res)
      
         #Effects
         beta0 <- alpha[1,1]
         for (i in 1:nPh) {
           beta.ph[i] <- alpha[i,1]-beta0
         }
         for (i in 1:nHealth) {
           beta.health[i] <- alpha[1,i]-beta0
         }
      
         for (i in 1:nPh) {
           for (j in 1:nHealth) {
      	   beta.int[i,j] <- alpha[i,j]-beta.ph[i]-beta.health[j]-beta0
      	 }
         }
         sd.Ph <- sd(beta.ph[])
         sd.Health <- sd(beta.health[])
         sd.Int <- sd(beta.int[,])
         sd.Resid <- sd(y.err[])
      
         #Pairwise comparisons of PH
         for (i in 1:nPh) {
           for (j in 1:nPh) {
      	   eff[i,j] <- (beta.ph[i] - beta.ph[j])
             eff2[i,j] <- ((alpha[i,1]+alpha[i,2])/2) - ((alpha[j,1]+alpha[j,2])/2)
      	 }
         }
       }
      "
      > ## write the model to a text file (I suggest you alter the path to somewhere more relevant 
      > ## to your system!)  
      > writeLines(modelString,con="../downloads/BUGSscripts/ws9.6bQ2.1b.txt")
      > stehman.list <- with(stehman,
      +                          list(brating=BRATING,
      +                          ph=as.numeric(as.factor(PH)), health=as.numeric(HEALTH),
      +                                  n=nrow(stehman), nPh=length(levels(as.factor(PH))),
      +                                  nHealth=length(levels(HEALTH))
      +                         )
      +            )
      > 
      > params <- c("beta.ph","beta.health","beta.int","sigma.res","sd.Ph","sd.Health","sd.Int","sd.Resid","eff","eff2")
      > adaptSteps = 1000
      > burnInSteps = 200
      > nChains = 3
      > numSavedSteps = 5000
      > thinSteps = 10
      > nIter = ceiling((numSavedSteps * thinSteps)/nChains)
      > 
      > library(R2jags)
      > ## 
      > stehman.r2jags2 <- jags(data=stehman.list,
      +           inits=NULL, #or inits=list(inits,inits,inits) # since there are three chains
      +           parameters.to.save=params,
      +           model.file="../downloads/BUGSscripts/ws9.6bQ2.1b.txt",
      +           n.chains=3,
      +           n.iter=nIter,
      +           n.burnin=burnInSteps,
      +       n.thin=thinSteps
      +           )
      
      Compiling model graph
         Resolving undeclared variables
         Allocating nodes
         Graph Size: 395
      
      Initializing model
      
      
      > print(stehman.r2jags2)
      
      Inference for Bugs model at "../downloads/BUGSscripts/ws9.6bQ2.1b.txt", fit using jags,
       3 chains, each with 16667 iterations (first 200 discarded), n.thin = 10
       n.sims = 4941 iterations saved
                     mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
      beta.health[1]   0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.health[2]   0.430   0.186   0.070   0.305   0.427   0.556   0.791 1.001  4900
      beta.int[1,1]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,1]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[3,1]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[1,2]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,2]   -0.008   0.288  -0.581  -0.195  -0.010   0.183   0.556 1.001  4900
      beta.int[3,2]   -0.215   0.267  -0.739  -0.390  -0.219  -0.034   0.294 1.001  4900
      beta.ph[1]       0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.ph[2]      -0.382   0.151  -0.677  -0.483  -0.382  -0.280  -0.083 1.001  4900
      beta.ph[3]      -0.067   0.154  -0.368  -0.169  -0.066   0.038   0.233 1.001  4900
      eff[1,1]         0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[2,1]        -0.382   0.151  -0.677  -0.483  -0.382  -0.280  -0.083 1.001  4900
      eff[3,1]        -0.067   0.154  -0.368  -0.169  -0.066   0.038   0.233 1.001  4900
      eff[1,2]         0.382   0.151   0.083   0.280   0.382   0.483   0.677 1.001  4900
      eff[2,2]         0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff[3,2]         0.315   0.153   0.017   0.211   0.316   0.417   0.617 1.001  4900
      eff[1,3]         0.067   0.154  -0.233  -0.038   0.066   0.169   0.368 1.001  4900
      eff[2,3]        -0.315   0.153  -0.617  -0.417  -0.316  -0.211  -0.017 1.001  4900
      eff[3,3]         0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff2[1,1]        0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff2[2,1]       -0.386   0.146  -0.674  -0.485  -0.384  -0.286  -0.101 1.001  4900
      eff2[3,1]       -0.174   0.136  -0.442  -0.266  -0.173  -0.083   0.087 1.001  3700
      eff2[1,2]        0.386   0.146   0.101   0.286   0.384   0.485   0.674 1.001  4900
      eff2[2,2]        0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      eff2[3,2]        0.212   0.147  -0.072   0.110   0.213   0.309   0.498 1.001  2600
      eff2[1,3]        0.174   0.136  -0.087   0.083   0.173   0.266   0.442 1.001  3700
      eff2[2,3]       -0.212   0.147  -0.498  -0.309  -0.213  -0.110   0.072 1.001  2600
      eff2[3,3]        0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      sd.Health        0.216   0.091   0.039   0.153   0.213   0.278   0.396 1.002  4700
      sd.Int           0.143   0.074   0.027   0.089   0.136   0.189   0.307 1.001  3200
      sd.Ph            0.179   0.059   0.067   0.138   0.179   0.218   0.296 1.001  4900
      sd.Resid         0.504   0.009   0.493   0.498   0.502   0.508   0.526 1.001  4900
      sigma.res        0.511   0.038   0.442   0.484   0.509   0.536   0.593 1.001  4900
      deviance       141.377   3.820 135.963 138.571 140.734 143.466 150.632 1.001  4900
      
      For each parameter, n.eff is a crude measure of effective sample size,
      and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
      
      DIC info (using the rule, pD = var(deviance)/2)
      pD = 7.3 and DIC = 148.7
      DIC is an estimate of expected predictive error (lower deviance is better).
      
  2. Before fully exploring the parameters, it is prudent to examine the convergence and mixing diagnostics. Chose either the effects or means parameterization model (they should yield much the same).
    • Trace plots
      View code
      > plot(as.mcmc(stehman.r2jags))
      
    • Raftery diagnostic
      View code
      > raftery.diag(as.mcmc(stehman.r2jags))
      
      [[1]]
      
      Quantile (q) = 0.025
      Accuracy (r) = +/- 0.005
      Probability (s) = 0.95 
      
      You need a sample size of at least 3746 with these values of q, r and s
      
      [[2]]
      
      Quantile (q) = 0.025
      Accuracy (r) = +/- 0.005
      Probability (s) = 0.95 
      
      You need a sample size of at least 3746 with these values of q, r and s
      
      [[3]]
      
      Quantile (q) = 0.025
      Accuracy (r) = +/- 0.005
      Probability (s) = 0.95 
      
      You need a sample size of at least 3746 with these values of q, r and s
      
      
    • Autocorrelation
      View code
      > autocorr.diag(as.mcmc(stehman.r2jags))
      
              Mean[1,1] Mean[2,1] Mean[3,1] Mean[1,2]  Mean[2,2] Mean[3,2] beta.health[1] beta.health[2]
      Lag 0    1.000000   1.00000  1.000000  1.000000  1.0000000  1.000000            NaN       1.000000
      Lag 10   0.015825  -0.02483  0.029052  0.001061 -0.0204654 -0.025140            NaN      -0.007364
      Lag 50  -0.014374  -0.02767  0.003399  0.035938 -0.0004691 -0.003156            NaN      -0.002253
      Lag 100  0.008493  -0.01021 -0.007150  0.018119 -0.0403369 -0.006234            NaN       0.003267
      Lag 500  0.016258   0.02617 -0.019063  0.022172 -0.0003144  0.010743            NaN       0.016014
              beta.int[1,1] beta.int[2,1] beta.int[3,1] beta.int[1,2] beta.int[2,2] beta.int[3,2]
      Lag 0             NaN           NaN           NaN           NaN      1.000000      1.000000
      Lag 10            NaN           NaN           NaN           NaN      0.006026     -0.002851
      Lag 50            NaN           NaN           NaN           NaN      0.007103      0.000095
      Lag 100           NaN           NaN           NaN           NaN     -0.006275     -0.012467
      Lag 500           NaN           NaN           NaN           NaN      0.007731      0.003171
              beta.ph[1] beta.ph[2] beta.ph[3]     beta0  deviance eff[1,1]  eff[2,1]  eff[3,1]  eff[1,2]
      Lag 0          NaN   1.000000   1.000000  1.000000  1.000000      NaN  1.000000  1.000000  1.000000
      Lag 10         NaN   0.005331   0.014305  0.015825 -0.006270      NaN  0.005331  0.014305  0.005331
      Lag 50         NaN  -0.041872  -0.023950 -0.014374 -0.005720      NaN -0.041872 -0.023950 -0.041872
      Lag 100        NaN  -0.004293   0.001759  0.008493 -0.001357      NaN -0.004293  0.001759 -0.004293
      Lag 500        NaN   0.032709   0.003853  0.016258 -0.025478      NaN  0.032709  0.003853  0.032709
              eff[2,2]  eff[3,2]  eff[1,3]  eff[2,3] eff[3,3] eff2[1,1] eff2[2,1] eff2[3,1] eff2[1,2]
      Lag 0        NaN  1.000000  1.000000  1.000000      NaN       NaN   1.00000  1.000000   1.00000
      Lag 10       NaN -0.001814  0.014305 -0.001814      NaN       NaN  -0.02276  0.008053  -0.02276
      Lag 50       NaN -0.022216 -0.023950 -0.022216      NaN       NaN   0.02488  0.008395   0.02488
      Lag 100      NaN -0.014981  0.001759 -0.014981      NaN       NaN  -0.02791  0.009600  -0.02791
      Lag 500      NaN  0.012182  0.003853  0.012182      NaN       NaN   0.01505  0.031804   0.01505
              eff2[2,2] eff2[3,2] eff2[1,3] eff2[2,3] eff2[3,3] sd.Health    sd.Int      sd.Ph  sd.Resid
      Lag 0         NaN   1.00000  1.000000   1.00000       NaN  1.000000  1.000000  1.000e+00  1.000000
      Lag 10        NaN  -0.03517  0.008053  -0.03517       NaN -0.007974 -0.007266  4.519e-05 -0.020409
      Lag 50        NaN  -0.01782  0.008395  -0.01782       NaN  0.001795  0.006918 -2.808e-02 -0.013002
      Lag 100       NaN  -0.02574  0.009600  -0.02574       NaN  0.001555 -0.025238 -1.349e-02  0.005629
      Lag 500       NaN   0.01089  0.031804   0.01089       NaN  0.018169 -0.003348  3.505e-02 -0.042728
              sigma.res
      Lag 0   1.0000000
      Lag 10  0.0120803
      Lag 50  0.0134666
      Lag 100 0.0007495
      Lag 500 0.0017084
      
  3. A summary figure would be nice
    View code
    > stehman.mcmc <- as.mcmc(stehman.r2jags$BUGSoutput$sims.matrix)
    > #get the column indexes for the parameters we are interested in
    > stehman.sum <- summary.mcmc(stehman.mcmc, terms=c("Mean")) #or "alpha"
    > stehman.sum <- cbind(stehman.sum, expand.grid(Ph=c("3","5.5","7"), Health=c("Diseased","Healthy")))
    > 
    > 
    > p1<-ggplot(stehman.sum,aes(y=mean, x=Ph, group=Health))+
    + geom_errorbar(aes(ymin=lower, ymax=upper),width=0,size=1.5,position=position_dodge(0.2))+
    + geom_errorbar(aes(ymin=lower.1, ymax=upper.1),width=0,position=position_dodge(0.2))+
    + geom_point(aes(shape=Health, fill=Health),size=3,position=position_dodge(0.2))+
    + scale_y_continuous("Bud emergence rating")+
    + scale_x_discrete("Ph")+
    + scale_shape_manual(name="Health",values=c(21,16), breaks=c("Diseased","Healthy"), labels=c("Diseased","Healthy"))+
    + scale_fill_manual(name="Health",values=c("white","black"), breaks=c("Diseased","Healthy"), labels=c("Diseased","Healthy"))+
    + murray_opts+opts(legend.position=c(1,0.99),legend.justification=c(1,1))
    > 
    > 
    > #Finite-population standard deviations
    > stehman.sd <- summary.mcmc(stehman.mcmc, terms=c("sd"))
    > rownames(stehman.sd) <- c("Health","Ph:Health","Ph", "Residuals")
    > stehman.sd$name <- factor(c("Health","Ph:Health","Ph", "Residuals"), levels=c("Residuals", "Ph:Health", "Health","Ph"))
    > 
    > p2<-ggplot(stehman.sd,aes(y=name, x=median))+
    + geom_vline(xintercept=0,linetype="dashed")+
    + geom_hline(xintercept=0)+
    + scale_x_continuous("Finite population \nvariance components (sd)")+
    + geom_errorbarh(aes(xmin=lower.1, xmax=upper.1), height=0, size=1)+
    + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0, size=1.5)+
    + geom_point(size=3)+
    + geom_text(aes(label=sprintf("(%4.1f%%)",Perc),vjust=-1))+
    + murray_opts+
    + opts(axis.line=theme_blank(),axis.title.y=theme_blank(),
    + axis.text.y=theme_text(size=12,hjust=1))
    > 
    > grid.newpage()
    > pushViewport(viewport(layout=grid.layout(1,2,5,5,respect=TRUE)))
    > pushViewport(viewport(layout.pos.col=1, layout.pos.row=1))
    > print(p1,newpage=FALSE)
    > popViewport(1)
    > pushViewport(viewport(layout.pos.col=2, layout.pos.row=1))
    > print(p2,newpage=FALSE)
    > popViewport(0)
    
    > 
    + 
    + #eff and eff2
    > stehman.eff <- summary.mcmc(stehman.mcmc, terms=c("eff\\["))
    > rows <- matrix(upper.tri(diag(c("3","5.5","7"))),ncol=1)[,1]
    > #rows <- c(matrix(upper.tri(diag(c("3","5.5","7"))),ncol=1)[,1],
    > #matrix(upper.tri(diag(c("3","5.5","7"))),ncol=1)[,1]
    > #)
    > nms <- c("3","5.5","7")
    > nms <- with(expand.grid(S1=nms, S2=nms),paste(S1,S2,sep=" vs "))
    > #nms <- c(nms,nms)
    > stehman.eff$name <- nms
    > stehman.eff <- stehman.eff[rows,]
    > 
    > p3 <- ggplot(stehman.eff, aes(y=name, x=mean))+
    + geom_vline(xintercept=0,linetype="dashed")+
    + geom_hline(xintercept=0)+
    + geom_errorbarh(aes(xmin=lower.1, xmax=upper.1), height=0, size=1)+
    + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0, size=1.5)+
    + geom_point(size=3)+murray_opts+
    + scale_x_continuous(name="Effect size")+
    + opts(axis.line=theme_blank(),axis.title.y=theme_blank(),
    + axis.text.y=theme_text(size=12,hjust=1))
    > p3
    
    > 
    + 
    + 
    + #cohen's D
    > eff<-NULL
    > eff.sd <- NULL
    > for (i in seq(1,7,b=2)) {
    +   eff <-cbind(eff,stehman.mcmc[,i] - stehman.mcmc[,i+1])
    +   eff.sd <- cbind(eff.sd,mean(apply(stehman.mcmc[,i:(i+1)],2,sd)))
    + }
    > cohen <- sweep(eff,2,eff.sd,"/")
    > adply(cohen,2,function(x) {
    +    data.frame(mean=mean(x), median=median(x),
    +                  HPDinterval(as.mcmc(x), p=0.68), HPDinterval(as.mcmc(x)))
    + })
    
    Error: unable to find an inherited method for function "HPDinterval", for signature "mcmc"

Two-factor ANOVA with substantial interactions

An ecologist studying a rocky shore at Phillip Island, in southeastern Australia, was interested in how clumps of intertidal mussels are maintained. In particular, he wanted to know how densities of adult mussels affected recruitment of young individuals from the plankton. As with most marine invertebrates, recruitment is highly patchy in time, so he expected to find seasonal variation, and the interaction between season and density - whether effects of adult mussel density vary across seasons - was the aspect of most interest.

The data were collected from four seasons, and with two densities of adult mussels. The experiment consisted of clumps of adult mussels attached to the rocks. These clumps were then brought back to the laboratory, and the number of baby mussels recorded. There were 3-6 replicate clumps for each density and season combination.

Download Quinn data set
Format of quinn.csv data files
SEASONDENSITYRECRUITSSQRTRECRUITSGROUP
SpringLow153.87SpringLow
..........
SpringHigh113.32SpringHigh
..........
SummerLow214.58SummerLow
..........
SummerHigh345.83SummerHigh
..........
AutumnLow143.74AutumnLow
..........
SEASONCategorical listing of Season in which mussel clumps were collected ­ independent variable
DENSITYCategorical listing of the density of mussels within mussel clump ­ independent variable
RECRUITSThe number of mussel recruits ­ response variable
SQRTRECRUITSSquare root transformation of RECRUITS - needed to meet the test assumptions
GROUPSCategorical listing of Season/Density combinations - used for checking ANOVA assumptions
Mussel

Open the quinn data file.

Show code
> quinn <- read.table("../downloads/data/quinn.csv", header = T, sep = ",",
+     strip.white = T)
> head(quinn)
  SEASON DENSITY RECRUITS SQRTRECRUITS      GROUP
1 Spring     Low       15        3.873  SpringLow
2 Spring     Low       10        3.162  SpringLow
3 Spring     Low       13        3.606  SpringLow
4 Spring     Low       13        3.606  SpringLow
5 Spring     Low        5        2.236  SpringLow
6 Spring    High       11        3.317 SpringHigh

Exploratory data analysis revealed non-normality and heterogeneous variances (probably due to non-normality). As the response is count data, a square-root transformation is one possible way of addressing these issues.

  1. Fit the JAGS model including all pairwise comparisons AND (just because we can), compare the bare surfaces to the surfaces with an algal biofilm. Use with either:
    • Effects parameterization
      recruitsi = β0 + β1j(i)*seasoni + β2j(i)*densityi + β3j(i)*seasoni*densityi + εi where ε ∼ N(0,σ2)
      View full code
      > modelString="
      data {
         srecruits <- sqrt(recruits)
      }
      model {
         #Likelihood
         for (i in 1:n) {
            srecruits[i]~dnorm(mean[i],tau.res)
            mean[i] <- beta0+beta.season[season[i]]+beta.density[density[i]]+
      	             beta.int[season[i],density[i]]
            y.err[i] <- srecruits[i] - mean[i]
         }
      
         #Priors and derivatives
         beta0 ~ dnorm(0,1.0E-6)
         beta.season[1] <- 0
         beta.int[1,1] <- 0
         for (i in 2:nSeason) {
           beta.season[i] ~ dnorm(0, 1.0E-6) #prior
           beta.int[i,1] <- 0
         }
         beta.density[1] <- 0
         for (i in 2:nDensity){
           beta.density[i] ~ dnorm(0, 1.0E-6) #prior cannot be 0, or else will have an 
           beta.int[1,i] <- 0
         }
         for (i in 2:nSeason) {
           for (j in 2:nDensity){
             beta.int[i,j] ~ dnorm(0, 1.0E-6)
           }
         }
      
         #Means
         for (i in 1:nSeason) {
           for (j in 1:nDensity) {
             Mean[i,j] <- beta0 + beta.season[i] + beta.density[j] + beta.int[i,j]
           }
         } 
      
         sigma.res ~ dgamma(0.001, 0.001) #prior on sd residuals
         tau.res <- 1 / (sigma.res * sigma.res)
         sd.Season <- sd(beta.season[])
         sd.Density <- sd(beta.density[])
         sd.Int <- sd(beta.int[,])
         sd.Resid <- sd(y.err[])
      
         #Pairwise comparisons of Density (in each Season)
         for (i in 1:nSeason) {
           eff[i] <- Mean[i,1] - Mean[i,2]
         }
         
       }
      "
      > ## write the model to a text file (I suggest you alter the path to somewhere more relevant 
      > ## to your system!)  
      > writeLines(modelString,con="../downloads/BUGSscripts/ws9.6bQ3.1a.txt")
      > quinn.list <- with(quinn,
      +                          list(recruits=RECRUITS,
      +                          season=as.numeric(SEASON), density=as.numeric(DENSITY),
      +                                  n=nrow(quinn), nSeason=length(levels(SEASON)),
      +                                  nDensity=length(levels(DENSITY))
      +                         )
      +            )
      > 
      > #params <- c("beta0","beta.season","beta.density","beta.int","sigma.res","sd.Season","sd.Density","sd.Int","sd.Resid", "Mean","eff","eff2")
      > params <- c("beta0","beta.season","beta.density","beta.int","sigma.res","Mean", "sd.Season","sd.Density","sd.Int","sd.Resid","eff")
      > 
      > adaptSteps = 1000
      > burnInSteps = 200
      > nChains = 3
      > numSavedSteps = 5000
      > thinSteps = 10
      > nIter = ceiling((numSavedSteps * thinSteps)/nChains)
      > 
      > library(R2jags)
      > ## 
      > quinn.r2jags <- jags(data=quinn.list,
      +           inits=NULL, #or inits=list(inits,inits,inits) # since there are three chains
      +           parameters.to.save=params,
      +           model.file="../downloads/BUGSscripts/ws9.6bQ3.1a.txt",
      +           n.chains=3,
      +           n.iter=nIter,
      +           n.burnin=burnInSteps,
      +       n.thin=thinSteps
      +           )
      
      Compiling data graph
         Resolving undeclared variables
         Allocating nodes
         Initializing
         Reading data back into data table
      Compiling model graph
         Resolving undeclared variables
         Allocating nodes
         Graph Size: 249
      
      Initializing model
      
      
      > print(quinn.r2jags)
      
      Inference for Bugs model at "../downloads/BUGSscripts/ws9.6bQ3.1a.txt", fit using jags,
       3 chains, each with 16667 iterations (first 200 discarded), n.thin = 10
       n.sims = 4941 iterations saved
                      mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
      Mean[1,1]         4.232   0.430   3.409   3.947   4.227   4.520   5.075 1.001  4900
      Mean[2,1]         3.023   0.423   2.182   2.752   3.034   3.296   3.857 1.001  4900
      Mean[3,1]         6.875   0.430   6.027   6.591   6.877   7.153   7.727 1.001  4900
      Mean[4,1]         2.272   0.427   1.437   1.987   2.267   2.559   3.114 1.002  1400
      Mean[1,2]         4.264   0.515   3.241   3.927   4.271   4.600   5.275 1.002  1100
      Mean[2,2]         3.309   0.465   2.384   2.994   3.308   3.616   4.255 1.001  4900
      Mean[3,2]         4.658   0.421   3.820   4.379   4.659   4.935   5.477 1.001  4900
      Mean[4,2]         0.925   0.613  -0.318   0.516   0.941   1.338   2.096 1.001  4900
      beta.density[1]   0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.density[2]   0.032   0.672  -1.260  -0.418   0.037   0.484   1.340 1.002  1700
      beta.int[1,1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[3,1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[4,1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[1,2]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,2]     0.254   0.914  -1.562  -0.348   0.250   0.863   2.062 1.002  2300
      beta.int[3,2]    -2.249   0.911  -4.056  -2.839  -2.253  -1.670  -0.413 1.001  4400
      beta.int[4,2]    -1.379   0.994  -3.298  -2.044  -1.380  -0.713   0.578 1.002  1200
      beta.season[1]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.season[2]   -1.209   0.603  -2.402  -1.620  -1.209  -0.814  -0.040 1.001  4900
      beta.season[3]    2.643   0.610   1.437   2.225   2.658   3.047   3.845 1.001  4900
      beta.season[4]   -1.960   0.604  -3.125  -2.372  -1.958  -1.550  -0.771 1.002  1300
      beta0             4.232   0.430   3.409   3.947   4.227   4.520   5.075 1.001  4900
      eff[1]           -0.032   0.672  -1.340  -0.484  -0.037   0.418   1.260 1.002  1700
      eff[2]           -0.286   0.620  -1.526  -0.695  -0.276   0.140   0.920 1.001  4900
      eff[3]            2.217   0.604   1.036   1.825   2.212   2.613   3.402 1.001  4900
      eff[4]            1.347   0.742  -0.134   0.848   1.351   1.845   2.821 1.002  1400
      sd.Density        0.269   0.202   0.012   0.111   0.224   0.381   0.758 1.002  4900
      sd.Int            0.935   0.270   0.454   0.745   0.922   1.107   1.504 1.001  4900
      sd.Resid          1.002   0.053   0.932   0.963   0.992   1.029   1.132 1.001  4900
      sd.Season         1.774   0.212   1.360   1.633   1.774   1.911   2.193 1.001  3200
      sigma.res         1.033   0.127   0.821   0.942   1.021   1.111   1.317 1.001  4900
      deviance        121.313   4.717 114.234 117.829 120.652 124.041 132.017 1.001  4900
      
      For each parameter, n.eff is a crude measure of effective sample size,
      and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
      
      DIC info (using the rule, pD = var(deviance)/2)
      pD = 11.1 and DIC = 132.4
      DIC is an estimate of expected predictive error (lower deviance is better).
      
    • Means parameterization
      srecruitsi = αj(i)*seasoni*densityi + εi where ε ∼ N(0,σ2)
      View full code
      > modelString="
      data {
         srecruits <- sqrt(recruits)
      }
      model {
         #Likelihood
         for (i in 1:n) {
            srecruits[i]~dnorm(mean[i],tau.res)
            mean[i] <- alpha[season[i],density[i]]
            y.err[i] <- srecruits[i] - mean[i]
         }
      
         #Priors and derivatives
         for (i in 1:nSeason) {
           for (j in 1:nDensity){
             alpha[i,j] ~ dnorm(0, 1.0E-6)
           }
         }
      
         sigma.res ~ dgamma(0.001, 0.001) #prior on sd residuals
         tau.res <- 1 / (sigma.res * sigma.res)
      
         #Effects
         beta0 <- alpha[1,1]
         for (i in 1:nSeason) {
           beta.season[i] <- alpha[i,1]-beta0
         }
         for (i in 1:nDensity) {
           beta.density[i] <- alpha[1,i]-beta0
         }
      
         for (i in 1:nSeason) {
           for (j in 1:nDensity) {
      	   beta.int[i,j] <- alpha[i,j]-beta.season[i]-beta.density[j]-beta0
      	 }
         }
         sd.Season <- sd(beta.season[])
         sd.Density <- sd(beta.density[])
         sd.Int <- sd(beta.int[,])
         sd.Resid <- sd(y.err[])
      
         #Pairwise comparisons of Density (in each Season)
         for (i in 1:nSeason) {
           eff[i] <- alpha[i,1] - alpha[i,2]
         }
      
       }
      "
      > ## write the model to a text file (I suggest you alter the path to somewhere more relevant 
      > ## to your system!)  
      > writeLines(modelString,con="../downloads/BUGSscripts/ws9.6bQ3.1b.txt")
      > quinn.list <- with(quinn,
      +                          list(recruits=RECRUITS,
      +                          season=as.numeric(SEASON), density=as.numeric(DENSITY),
      +                                  n=nrow(quinn), nSeason=length(levels(SEASON)),
      +                                  nDensity=length(levels(DENSITY))
      +                         )
      +            )
      > 
      > params <- c("alpha","beta0","beta.season","beta.density","beta.int","sigma.res","sd.Season","sd.Density","sd.Int","sd.Resid","eff")
      > adaptSteps = 1000
      > burnInSteps = 200
      > nChains = 3
      > numSavedSteps = 5000
      > thinSteps = 10
      > nIter = ceiling((numSavedSteps * thinSteps)/nChains)
      > 
      > library(R2jags)
      > ## 
      > quinn.r2jags2 <- jags(data=quinn.list,
      +           inits=NULL, #or inits=list(inits,inits,inits) # since there are three chains
      +           parameters.to.save=params,
      +           model.file="../downloads/BUGSscripts/ws9.6bQ3.1b.txt",
      +           n.chains=3,
      +           n.iter=nIter,
      +           n.burnin=burnInSteps,
      +       n.thin=thinSteps
      +           )
      
      Compiling data graph
         Resolving undeclared variables
         Allocating nodes
         Initializing
         Reading data back into data table
      Compiling model graph
         Resolving undeclared variables
         Allocating nodes
         Graph Size: 263
      
      Initializing model
      
      
      > print(quinn.r2jags2)
      
      Inference for Bugs model at "../downloads/BUGSscripts/ws9.6bQ3.1b.txt", fit using jags,
       3 chains, each with 16667 iterations (first 200 discarded), n.thin = 10
       n.sims = 4941 iterations saved
                      mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
      alpha[1,1]        4.225   0.426   3.392   3.940   4.228   4.510   5.070 1.001  4900
      alpha[2,1]        3.026   0.430   2.200   2.741   3.022   3.317   3.883 1.001  4900
      alpha[3,1]        6.859   0.424   6.018   6.581   6.857   7.140   7.703 1.002  2400
      alpha[4,1]        2.280   0.428   1.441   1.995   2.283   2.552   3.115 1.001  4700
      alpha[1,2]        4.246   0.525   3.212   3.909   4.248   4.600   5.264 1.002  1800
      alpha[2,2]        3.285   0.462   2.386   2.974   3.277   3.590   4.206 1.001  4900
      alpha[3,2]        4.648   0.422   3.811   4.372   4.649   4.921   5.502 1.001  4900
      alpha[4,2]        0.941   0.609  -0.267   0.535   0.944   1.346   2.122 1.001  4900
      beta.density[1]   0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.density[2]   0.021   0.674  -1.317  -0.412   0.024   0.465   1.337 1.001  3400
      beta.int[1,1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[3,1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[4,1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[1,2]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.int[2,2]     0.237   0.929  -1.582  -0.378   0.232   0.841   2.062 1.001  4900
      beta.int[3,2]    -2.232   0.895  -3.975  -2.815  -2.239  -1.666  -0.393 1.001  4900
      beta.int[4,2]    -1.361   0.995  -3.308  -2.014  -1.361  -0.698   0.564 1.002  2100
      beta.season[1]    0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
      beta.season[2]   -1.198   0.603  -2.367  -1.595  -1.208  -0.799  -0.003 1.001  4900
      beta.season[3]    2.635   0.599   1.440   2.261   2.630   3.021   3.816 1.002  3600
      beta.season[4]   -1.945   0.600  -3.130  -2.353  -1.949  -1.553  -0.787 1.001  4900
      beta0             4.225   0.426   3.392   3.940   4.228   4.510   5.070 1.001  4900
      eff[1]           -0.021   0.674  -1.337  -0.465  -0.024   0.412   1.317 1.001  3400
      eff[2]           -0.259   0.627  -1.481  -0.679  -0.263   0.156   0.972 1.001  4900
      eff[3]            2.211   0.590   1.025   1.820   2.210   2.608   3.370 1.004  1200
      eff[4]            1.339   0.735  -0.097   0.862   1.338   1.811   2.785 1.001  4300
      sd.Density        0.266   0.207   0.010   0.104   0.220   0.382   0.760 1.001  4900
      sd.Int            0.928   0.268   0.451   0.745   0.913   1.090   1.491 1.001  4900
      sd.Resid          1.002   0.054   0.929   0.963   0.992   1.029   1.140 1.001  4900
      sd.Season         1.765   0.212   1.349   1.628   1.762   1.904   2.193 1.002  1900
      sigma.res         1.032   0.129   0.818   0.940   1.019   1.111   1.316 1.001  4900
      deviance        121.319   4.838 114.256 117.725 120.497 124.008 133.096 1.001  3300
      
      For each parameter, n.eff is a crude measure of effective sample size,
      and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
      
      DIC info (using the rule, pD = var(deviance)/2)
      pD = 11.7 and DIC = 133.0
      DIC is an estimate of expected predictive error (lower deviance is better).
      
  2. Before fully exploring the parameters, it is prudent to examine the convergence and mixing diagnostics. Chose either the effects or means parameterization model (they should yield much the same).
    • Trace plots
      View code
      > plot(as.mcmc(quinn.r2jags))
      
    • Raftery diagnostic
      View code
      > raftery.diag(as.mcmc(quinn.r2jags))
      
      [[1]]
      
      Quantile (q) = 0.025
      Accuracy (r) = +/- 0.005
      Probability (s) = 0.95 
      
      You need a sample size of at least 3746 with these values of q, r and s
      
      [[2]]
      
      Quantile (q) = 0.025
      Accuracy (r) = +/- 0.005
      Probability (s) = 0.95 
      
      You need a sample size of at least 3746 with these values of q, r and s
      
      [[3]]
      
      Quantile (q) = 0.025
      Accuracy (r) = +/- 0.005
      Probability (s) = 0.95 
      
      You need a sample size of at least 3746 with these values of q, r and s
      
      
    • Autocorrelation
      View code
      > autocorr.diag(as.mcmc(quinn.r2jags))
      
               Mean[1,1] Mean[2,1] Mean[3,1]  Mean[4,1] Mean[1,2] Mean[2,2]  Mean[3,2] Mean[4,2]
      Lag 0    1.0000000  1.000000  1.000000  1.0000000  1.000000  1.000000  1.000e+00  1.000000
      Lag 10  -0.0082297  0.004866 -0.014512  0.0004189  0.019037  0.013411  4.773e-03 -0.008294
      Lag 50   0.0008819 -0.003041 -0.005565  0.0035012  0.007405  0.013454 -1.337e-06  0.016589
      Lag 100 -0.0060674 -0.009696 -0.006897 -0.0111460 -0.001689  0.014131 -2.677e-03  0.013051
      Lag 500 -0.0119590 -0.001043 -0.019341  0.0072319 -0.008894 -0.009759 -2.027e-02  0.002028
              beta.density[1] beta.density[2] beta.int[1,1] beta.int[2,1] beta.int[3,1] beta.int[4,1]
      Lag 0               NaN        1.000000           NaN           NaN           NaN           NaN
      Lag 10              NaN        0.003179           NaN           NaN           NaN           NaN
      Lag 50              NaN        0.007616           NaN           NaN           NaN           NaN
      Lag 100             NaN       -0.010429           NaN           NaN           NaN           NaN
      Lag 500             NaN       -0.011344           NaN           NaN           NaN           NaN
              beta.int[1,2] beta.int[2,2] beta.int[3,2] beta.int[4,2] beta.season[1] beta.season[2]
      Lag 0             NaN     1.0000000      1.000000      1.000000            NaN       1.000000
      Lag 10            NaN     0.0299663     -0.007553      0.005904            NaN      -0.011382
      Lag 50            NaN     0.0000938      0.003024     -0.005636            NaN      -0.006023
      Lag 100           NaN     0.0022546     -0.012278     -0.008113            NaN      -0.010964
      Lag 500           NaN    -0.0093048     -0.002738     -0.025463            NaN      -0.010620
              beta.season[3] beta.season[4]      beta0  deviance    eff[1]    eff[2]    eff[3]    eff[4]
      Lag 0          1.00000       1.000000  1.0000000  1.000000  1.000000  1.000000  1.000000  1.000000
      Lag 10        -0.02272      -0.002790 -0.0082297 -0.004244  0.003179  0.028275 -0.001557 -0.009585
      Lag 50         0.01592       0.004193  0.0008819  0.001471  0.007616 -0.004404 -0.015173  0.006090
      Lag 100       -0.01495      -0.017531 -0.0060674  0.008789 -0.010429  0.015171 -0.019009 -0.009054
      Lag 500       -0.01460      -0.019725 -0.0119590  0.008385 -0.011344 -0.007830 -0.026246 -0.003047
              sd.Density    sd.Int  sd.Resid  sd.Season sigma.res
      Lag 0     1.000000  1.000000  1.000000  1.0000000  1.000000
      Lag 10   -0.016133 -0.017361 -0.013893  0.0028862  0.008276
      Lag 50    0.006714 -0.005050  0.017671 -0.0029169 -0.009613
      Lag 100   0.019943  0.004759  0.008819 -0.0102643 -0.025626
      Lag 500   0.014718 -0.025343  0.013722  0.0006947 -0.016266
      
  3. A summary figure would be nice
    View code
    > quinn.mcmc <- as.mcmc(quinn.r2jags$BUGSoutput$sims.matrix)
    > #get the column indexes for the parameters we are interested in
    > quinn.sum <- summary.mcmc(quinn.mcmc, terms=c("Mean")) #or "alpha"
    > quinn.sum <- cbind(quinn.sum, expand.grid(Season=c("Autumn","Spring","Summer","Winter"), Density=c("High","Low")))
    > 
    > quinn.sum$Season <- factor(quinn.sum$Season, levels=c("Spring","Summer","Autumn","Winter"))
    > library(scales)
    > 
    > p1<-ggplot(quinn.sum,aes(y=mean, x=Season, group=Density))+
    + geom_errorbar(aes(ymin=lower, ymax=upper),width=0,size=1.5,position=position_dodge(0.2))+
    + geom_errorbar(aes(ymin=lower.1, ymax=upper.1),width=0,position=position_dodge(0.2))+
    + geom_point(aes(shape=Density, fill=Density),size=3,position=position_dodge(0.2))+
    + scale_y_continuous("Number of mussel recruits",
    +   trans=trans_new("x2",function(x) x^2,"sqrt"),
    +   breaks=sqrt(pretty(c(0,60))),
    +   labels=pretty(c(0,60)))+
    + scale_x_discrete("Season")+
    + scale_shape_manual(name="Density",values=c(21,16), breaks=c("High","Low"), labels=c("High","Low"))+
    + scale_fill_manual(name="Density",values=c("white","black"), breaks=c("High","Low"), labels=c("High","Low"))+
    + murray_opts+opts(legend.position=c(1,0.99),legend.justification=c(1,1))
    > p1
    
    > 
    + 
    + #Finite-population standard deviations
    > quinn.sd <- summary.mcmc(quinn.mcmc, terms=c("sd"))
    > rownames(quinn.sd) <- c("Density","Season:Density","Residuals","Season")
    > quinn.sd$name <- factor(c("Density","Season:Density","Residuals","Season"), levels=c("Residuals", "Season:Density", "Density","Season"))
    > 
    > p2<-ggplot(quinn.sd,aes(y=name, x=median))+
    + geom_vline(xintercept=0,linetype="dashed")+
    + geom_hline(xintercept=0)+
    + scale_x_continuous("Finite population \nvariance components (sd)")+
    + geom_errorbarh(aes(xmin=lower.1, xmax=upper.1), height=0, size=1)+
    + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0, size=1.5)+
    + geom_point(size=3)+
    + geom_text(aes(label=sprintf("(%4.1f%%)",Perc),vjust=-1))+
    + murray_opts+
    + opts(axis.line=theme_blank(),axis.title.y=theme_blank(),
    + axis.text.y=theme_text(size=12,hjust=1))
    > 
    > grid.newpage()
    > pushViewport(viewport(layout=grid.layout(1,2,5,5,respect=TRUE)))
    > pushViewport(viewport(layout.pos.col=1, layout.pos.row=1))
    > print(p1,newpage=FALSE)
    > popViewport(1)
    > pushViewport(viewport(layout.pos.col=2, layout.pos.row=1))
    > print(p2,newpage=FALSE)
    > popViewport(0)
    
    > 
    + 
    + #eff
    > quinn.eff <- summary.mcmc(quinn.mcmc, terms=c("eff\\["))
    > nms <- c("Autumn","Spring","Summer","Winter")
    > quinn.eff$name <- nms
    > 
    > p3 <- ggplot(quinn.eff, aes(y=name, x=mean))+
    + geom_vline(xintercept=0,linetype="dashed")+
    + geom_hline(xintercept=0)+
    + geom_errorbarh(aes(xmin=lower.1, xmax=upper.1), height=0, size=1)+
    + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0, size=1.5)+
    + geom_point(size=3)+murray_opts+
    + scale_x_continuous(name="Effect size")+
    + opts(axis.line=theme_blank(),axis.title.y=theme_blank(),
    + axis.text.y=theme_text(size=12,hjust=1))
    > p3
    
    > 
    + 
    + 
    + #cohen's D
    > #eff<-NULL
    > #eff.sd <- NULL
    > #for (i in seq(1,7,b=2)) {
    > #  eff <-cbind(eff,quinn.mcmc[,i] - quinn.mcmc[,i+1])
    > #  eff.sd <- cbind(eff.sd,mean(apply(quinn.mcmc[,i:(i+1)],2,sd)))
    > #}
    > #cohen <- sweep(eff,2,eff.sd,"/")
    > #adply(cohen,2,function(x) {
    > #   data.frame(mean=mean(x), median=median(x),
    > #		 HPDinterval(as.mcmc(x), p=0.68), HPDinterval(as.mcmc(x)))
    > #})
    

  Frequentist pooled variances t-test

End of instructions