Tukey (Tukey 1949) proposed a test for non-additvity among a pair of effects in a two-way model, without interaction. (See ‘Some Consequences …’, Table VII for examples). This test is sometimes recommended to test the assumption of non-iteraction between treatments and blocks.
We’ll use, for this discussion, data from (Littell et al. 1996). These data have been used extensively for spatial model analysis. See the associated ARM files
More important, we wish to see the trial map.
We see that the replicates are not uniform in shape, nor are blocks orthogonal to the apparent spatial fertility pattern. We should consider if there is some interaction between treatment estimates and blocks, and if spatial models reduce this interaction.
If we consider the analysis of RCB designs as a classical two-way model with interaction, we might write the model as
\[ y_{ij} = \mu + \tau _i + \rho _j + (\tau\rho)_{ij} + \epsilon _{ij} \]
where \(\tau _i\) is the treatment effect, \(\rho _j\) is the bloc (replicate) and \((\tau\rho)_{ij}\) is the term for treatment-replicate interaction.
We can’t test for the interaction effects directly, but we can test regression of the main effects. We write
\[ y_{ij} = \mu + \tau _i + \rho _j + \lambda \tau_i \rho _j + \epsilon _{ij} \]
In the general 2-factor model, tests of interaction (A \(\times\) B) from the analysis of variance requires sufficient replication to allow for an error term with non-0 degrees of freedom. In the unreplicated case, the error term consumes all degrees of freedom available for interaction (or vice-verse).
summary(aov(assessment1 ~ treatment + replicate,data=arm.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 55 2387 43.4 0.875 0.712
## replicate 3 1809 603.0 12.162 3.13e-07 ***
## Residuals 165 8181 49.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(assessment1 ~ treatment + replicate + treatment:replicate, data=arm.dat))## Df Sum Sq Mean Sq
## treatment 55 2387 43.4
## replicate 3 1809 603.0
## treatment:replicate 165 8181 49.6summary(aov(assessment1 ~ treatment + replicate + Error(treatment:replicate), data=arm.dat))## Warning in aov(assessment1 ~ treatment + replicate +
## Error(treatment:replicate), : Error() model is singular##
## Error: treatment:replicate
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 55 2387 43.4 0.875 0.712
## replicate 3 1809 603.0 12.162 3.13e-07 ***
## Residuals 165 8181 49.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1However, a test of additivity of the two factors A and B can be determined without computing a SS component for the AB interaction.
First, find \(\mathbf{\alpha}\) and \(\mathbf{\beta}\) for the model \[ y_{ij} = \mu + \alpha _i + \beta _j + \epsilon _{ij} \]
We then find \[ \begin{aligned} \hat{\mu} ^{\star} & = & \hat{\mu} + \bar{\hat{\alpha}}_{.} + \bar{\hat{\beta}}_{.}\\ \hat{\alpha_i} ^{\star} & = & \hat{\alpha_i} - \bar{\hat{\alpha}}_{.}\\ \hat{\beta_i} ^{\star} & = & \hat{\beta_i} - \bar{\hat{\beta}}_{.} \end{aligned} \]
We get a 1 d.f. F-test of \(\hat{\alpha}\hat{\beta}\). We encapsulate this in a function:
tukey.1df <- function(data,
response = 'assessment1',
AName = 'treatment',
BName = 'replicate') {
#first column is variable name (A)
table <- tapply(get(response,data), list(get(AName,data),get(BName,data)),mean)
#fit crd model
ab.lm <- NULL
modelString <- paste(response,' ~ ',AName,' + ',BName)
ab.lm <- lm(as.formula(modelString),data=data)
#extract coefficients
a.levels <- levels(get(AName,data))
b.levels <- levels(get(BName,data))
a.max <- length(a.levels)
b.max <- length(b.levels)
m <- coef(ab.lm)[1]
a <- as.vector(c(0,coef(ab.lm)[2:(a.max)]))
b <- as.vector(c(0,coef(ab.lm)[(a.max+1):(a.max+b.max-1)]))
a.mean <- mean(a)
b.mean <- mean(b)
a.star <- a-a.mean
b.star <- b-b.mean
m.star = m + a.mean + b.mean
data$a.star <- a.star[get(AName,data)]
data$b.star <- b.star[get(BName,data)]
#analyze
modelString <- paste(response,' ~ ',AName,' + ',BName,'+ a.star: b.star')
ab.lm.2 <- lm(as.formula(modelString),data=data)
tbl <- anova(ab.lm.2)
return(list(F=tbl[3, 4], pF = tbl[3, 5], lm = ab.lm.2, a.star = a.star, b.star = b.star, m.star = m.star, table = table))
}The test for these data follows
tdf <- tukey.1df(data=arm.dat,response="assessment1")
tdf$F## [1] 3.305806tdf$pF## [1] 0.07086016In ARM Diagnostics, we only report the \(p(>F)\), which in this case is not signficant to \(\alpha=0.05\), but is larger than we might be comfortable with. If we examine the AOV
summary(aov(tdf$lm))## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 55 2387 43.4 0.888 0.6904
## replicate 3 1809 603.0 12.332 2.58e-07 ***
## a.star:b.star 1 162 161.7 3.306 0.0709 .
## Residuals 164 8019 48.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The residual sums of squares for this AOV is 48.9, while the SS associated with interaction is 161.7. This suggests that the residual error term in the origina (RCB) AOV may be too conservative (consider the analogy of using residual error instead of Replicate:A error in testing A means in a split-plot experiment). We might instead prefer a spatial analysis that might provide a better estimate of the confounding of treatment effects with spatial effects.
Note that the blocks are not uniformly shaped This may contribute to the interaction between blocks and treatments; we may not be sampling over the same spatial error in each replicate. See (Schabenberger and Pierce 2001) for more discussion.
summary(tdf$lm)##
## Call:
## lm(formula = as.formula(modelString), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.0668 -2.9097 -0.1218 3.7232 15.7267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.48549 3.58882 8.773 2.15e-15 ***
## treatment2 -3.36250 4.94464 -0.680 0.4974
## treatment3 -3.87500 4.94464 -0.784 0.4344
## treatment4 -7.78750 4.94464 -1.575 0.1172
## treatment5 0.86250 4.94464 0.174 0.8617
## treatment6 -1.37500 4.94464 -0.278 0.7813
## treatment7 -8.22500 4.94464 -1.663 0.0981 .
## treatment8 -2.43750 4.94464 -0.493 0.6227
## treatment9 -4.92500 4.94464 -0.996 0.3207
## treatment10 -1.80000 4.94464 -0.364 0.7163
## treatment11 -5.31250 4.94464 -1.074 0.2842
## treatment12 -0.87500 4.94464 -0.177 0.8598
## treatment13 -2.88750 4.94464 -0.584 0.5600
## treatment14 -2.05000 4.94464 -0.415 0.6790
## treatment15 -5.16250 4.94464 -1.044 0.2980
## treatment16 -6.75000 4.94464 -1.365 0.1741
## treatment17 -9.71250 4.94464 -1.964 0.0512 .
## treatment18 0.68750 4.94464 0.139 0.8896
## treatment19 -7.87500 4.94464 -1.593 0.1132
## treatment20 -8.91250 4.94464 -1.802 0.0733 .
## treatment21 -3.05000 4.94464 -0.617 0.5382
## treatment22 -7.71250 4.94464 -1.560 0.1207
## treatment23 -5.15000 4.94464 -1.042 0.2992
## treatment24 1.50000 4.94464 0.303 0.7620
## treatment25 3.21250 4.94464 0.650 0.5168
## treatment26 -5.65000 4.94464 -1.143 0.2548
## treatment27 -2.58750 4.94464 -0.523 0.6015
## treatment28 -7.42500 4.94464 -1.502 0.1351
## treatment29 -4.90000 4.94464 -0.991 0.3232
## treatment30 0.32500 4.94464 0.066 0.9477
## treatment31 -0.11250 4.94464 -0.023 0.9819
## treatment32 -7.90000 4.94464 -1.598 0.1120
## treatment33 -4.31250 4.94464 -0.872 0.3844
## treatment34 -3.13750 4.94464 -0.635 0.5266
## treatment35 -8.06250 4.94464 -1.631 0.1049
## treatment36 -1.76250 4.94464 -0.356 0.7220
## treatment37 -4.82500 4.94464 -0.976 0.3306
## treatment38 -5.52500 4.94464 -1.117 0.2655
## treatment39 -3.52500 4.94464 -0.713 0.4769
## treatment40 -9.02500 4.94464 -1.825 0.0698 .
## treatment41 -6.18750 4.94464 -1.251 0.2126
## treatment42 -2.62500 4.94464 -0.531 0.5962
## treatment43 -4.43750 4.94464 -0.897 0.3708
## treatment44 -7.63750 4.94464 -1.545 0.1244
## treatment45 -0.03750 4.94464 -0.008 0.9940
## treatment46 -3.75000 4.94464 -0.758 0.4493
## treatment47 1.82500 4.94464 0.369 0.7125
## treatment48 -6.21250 4.94464 -1.256 0.2108
## treatment49 -5.02500 4.94464 -1.016 0.3110
## treatment50 1.06250 4.94464 0.215 0.8301
## treatment51 -8.25000 4.94464 -1.668 0.0971 .
## treatment52 -1.91250 4.94464 -0.387 0.6994
## treatment53 0.67500 4.94464 0.137 0.8916
## treatment54 -1.03750 4.94464 -0.210 0.8341
## treatment55 -8.20000 4.94464 -1.658 0.0992 .
## treatment56 -5.83750 4.94464 -1.181 0.2395
## replicate2 1.05000 1.32151 0.795 0.4280
## replicate3 -2.99732 1.32151 -2.268 0.0246 *
## replicate4 -6.24464 1.32151 -4.725 4.91e-06 ***
## a.star:b.star -0.09156 0.05036 -1.818 0.0709 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.993 on 164 degrees of freedom
## Multiple R-squared: 0.3521, Adjusted R-squared: 0.119
## F-statistic: 1.511 on 59 and 164 DF, p-value: 0.0223For reference, we see a small negative correlation between a.star (treatment effect) and b.star (replicate effect). Perhaps varieties have different yield stability w.r.t. micro-environment?
Below are my notes to myself as I worked through Tukey’s example.
tukey.dat <- data.frame(
value=c(12,4,2,5,4,-2,-4,-5,4,-3,-7,-2),
squares=c(14,2,1,2,2,0,2,2,2,1,5,0),
row=as.factor(c(1,1,1,1,2,2,2,2,3,3,3,3)),
col=as.factor(c(1,2,3,4,1,2,3,4,1,2,3,4)))
#tukey.dat$squares <- round((tukey.dat$value*tukey.dat$value)/10)summary(aov(value~row+col,data=tukey.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## row 2 155.17 77.58 42.97 0.000278 ***
## col 3 156.67 52.22 28.92 0.000577 ***
## Residuals 6 10.83 1.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(value~col+row,data=tukey.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## col 3 156.67 52.22 28.92 0.000577 ***
## row 2 155.17 77.58 42.97 0.000278 ***
## Residuals 6 10.83 1.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(squares~row*col,data=tukey.dat))## Df Sum Sq Mean Sq
## row 2 24.50 12.25
## col 3 46.92 15.64
## row:col 6 84.83 14.14tapply(tukey.dat$value,list(tukey.dat$row),sum)## 1 2 3
## 23 -7 -8tapply(tukey.dat$value,list(tukey.dat$col),sum)## 1 2 3 4
## 20 -1 -9 -2row.means <- tapply(tukey.dat$value, list(tukey.dat$row),mean)
tukey.dat$row.means <- row.means[tukey.dat$row]
tukey.dat$row.means## [1] 5.75 5.75 5.75 5.75 -1.75 -1.75 -1.75 -1.75 -2.00 -2.00 -2.00
## [12] -2.00(tukey.dat$value-tukey.dat$row.means)^2## [1] 39.0625 3.0625 14.0625 0.5625 33.0625 0.0625 5.0625 10.5625
## [9] 36.0000 1.0000 25.0000 0.0000sum((tukey.dat$value-tukey.dat$row.means)^2)## [1] 167.5col.means <- tapply(tukey.dat$value, list(tukey.dat$col),mean)
tukey.dat$col.means <- col.means[tukey.dat$col]
tukey.dat$col.means## [1] 6.6666667 -0.3333333 -3.0000000 -0.6666667 6.6666667 -0.3333333
## [7] -3.0000000 -0.6666667 6.6666667 -0.3333333 -3.0000000 -0.6666667(tukey.dat$value-tukey.dat$col.means)^2## [1] 28.444444 18.777778 25.000000 32.111111 7.111111 2.777778 1.000000
## [8] 18.777778 7.111111 7.111111 16.000000 1.777778sum((tukey.dat$value-tukey.dat$col.means)^2)## [1] 166summary(aov(value~row,data=tukey.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## row 2 155.2 77.58 4.169 0.0523 .
## Residuals 9 167.5 18.61
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(value~col,data=tukey.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## col 3 156.7 52.22 2.517 0.132
## Residuals 8 166.0 20.75The first doesn’t match Tukey’s AOV, but the second does.
tdf <- tukey.1df(data=tukey.dat,AName="col",BName="row",response="squares")
tdf## $F
## [1] 7.898852
##
## $pF
## [1] 0.03752657
##
## $lm
##
## Call:
## lm(formula = as.formula(modelString), data = data)
##
## Coefficients:
## (Intercept) col2 col3 col4 row2
## 8.0000 -5.0000 -3.3333 -4.6667 -3.2500
## row3 a.star:b.star
## -2.7500 0.7364
##
##
## $a.star
## [1] 3.25000000 -1.75000000 -0.08333333 -1.41666667
##
## $b.star
## [1] 2.00 -1.25 -0.75
##
## $m.star
## (Intercept)
## 2.75
##
## $table
## 1 2 3
## 1 14 2 2
## 2 2 0 1
## 3 1 2 5
## 4 2 2 0Duplicating Tukey’s table 3
| . | Sums | Means | Deviations | Sums of x products | ||||
|---|---|---|---|---|---|---|---|---|
| . | 14 | 2 | 1 | 2 | 19 | 4.75 | 2.0 | 38.4 |
| . | 2 | 0 | 2 | 2 | 6 | 1.50 | -1.2 | 3.6 |
| . | 2 | 1 | 5 | 0 | 8 | 2.00 | -0.8 | 4.6 |
| Sums | 18 | 3 | 8 | 4 | 33 | . | 0.0 | 68.8 |
| Means | 6.00 | 1.00 | 2.67 | 1.33 | . | 2.75 | 6.08 | |
| Deviations | 3.2 | -1.8 | 0.0 | -1.4 | 0.0 | 15.44 | 50.9 |
These results don’t match exactly Tukey’s reported MS for non-additivity, but we can check that this may be due to rounding error
Tukey’s original
tukey.table <- matrix(tukey.dat$squares,nrow=3,byrow=TRUE)
a.deviations <- c(3.2,-1.8,0.0,-1.4)
b.deviations <- c(2.0,-1.2,-0.8)
sum(tukey.table[1,]*a.deviations)## [1] 38.4sum(tukey.table[2,]*a.deviations)## [1] 3.6sum(tukey.table[3,]*a.deviations)## [1] 4.6sum(c(sum(tukey.table[1,]*a.deviations),sum(tukey.table[2,]*a.deviations),sum(tukey.table[3,]*a.deviations))*b.deviations)## [1] 68.8sum(a.deviations*a.deviations)## [1] 15.44sum(b.deviations*b.deviations)## [1] 6.08(sum(a.deviations*a.deviations))*(sum(b.deviations*b.deviations))## [1] 93.8752(sum(c(sum(tukey.table[1,]*a.deviations),sum(tukey.table[2,]*a.deviations),sum(tukey.table[3,]*a.deviations))*b.deviations)^2)/((sum(a.deviations*a.deviations))*(sum(b.deviations*b.deviations)))## [1] 50.42269(68.8*68.8)/(15.44*6.08)## [1] 50.42269So, there’s a typo
a.deviations <- tdf$a.star
b.deviations <- tdf$b.star
sum(tukey.table[1,]*a.deviations)## [1] 39.08333sum(tukey.table[2,]*a.deviations)## [1] 3.5sum(tukey.table[3,]*a.deviations)## [1] 4.333333sum(c(sum(tukey.table[1,]*a.deviations),sum(tukey.table[2,]*a.deviations),sum(tukey.table[3,]*a.deviations))*b.deviations)## [1] 70.54167sum(a.deviations*a.deviations)## [1] 15.63889sum(b.deviations*b.deviations)## [1] 6.125(sum(a.deviations*a.deviations))*(sum(b.deviations*b.deviations))## [1] 95.78819(sum(c(sum(tukey.table[1,]*a.deviations),sum(tukey.table[2,]*a.deviations),sum(tukey.table[3,]*a.deviations))*b.deviations)^2)/((sum(a.deviations*a.deviations))*(sum(b.deviations*b.deviations)))## [1] 51.94927Littell, R C, G A Milliken, W W Stroup, and R D Wolfinger. 1996. SAS System for Mixed Models. SAS Institute, Cary, NC.
Schabenberger, Oliver, and Francis J. Pierce. 2001. Contemporary Statistical Models for the Plant and Soil Sciences [Hardcover]. CRC Press.
Tukey, John W. 1949. “One Degree of Freedom for Non-Additivity.” Biometrics 5 (3): 232–42.