43  ANOVA - M3L8

Bayesian Statistics: Techniques and Models

An overview of ANOVA in the context of Bayesian statistics.
Monte Carlo Estimation
Author

Oren Bochman

Keywords

ANOVA, Bayesian statistics, R programming, statistical modeling, Analysis of Variance

43.1 Introduction to ANOVA 🎥

Introduction to ANOVA

Introduction to ANOVA

43.2 One way ANOVA model using JAGS

43.2.1 Data & EDA

As an example of a one-way ANOVA, we’ll look at the Plant Growth data in R.

Listing 43.1: Plant Growth Query
data("PlantGrowth")
#?PlantGrowth
head(PlantGrowth)
  weight group
1   4.17  ctrl
2   5.58  ctrl
3   5.18  ctrl
4   6.11  ctrl
5   4.50  ctrl
6   4.61  ctrl

We first load the dataset (Listing 43.1)

Because the explanatory variable group is a factor and not continuous, we choose to visualize the data with box plots rather than scatter plots.

boxplot(weight ~ group, data=PlantGrowth)
Figure 43.1: PlantGrowth boxplot

The box plots summarize the distribution of the data for each of the three groups. It appears that treatment 2 has the highest mean yield. It might be questionable whether each group has the same variance, but we’ll assume that is the case.

43.2.2 Modeling

Again, we can start with the reference analysis (with a noninformative prior) with a linear model in R.

lmod = lm(weight ~ group, data=PlantGrowth)
summary(lmod)

Call:
lm(formula = weight ~ group, data = PlantGrowth)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0710 -0.4180 -0.0060  0.2627  1.3690 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.0320     0.1971  25.527   <2e-16 ***
grouptrt1    -0.3710     0.2788  -1.331   0.1944    
grouptrt2     0.4940     0.2788   1.772   0.0877 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6234 on 27 degrees of freedom
Multiple R-squared:  0.2641,    Adjusted R-squared:  0.2096 
F-statistic: 4.846 on 2 and 27 DF,  p-value: 0.01591
anova(lmod)
Analysis of Variance Table

Response: weight
          Df  Sum Sq Mean Sq F value  Pr(>F)  
group      2  3.7663  1.8832  4.8461 0.01591 *
Residuals 27 10.4921  0.3886                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(lmod) # for graphical residual analysis
Figure 43.2: Graphical residual analysis
Figure 43.3: Graphical residual analysis
Figure 43.4: Graphical residual analysis
Figure 43.5: Graphical residual analysis

The default model structure in R is the linear model with dummy indicator variables. Hence, the “intercept” in this model is the mean yield for the control group. The two other parameters are the estimated effects of treatments 1 and 2. To recover the mean yield in treatment group 1, you would add the intercept term and the treatment 1 effect. To see how R sets the model up, use the model.matrix(lmod) function to extract the X matrix.

The anova() function in R compares variability of observations between the treatment groups to variability within the treatment groups to test whether all means are equal or whether at least one is different. The small p-value here suggests that the means are not all equal.

Let’s fit the cell means model in JAGS.

library("rjags")
mod_string = " model {
    for (i in 1:length(y)) {
        y[i] ~ dnorm(mu[grp[i]], prec)
    }
    
    for (j in 1:3) {
        mu[j] ~ dnorm(0.0, 1.0/1.0e6)
    }
    
    prec ~ dgamma(5/2.0, 5*1.0/2.0)
    sig = sqrt( 1.0 / prec )
} "

set.seed(82)
str(PlantGrowth)
'data.frame':   30 obs. of  2 variables:
 $ weight: num  4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
 $ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...
data_jags = list(y=PlantGrowth$weight, 
              grp=as.numeric(PlantGrowth$group))

params = c("mu", "sig")

inits = function() {
    inits = list("mu"=rnorm(3,0.0,100.0), "prec"=rgamma(1,1.0,1.0))
}

mod = jags.model(textConnection(mod_string), data=data_jags, inits=inits, n.chains=3)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 30
   Unobserved stochastic nodes: 4
   Total graph size: 74

Initializing model
update(mod, 1e3)

mod_sim = coda.samples(model=mod,
                        variable.names=params,
                        n.iter=5e3)
mod_csim = as.mcmc(do.call(rbind, mod_sim)) # combined chains

43.2.3 Model checking

As usual, we check for convergence of our MCMC.

par(mar = c(2.5, 1, 2.5, 1))
plot(mod_sim)

gelman.diag(mod_sim)
Potential scale reduction factors:

      Point est. Upper C.I.
mu[1]          1          1
mu[2]          1          1
mu[3]          1          1
sig            1          1

Multivariate psrf

1
autocorr.diag(mod_sim)
              mu[1]         mu[2]        mu[3]          sig
Lag 0   1.000000000  1.0000000000  1.000000000  1.000000000
Lag 1  -0.003295982  0.0032550285  0.014229944  0.093790272
Lag 5  -0.002099076  0.0004635964 -0.015080105 -0.004950599
Lag 10 -0.005994920 -0.0016405035  0.002439778  0.003699796
Lag 50  0.007647856  0.0027188925 -0.008915493 -0.008770126
effectiveSize(mod_sim)
   mu[1]    mu[2]    mu[3]      sig 
15000.00 15152.88 14977.57 11527.49 
Figure 43.6: MCMC convergence diagnostics

We can also look at the residuals to see if there are any obvious problems with our model choice.

(pm_params = colMeans(mod_csim))
    mu[1]     mu[2]     mu[3]       sig 
5.0311341 4.6644016 5.5220571 0.7139621 
yhat = pm_params[1:3][data_jags$grp]
resid = data_jags$y - yhat
plot(resid)
Figure 43.7: Residuals vs Index
plot(yhat, resid)
Figure 43.8: Residuals vs Fitted values for PlantGrowth model

Again, it might be appropriate to have a separate variance for each group. We will have you do that as an exercise.

43.2.4 Results

Let’s look at the posterior summary of the parameters.

summary(mod_sim)

Iterations = 1001:6000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

       Mean      SD  Naive SE Time-series SE
mu[1] 5.031 0.22776 0.0018596      0.0018596
mu[2] 4.664 0.22758 0.0018581      0.0018490
mu[3] 5.522 0.22782 0.0018602      0.0018615
sig   0.714 0.09329 0.0007617      0.0008709

2. Quantiles for each variable:

        2.5%    25%    50%    75%  97.5%
mu[1] 4.5813 4.8780 5.0304 5.1825 5.4792
mu[2] 4.2164 4.5151 4.6630 4.8149 5.1154
mu[3] 5.0677 5.3731 5.5225 5.6728 5.9738
sig   0.5592 0.6484 0.7042 0.7684 0.9269
HPDinterval(mod_csim)
          lower     upper
mu[1] 4.5751119 5.4713908
mu[2] 4.1972351 5.0959196
mu[3] 5.0870461 5.9868314
sig   0.5415318 0.8962747
attr(,"Probability")
[1] 0.95

The HPDinterval() function in the coda package calculates intervals of highest posterior density for each parameter.

We are interested to know if one of the treatments increases mean yield. It is clear that treatment 1 does not. What about treatment 2?

mean(mod_csim[,3] > mod_csim[,1])
[1] 0.9368667

There is a high posterior probability that the mean yield for treatment 2 is greater than the mean yield for the control group.

It may be the case that treatment 2 would be costly to put into production. Suppose that to be worthwhile, this treatment must increase mean yield by 10%. What is the posterior probability that the increase is at least that?

mean(mod_csim[,3] > 1.1*mod_csim[,1])
[1] 0.4848667

We have about 50/50 odds that adopting treatment 2 would increase mean yield by at least 10%.

43.3 Two Factor ANOVA

43.3.1 Data

Let’s explore an example with two factors. We’ll use the Warpbreaks data set in R. Check the documentation for a description of the data by typing ?warpbreaks.

data("warpbreaks")
#?warpbreaks
head(warpbreaks)
  breaks wool tension
1     26    A       L
2     30    A       L
3     54    A       L
4     25    A       L
5     70    A       L
6     52    A       L
Table 43.1: Preview of first few rows of warpbreaks data
# This chunk is for displaying the output that was previously static.
# If the static output below is preferred, this chunk can be removed 
# and the static output remains unlabelled as it's not a code cell.
# For a labeled table, this chunk should generate it.
# The original file had static output here:
##   breaks wool tension
## 1     26    A       L
## 2     30    A       L
## 3     54    A       L
## 4     25    A       L
## 5     70    A       L
## 6     52    A       L
# To make this a labeled table from code:
head(warpbreaks)
  breaks wool tension
1     26    A       L
2     30    A       L
3     54    A       L
4     25    A       L
5     70    A       L
6     52    A       L
Table 43.2: Contingency table of wool type vs tension level
table(warpbreaks$wool, warpbreaks$tension)
   
    L M H
  A 9 9 9
  B 9 9 9

Again, we visualize the data with box plots.

boxplot(breaks ~ wool + tension, data=warpbreaks)
Figure 43.9: Warpbreaks boxplot
boxplot(log(breaks) ~ wool + tension, data=warpbreaks)
Figure 43.10: Warpbreaks boxplot with log-transformed breaks

The different groups have more similar variance if we use the logarithm of breaks. From this visualization, it looks like both factors may play a role in the number of breaks. It appears that there is a general decrease in breaks as we move from low to medium to high tension. Let’s start with a one-way model using tension only.

43.3.2 One-way model

mod1_string = " model {
    for( i in 1:length(y)) {
        y[i] ~ dnorm(mu[tensGrp[i]], prec)
    }
    
    for (j in 1:3) {
        mu[j] ~ dnorm(0.0, 1.0/1.0e6)
    }
    
    prec ~ dgamma(5/2.0, 5*2.0/2.0)
    sig = sqrt(1.0 / prec)
} "

set.seed(83)
str(warpbreaks)
'data.frame':   54 obs. of  3 variables:
 $ breaks : num  26 30 54 25 70 52 51 26 67 18 ...
 $ wool   : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ tension: Factor w/ 3 levels "L","M","H": 1 1 1 1 1 1 1 1 1 2 ...
data1_jags = list(y=log(warpbreaks$breaks), tensGrp=as.numeric(warpbreaks$tension))

params1 = c("mu", "sig")

mod1 = jags.model(textConnection(mod1_string), data=data1_jags, n.chains=3)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 54
   Unobserved stochastic nodes: 4
   Total graph size: 123

Initializing model
update(mod1, 1e3)

mod1_sim = coda.samples(model=mod1,
                        variable.names=params1,
                        n.iter=5e3)
## convergence diagnostics
plot(mod1_sim)
Figure 43.11: MCMC convergence diagnostics for one-way tension model
gelman.diag(mod1_sim)
Potential scale reduction factors:

      Point est. Upper C.I.
mu[1]          1          1
mu[2]          1          1
mu[3]          1          1
sig            1          1

Multivariate psrf

1
autocorr.diag(mod1_sim)
              mu[1]        mu[2]         mu[3]          sig
Lag 0   1.000000000  1.000000000  1.0000000000  1.000000000
Lag 1  -0.001543382 -0.011482925 -0.0048970151  0.051048754
Lag 5  -0.006377853 -0.017028742 -0.0054862444 -0.021399059
Lag 10  0.001367067  0.009851261  0.0065758333 -0.001345930
Lag 50  0.021506593  0.012051621 -0.0006951465  0.001122476
effectiveSize(mod1_sim)
   mu[1]    mu[2]    mu[3]      sig 
14721.39 14972.50 15899.17 13717.48 

The 95% posterior interval for the mean of group 2 (medium tension) overlaps with both the low and high groups, but the intervals for low and high group only slightly overlap. That is a pretty strong indication that the means for low and high tension are different. Let’s collect the DIC for this model and move on to the two-way model.

dic1 = dic.samples(mod1, n.iter=1e3)

43.3.3 Two-way additive model

With two factors, one with two levels and the other with three, we have six treatment groups, which is the same situation we discussed when introducing multiple factor ANOVA. We will first fit the additive model which treats the two factors separately with no interaction. To get the X matrix (or design matrix) for this model, we can create it in R.

Table 43.3: Head of the design matrix for the additive model
X = model.matrix( ~ wool + tension, data=warpbreaks)
head(X)
  (Intercept) woolB tensionM tensionH
1           1     0        0        0
2           1     0        0        0
3           1     0        0        0
4           1     0        0        0
5           1     0        0        0
6           1     0        0        0
Table 43.4: Tail of the design matrix for the additive model
tail(X)
   (Intercept) woolB tensionM tensionH
49           1     1        0        1
50           1     1        0        1
51           1     1        0        1
52           1     1        0        1
53           1     1        0        1
54           1     1        0        1

By default, R has chosen the mean for wool A and low tension to be the intercept. Then, there is an effect for wool B, and effects for medium tension and high tension, each associated with dummy indicator variables.

mod2_string = " model {
    for( i in 1:length(y)) {
        y[i] ~ dnorm(mu[i], prec)
        mu[i] = int + alpha*isWoolB[i] + beta[1]*isTensionM[i] + beta[2]*isTensionH[i]
    }
    
    int ~ dnorm(0.0, 1.0/1.0e6)
    alpha ~ dnorm(0.0, 1.0/1.0e6)
    for (j in 1:2) {
        beta[j] ~ dnorm(0.0, 1.0/1.0e6)
    }
    
    prec ~ dgamma(3/2.0, 3*1.0/2.0)
    sig = sqrt(1.0 / prec)
} "

data2_jags = list(y=log(warpbreaks$breaks), isWoolB=X[,"woolB"], isTensionM=X[,"tensionM"], isTensionH=X[,"tensionH"])

params2 = c("int", "alpha", "beta", "sig")

mod2 = jags.model(textConnection(mod2_string), data=data2_jags, n.chains=3)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 54
   Unobserved stochastic nodes: 5
   Total graph size: 243

Initializing model
update(mod2, 1e3)

mod2_sim = coda.samples(model=mod2,
                        variable.names=params2,
                        n.iter=5e3)
## convergence diagnostics
plot(mod2_sim)

gelman.diag(mod2_sim)    # Corrected from mod1_sim
Potential scale reduction factors:

        Point est. Upper C.I.
alpha            1          1
beta[1]          1          1
beta[2]          1          1
int              1          1
sig              1          1

Multivariate psrf

1
autocorr.diag(mod2_sim)  # Corrected from mod1_sim
               alpha      beta[1]      beta[2]         int         sig
Lag 0   1.000000e+00  1.000000000  1.000000000 1.000000000 1.000000000
Lag 1   4.949732e-01  0.493664011  0.488555715 0.740486160 0.069036596
Lag 5   3.798922e-02  0.105554167  0.095890680 0.166494047 0.014138106
Lag 10 -4.884794e-05 -0.008619817 -0.010211226 0.001239427 0.015050523
Lag 50 -2.346203e-04 -0.003652400 -0.001246754 0.006844612 0.002058968
effectiveSize(mod2_sim) # Corrected from mod1_sim
    alpha   beta[1]   beta[2]       int       sig 
 4996.075  4095.918  4037.636  2551.051 12259.928 
Figure 43.12: Convergence and diagnostics for the additive two-way ANOVA model
Figure 43.13: Convergence and diagnostics for the additive two-way ANOVA model

Let’s summarize the results, collect the DIC for this model, and compare it to the first one-way model.

summary(mod2_sim)

Iterations = 1001:6000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean      SD  Naive SE Time-series SE
alpha   -0.1521 0.12461 0.0010174      0.0017642
beta[1] -0.2886 0.15075 0.0012309      0.0023670
beta[2] -0.4901 0.15036 0.0012277      0.0023658
int      3.5769 0.12139 0.0009912      0.0024042
sig      0.4536 0.04517 0.0003688      0.0004085

2. Quantiles for each variable:

           2.5%     25%     50%      75%     97.5%
alpha   -0.3959 -0.2344 -0.1512 -0.07018  0.094806
beta[1] -0.5850 -0.3891 -0.2905 -0.18832  0.007825
beta[2] -0.7830 -0.5905 -0.4903 -0.38986 -0.192312
int      3.3352  3.4954  3.5765  3.65790  3.817546
sig      0.3758  0.4216  0.4502  0.48164  0.551043
(dic2 = dic.samples(mod2, n.iter=1e3))
Mean deviance:  55.62 
penalty 5.221 
Penalized deviance: 60.84 
dic1
Mean deviance:  66.45 
penalty 4.064 
Penalized deviance: 70.52 

This suggests there is much to be gained adding the wool factor to the model. Before we settle on this model however, we should consider whether there is an interaction. Let’s look again at the box plot with all six treatment groups.

boxplot(log(breaks) ~ wool + tension, data=warpbreaks)
Figure 43.14: Re-examining boxplot of log(breaks) by wool and tension for interaction effects

Our two-way model has a single effect for wool B and the estimate is negative. If this is true, then we would expect wool B to be associated with fewer breaks than its wool A counterpart on average. This is true for low and high tension, but it appears that breaks are higher for wool B when there is medium tension. That is, the effect for wool B is not consistent across tension levels, so it may appropriate to add an interaction term. In R, this would look like:

lmod2 = lm(log(breaks) ~ .^2, data=warpbreaks)
summary(lmod2)

Call:
lm(formula = log(breaks) ~ .^2, data = warpbreaks)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.81504 -0.27885  0.04042  0.27319  0.64358 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      3.7179     0.1247  29.824  < 2e-16 ***
woolB           -0.4356     0.1763  -2.471  0.01709 *  
tensionM        -0.6012     0.1763  -3.410  0.00133 ** 
tensionH        -0.6003     0.1763  -3.405  0.00134 ** 
woolB:tensionM   0.6281     0.2493   2.519  0.01514 *  
woolB:tensionH   0.2221     0.2493   0.891  0.37749    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.374 on 48 degrees of freedom
Multiple R-squared:  0.3363,    Adjusted R-squared:  0.2672 
F-statistic: 4.864 on 5 and 48 DF,  p-value: 0.001116

Adding the interaction, we get an effect for being in wool B and medium tension, as well as for being in wool B and high tension. There are now six parameters for the mean, one for each treatment group, so this model is equivalent to the full cell means model. Let’s use that.

43.3.4 Two-way cell means model

In this new model, \mu will be a matrix with six entries, each corresponding to a treatment group.

mod3_string = " model {
    for( i in 1:length(y)) {
        y[i] ~ dnorm(mu[woolGrp[i], tensGrp[i]], prec)
    }
    
    for (j in 1:max(woolGrp)) {
        for (k in 1:max(tensGrp)) {
            mu[j,k] ~ dnorm(0.0, 1.0/1.0e6)
        }
    }
    
    prec ~ dgamma(3/2.0, 3*1.0/2.0)
    sig = sqrt(1.0 / prec)
} "

str(warpbreaks)
'data.frame':   54 obs. of  3 variables:
 $ breaks : num  26 30 54 25 70 52 51 26 67 18 ...
 $ wool   : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ tension: Factor w/ 3 levels "L","M","H": 1 1 1 1 1 1 1 1 1 2 ...
data3_jags = list(y=log(warpbreaks$breaks), woolGrp=as.numeric(warpbreaks$wool), tensGrp=as.numeric(warpbreaks$tension))

params3 = c("mu", "sig")

mod3 = jags.model(textConnection(mod3_string), data=data3_jags, n.chains=3)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 54
   Unobserved stochastic nodes: 7
   Total graph size: 179

Initializing model
update(mod3, 1e3)

mod3_sim = coda.samples(model=mod3,
                        variable.names=params3,
                        n.iter=5e3)
mod3_csim = as.mcmc(do.call(rbind, mod3_sim))
plot(mod3_sim)
Figure 43.15: Traceplots for the cell means model
Figure 43.16: Traceplots for the cell means model
## convergence diagnostics
gelman.diag(mod3_sim)
Potential scale reduction factors:

        Point est. Upper C.I.
mu[1,1]          1          1
mu[2,1]          1          1
mu[1,2]          1          1
mu[2,2]          1          1
mu[1,3]          1          1
mu[2,3]          1          1
sig              1          1

Multivariate psrf

1
autocorr.diag(mod3_sim)
             mu[1,1]       mu[2,1]      mu[1,2]      mu[2,2]      mu[1,3]
Lag 0   1.0000000000  1.0000000000  1.000000000  1.000000000  1.000000000
Lag 1  -0.0031449739  0.0088287336  0.006110741 -0.006245016  0.007298649
Lag 5   0.0026387814 -0.0053453974  0.008388578 -0.007357227 -0.003745448
Lag 10  0.0009629569 -0.0032304891  0.015507145 -0.011358803 -0.007156927
Lag 50 -0.0037293287 -0.0007610572 -0.009355279 -0.008208487  0.010905305
             mu[2,3]          sig
Lag 0   1.000000e+00  1.000000000
Lag 1  -1.290280e-02  0.114523745
Lag 5   2.342043e-05  0.006218406
Lag 10  6.727064e-03  0.010898163
Lag 50  2.500647e-03 -0.004887363
effectiveSize(mod3_sim)
 mu[1,1]  mu[2,1]  mu[1,2]  mu[2,2]  mu[1,3]  mu[2,3]      sig 
15000.00 14774.09 14325.64 14967.22 14774.57 15350.80 11208.04 
raftery.diag(mod3_sim)
[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                               
         Burn-in  Total Lower bound  Dependence
         (M)      (N)   (Nmin)       factor (I)
 mu[1,1] 2        3741  3746         0.999     
 mu[2,1] 2        3680  3746         0.982     
 mu[1,2] 2        3741  3746         0.999     
 mu[2,2] 2        3741  3746         0.999     
 mu[1,3] 2        3866  3746         1.030     
 mu[2,3] 2        3680  3746         0.982     
 sig     2        3680  3746         0.982     


[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                               
         Burn-in  Total Lower bound  Dependence
         (M)      (N)   (Nmin)       factor (I)
 mu[1,1] 2        3866  3746         1.030     
 mu[2,1] 2        3930  3746         1.050     
 mu[1,2] 2        3803  3746         1.020     
 mu[2,2] 2        3803  3746         1.020     
 mu[1,3] 3        4129  3746         1.100     
 mu[2,3] 2        3680  3746         0.982     
 sig     1        3712  3746         0.991     


[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                               
         Burn-in  Total Lower bound  Dependence
         (M)      (N)   (Nmin)       factor (I)
 mu[1,1] 2        3620  3746         0.966     
 mu[2,1] 2        3803  3746         1.020     
 mu[1,2] 2        3803  3746         1.020     
 mu[2,2] 2        3803  3746         1.020     
 mu[1,3] 2        3803  3746         1.020     
 mu[2,3] 2        3741  3746         0.999     
 sig     2        3680  3746         0.982     

Let’s compute the DIC and compare with our previous models.

(dic3 = dic.samples(mod3, n.iter=1e3))
Mean deviance:  52.05 
penalty 7.185 
Penalized deviance: 59.24 
dic2
Mean deviance:  55.62 
penalty 5.221 
Penalized deviance: 60.84 
dic1
Mean deviance:  66.45 
penalty 4.064 
Penalized deviance: 70.52 

This suggests that the full model with interaction between wool and tension (which is equivalent to the cell means model) is the best for explaining/predicting warp breaks.

43.3.5 Results

summary(mod3_sim)

Iterations = 1001:6000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean      SD  Naive SE Time-series SE
mu[1,1] 3.7179 0.14854 0.0012128      0.0012129
mu[2,1] 3.2821 0.14809 0.0012091      0.0012186
mu[1,2] 3.1165 0.14918 0.0012180      0.0012488
mu[2,2] 3.3081 0.15043 0.0012283      0.0012297
mu[1,3] 3.1174 0.14833 0.0012111      0.0012210
mu[2,3] 2.9058 0.14723 0.0012021      0.0011887
sig     0.4427 0.04511 0.0003683      0.0004282

2. Quantiles for each variable:

          2.5%   25%    50%    75%  97.5%
mu[1,1] 3.4207 3.620 3.7179 3.8156 4.0118
mu[2,1] 2.9935 3.183 3.2822 3.3803 3.5780
mu[1,2] 2.8176 3.018 3.1158 3.2159 3.4132
mu[2,2] 3.0088 3.209 3.3083 3.4075 3.6034
mu[1,3] 2.8264 3.019 3.1170 3.2164 3.4073
mu[2,3] 2.6189 2.806 2.9058 3.0038 3.1987
sig     0.3636 0.411 0.4391 0.4708 0.5421
HPDinterval(mod3_csim)
            lower     upper
mu[1,1] 3.4211281 4.0119744
mu[2,1] 2.9850428 3.5675903
mu[1,2] 2.8162798 3.4112760
mu[2,2] 3.0137257 3.6068973
mu[1,3] 2.8229589 3.4028172
mu[2,3] 2.6252050 3.2036412
sig     0.3589913 0.5347668
attr(,"Probability")
[1] 0.95
par(mfrow=c(3,2)) # arrange frame for plots
densplot(mod3_csim[,1:6], xlim=c(2.0, 4.5))
Figure 43.17: Posterior densities for cell means

It might be tempting to look at comparisons between each combination of treatments, but we warn that this could yield spurious results. When we discussed the statistical modeling cycle, we said it is best not to search your results for interesting hypotheses, because if there are many hypotheses, some will appear to show “effects” or “associations” simply due to chance. Results are most reliable when we determine a relatively small number of hypotheses we are interested in beforehand, collect the data, and statistically evaluate the evidence for them.

One question we might be interested in with these data is finding the treatment combination that produces the fewest breaks. To calculate this, we can go through our posterior samples and for each sample, find out which group has the smallest mean. These counts help us determine the posterior probability that each of the treatment groups has the smallest mean.

Table 43.5: Posterior probabilities of each treatment group having the smallest mean break rate
prop.table( table( apply(mod3_csim[,1:6], 1, which.min) ) )

         2          3          4          5          6 
0.01653333 0.11993333 0.01173333 0.11546667 0.73633333 

The evidence supports wool B with high tension as the treatment that produces the fewest breaks.