Regression Adjustments

Study: Meadowfoam

Study: Meadowfoam

Meadowfoam is a flowering California native plant.

Can we effect the number of flowers produced by plants by providing extra light treatment, which varies in its intensity as well as it’s timing?

A model for a categorical and continuous effect (ANCOVA):

\[ y_{ij} = \mu + \alpha_i + \beta_1 x_{ij} + \epsilon_{ij} \]

Each group has a linear model with the same effect on \(x\) and different effect on the categorical effect (two intercepts).

\[\begin{align} y_{i1} &= \mu + \alpha_1 + \beta_1 x_{i1} + \epsilon_{i1} = \mu_1 + \beta_1 x_{i1} + \epsilon_{i1} \\ y_{i2} &= \mu + \alpha_2 + \beta_1 x_{i2} + \epsilon_{i2} = \mu_2 + \beta_1 x_{i2} + \epsilon_{i2} \end{align}\]

Regression Adjustments

Meadowfoam

meadowfoam
# A tibble: 24 × 3
   flowers time  intensity
     <dbl> <chr>     <dbl>
 1    62.3 late        150
 2    77.4 late        150
 3    55.3 late        300
 4    54.2 late        300
 5    49.6 late        450
 6    61.9 late        450
 7    39.4 late        600
 8    45.7 late        600
 9    31.3 late        750
10    44.9 late        750
# ℹ 14 more rows

Let’s say that only time was randomly assigned and intensity was another variable that was recorded that we know is a good predictor of flowers.

What we’ve been doing

Since time is the only experimental factor, we would fit this model:

options(contrasts=c("contr.sum","contr.poly"))
m1 <- lm(flowers ~ time, data = meadowfoam)
summary(m1)$coef
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 56.137500   2.556552 21.958286 1.874990e-16
time1        6.079166   2.556552  2.377877 2.652621e-02

Adjust for intensity

m_full <- lm(flowers ~ intensity + time, data = meadowfoam)
m_no_time <- lm(flowers ~ intensity, data = meadowfoam)
summary(m_full)$coef
               Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 77.38500017 2.99815637 25.810862 2.166375e-17
intensity   -0.04047143 0.00513237 -7.885525 1.036787e-07
time1        6.07916641 1.31477854  4.623719 1.463777e-04

summary(m_no_time)$coef
               Estimate  Std. Error   t value     Pr(>|t|)
(Intercept) 77.38500017 4.161186173 18.596861 6.059011e-15
intensity   -0.04047143 0.007123293 -5.681562 1.029503e-05

anova(m_no_time, m_full)
Analysis of Variance Table

Model 1: flowers ~ intensity
Model 2: flowers ~ intensity + time
  Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
1     22 1758.19                                  
2     21  871.24  1    886.95 21.379 0.0001464 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cautions

When fitting regression models to experimental data:

  1. Non-experimental variables have no causal interpretation.

  2. Regression adjustments generally should not be used in very small datasets (\(N < 20\)) because df are better spent estimating causal effects.

Collinearity

Simulated data with two explanatory variables.

# A tibble: 100 × 3
        y      x1      x2
    <dbl>   <dbl>   <dbl>
 1  3.58  -0.349  -0.708 
 2  2.03  -0.756   1.12  
 3  8.59  -1.04   -1.07  
 4 -7.73   1.21    1.41  
 5 -2.51   1.24   -1.90  
 6 -7.48   1.74   -0.983 
 7  0.319  0.0977 -0.179 
 8 17.2   -2.52   -1.93  
 9  4.07  -0.789   0.467 
10 -3.49   0.576  -0.0249
# ℹ 90 more rows

pairs(weak_collinearity)

cor(weak_collinearity)
            y         x1         x2
y   1.0000000 -0.9119516 -0.3798014
x1 -0.9119516  1.0000000  0.0144365
x2 -0.3798014  0.0144365  1.0000000

Collinearity: describes the presence of linear relationships between the explanatory variables.

Results from Linear Models

\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \]

\[ Var(\hat{\boldsymbol{\beta}}) = \sigma^2 (\mathbf{X}^\top \mathbf{X})^{-1} \]

\[ Var(\hat{\beta}_j) = \frac{\sigma^2}{(n-1) Var(X_j)} \cdot \frac{1}{1 - R^2_j} \]

where \(R^2_j\) is the \(R^2\) from regressing \(X_j\) on other \(X\). The second term is called the Variance Inflation Factor (VIF).

Weak Collinearity

m1 <- lm(y ~ x1 + x2, data = weak_collinearity)
summary(m1)$coef
               Estimate Std. Error     t value     Pr(>|t|)
(Intercept)  0.05407907 0.09898302   0.5463469 5.860820e-01
x1          -5.04493666 0.10402428 -48.4976860 8.096236e-70
x2          -1.95018390 0.09941965 -19.6156791 1.569072e-35

This linear model corresponds to a plane in 3D.

Strong Collinearity

# A tibble: 100 × 3
        y     x1      x2
    <dbl>  <dbl>   <dbl>
 1  5.93  -0.766 -0.806 
 2  0.846 -0.200  0.0721
 3  2.33  -0.393 -0.427 
 4  8.43  -1.31  -1.24  
 5  1.23  -0.204 -0.133 
 6 -5.82   0.687  0.767 
 7 -2.19   0.350  0.120 
 8  3.68  -0.339 -0.255 
 9 -1.48   0.297 -0.0168
10 -5.03   0.722  0.616 
# ℹ 90 more rows

Strong Collinearity

pairs(strong_collinearity)

Strong collinearity

cor(strong_collinearity)
            y         x1         x2
y   1.0000000 -0.9875549 -0.9791803
x1 -0.9875549  1.0000000  0.9841485
x2 -0.9791803  0.9841485  1.0000000

High collinearity leads to high VIF.

Designed Experiment

# A tibble: 100 × 3
         y    x1    x2
     <dbl> <dbl> <dbl>
 1  20.6      -3    -3
 2  12.0      -1    -3
 3   6.52      0    -3
 4   1.19      1    -3
 5  -9.76      3    -3
 6  18.9      -3    -1
 7   6.23     -1    -1
 8   0.981     0    -1
 9  -3.15      1    -1
10 -12.3       3    -1
# ℹ 90 more rows

Designed Experiment

pairs(designed_exp)

Design Experiment

cor(designed_exp)
            y         x1        x2
y   1.0000000 -0.9256147 -0.369353
x1 -0.9256147  1.0000000  0.000000
x2 -0.3693530  0.0000000  1.000000

No collinearity leads to . . .