Part (a)

Explain the difference between crossed factors and nested factors.

Factors \(A\) and \(B\) are crossed in an experimental design if the levels of \(B\) are the same at each level of \(A\).

Factor \(B\) is nested within factor \(A\) if the levels of \(B\) are different for each of the levels of factor \(A\).

Part (b)

Write the model for a nested design. Provide algebraic formulas for the estimates \(\hat{\tau}_i\) and \(\hat{\beta})_{j(i)}\).

Given that \(A\) and \(B\) are fixed factors, \[ Y_{i j k} = \mu + \tau_i + \beta_{j(i)} + \epsilon_{i j k} \begin{cases} i = 1,\ldots,a\\ j = 1,\ldots,b\\ k = 1,\ldots,n \end{cases} \] where \(\{\tau_i\}\) (\(\sum_i \tau_i = 0\)) with \(\operatorname{df_{A}}=a-1\) and \(\{\beta_(j(i))\}\) (\(\sum_j \beta_{j(i)} = 0\) for \(i=1,\ldots,a\)) with \(\operatorname{df_{B}} = a(b-1)\) are fixed effects.

Estimators are given by \[ \hat{\tau}_i = \bar{Y}_{i\cdot\cdot} - \bar{Y}_{\cdot\cdot\cdot} \] and \[ \hat{\beta}_{j(i)} = \bar{Y}_{i j \cdot} - \bar{Y}_{i \cdot\cdot}. \]

We see that level factors of \(A\) are compared with \(\{\hat{\tau}_i\}\) while level factors of \(B\) are compared only with the same level of \(A\), e.g., at the \(i\)-th level of \(A\), level factors of \(B\) are compared with \(\{\beta_j(i) \}\).

Part (c)

Compute the \(F_A\) statistic for testing factor \(A\) effects, and the \(F_{B(A)}\) statistic for testing nested factor \(B\) effects. Compute the \(p\)-values. Provide an overall interpretation, stated in the context of the problem.

For fixed factor \(A\) and fixed factor \(B\) nested in \(A\), \[ F_A = \operatorname{MS_{A}} / \operatorname{MS_{E}} \] and \[ F_{B(A)} = \operatorname{MS_{B(A)}} / \operatorname{MS_{E}} \] where \[\begin{align*} \operatorname{MS_{A}} &= \frac{b n \sum_{i=1}^{a} \hat{\tau}_i^2 }{a-1},\\ \operatorname{MS_{B(A)}} &= \frac{n \sum_{i=1}^{a} \sum_{j=1}^{b} \hat{\beta}_{j(i)}^2}{a(b-1)},\\ \operatorname{MS_{E}} &= \frac{\sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} (Y_{i j k} - \bar{Y}_{i j \cdot})}{a b (n-1)}. \end{align*}\]

Under the model \(\tau_1 = \cdots \tau_a = 0\), \[ F_A \sim F(\operatorname{df_{A}},\operatorname{df_{E}}) = F(a-1,a b (n-1)) \] and under \(\beta_{j(i)} = 0\) for all \(i,j\), \[ F_{B(A)} \sim F(\operatorname{df_{B(A)}},\operatorname{df_{E}}) = F(a(b-1),a b (n-1)). \]

We perform these computations in R with:

library("readxl")
data = read_excel("handout9data.xlsx")
A = as.factor(na.omit(data$mchine))
B = as.factor(na.omit(data$operator))
y = na.omit(data$surface)

# use contrasts to define parameter restrictions for the fixed effects in the
# model
contrasts(A)=contr.sum
contrasts(B)=contr.sum

# to fit a model with a nested fixed effect, we use the / notation within the
# aov command.
nested.mod = aov(y ~ A/B)
summary(nested.mod) # compute F statistics for factor effects

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            3   3618  1205.9  14.271 0.000291 ***
## A:B          8   2818   352.2   4.168 0.013408 *  
## Residuals   12   1014    84.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We see that \(F_A = 14.271\) (\(p\)-value \(= .000\)) and \(F_{B(A)} = 4.168\) (\(p\)-value \(= .013\)).

Part (d)

Compute estimates for each of the effect parameters. Identify which machine performs best, and which operator performs best on each machine. (Higher scores of response are preferred.)

# compute parameter estimates
estimates = dummy.coef(nested.mod)
estimates$A

##          1          2          3          4 
##   2.833333  17.333333  -3.333333 -16.833333

library("multcomp")
# The following commands are used to perform pairwise comparisons between machines.
# The glht command should be familiar from previous sections on multiple comparisons.
compare.A = glht(nested.mod,linfct = mcp(A="Tukey"))
c.m = summary(compare.A,test=adjusted("none"))
c.m

## 
## 	 Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = y ~ A/B)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)    
## 2 - 1 == 0   14.500      5.307   2.732  0.01819 *  
## 3 - 1 == 0   -6.167      5.307  -1.162  0.26785    
## 4 - 1 == 0  -19.667      5.307  -3.706  0.00300 ** 
## 3 - 2 == 0  -20.667      5.307  -3.894  0.00213 ** 
## 4 - 2 == 0  -34.167      5.307  -6.438 3.22e-05 ***
## 4 - 3 == 0  -13.500      5.307  -2.544  0.02576 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)

cld(c.m)

##   1   2   3   4 
## "a" "b" "a" "c"

library("lsmeans")
# comparisons only involving operators (B) within machines (A)
x.mod = aov(y~A*B)
lsmeans(x.mod, pairwise ~ A:B, adjust="none")

Part (e)

Explain why it is not possible to directly compare operators across machines in the above design.

We cannot directly compare operators (factor \(B\) levels) across machines (factor \(A\) levels) since we only have data for operator performance with respect to a particular machine (\(A\) level).

However, we can compare operators within the same machine, or compare operator \(\times\) machine combinations.

Part (a)

Provide the algebraic formulas for \(\operatorname{MS_{A}}\), \(\operatorname{MS_{B(A)}}\), and \(\operatorname{MS_{E}}\).

\[\begin{align*} \operatorname{MS_{A}} &= \frac{b n \sum_{i=1}^{a} \hat{\tau}_i }{a-1},\\ \operatorname{MS_{B(A)}} &= \frac{n \sum_{i=1}^{a} \sum_{j=1}^{b} \hat{\beta}_(j(i)}{a(b-1)},\\ \operatorname{MS_{E}} &= \frac{\sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} (Y_{i j k} - \bar{Y}_{i j \cdot})}{a b (n-1)}. \end{align*}\]

Part (b)

State the expected value for each of the mean squares.

\[\begin{align*} E(\operatorname{MS_{A}}) &= \sigma^2 + n \sigma_{\beta}^2 + \frac{b n}{a-1} \sum_{i=1}^{a} \tau_i^2,\\ E(\operatorname{MS_{B(A)}}) &= \sigma^2 + n \sigma_{\beta}^2,\\ E(\operatorname{MS_{E}}) &= \sigma^2. \end{align*}\]

Part (c)

Test for differences between the production processes. Write the \(F_A\) statistic as a ratio of mean squares. Compute \(F_A\) and the \(p\)-value. Provide an interpretation, stated in the context of the problem.

Under the null model \(H_0 : \tau_1 = \cdots = \tau_a = 0\), \(\operatorname{MS_{A}}\) and \(\operatorname{MS_{B(A)}}\) have the same expected value, and thus an appropriate test statistic for testing evidence of \(H_0\) is \[ F_A = \frac{\operatorname{MS_{A}}}{\operatorname{MS_{B(A)}}}. \]

We load the data in R with:

A = as.factor(na.omit(data$proc))
B = as.factor(na.omit(data$batc))
y = na.omit(data$burn_rate)

We compute \(F_A\) in R with:

nested.test(A,B,y)

##                           SS df        MS
## Fixed Effect A      676.0556  2 338.02778
## Random Effect B(A) 2077.5833  9 230.84259
## Error               454.0000 24  18.91667
##  F-test for fixed effect   p-value
##                 1.464322 0.2814697
##  error.var    B.var
##   18.91667 70.64198

We see that \(F_A = 1.464\) (\(p\)-value \(= .281\)).

Part (d)

Explain why \(\operatorname{MS_{E}}\) is the incorrect error term to use when the nested factor is random. In particular, comment on the pertinent sample size.

I touched on this in part (c) (\(\operatorname{MS_{A}}\) and \(\operatorname{MS_{B(A)}}\) have the same expected value under the null model \(\tau_1 = \cdots = \tau_a = 0\)).

In addition, we may also think of the batch as the experimental unit and thus the appropriate error term is then a measure of batch variance.

Finally, taking repeat measurements on each random factor level does not increase the pertinent sample size.

Part (e)

Illustrate how evidence in favor of a process effect would be overstated if \(\operatorname{MS_{E}}\) is used when computing the test statistic.

Under the null model, \(H_0 : \tau_1 = \cdots = \tau_a = 0\), \(E(\operatorname{MS_{A}}) = \sigma^2 + n \sigma_{\beta}^2\), and \(E(\operatorname{MS_{E}}) = \sigma^2\). The reference distribution \(F(\operatorname{df_{A}},\operatorname{df_{E}})\) assumes the numerator and denominator have the same expected value under \(H_0\), thus the test statistic \[ \frac{\operatorname{MS_{A}}}{\operatorname{MS_{E}}} \] will be inflated and consequently overstates the effect of the process (factor \(A\)).

More generally, since \(n \sigma_{\beta}^2 > 0\), even if \(H_0\) is not the case, the significance of the evidence in favor of a process effect will necessarily be overstated.

Part (f)

Compute estimates of the batch variance and the measurement variance.

An estimator for the measurement variance (within batches) is \(\hat{\sigma}^2 = \operatorname{MS_{E}}\) and an estimator for the batch variance (between batches) is given by \(\hat{\sigma}_{\beta}^2 = \frac{\operatorname{MS_{B(A)}} - \operatorname{MS_{E}}}{n}\).

In part (c), nested.test showed the outputs \[ \hat{\sigma}^2 = 18.917 \] and \[ \hat{\sigma}_{\beta}^2 = 70.642. \]