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 ctrl43.1 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.
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)
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.01591anova(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 ' ' 1plot(lmod) # for 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 modelupdate(mod, 1e3)
mod_sim = coda.samples(model=mod,
variable.names=params,
n.iter=5e3)
mod_csim = as.mcmc(do.call(rbind, mod_sim)) # combined chains43.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
1autocorr.diag(mod_sim) mu[1] mu[2] mu[3] sig
Lag 0 1.000000000 1.000000000 1.000000000 1.000000000
Lag 1 0.010128434 -0.017719050 0.019278559 0.087841309
Lag 5 0.002314983 0.001378619 0.005718167 0.004071097
Lag 10 0.019904058 -0.012208490 0.005288907 0.001012025
Lag 50 -0.004881143 -0.016100169 0.007174652 0.002744497effectiveSize(mod_sim) mu[1] mu[2] mu[3] sig
14759.76 15414.04 14493.54 12577.46 
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.0334806 4.6590470 5.5287783 0.7123182 yhat = pm_params[1:3][data_jags$grp]
resid = data_jags$y - yhat
plot(resid)
plot(yhat, resid)
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.0335 0.22598 0.0018451 0.0018605
mu[2] 4.6590 0.22839 0.0018648 0.0018411
mu[3] 5.5288 0.22844 0.0018652 0.0019004
sig 0.7123 0.09213 0.0007522 0.0008219
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
mu[1] 4.5940 4.883 5.0332 5.1828 5.4783
mu[2] 4.2114 4.508 4.6572 4.8109 5.1061
mu[3] 5.0756 5.378 5.5284 5.6793 5.9748
sig 0.5573 0.647 0.7028 0.7685 0.9171HPDinterval(mod_csim) lower upper
mu[1] 4.5934199 5.4775181
mu[2] 4.1935421 5.0850689
mu[3] 5.0739542 5.9725605
sig 0.5473127 0.8999281
attr(,"Probability")
[1] 0.95The 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.9382There 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.4906667We 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# 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 Ltable(warpbreaks$wool, warpbreaks$tension)
L M H
A 9 9 9
B 9 9 9Again, we visualize the data with box plots.
boxplot(breaks ~ wool + tension, data=warpbreaks)
boxplot(log(breaks) ~ wool + tension, data=warpbreaks)
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 modelupdate(mod1, 1e3)
mod1_sim = coda.samples(model=mod1,
variable.names=params1,
n.iter=5e3)## convergence diagnostics
plot(mod1_sim)
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
1autocorr.diag(mod1_sim) mu[1] mu[2] mu[3] sig
Lag 0 1.0000000000 1.000000000 1.000000000 1.000000000
Lag 1 0.0001393307 0.020800782 0.009341394 0.041307730
Lag 5 -0.0038219169 -0.007453829 0.001548361 -0.004523855
Lag 10 -0.0029031913 0.008130131 0.005825012 0.010440068
Lag 50 -0.0011783744 0.005062078 -0.008345652 -0.008374586effectiveSize(mod1_sim) mu[1] mu[2] mu[3] sig
15000.00 14680.53 15000.00 14889.24 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.
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 0tail(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 1By 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 modelupdate(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_simPotential 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
1autocorr.diag(mod2_sim) # Corrected from mod1_sim alpha beta[1] beta[2] int sig
Lag 0 1.000000000 1.0000000000 1.00000000 1.00000000 1.000000000
Lag 1 0.493282226 0.5017317196 0.50366787 0.75149396 0.071671111
Lag 5 0.008396312 0.0913447330 0.09017855 0.15572751 0.019419206
Lag 10 -0.004195011 0.0042637337 0.01992147 0.01326364 -0.009419705
Lag 50 -0.002496894 -0.0001677716 0.01443150 0.02241556 -0.002429448effectiveSize(mod2_sim) # Corrected from mod1_sim alpha beta[1] beta[2] int sig
5289.301 3858.654 3898.090 2636.393 11705.833 
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.1483 0.12499 0.0010205 0.0017246
beta[1] -0.2838 0.15263 0.0012462 0.0024574
beta[2] -0.4854 0.15410 0.0012582 0.0024702
int 3.5714 0.12519 0.0010222 0.0024378
sig 0.4545 0.04503 0.0003677 0.0004188
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha -0.3925 -0.2303 -0.1490 -0.06458 0.09815
beta[1] -0.5883 -0.3835 -0.2839 -0.18240 0.01860
beta[2] -0.7893 -0.5893 -0.4859 -0.38190 -0.18116
int 3.3225 3.4901 3.5710 3.65366 3.82012
sig 0.3763 0.4222 0.4512 0.48276 0.55174(dic2 = dic.samples(mod2, n.iter=1e3))Mean deviance: 55.66
penalty 5.309
Penalized deviance: 60.97 dic1Mean deviance: 66.62
penalty 4.107
Penalized deviance: 70.73 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)
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.001116Adding 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 modelupdate(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)
## 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
1autocorr.diag(mod3_sim) mu[1,1] mu[2,1] mu[1,2] mu[2,2] mu[1,3]
Lag 0 1.000000000 1.000000000 1.000000000 1.000000e+00 1.000000000
Lag 1 -0.006689582 -0.014285293 0.005405337 -1.853014e-02 -0.010255203
Lag 5 -0.006500033 0.006569306 -0.013855270 8.833668e-03 0.003399051
Lag 10 -0.011961976 0.008401321 0.007407778 2.288671e-05 -0.001367783
Lag 50 -0.005116613 0.009546109 -0.002036574 6.815046e-04 0.001483592
mu[2,3] sig
Lag 0 1.000000000 1.000000000
Lag 1 -0.010308626 0.112982183
Lag 5 0.004423700 0.001134386
Lag 10 -0.010499616 -0.001221645
Lag 50 -0.005008636 0.008543564effectiveSize(mod3_sim) mu[1,1] mu[2,1] mu[1,2] mu[2,2] mu[1,3] mu[2,3] sig
15215.94 16673.80 14603.70 15212.54 14990.01 15000.00 11632.76 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 3803 3746 1.020
mu[2,1] 1 3741 3746 0.999
mu[1,2] 2 3741 3746 0.999
mu[2,2] 2 3741 3746 0.999
mu[1,3] 2 3803 3746 1.020
mu[2,3] 2 3741 3746 0.999
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 3741 3746 0.999
mu[2,1] 2 3741 3746 0.999
mu[1,2] 2 3930 3746 1.050
mu[2,2] 2 3680 3746 0.982
mu[1,3] 2 3803 3746 1.020
mu[2,3] 2 3930 3746 1.050
sig 2 3930 3746 1.050
[[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 3866 3746 1.030
mu[1,2] 2 3741 3746 0.999
mu[2,2] 2 3995 3746 1.070
mu[1,3] 2 3561 3746 0.951
mu[2,3] 2 3680 3746 0.982
sig 2 3930 3746 1.050 Let’s compute the DIC and compare with our previous models.
(dic3 = dic.samples(mod3, n.iter=1e3))Mean deviance: 52.09
penalty 7.192
Penalized deviance: 59.29 dic2Mean deviance: 55.66
penalty 5.309
Penalized deviance: 60.97 dic1Mean deviance: 66.62
penalty 4.107
Penalized deviance: 70.73 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.7168 0.14987 0.0012237 0.0012153
mu[2,1] 3.2820 0.14851 0.0012126 0.0011539
mu[1,2] 3.1181 0.14837 0.0012114 0.0012292
mu[2,2] 3.3081 0.14910 0.0012174 0.0012092
mu[1,3] 3.1170 0.14847 0.0012123 0.0012134
mu[2,3] 2.9049 0.14789 0.0012075 0.0012075
sig 0.4433 0.04501 0.0003675 0.0004182
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
mu[1,1] 3.4221 3.6179 3.718 3.8159 4.0119
mu[2,1] 2.9901 3.1824 3.282 3.3822 3.5732
mu[1,2] 2.8234 3.0195 3.119 3.2156 3.4077
mu[2,2] 3.0156 3.2093 3.306 3.4072 3.6010
mu[1,3] 2.8234 3.0181 3.117 3.2158 3.4112
mu[2,3] 2.6160 2.8053 2.905 3.0035 3.1987
sig 0.3666 0.4113 0.440 0.4709 0.5425HPDinterval(mod3_csim) lower upper
mu[1,1] 3.4219871 4.0119276
mu[2,1] 2.9993067 3.5794475
mu[1,2] 2.8316098 3.4156071
mu[2,2] 3.0203452 3.6047627
mu[1,3] 2.8220146 3.4090290
mu[2,3] 2.6127096 3.1935279
sig 0.3602622 0.5338967
attr(,"Probability")
[1] 0.95par(mfrow=c(3,2)) # arrange frame for plots
densplot(mod3_csim[,1:6], xlim=c(2.0, 4.5))
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.
prop.table( table( apply(mod3_csim[,1:6], 1, which.min) ) )
1 2 3 4 5 6
6.666667e-05 1.620000e-02 1.134667e-01 1.233333e-02 1.197333e-01 7.382000e-01 The evidence supports wool B with high tension as the treatment that produces the fewest breaks.