Chapter 7 Unbalanced Two-Factor Analysis of Variance
When sample sizes are not all equal, the sums of squares cannot be obtained simply among the cell, marginal, and overall means. Even in well planned controlled experiments, observations may not be used due to malfunction or subjects dropping out. In observational studies, it may not be feasible to get equal number of units from the various sub-populations.
The model can be fit based on regression models in scalar and matrix form. In this chapter, we describe the process based on the treatment effects model. The balanced case can be formed in this manner as well.
7.1 Statistical Model and the Analysis of Variance
\[ Y_{ijk}=\mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\epsilon_{ijk}\qquad i=1,\ldots,a; \quad j=1,\ldots,b; \quad k=1,\ldots,n_{ij} \qquad \epsilon_{ijk} \sim NID\left(0,\sigma^2\right) \]
\[ n_{i\bullet} = \sum_{j=1}^b n_{ij} \quad n_{\bullet j} = \sum_{i=1}^a n_{ij} \quad n_T=\sum_{i=1}^a\sum_{j=1}^b n_{ij} \qquad Y_{ij\bullet}=\sum_{k=1}^{n_{ij}}Y_{ijk} \quad \overline{Y}_{ij\bullet} = \frac{Y_{ij\bullet}}{n_{ij}} \]
\[ \sum_{i=1}^a\alpha_i = \sum_{j=1}^b\beta_j = \sum_{i=1}^a(\alpha\beta)_{ij} = \sum_{j=1}^b(\alpha\beta)_{ij} = 0 \] \[ \alpha_a=-\sum_{i=1}^{a-1}\alpha_i \quad \beta_b = -\sum_{j=1}^{b-1}\beta_j \quad (\alpha\beta)_{aj}=-\sum_{i=1}^{a-1}(\alpha\beta)_{ij} \quad j=1,\ldots,b \quad (\alpha\beta)_{ib}=-\sum_{j=1}^{b-1}(\alpha\beta)_{ij} \quad i=1,\ldots,a \]
To fit this model, construct \(a-1\) \(X\) variables for the factor A effects, \(b-1\) \(X\) variables for the factor B effects, and then obtain the \((a-1)(b-1)\) cross-products of these variables for the interaction effects.
\[ X_{ijk1}^A=\left\{ \begin{array}{r@{\quad:\quad}l} 1 & i=1 \\ 0 & i=2,\ldots,a-1\\ -1 & i=a \end{array} \right. \qquad \cdots \qquad X_{ijk,a-1}^A=\left\{ \begin{array}{r@{\quad:\quad}l} 1 & i=a-1 \\ 0 & i=1,\ldots,a-2\\ -1 & i=a \end{array} \right. \]
\[ X_{ijk1}^B=\left\{ \begin{array}{r@{\quad:\quad}l} 1 & j=1 \\ 0 & j=2,\ldots,b-1\\ -1 & j=b \end{array} \right. \qquad \cdots \qquad X_{ijk,b-1}^B=\left\{ \begin{array}{r@{\quad:\quad}l} 1 & j=b-1 \\ 0 & j=1,\ldots,b-2\\ -1 & j=b \end{array} \right. \] \[ Y_{ijk}=\mu_{\bullet\bullet} + \alpha_1X_{ijk1}^A + \cdots +\alpha_{a-1}X_{ijk,a-1}^A + \beta_1X_{ijk1}^B + \cdots + \beta_{b-1}X_{ijk,b-1}^B + \] \[ + (\alpha\beta)_{11}X_{ijk1}^AX_{ijk1}^B + \cdots + (\alpha\beta)_{a-1,b-1}X_{ijk,a-1}^AX_{ijk,b-1}^B + \epsilon_{ijk} \]
The testing strategy is done as it was in the balanced case. Fit various models and compare them using the general linear test procedure.
- Model 1 - Include main effects for factors A and B as well as interaction effects.
- Model 2 - Include main effects for factors A and B, but no interaction effects
- Model 3 - Include main effects for factor B and interaction effects, but no factor A effects
- Model 4 - Include main effects for factor A and interaction effects, but no factor B effects
To test for interaction effects, controlling for main effects for factors A and B, we compare Model 1 (Complete) and Model 2 (Reduced).
To test for factor A effects, controlling for factor B effects and interaction effects, we compare Model 1 (Complete) and Model 3 (Reduced).
To test for factor B effects, controlling for factor A effects and interaction effects, we compare Model 1 (Complete) and Model 4 (Reduced).
When the test for interaction is “very non-significant” (high \(P\)-value), often tests for factor A and B effects use Model 2 as the Complete Model.
7.1.1 Example - Age at Peak for Famous Writers
A study considered the age at career peak for famous writers (Simonson 2007). The factors were Style (Factor A: Conceptualists \((i=1)\), Experimentalists \((i=2)\)) and Writer Type (Factor B:Poets \((j=1)\), Novelists \((j=2)\)). Thus, there are \(a=b=2\) levels for the factors, there will only be \(a-1=1\) \(X\) variable for factor A and \(b-1=1\) \(X\) variable for factor B, and \((a-1)(b-1)=1\) cross-product term for the AB interaction. The numbers of writers for the \(ab=4\) conditions are \(n_{11}=n_{12}=5,n_{21}=6,n_{22}=7\). The data and the generated X variables are given in the following table. The regression model is fit directly, then with the aov function.
| Writer | Style | Type | yearPeak | X1A | X1B | X1AX1B |
|---|---|---|---|---|---|---|
| Eliot | 1 | 1 | 23 | 1 | 1 | 1 |
| Cummings | 1 | 1 | 26 | 1 | 1 | 1 |
| Plath | 1 | 1 | 30 | 1 | 1 | 1 |
| Pound | 1 | 1 | 30 | 1 | 1 | 1 |
| Wilbur | 1 | 1 | 34 | 1 | 1 | 1 |
| Fitzgerald | 1 | 2 | 29 | 1 | -1 | -1 |
| Hemingway | 1 | 2 | 30 | 1 | -1 | -1 |
| Melville | 1 | 2 | 32 | 1 | -1 | -1 |
| Lawrence | 1 | 2 | 35 | 1 | -1 | -1 |
| Joyce | 1 | 2 | 40 | 1 | -1 | -1 |
| Bishop | 2 | 1 | 29 | -1 | 1 | -1 |
| Moore | 2 | 1 | 32 | -1 | 1 | -1 |
| Williams | 2 | 1 | 40 | -1 | 1 | -1 |
| Lowell | 2 | 1 | 41 | -1 | 1 | -1 |
| Stevens | 2 | 1 | 42 | -1 | 1 | -1 |
| Frost | 2 | 1 | 48 | -1 | 1 | -1 |
| James | 2 | 2 | 38 | -1 | -1 | 1 |
| Faulkner | 2 | 2 | 39 | -1 | -1 | 1 |
| Dickens | 2 | 2 | 41 | -1 | -1 | 1 |
| Woolf | 2 | 2 | 45 | -1 | -1 | 1 |
| Conrad | 2 | 2 | 47 | -1 | -1 | 1 |
| Twain | 2 | 2 | 50 | -1 | -1 | 1 |
| Hardy | 2 | 2 | 51 | -1 | -1 | 1 |
##
## Call:
## lm(formula = yearPeak ~ X1A + X1B + X1AX1B, data = ap1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.667 -3.814 1.333 2.952 9.333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.2238 1.1402 31.769 < 2e-16 ***
## X1A -5.3238 1.1402 -4.669 0.000167 ***
## X1B -2.5905 1.1402 -2.272 0.034902 *
## X1AX1B 0.2905 1.1402 0.255 0.801652
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.415 on 19 degrees of freedom
## Multiple R-squared: 0.5978, Adjusted R-squared: 0.5343
## F-statistic: 9.413 on 3 and 19 DF, p-value: 0.0005033## ANOVA model (Treatment effects)
options(contrasts=c("contr.sum", "contr.poly"))
ap.aov1 <- aov(yearPeak ~ factor(Style) + factor(Type) + factor(Style)*factor(Type), data=ap1)
summary.lm(ap.aov1)##
## Call:
## aov(formula = yearPeak ~ factor(Style) + factor(Type) + factor(Style) *
## factor(Type), data = ap1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.667 -3.814 1.333 2.952 9.333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.2238 1.1402 31.769 < 2e-16 ***
## factor(Style)1 -5.3238 1.1402 -4.669 0.000167 ***
## factor(Type)1 -2.5905 1.1402 -2.272 0.034902 *
## factor(Style)1:factor(Type)1 0.2905 1.1402 0.255 0.801652
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.415 on 19 degrees of freedom
## Multiple R-squared: 0.5978, Adjusted R-squared: 0.5343
## F-statistic: 9.413 on 3 and 19 DF, p-value: 0.0005033## Analysis of Variance Table
##
## Response: yearPeak
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Style) 1 667.75 667.75 22.7758 0.0001324 ***
## factor(Type) 1 158.26 158.26 5.3979 0.0314109 *
## factor(Style):factor(Type) 1 1.90 1.90 0.0649 0.8016516
## Residuals 19 557.05 29.32
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We see from the regression coefficients (that control for all other factors) that there is no evidence of an interaction, but there is evidence of main effects for Style and Writer Type. Formal tests are given below. Given that each test has 1 degree of freedom, we use \(t\)-tests.
\[ H_0^{AB}: (\alpha\beta)_{11}=(\alpha\beta)_{12}=(\alpha\beta)_{21}=(\alpha\beta)_{22}=0 \qquad t_{AB}^*=\frac{0.2905}{1.1402}=0.255 \qquad P_{AB}=.8017 \] \[ H_0^A:\alpha_1=\alpha_2=0 \qquad t_A^*=\frac{-5.3238}{1.1402}=-4.669 \qquad P_A=.0002 \] \[ H_0^B: \beta_1=\beta_2=0 \qquad t_B^*=\frac{-2.5905}{1.1402}=-2.272 \qquad P_B=.0349 \] As there are strong evidence of main effects and no suggestion of an interaction, we will fit Model 2 for making comparisons, but not fit Models 3 or 4.
##
## Call:
## lm(formula = yearPeak ~ X1A + X1B, data = ap1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.940 -4.027 1.060 2.933 9.060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.234 1.113 32.566 < 2e-16 ***
## X1A -5.334 1.113 -4.794 0.000111 ***
## X1B -2.628 1.104 -2.380 0.027397 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.287 on 20 degrees of freedom
## Multiple R-squared: 0.5964, Adjusted R-squared: 0.5561
## F-statistic: 14.78 on 2 and 20 DF, p-value: 0.0001146##
## Call:
## aov(formula = yearPeak ~ factor(Style) + factor(Type), data = ap1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.940 -4.027 1.060 2.933 9.060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.234 1.113 32.566 < 2e-16 ***
## factor(Style)1 -5.334 1.113 -4.794 0.000111 ***
## factor(Type)1 -2.628 1.104 -2.380 0.027397 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.287 on 20 degrees of freedom
## Multiple R-squared: 0.5964, Adjusted R-squared: 0.5561
## F-statistic: 14.78 on 2 and 20 DF, p-value: 0.0001146## Analysis of Variance Table
##
## Response: yearPeak
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Style) 1 667.75 667.75 23.8930 8.894e-05 ***
## factor(Type) 1 158.26 158.26 5.6627 0.0274 *
## Residuals 20 558.95 27.95
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Based on the additive model (Model 2), we obtain the following estimates and predicted ages at peak.
\[ \hat{\mu}_{\bullet\bullet}=36.234 \qquad \hat{\alpha}_1=-5.334 \quad \hat{\alpha}_2=5.334 \qquad \hat{\beta}_1=-2.628 \quad \hat{\beta}_2=2.628 \] \[ \mbox{Conceptualists/Poets:} \quad \hat{Y}_{11}=\hat{\mu}_{\bullet\bullet}+\hat{\alpha}_1+\hat{\beta}_1= 36.234-5.334-2.628=28.272 \]
\[ \mbox{Conceptualists/Novelists:} \quad \hat{Y}_{12}=\hat{\mu}_{\bullet\bullet}+\hat{\alpha}_1+\hat{\beta}_2= 36.234-5.334+2.628=33.528 \]
\[ \mbox{Experimantalists/Poets:} \quad \hat{Y}_{21}=\hat{\mu}_{\bullet\bullet}+\hat{\alpha}_2+\hat{\beta}_1= 36.234+5.334-2.628=38.940 \]
\[ \mbox{Experimantalists/Novelists:} \quad \hat{Y}_{21}=\hat{\mu}_{\bullet\bullet}+\hat{\alpha}_2+\hat{\beta}_2= 36.234+5.334+2.628=44.196 \]
Experimentalists take longer on average to reach their peak than conceptualists, and novelists take longer than poets.
\[\nabla \]
A second example with more levels of factor B is included here.
7.1.2 Example - Makiwara Punching Boards
An experiment was conducted to compare two types of karate boards, made from 4 wood types (Smith, Niiler, and McCullough 2010). The board types were stacked \((i=1)\) and tapered \((i=2)\), the wood types were: cherry \((j=1)\), ash \((j=2)\), fir \((j=3)\), and oak \((j=4)\). Apparently due to breakage, there were unequal numbers of replicates among the \(ab=2(4)=8\) treatments. We use the aov function in R with the treatment effects form to fit Models 1 (Interaction) and 2 (Additive) to test whether there is an interaction between board and wood types. The response \(Y\) is the deflection (millimeters) for the Makiwara board, with a total of \(n_T=336\) measurements.
mb1 <- read.table("https://www.stat.ufl.edu/~winner/data/karate_board.dat",
col.names=c("trt.id", "wood", "board", "id", "deflect"))
head(mb1)## trt.id wood board id deflect
## 1 1 1 1 1 144.3
## 2 1 1 1 2 125.9
## 3 1 1 1 3 263.2
## 4 1 1 1 4 114.6
## 5 1 1 1 5 242.5
## 6 1 1 1 6 141.9## trt.id wood board id deflect
## 331 8 4 2 40 56.6
## 332 8 4 2 41 123.5
## 333 8 4 2 42 12.0
## 334 8 4 2 43 62.0
## 335 8 4 2 44 73.3
## 336 8 4 2 45 44.9mb1$wood.f <- factor(mb1$wood, levels=1:4, labels=c("cherry", "ash", "fir", "oak"))
mb1$board.f <- factor(mb1$board, levels=1:2, labels=c("stacked", "tapered"))
tapply(mb1$deflect, list(mb1$board.f,mb1$wood.f), mean)## cherry ash fir oak
## stacked 118.00000 105.9944 93.99706 89.99556
## tapered 77.99556 76.0000 55.00000 49.99333##
## Call:
## aov(formula = deflect ~ board.f * wood.f, data = mb1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -97.600 -35.775 -7.246 25.705 281.104
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 83.3720 3.0838 27.035 < 2e-16 ***
## board.f1 18.6248 3.0838 6.039 4.19e-09 ***
## wood.f1 14.6258 5.2834 2.768 0.00596 **
## wood.f2 7.6252 5.4085 1.410 0.15953
## wood.f3 -8.8735 5.4678 -1.623 0.10558
## board.f1:wood.f1 1.3775 5.2834 0.261 0.79448
## board.f1:wood.f2 -3.6275 5.4085 -0.671 0.50288
## board.f1:wood.f3 0.8738 5.4678 0.160 0.87314
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 56.2 on 328 degrees of freedom
## Multiple R-squared: 0.1359, Adjusted R-squared: 0.1174
## F-statistic: 7.366 on 7 and 328 DF, p-value: 3.178e-08## Analysis of Variance Table
##
## Response: deflect
## Df Sum Sq Mean Sq F value Pr(>F)
## board.f 1 115479 115479 36.5595 4.03e-09 ***
## wood.f 3 45955 15318 4.8496 0.002575 **
## board.f:wood.f 3 1444 481 0.1524 0.928115
## Residuals 328 1036045 3159
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Call:
## aov(formula = deflect ~ board.f + wood.f, data = mb1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -96.220 -35.852 -7.682 26.736 282.422
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 83.431 3.067 27.205 < 2e-16 ***
## board.f1 18.683 3.067 6.092 3.09e-09 ***
## wood.f1 14.506 5.252 2.762 0.00607 **
## wood.f2 7.976 5.359 1.488 0.13757
## wood.f3 -9.046 5.407 -1.673 0.09525 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 55.99 on 331 degrees of freedom
## Multiple R-squared: 0.1346, Adjusted R-squared: 0.1242
## F-statistic: 12.88 on 4 and 331 DF, p-value: 9.383e-10## Analysis of Variance Table
##
## Response: deflect
## Df Sum Sq Mean Sq F value Pr(>F)
## board.f 1 115479 115479 36.8425 3.506e-09 ***
## wood.f 3 45955 15318 4.8872 0.002445 **
## Residuals 331 1037489 3134
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1When comparing the mean deflection differences (stacked-tapered) within the 4 wood types, based on the tapply function, we observe: 40 (cherry), 30 (ash), 39 (fir), and 41 (oak). While these differences are not exactly the same, they are very consistent. To test for interaction, we directly use the ANOVA table from Model 1.
\[ H_0: (\alpha\beta)_{11}= \cdots = (\alpha\beta)_{42}=0 \qquad H_A:\mbox{ Not all } (\alpha\beta)_{ij}= 0 \]
\[ \mbox{Test Statistic:} \quad F_{AB}^*=\frac{MSAB}{MSE}=\frac{481}{3159}=0.1524 \qquad \mbox{Rejection Region:}\quad F_{AB}^*\ge F_{.95;3,328}=2.632 \]
The \(P\)-value is large, \(P_{AB}=.9281\), giving no evidence of any interaction effects.
Based on this, we test for main effects for Board and Wood Effects, based on Model 2, the additive model.
For Board effects, \(H_0^A:\alpha_1=\alpha_2=0\), we obtain the following test.
\[ \mbox{Test Statistic:} \quad F_{A}^*=\frac{MSA}{MSE}=\frac{115479}{3134}=36.8425 \qquad \mbox{Rejection Region:}\quad F_{A}^*\ge F_{.95;1,331}=3.870 \] For Wood effects, \(H_0^B:\beta_1=\cdots = \beta_4=0\), we obtain the following test.
\[ \mbox{Test Statistic:} \quad F_{B}^*=\frac{MSB}{MSE}=\frac{15318}{3134}=4.8872 \qquad \mbox{Rejection Region:}\quad F_{B}^*\ge F_{.95;3,331}=2.632 \] Both results are significant, with \(P_A<.0001\) and \(P_B=.0024\).
\[\nabla\]
7.2 Least Squares Estimators and Contrasts/Linear Functions Among Means
In this section, the formulas for estimators and estimated standard errors for means and linear functions are given. For treatment (cell) means, we have the following results.
\[ \mbox{Parameter:} \quad \mu_{ij} \qquad \mbox{Estimator:} \quad \hat{\mu}_{ij} = \frac{\sum_{k=1}^{n_{ij}}Y_{ijk}}{n_{ij}}= \overline{Y}_{ij\bullet} \] \[ \mbox{Estimated Standard Error:} \quad s\left\{\hat{\mu}_{ij}\right\}=\sqrt{\frac{MSE}{n_{ij}}} \] For factor A level means, we have the following results. \[ \mbox{Parameter:} \quad \mu_{i\bullet} = \frac{\sum_{j=1}^b\mu_{ij}}{b} \qquad \mbox{Estimator:} \quad \hat{\mu}_{i\bullet} = \frac{\sum_{j=1}^b\overline{Y}_{ij\bullet}}{b} \] \[ \mbox{Estimated Standard Error:} \quad s\left\{\hat{\mu}_{i\bullet}\right\}= \sqrt{\frac{MSE}{b^2}\sum_{j=1}^b\frac{1}{n_{ij}}} \] For factor B level means, we have the following results. \[ \mbox{Parameter:} \quad \mu_{\bullet j} = \frac{\sum_{i=1}^a\mu_{ij}}{a} \qquad \mbox{Estimator:} \quad \hat{\mu}_{\bullet j} = \frac{\sum_{i=1}^a\overline{Y}_{ij\bullet}}{a} \] \[ \mbox{Estimated Standard Error:} \quad s\left\{\hat{\mu}_{\bullet j}\right\}= \sqrt{\frac{MSE}{a^2}\sum_{i=1}^a\frac{1}{n_{ij}}} \]
For contrasts or linear functions among factor A and B levels, we have the following parameters, estimators, and estimated standard errors.
\[ L_A = \sum_{i=1}^a c_i\mu_{i\bullet} \qquad \hat{L}_A = \sum_{i=1}^a c_i\hat{\mu}_{i\bullet} \qquad s\left\{\hat{L}_A\right\} = \sqrt{\frac{MSE}{b^2}\sum_{i=1}^ac_i^2\sum_{j=1}^b\frac{1}{n_{ij}}} \]
\[ L_B = \sum_{j=1}^b c_j\mu_{\bullet j} \qquad \hat{L}_B = \sum_{j=1}^b c_j\hat{\mu}_{\bullet j} \qquad s\left\{\hat{L}_B\right\} = \sqrt{\frac{MSE}{a^2}\sum_{j=1}^bc_j^2\sum_{i=1}^a\frac{1}{n_{ij}}} \] Finally, the results are given for contrasts and linear functions among treatment (cell) means.
\[ L_{AB} = \sum_{i=1}^a \sum_{j=1}^b c_{ij}\mu_{ij} \qquad \hat{L}_{AB} = \sum_{i=1}^a\sum_{j=1}^b c_{ij}\hat{\mu}_{ij} \qquad s\left\{\hat{L}_{AB}\right\} = \sqrt{MSE\sum_{i=1}^a\sum_{j=1}^b\frac{c_{ij}^2}{n_{ij}}} \]
The standard error multipliers for contrasts/linear functions based on the Scheffe (all comparisons), Bonferroni (\(g\) pre-planned comparisons), and Tukey (all pairwise comparisons) methods are given here.
For single comparisons, use \(t_{1-\alpha/2;n_T-ab}\).
For comparisons among treatment (cell) means, use the following multipliers.
\[ \mbox{Scheffe:} \quad S=\sqrt{(ab-1)F_{1-\alpha;ab-1,n_T-ab}} \] \[ \mbox{Bonferroni:} \quad B=t_{1-\alpha/(2g); n_T-ab} \] \[ \mbox{Tukey:} \quad T=\frac{1}{\sqrt{2}}q_{1-\alpha;ab,n_T-ab} \]
Standard error multipliers for contrasts and linear functions among levels of factors A and B are given here.
\[ \mbox{Scheffe:} \qquad \mbox{A:} \quad S_A=\sqrt{(a-1)F_{1-\alpha;a-1,n_T-ab}} \qquad \quad \mbox{B:} \quad S_B=\sqrt{(b-1)F_{1-\alpha;b-1,n_T-ab}} \] \[ \mbox{Bonferroni (Factors A and B):} \quad B=t_{1-\alpha/(2g); n_T-ab} \] \[ \mbox{Tukey:} \qquad \mbox{A:} \quad T_A=\frac{1}{\sqrt{2}}q_{1-\alpha;a,n_T-ab} \qquad \mbox{B:} \quad T_B=\frac{1}{\sqrt{2}}q_{1-\alpha;b,n_T-ab} \]
7.2.1 Example - Age at Peak for Famous Writers
Although the interaction was not significant for the Age at Peak of the writers, we will use the Tukey method to compare all pairs of Style/Type means. The Mean Square Error for the interaction model was \(MSE=29.32\) with degrees of freedom \(n_T-ab=23-2(2)=19\).
\[ \overline{y}_{11\bullet} = 28.60 \qquad n_{11}=5 \qquad s\left\{\overline{Y}_{11\bullet}\right\}=\sqrt{\frac{MSE}{n_{11}}}=\sqrt{\frac{29.32}{5}}=2.42 \] \[ \overline{y}_{12\bullet} = 33.20 \qquad n_{12}=5 \qquad s\left\{\overline{Y}_{12\bullet}\right\}=\sqrt{\frac{MSE}{n_{12}}}=\sqrt{\frac{29.32}{5}}=2.42 \] \[ \overline{y}_{21\bullet} = 38.67 \qquad n_{21}=6 \qquad s\left\{\overline{Y}_{21\bullet}\right\}=\sqrt{\frac{MSE}{n_{21}}}=\sqrt{\frac{29.32}{6}}=2.21 \] \[ \overline{y}_{22\bullet} = 44.43 \qquad n_{22}=7 \qquad s\left\{\overline{Y}_{22\bullet}\right\}=\sqrt{\frac{MSE}{n_{22}}}=\sqrt{\frac{29.32}{7}}=2.05 \] \[ T = \frac{1}{\sqrt{2}}q_{.95;4,19}=\frac{1}{\sqrt{2}}(3.977)=2.812 \qquad s\left\{\overline{Y}_{ij\bullet}-\overline{Y}_{i'j'\bullet}\right\}= \sqrt{MSE\left(\frac{1}{n_{ij}}+\frac{1}{n_{i'j'}}\right)} \] \[ \overline{y}_{11\bullet}-\overline{y}_{12\bullet} = 28.60-33.20=-4.60 \qquad s\left\{\overline{Y}_{11\bullet}-\overline{Y}_{12\bullet}\right\} = \sqrt{29.32\left(\frac{1}{5}+\frac{1}{5}\right)}=3.43 \] \[ \overline{y}_{11\bullet}-\overline{y}_{21\bullet} = 28.60-38.67=-10.07 \qquad s\left\{\overline{Y}_{11\bullet}-\overline{Y}_{21\bullet}\right\} = \sqrt{29.32\left(\frac{1}{5}+\frac{1}{6}\right)}=3.28 \] \[ \overline{y}_{11\bullet}-\overline{y}_{22\bullet} = 28.60-44.43=-15.83 \qquad s\left\{\overline{Y}_{11\bullet}-\overline{Y}_{22\bullet}\right\} = \sqrt{29.32\left(\frac{1}{5}+\frac{1}{7}\right)}=3.17 \] \[ \overline{y}_{12\bullet}-\overline{y}_{21\bullet} = 33.20-38.67=-5.47 \qquad s\left\{\overline{Y}_{12\bullet}-\overline{Y}_{21\bullet}\right\} = \sqrt{29.32\left(\frac{1}{5}+\frac{1}{6}\right)}=3.28 \] \[ \overline{y}_{12\bullet}-\overline{y}_{22\bullet} = 33.20-44.43=-11.23 \qquad s\left\{\overline{Y}_{12\bullet}-\overline{Y}_{22\bullet}\right\} = \sqrt{29.32\left(\frac{1}{5}+\frac{1}{7}\right)}=3.17 \] \[ \overline{y}_{21\bullet}-\overline{y}_{22\bullet} = 38.67-44.43=-5.76 \qquad s\left\{\overline{Y}_{21\bullet}-\overline{Y}_{22\bullet}\right\} = \sqrt{29.32\left(\frac{1}{6}+\frac{1}{7}\right)}=3.01 \]
\[ \mbox{Conc/Poet-Conc/Novel:} \quad HSD=2.812(3.43)=9.65 \qquad \overline{y}_{11\bullet}-\overline{y}_{12\bullet} =-4.60 \] \[ \mbox{Conc/Poet-Exp/Poet:} \quad HSD=2.812(3.28)=9.22 \qquad \overline{y}_{11\bullet}-\overline{y}_{21\bullet} = 28.60-38.67=-10.07 \] \[ \mbox{Conc/Poet-Exp/Novel:} \quad HSD=2.812(3.17)=8.92 \qquad \overline{y}_{11\bullet}-\overline{y}_{22\bullet} =-15.83 \] \[ \mbox{Conc/Nov-Exp/Poet:} \quad HSD=2.812(3.28)=9.22 \qquad \overline{y}_{12\bullet}-\overline{y}_{21\bullet} = -5.47 \] \[ \mbox{Conc/Nov-Exp/Nov:} \quad HSD=2.812(3.17)=8.92 \qquad \overline{y}_{12\bullet}-\overline{y}_{22\bullet} = -11.23 \] \[ \mbox{Exp/Poet-Exp/Nov:} \quad HSD=2.812(3.01)=8.47 \qquad \overline{y}_{21\bullet}-\overline{y}_{22\bullet} = -5.76 \]
The groups that are significantly different, based on Tukey’s HSD are as follow.
- Conceptualist/Poet is significantly lower than Experimentalist/Poet
- Conceptualist/Poet is significantly lower than Experimentalist/Novelist
- Conceptualist/Novelist is significantly lower than Experimentalist/Novelist
The following R program runs Tukey’s HSD on the marginal means for factors A and B, and on the cell means as well.
ap1 <- read.csv("author_peak1.csv")
options(contrasts=c("contr.sum", "contr.poly"))
ap.aov1 <- aov(yearPeak ~ factor(Style) + factor(Type) + factor(Style)*factor(Type), data=ap1)
TukeyHSD(ap.aov1)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = yearPeak ~ factor(Style) + factor(Type) + factor(Style) * factor(Type), data = ap1)
##
## $`factor(Style)`
## diff lwr upr p adj
## 2-1 10.86923 6.102333 15.63613 0.0001324
##
## $`factor(Type)`
## diff lwr upr p adj
## 2-1 5.247378 0.5167309 9.978024 0.0315205
##
## $`factor(Style):factor(Type)`
## diff lwr upr p adj
## 2:1-1:1 10.066667 0.8473951 19.285938 0.0293965
## 1:2-1:1 4.600000 -5.0292152 14.229215 0.5481270
## 2:2-1:1 15.828571 6.9136505 24.743492 0.0004327
## 1:2-2:1 -5.466667 -14.6859383 3.752605 0.3672669
## 2:2-2:1 5.761905 -2.7085734 14.232383 0.2559999
## 2:2-1:2 11.228571 2.3136505 20.143492 0.0107243\[\nabla\]
7.2.2 Example - Makiwara Punching Boards
For the Karate board experiment, we will conduct tests among levels of factors A and B, respectively. The Mean Square Error was \(MSE=3159\) with degrees of freedom \(n_T-ab=336-2(4)=328\) for Model 1 (with interaction).
For comparing the Board Types (Stacked-Tapered), we obtain the following cell means from the tapply function previously.
\[ \overline{y}_{11\bullet}=118 \quad \overline{y}_{12\bullet}=106 \quad \overline{y}_{13\bullet}=94 \quad \overline{y}_{14\bullet}=90 \quad \Rightarrow \quad \hat{\mu}_{1\bullet}=\frac{118+106+94+90}{4}=102 \] \[ \overline{y}_{21\bullet}=78 \quad \overline{y}_{22\bullet}=76 \quad \overline{y}_{23\bullet}=55 \quad \overline{y}_{24\bullet}=50 \quad \Rightarrow \quad \hat{\mu}_{2\bullet}=\frac{78+76+55+50}{4}=64.75 \]
The standard errors are computed here, where the sample sizes for each cell are given first.
\[ n_{11}=41 \quad n_{12}=36 \quad n_{13}=34 \quad n_{14}=n_{21}=n_{22}=n_{23}=n_{24}=45 \] \[ s\left\{\hat{\mu}_{1\bullet}\right\}= \sqrt{\frac{3159}{4^2}\left(\frac{1}{41}+\frac{1}{36}+\frac{1}{34}+\frac{1}{45}\right)} =4.53 \] \[ s\left\{\hat{\mu}_{2\bullet}\right\}= \sqrt{\frac{3159}{4^2}\left(\frac{1}{45}+\frac{1}{45}+\frac{1}{45}+\frac{1}{45}\right)} =4.19 \] The standard error of the difference \(\hat{\mu}_{1\bullet}-\hat{\mu}_{2\bullet}\) is the square root of the sum of their variances, as these are independent samples. Also, we obtain the critical \(t\)-value, and the 95% Confidence Interval for \(\mu_{1\bullet}-\mu_{2\bullet}\).
\[ s\left\{\hat{\mu}_{1\bullet}-\hat{\mu}_{2\bullet}\right\} = \sqrt{(4.53)^2+(4.19)^2}=6.17 \qquad t_{.975;328}=1.967 \] \[ \mbox{95% Confidence Interval:} \quad (102-64.75) \pm 1.967(6.17) \quad \equiv \quad 37.25 \pm 12.14 \quad \equiv \quad (25.11,49.39) \] Thus, the mean deflection is higher for the stacked board than the tapered on average (likely by somewhere between 25-50 millimeters). We could also compare all pairs of Wood Types. The computations are similar, but there are 4(3)/2 = 6 possible pairs. We will use the TukeyHSD function along with the aov object in R below.
mb1 <- read.table("https://www.stat.ufl.edu/~winner/data/karate_board.dat",
col.names=c("trt.id", "wood", "board", "id", "deflect"))
head(mb1)## trt.id wood board id deflect
## 1 1 1 1 1 144.3
## 2 1 1 1 2 125.9
## 3 1 1 1 3 263.2
## 4 1 1 1 4 114.6
## 5 1 1 1 5 242.5
## 6 1 1 1 6 141.9## trt.id wood board id deflect
## 331 8 4 2 40 56.6
## 332 8 4 2 41 123.5
## 333 8 4 2 42 12.0
## 334 8 4 2 43 62.0
## 335 8 4 2 44 73.3
## 336 8 4 2 45 44.9mb1$wood.f <- factor(mb1$wood, levels=1:4, labels=c("cherry", "ash", "fir", "oak"))
mb1$board.f <- factor(mb1$board, levels=1:2, labels=c("stacked", "tapered"))
tapply(mb1$deflect, list(mb1$board.f,mb1$wood.f), mean)## cherry ash fir oak
## stacked 118.00000 105.9944 93.99706 89.99556
## tapered 77.99556 76.0000 55.00000 49.99333## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = deflect ~ board.f * wood.f, data = mb1)
##
## $board.f
## diff lwr upr p adj
## tapered-stacked -37.17265 -49.26685 -25.07845 0## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = deflect ~ board.f * wood.f, data = mb1)
##
## $wood.f
## diff lwr upr p adj
## ash-cherry -6.535911 -29.00664 15.9348190 0.8761668
## fir-cherry -23.560409 -46.17714 -0.9436753 0.0375052
## oak-cherry -27.937478 -49.82197 -6.0529827 0.0059563
## fir-ash -17.024498 -39.97301 5.9240184 0.2234297
## oak-ash -21.401567 -43.62878 0.8256437 0.0639197
## oak-fir -4.377069 -26.75187 17.9977342 0.9578232In terms of the Wood Types, fir has significantly lower deflection than cherry, and oak has significantly lower deflection than cherry.
\[\nabla\]