RCBD Diagnostics

Multiple comparison and model checks for RCBDs

Kris Sankaran true
10-05-2021

Readings 4.1, Rmarkdown

Multiple Comparisons

  1. As before, we may want to use contrasts to decide on which treatment effects are different We can continue to use the same multiple comparisons procedures as before, but have to account for a few differences,
    • The number of samples per treatment \(n\) is replaced by the number of blocks \(b\).
    • The d.f. for \(MS_{E}\) is now \((a - 1)(b - 1)\), not \(N - a\).

  1. For example, the cutoff in Fisher’s LSD becomes \[ t_{\frac{\alpha}{2}, \left(a - 1\right)\left(b - 1\right)}\sqrt{\frac{2 MS_{E}}{b}} \]

Model Diagnostics

  1. There are two key assumptions when we use RCBDs,
    • \(\epsilon \sim \mathcal{N}\left(0, \sigma^2\right)\) are independently distributed. Note that this also implies homoskedasticity (the variances are not changing from sample to sample).
    • Additivity (i.e. “no interaction”). The treatment effects \(\tau_i\) need to be the same across all blocks.
An example where the additivity assumption fails.

Figure 1: An example where the additivity assumption fails.

  1. We’ve seen the first assumption before, and can continue to use normal probability plots and residual analysis to check it.

  2. One way to check additivity is to look at residuals, and see whether they are consistently lower / higher in some blocks.

  3. What can we do if we find an interaction effect?

    • Sometimes the interaction effects will disappear after transforming the response (e.g., using \(\log\) or \(\sqrt{}\))
    • Otherwise, another design may be necessary. Factorial designs (to be discussed soon) allow for inference even in the presence of interactions.

Code Example

  1. We will continue the graft example from the previous notes. First, let’s refit the RCBD model.
library(readr)
library(tidyr)
graft <- read_table2("https://uwmadison.box.com/shared/static/0ciblk4z2f3k6zizbj4plg3q33w1d0n6.txt") %>%
  pivot_longer(-Pressure, names_to = "batch", values_to = "yield")
graft$Pressure <- as.factor(graft$Pressure)
fit <- lm(yield ~ Pressure + batch, graft)
aov_table <- aov(fit)
  1. We can fit contrasts and correct for multiple comparisons using the same type of code as for ANOVA without batch effects.
library(gmodels)
library(DescTools)

contrasts <- matrix(
  c(1, 1, -1, -1,
    1, 0, 0, -1),
  nrow = 2, byrow=TRUE
)

fit.contrast(aov_table, "Pressure", contrasts)

                         Estimate Std. Error  t value     Pr(>|t|)
Pressure c=( 1 1 -1 -1 ) 9.816667   2.209940 4.442052 0.0004751988
Pressure c=( 1 0 0 -1 )  7.050000   1.562663 4.511528 0.0004136854
attr(,"class")
[1] "fit_contrast"
PostHocTest(aov_table, method = "scheffe", contrasts = contrasts)


  Posthoc multiple comparisons of means: Scheffe Test 
    95% family-wise confidence level

$Pressure
                           diff     lwr.ci     upr.ci   pval    
8500,8700,8900,9100-  184.50000  177.31746  191.68254 <2e-16 ***
8500,8900-             88.91667   83.83785   93.99549 <2e-16 ***
-8500,8900            -92.81667  -97.89549  -87.73785 <2e-16 ***
-8500,8700,8900,9100 -174.68333 -181.86587 -167.50080 <2e-16 ***

$batch
                                               diff     lwr.ci
Batch1,Batch2,Batch3,Batch4,Batch5,Batch6-  177.450  168.65322
Batch1,Batch3,Batch5-                        91.000   84.77974
-Batch1,Batch3,Batch5                       -85.325  -91.54526
-Batch1,Batch2,Batch3,Batch4,Batch5,Batch6 -177.450 -186.24678
                                               upr.ci   pval    
Batch1,Batch2,Batch3,Batch4,Batch5,Batch6-  186.24678 <2e-16 ***
Batch1,Batch3,Batch5-                        97.22026 <2e-16 ***
-Batch1,Batch3,Batch5                       -79.10474 <2e-16 ***
-Batch1,Batch2,Batch3,Batch4,Batch5,Batch6 -168.65322 <2e-16 ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov_table)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = fit)

$Pressure
               diff        lwr       upr     p adj
8700-8500 -1.133333  -5.637161  3.370495 0.8854831
8900-8500 -3.900000  -8.403828  0.603828 0.1013084
9100-8500 -7.050000 -11.553828 -2.546172 0.0020883
8900-8700 -2.766667  -7.270495  1.737161 0.3245644
9100-8700 -5.916667 -10.420495 -1.412839 0.0086667
9100-8900 -3.150000  -7.653828  1.353828 0.2257674

$batch
                diff         lwr        upr     p adj
Batch2-Batch1  2.050  -4.1680828  8.2680828 0.8853016
Batch3-Batch1  3.300  -2.9180828  9.5180828 0.5376297
Batch4-Batch1  2.850  -3.3680828  9.0680828 0.6757699
Batch5-Batch1 -2.375  -8.5930828  3.8430828 0.8105903
Batch6-Batch1  6.750   0.5319172 12.9680828 0.0297368
Batch3-Batch2  1.250  -4.9680828  7.4680828 0.9845521
Batch4-Batch2  0.800  -5.4180828  7.0180828 0.9980198
Batch5-Batch2 -4.425 -10.6430828  1.7930828 0.2483499
Batch6-Batch2  4.700  -1.5180828 10.9180828 0.1986961
Batch4-Batch3 -0.450  -6.6680828  5.7680828 0.9998784
Batch5-Batch3 -5.675 -11.8930828  0.5430828 0.0837504
Batch6-Batch3  3.450  -2.7680828  9.6680828 0.4925715
Batch5-Batch4 -5.225 -11.4430828  0.9930828 0.1263042
Batch6-Batch4  3.900  -2.3180828 10.1180828 0.3674672
Batch6-Batch5  9.125   2.9069172 15.3430828 0.0027838
  1. To check the normality and additivity assumptions, we can plot the residuals against the batch. Even though additivity doesn’t seem to be a problem (no group of residuals is consistently \(> 0\) or \(< 0\)), equality of variance across batches seems to be violated. The variance of Batch 6’s residuals is too small.
library(ggplot2)
graft$residual <- resid(fit)
ggplot(graft) +
  geom_point(aes(x = batch, y = residual, col = Pressure)) +
  scale_color_brewer(palette = "Set2")