Let’s think more about the two-way ANOVA. The most general model is:

\[y_{ij} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ij}\]

for two factors, \(\alpha\) with \(I\) levels and \(\beta\) with \(J\) levels.

• Note: not all parameters are identifiable.

• Interaction \((\alpha\beta)_{ij}\) is the part of the mean response not attributable to the additive effects of the factors.
  • Ex: you might not want to eat straight sugar or drink heavy cream, but whipped cream (a combination of the two) is much more pleasant.
• Interactions make models harder to interpret, since we can’t examine \(\alpha\) and \(\beta\) independently of each other.

Consider the following 2x2 examples

The response is some measure of drug performance.

• Left has no interaction, right does.

• When there’s an interaction, not a clear way to define the individual effects.
  • Is the gender effect the effect for F? For M? For their average?
• We can also create a model \(y_{ijk}=\mu_{ijk} + \epsilon_{ijk}\), treating it as a one-way ANOVA with \(I\times J\) total levels.

data("composite")
composite

##   strength laser   tape
## 1    25.66   40W   slow
## 2    29.15   50W   slow
## 3    35.73   60W   slow
## 4    28.00   40W medium
## 5    35.09   50W medium
## 6    39.56   60W medium
## 7    20.65   40W   fast
## 8    29.79   50W   fast
## 9    35.66   60W   fast

• One observation for each combination of laser and tape.

We can fit a main effects model:

mod1 <- lm(strength ~ laser + tape, data=composite)

but this loses sight of possible interactions.

However, if we try to fit an interaction model, parameters = observations

par(mfrow=c(1,2))
with (composite, interaction.plot (laser, tape, strength, legend=F))
with (composite, interaction.plot (tape, laser, strength, legend=F))

Those don’t look like significant interactions, but here’s another way to test. We fit the model

\[y_{ij} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ij}\]

and then test if the interaction term is zero.

• Note that this does assume the interaction is multiplicative, which may not always be the case.

coef (mod1)

## (Intercept)    laser50W    laser60W  tapemedium    tapefast 
##   23.917778    6.573333   12.213333    4.036667   -1.480000

lasercoefs <- rep(c(0,6.5733, 12.2133),3)
tapecoefs <- rep(c(0,4.0367,-1.4800), each=3)
h <- update(mod1, . ~ . + I(lasercoefs*tapecoefs))

anova(h)

## Analysis of Variance Table
## 
## Response: strength
##                           Df  Sum Sq Mean Sq F value   Pr(>F)   
## laser                      2 224.184 112.092 36.8201 0.007745 **
## tape                       2  48.919  24.459  8.0344 0.062401 . 
## I(lasercoefs * tapecoefs)  1   1.370   1.370  0.4501 0.550340   
## Residuals                  3   9.133   3.044                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Looks like the interaction is not significant

anova(mod1)

## Analysis of Variance Table
## 
## Response: strength
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## laser      2 224.184 112.092 42.6893 0.002003 **
## tape       2  48.919  24.459  9.3151 0.031242 * 
## Residuals  4  10.503   2.626                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can then look at the plots and coefficients to see how strength changes with laser power and tape speed.

One thing we’re not taking into account is that these factors are ordered.

composite$laser <- as.ordered(composite$laser)
composite$tape <- as.ordered(composite$tape)
g <- lm(strength ~ laser + tape, composite)
summary(g) #note the anova itself will look the same

Estimate Std. Error t value Pr(>|t|)

(Intercept) 31.0322 0.5401 57.452 5.5e-07 ***

laser.L 8.6361 0.9355 9.231 0.000765 ***

laser.Q -0.3810 0.9355 -0.407 0.704655

tape.L -1.0465 0.9355 -1.119 0.325944

tape.Q -3.9001 0.9355 -4.169 0.014045 *

Multiple R-squared: 0.963, Adjusted R-squared: 0.9259

F-statistic: 26 on 4 and 4 DF, p-value: 0.004013

Two-way ANOVA with more than one observation per cell:

Consider the setting where \(n_{ij}=n>1\) for all \(i,j\).

• We call this design orthogonal

Our model \[y_{ij} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ij}\] is now identifiable and can be examined in the usual ways.

data(pvc); attach(pvc)
par(mfrow=c(1,2))
stripchart (psize ~ resin, xlab="Particle size", ylab="Resin railcar")
stripchart (psize ~ operator, xlab="Particle size", ylab="0perator")

Can be hard to tell random noise vs. actual interactions.

mod2 <- lm(psize ~ operator*resin, pvc)
anova(mod2)

## Analysis of Variance Table
## 
## Response: psize
##                Df  Sum Sq Mean Sq F value    Pr(>F)    
## operator        2  20.718  10.359  7.0072   0.00401 ** 
## resin           7 283.946  40.564 27.4388 5.661e-10 ***
## operator:resin 14  14.335   1.024  0.6926   0.75987    
## Residuals      24  35.480   1.478                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(lm(psize ~ operator+resin, pvc))

## Analysis of Variance Table
## 
## Response: psize
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## operator   2  20.718  10.359   7.902   0.00135 ** 
## resin      7 283.946  40.564  30.943 8.111e-14 ***
## Residuals 38  49.815   1.311                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

par(mfrow=c(1,2))
plot(mod2, which=c(2,1))

• Some evidence of outliers.
• Symmetry is because, for the full model, the fitted value for each combination of factors is the mean of the two replicates.

mod4 <- lm(psize ~ operator+resin, pvc[-45,]) # remove outlier
anova(mod4)

## Analysis of Variance Table
## 
## Response: psize
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## operator   2  28.215  14.108  15.973 9.982e-06 ***
## resin      7 283.479  40.497  45.853 2.497e-16 ***
## Residuals 37  32.678   0.883                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov(mod4))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mod4)
## 
## $operator
##          diff       lwr        upr     p adj
## 2-1 -0.262500 -1.073718  0.5487178 0.7114176
## 3-1 -1.777083 -2.601711 -0.9524561 0.0000184
## 3-2 -1.514583 -2.339211 -0.6899561 0.0001985
## 
## $resin
##           diff         lwr        upr     p adj
## 2-1 -1.0333333 -2.77520402  0.7085374 0.5566387
## 3-1 -5.8000000 -7.54187069 -4.0581293 0.0000000
## 4-1 -6.1833333 -7.92520402 -4.4414626 0.0000000
## 5-1 -4.8000000 -6.54187069 -3.0581293 0.0000000
## 6-1 -5.4500000 -7.19187069 -3.7081293 0.0000000
## 7-1 -3.6894444 -5.51633383 -1.8625551 0.0000037
## 8-1 -0.1833333 -1.92520402  1.5585374 0.9999706
## 3-2 -4.7666667 -6.50853735 -3.0247960 0.0000000
## 4-2 -5.1500000 -6.89187069 -3.4081293 0.0000000
## 5-2 -3.7666667 -5.50853735 -2.0247960 0.0000009
## 6-2 -4.4166667 -6.15853735 -2.6747960 0.0000000
## 7-2 -2.6561111 -4.48300050 -0.8292217 0.0009360
## 8-2  0.8500000 -0.89187069  2.5918707 0.7663296
## 4-3 -0.3833333 -2.12520402  1.3585374 0.9962814
## 5-3  1.0000000 -0.74187069  2.7418707 0.5962925
## 6-3  0.3500000 -1.39187069  2.0918707 0.9978942
## 7-3  2.1105556  0.28366617  3.9374449 0.0141248
## 8-3  5.6166667  3.87479598  7.3585374 0.0000000
## 5-4  1.3833333 -0.35853735  3.1252040 0.2072051
## 6-4  0.7333333 -1.00853735  2.4752040 0.8725789
## 7-4  2.4938889  0.66699950  4.3207783 0.0021575
## 8-4  6.0000000  4.25812931  7.7418707 0.0000000
## 6-5 -0.6500000 -2.39187069  1.0918707 0.9274222
## 7-5  1.1105556 -0.71633383  2.9374449 0.5263790
## 8-5  4.6166667  2.87479598  6.3585374 0.0000000
## 7-6  1.7605556 -0.06633383  3.5874449 0.0658042
## 8-6  5.2666667  3.52479598  7.0085374 0.0000000
## 8-7  3.5061111  1.67922172  5.3330005 0.0000099

• The factors in these models are fixed effects.
  • If factors were randomly selected from larger population, they would be random effects, which are a little more complicated to deal with.
• Also important that observations be genuine replications.
  • If they are correlated, the analysis will need to be adjusted.
  • Repeated measurements are not replications and will be correlated.
    • Repeated measures ANOVA, etc.

Suppose we have three or more factors \(\alpha, \beta, \gamma, \dots\) each with some number of levels \(l_{\alpha}, l_{\beta}, l_{\gamma}, \dots\).

• A full factorial experiment has at least one run for each combination of factor levels (equal to \(l_{\alpha}\times l_{\beta}\times l_{\gamma}\times\dots\))
  • This can easily be a very large number of treatments!
  • Could also include all possible interactions, higher order terms, etc.
  • Usually not performed for more than three or four factors.

• If no interactions significant, we get several one-way experiments for the cost of one.
• Efficient with experimental resources.
• Often done without replication due to cost concerns.
  • May need to assume some higher order interactions are zero in order to have enough degrees of freedom to examine lower level effects.

• Use only a fraction of the number of runs done in a full factorial experiment.

Consider an experiment with 7 factors, each with 2 levels:

• Assuming higher order interactions are negligible will save us from doing \(2^7=128\) runs.
• 8 runs would allow us to estimate main effects.
• Hard to find a good design to examine all 2-way effects without a lot of runs.

Goal: estimate as many parameters as possible with as few data points as possible.

• May not have many degrees of freedom left over.
• \(\sigma^2\) must be small for this to work (to distinguish significant effects).

data(speedo)
speedo # two levels of each factor: +/-

##    h d l b j f n a i e m c k g o      y
## 1  - - + - + + - - + + - + - - + 0.4850
## 2  + - - - - + + - - + + + + - - 0.5750
## 3  - + - - + - + - + - + + - + - 0.0875
## 4  + + + - - - - - - - - + + + + 0.1750
## 5  - - + + - - + - + + - - + + - 0.1950
## 6  + - - + + - - - - + + - - + + 0.1450
## 7  - + - + - + - - + - + - + - + 0.2250
## 8  + + + + + + + - - - - - - - - 0.1750
## 9  - - + - + + - + - - + - + + - 0.1250
## 10 + - - - - + + + + - - - - + + 0.1200
## 11 - + - - + - + + - + - - + - + 0.4550
## 12 + + + - - - - + + + + - - - - 0.5350
## 13 - - + + - - + + - - + + - - + 0.1700
## 14 + - - + + - - + + - - + + - - 0.2750
## 15 - + - + - + - + - + - + - + - 0.3425
## 16 + + + + + + + + + + + + + + + 0.5825

mod5 <- lm(y ~ ., speedo)
summary(mod5) # all residuals are 0!

## 
## Call:
## lm(formula = y ~ ., data = speedo)
## 
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.5825000        NaN     NaN      NaN
## h-          -0.0621875        NaN     NaN      NaN
## d-          -0.0609375        NaN     NaN      NaN
## l-          -0.0271875        NaN     NaN      NaN
## b-           0.0559375        NaN     NaN      NaN
## j-           0.0009375        NaN     NaN      NaN
## f-          -0.0740625        NaN     NaN      NaN
## n-          -0.0065625        NaN     NaN      NaN
## a-          -0.0678125        NaN     NaN      NaN
## i-          -0.0428125        NaN     NaN      NaN
## e-          -0.2453125        NaN     NaN      NaN
## m-          -0.0278125        NaN     NaN      NaN
## c-          -0.0896875        NaN     NaN      NaN
## k-          -0.0684375        NaN     NaN      NaN
## g-           0.1403125        NaN     NaN      NaN
## o-          -0.0059375        NaN     NaN      NaN
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 15 and 0 DF,  p-value: NA

We have no degrees of freedom, as number cases = number parameters.

model.matrix(mod5) # + is coded as 0; - as 1 bx of ASCII alphabet

##    (Intercept) h- d- l- b- j- f- n- a- i- e- m- c- k- g- o-
## 1            1  1  1  0  1  0  0  1  1  0  0  1  0  1  1  0
## 2            1  0  1  1  1  1  0  0  1  1  0  0  0  0  1  1
## 3            1  1  0  1  1  0  1  0  1  0  1  0  0  1  0  1
## 4            1  0  0  0  1  1  1  1  1  1  1  1  0  0  0  0
## 5            1  1  1  0  0  1  1  0  1  0  0  1  1  0  0  1
## 6            1  0  1  1  0  0  1  1  1  1  0  0  1  1  0  0
## 7            1  1  0  1  0  1  0  1  1  0  1  0  1  0  1  0
## 8            1  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1
## 9            1  1  1  0  1  0  0  1  0  1  1  0  1  0  0  1
## 10           1  0  1  1  1  1  0  0  0  0  1  1  1  1  0  0
## 11           1  1  0  1  1  0  1  0  0  1  0  1  1  0  1  0
## 12           1  0  0  0  1  1  1  1  0  0  0  0  1  1  1  1
## 13           1  1  1  0  0  1  1  0  0  1  1  0  0  1  1  0
## 14           1  0  1  1  0  0  1  1  0  0  1  1  0  0  1  1
## 15           1  1  0  1  0  1  0  1  0  1  0  1  0  1  0  1
## 16           1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## attr(,"assign")
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## attr(,"contrasts")
## attr(,"contrasts")$h
## [1] "contr.treatment"
## 
## attr(,"contrasts")$d
## [1] "contr.treatment"
## 
## attr(,"contrasts")$l
## [1] "contr.treatment"
## 
## attr(,"contrasts")$b
## [1] "contr.treatment"
## 
## attr(,"contrasts")$j
## [1] "contr.treatment"
## 
## attr(,"contrasts")$f
## [1] "contr.treatment"
## 
## attr(,"contrasts")$n
## [1] "contr.treatment"
## 
## attr(,"contrasts")$a
## [1] "contr.treatment"
## 
## attr(,"contrasts")$i
## [1] "contr.treatment"
## 
## attr(,"contrasts")$e
## [1] "contr.treatment"
## 
## attr(,"contrasts")$m
## [1] "contr.treatment"
## 
## attr(,"contrasts")$c
## [1] "contr.treatment"
## 
## attr(,"contrasts")$k
## [1] "contr.treatment"
## 
## attr(,"contrasts")$g
## [1] "contr.treatment"
## 
## attr(,"contrasts")$o
## [1] "contr.treatment"

Can’t do usual F-tests, so need different approach.

• Suppose no significant effects and errors normally distributed.
  • Estimated effects would be linear combos of errors and so also normal.
  • We can make a normal QQ plot of main effects.
    • Idea: outliers represent something interesting happening (significant effects).

qqnorm(coef(mod5)[-1], pch=names(coef(mod5)[-1]))

So we might further examine effects of e and g.