2024-02-01

Scenario 1: Little evidence for interaction

Interaction is “non-significant”

Should we remove from model?

Depends:

  • Traditionally: cull anything that’s not “significant” — bad.
  • Observational study, many EVs: Yes if \(P\) is very weak — see later: model simplification.
  • Planned Experiment: No.

Reporting

  • Keep the full model: blocking, covariate and interaction terms stay!
  • Give the ANOVA for the full model.
  • If interaction can be ignored: two simple stories

Back to Wheat example

m.2 <- lm(yield ~ block + sow.dens + variety + sow.dens:variety, Wheat)
anova(m.2)
## Analysis of Variance Table
## 
## Response: yield
##                  Df Sum Sq Mean Sq F value  Pr(>F)
## block             2 0.3937 0.19683  0.5193 0.60599
## sow.dens          3 5.8736 1.95786  5.1650 0.01303
## variety           1 2.1474 2.14742  5.6651 0.03206
## sow.dens:variety  3 0.0566 0.01887  0.0498 0.98470
## Residuals        14 5.3069 0.37906

Clearly no evidence for sow.dens:variety \((P=0.985)\), but we keep the model and use it to tell two simple stories.

  • Story 1: effect of variety
  • Story 2: effect of sow.dens

Effect of variety

cf <- coefficients( summary(m.2) )  # save table of coefs
est <- cf[, "Estimate"]             # save `Estimate` column
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.5400 0.126 67.900 0.0000
block1 -0.0847 0.178 -0.477 0.6410
block2 0.1810 0.178 1.020 0.3260
sow.dens1 -0.2800 0.218 -1.290 0.2190
sow.dens2 -0.4470 0.218 -2.050 0.0593
sow.dens3 -0.1030 0.218 -0.475 0.6420
variety1 0.2990 0.126 2.380 0.0321
sow.dens1:variety1 0.0442 0.218 0.203 0.8420
sow.dens2:variety1 0.0500 0.218 0.230 0.8220
sow.dens3:variety1 -0.0635 0.218 -0.292 0.7750

Expected Marginal Means: Easy with sum contrasts:

  • EMM of Variety 1: est["(Intercept)"] + est["variety1"] \(= 8.84\)
  • EMM of Variety 2: est["(Intercept)"] - est["variety1"] \(= 8.24\)

Difference between Variety 1 and 2: est["variety1"] - -est["variety1"]
or 2 * est["variety1"] \(= 0.598\)

Effect of sow.dens

Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.5400 0.126 67.900 0.0000
block1 -0.0847 0.178 -0.477 0.6410
block2 0.1810 0.178 1.020 0.3260
sow.dens1 -0.2800 0.218 -1.290 0.2190
sow.dens2 -0.4470 0.218 -2.050 0.0593
sow.dens3 -0.1030 0.218 -0.475 0.6420
variety1 0.2990 0.126 2.380 0.0321
sow.dens1:variety1 0.0442 0.218 0.203 0.8420
sow.dens2:variety1 0.0500 0.218 0.230 0.8220
sow.dens3:variety1 -0.0635 0.218 -0.292 0.7750
  • EMM of Density 1: est["(Intercept)"] + est["sow.dens1"] \(= 8.26\)
  • EMM of Density 2: est["(Intercept)"] + est["sow.dens2"] \(= 8.09\)
  • EMM of Density 3: est["(Intercept)"] + est["sow.dens3"] \(= 8.43\)
  • EMM of Density 4: est["(Intercept)"] - (est["sow.dens1"]+est["sow.dens2"]+est["sow.dens3"]) \(= 9.37\)
  • Difference between Density 1 and 2:
    est["sow.dens1"] - est["sow.dens2"] \(= 0.167\)

Easy EMMs with emmeans package

library(emmeans)
emmeans(m.2, ~ variety) %>%  # the rest is just making pretty table
  kbl %>% kable_classic(full_width = F, html_font = "Cambria")
## NOTE: Results may be misleading due to involvement in interactions
variety emmean SE df lower.CL upper.CL
1 8.83550 0.1777319 14 8.454303 9.216697
2 8.23725 0.1777319 14 7.856053 8.618447 
emmeans(m.2, ~ sow.dens) %>%  # the rest is just making pretty table
  kbl %>% kable_classic(full_width = F, html_font = "Cambria")
## NOTE: Results may be misleading due to involvement in interactions
sow.dens emmean SE df lower.CL upper.CL
1 8.2560 0.2513509 14 7.716906 8.795094
2 8.0895 0.2513509 14 7.550406 8.628594
3 8.4330 0.2513509 14 7.893906 8.972094
4 9.3670 0.2513509 14 8.827906 9.906094

Plotting the effects

emmip(m.2, ~ variety, CI=TRUE)

emmip(m.2, ~ sow.dens, CI=TRUE)

Scenario 2: Evidence for interaction

Interaction is “significant:” Wheat2 data

ANOVA

m.3 <- lm(yield ~ block + sow.dens + variety + sow.dens:variety, Wheat2)
( a <- anova(m.3) )  # "save and print" trick
## Analysis of Variance Table
## 
## Response: yield
##                  Df Sum Sq Mean Sq F value   Pr(>F)
## block             2 0.3937  0.1968  0.5193 0.605986
## sow.dens          3 1.7790  0.5930  1.5644 0.242244
## variety           1 5.6833  5.6833 14.9930 0.001692
## sow.dens:variety  3 5.8422  1.9474  5.1374 0.013275
## Residuals        14 5.3069  0.3791
  • Interaction sow.dens:variety \(P=0.0133\).
  • Note that here, sow.dens \(P=0.242\),
    but sow.dens is clearly important!

Interaction Plot

emmip(m.3, variety ~ sow.dens, CI=TRUE) +
  scale_color_manual(values=c("deepskyblue", "darkorange"))

Estimated Means (marginal on block!)

emmeans(m.3, ~ sow.dens + variety)
##  sow.dens variety emmean    SE df lower.CL upper.CL
##  1        1         8.60 0.355 14     7.84     9.36
##  2        1         9.94 0.355 14     9.18    10.70
##  3        1         9.67 0.355 14     8.91    10.43
##  4        1         8.64 0.355 14     7.87     9.40
##  1        2         7.91 0.355 14     7.15     8.68
##  2        2         7.74 0.355 14     6.98     8.50
##  3        2         8.20 0.355 14     7.43     8.96
##  4        2         9.10 0.355 14     8.34     9.86
## 
## Results are averaged over the levels of: block 
## Confidence level used: 0.95

Coefficients

cf <- coefficients( summary(m.3) )  # save table of coefs
est <- cf[, "Estimate"]             # save `Estimate` column
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.7200 0.126 69.400 0.00000
block1 -0.0847 0.178 -0.477 0.64100
block2 0.1810 0.178 1.020 0.32600
sow.dens1 -0.4680 0.218 -2.150 0.04960
sow.dens2 0.1160 0.218 0.531 0.60400
sow.dens3 0.2090 0.218 0.961 0.35300
variety1 0.4870 0.126 3.870 0.00169
sow.dens1:variety1 -0.1430 0.218 -0.658 0.52100
sow.dens2:variety1 0.6130 0.218 2.810 0.01380
sow.dens3:variety1 0.2490 0.218 1.140 0.27200

Grand Mean

  • Grand mean: est["(Intercept)"] \(=8.72\)

EMM for Variety 1

  • EMM for Variety 1 (across all sow densities):
    est["(Intercept)"] + est["variety1"] \(=9.21\)

EMM for sow density 2

  • EMM for sow density 2 (across both varieties):
    est["(Intercept)"] + est["sow.dens2"] \(=8.84\)

Yield for Variety 1 at sow density 2?

  • Sow density and Variety are non-additive:
    est["(Intercept)"] + est["sow.dens2"] + est["variety1"] + est["sow.dens2:variety1"] \(=9.94\)

Yield for Variety 2 at sow density 3?

  • Sow density and Variety are non-additive:
    est["(Intercept)"] + est["sow.dens3"] - est["variety1"] - est["sow.dens3:variety1"] \(=8.2\)