LECTURE 9: linear models part 2: categorical predictor w/ >2 levels

Outline

1) Motivating example

2) Extending the linear model

3) Dummy variables with > 2 levels

4) F-tests

5) Comparing linear model and ANOVA approaches

Categorical predictors

In lecture 4, we learned about fitting models with categorical predictors with 1 (measurement error) or 2 (frog) levels

• Linear models with these categorical predictors correspond to one-sample or two-sample t-tests

• Binary categorical predictors are quite common in ecological studies (e.g., male vs female, juvenile vs adult, disturbed vs undisturbed)

However, you will often encounter categorical predictors with >2 levels

• What are some examples?

Motivating example

Foresters are studying the effect of 4 different fertilizers (treatments) on the growth of loblolly pine, which are grown on 3 plots (replicates) receiving each treatment. Data are average height per plot after 5 years:

	Treatment
Replicate	A	B	C	D
1	11	9	7	5
2	10	8	6	4
3	9	7	5	3

Notation

• The number of groups (treatments) is \(a=4\)

• The number of observations within each group (replicates) is \(n=3\)

• \(y_{ij}\) denotes the \(j\)th observation from the \(i\)th group

a brief tangent

What counts as an observation?

Experimental unit

the physical unit that receives a particular treatment

Observational unit

the physical unit on which measurements are taken

These are not always the same!

Examples

• Agricultural fields given different fertilizer, crop yield measured

• Rats given different diets, disease state measured

• Microcosm given different predator abundance, tadpole growth measured

Motivating example

Question: Is there a difference in growth among the four treatment groups?

Motivating example

Hypotheses:

• \(\large H_0 : \mu_A = \mu_B = \mu_C = \mu_D\)

• \(\large H_a :\) At least one inequality

How should we test the null?

• We could do this using 6 t-tests

• But this would alter the overall (experiment-wise) \(\alpha\) level because each individual test has a chance (usually \(\alpha = 0.05\)) of incorrectly rejecting a true null hypothesis, and this is multiplied when multiple tests are used¹
  

Extending the linear model

Fortunately, extending the previous linear model to include additional predictor levels is “straightforward”

• The model looks familiar

\[\large y_{ij} = \beta_0 + \beta_ix_{ij} + \epsilon_{ij}\]

\[\large \epsilon_{ij} \sim normal(0, \sigma)\]

• What’s tricky is how we code \(\large x_{ij}\) when it is more than just 0 or 1

• Fortunately, what we learned about how R codes dummy variables will help

Dummy variable coding

View a subset (rows 1, 4, 7, & 10) of the data:

data(loblollydata)
loblollydata[c(1,4,7,10),]

   Treatment Replicate Height
1          A         1     11
4          B         1      9
7          C         1      7
10         D         1      5

View the corresponding rows of the model matrix²

fit1 <- lm(Height ~ Treatment, data = loblollydata)

model.matrix(fit1)[c(1,4,7,10),]

   (Intercept) TreatmentB TreatmentC TreatmentD
1            1          0          0          0
4            1          1          0          0
7            1          0          1          0
10           1          0          0          1

Dummy variable coding

model.matrix(fit1)[c(1,4,7,10),]

   (Intercept) TreatmentB TreatmentC TreatmentD
1            1          0          0          0
4            1          1          0          0
7            1          0          1          0
10           1          0          0          1

The model matrix has one row for each observation and one column for every parameter in the model

• How many parameters is this model estimating?

Multiplying the model matrix by the parameter matrix gives the expected value of each observation³

\[\large E[y_{ij}] = \beta_0 + \beta_1 I(B)_{ij} + \beta_2 I(C)_{ij} + \beta_3 I(D)_{ij}\]

where \(I(B/C/D)_{ij}\) are dummy variables (0/1 depending on treatment)

• What is the interpretation of each parameter?

Motivating example

broom::tidy(fit1)

# A tibble: 4 × 5
  term        estimate std.error statistic     p.value
  <chr>          <dbl>     <dbl>     <dbl>       <dbl>
1 (Intercept)    10.0      0.577     17.3  0.000000126
2 TreatmentB     -2.00     0.816     -2.45 0.0400     
3 TreatmentC     -4.00     0.816     -4.90 0.00120    
4 TreatmentD     -6.00     0.816     -7.35 0.0000801  

• Be sure you understand what each parameter means!

• t-statistics and p-values are calculated and interpreted in the same way as the frog example

• But wait, the null hypothesis was \(\mu_A = \mu_B = \mu_C = \mu_D\)

Regression f-statistic

summary(fit1)

Call:
lm(formula = Height ~ Treatment, data = loblollydata)

Residuals:
   Min     1Q Median     3Q    Max 
    -1     -1      0      1      1 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   10.000      0.577   17.32  1.3e-07 ***
TreatmentB    -2.000      0.816   -2.45   0.0400 *  
TreatmentC    -4.000      0.816   -4.90   0.0012 ** 
TreatmentD    -6.000      0.816   -7.35  8.0e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1 on 8 degrees of freedom
Multiple R-squared:  0.882, Adjusted R-squared:  0.838 
F-statistic:   20 on 3 and 8 DF,  p-value: 0.000449

Notice the last line of this output. Where did this p-value come from and what does it mean?

Regression f-statistic

The F-statistic tests whether all regression coefficients (other than the intercept) are all 0

• \(H_0\): All regression coefficients are 0

• \(H_A\): At least one regression coefficient is not 0

F-statistics are the ratio of sample variances

\[F = \frac{s^2_A}{s^2_B}\]

• Null hypothesis is that population variances are equal ( \(\sigma^2_A = \sigma^2_B\) )

• Always positive (variances cannot be negative)

• Usually > 1 (by putting larger variance in the numerator)

Where do the variances come from in the linear model and what do they tell us about differences among groups?

Regression f-statistic

To understand why the test is based on variance, it is helpful to consider several types of means:

• Grand mean

\[\large \bar{y}. = \frac{\sum_i \sum_jy_{ij}}{a \times n}\]

Regression f-statistic

To understand why the test is based on variance, it is helpful to consider several types of means:

• Grand mean

\[\large \bar{y}. = \frac{\sum_i \sum_jy_{ij}}{a \times n}\]

• Group means

\[\large \bar{y}_i = \frac{\sum_j y_{ij}}{n}\]

Regression f-statistic

To understand why the test is based on variance, it is helpful to consider several types of means:

• Grand mean

\[\large \bar{y}. = \frac{\sum_i \sum_jy_{ij}}{a \times n}\]

• Group means

\[\large \bar{y}_i = \frac{\sum_j y_{ij}}{n}\]

We can now decompose the observations as⁴:

\[y_{ij} = \color{#446E9B}{\bar{y}.} + \color{#D47500}{(\bar{y}_i - \bar{y}.)} + \color{#3CB521}{(y_{ij} - \bar{y}_i)}\]

The additive model

The decomposition:

\[\large y_{ij} = \color{#446E9B}{\bar{y}.} + \color{#D47500}{(\bar{y}_i - \bar{y}.)} + \color{#3CB521}{(y_{ij} - \bar{y}_i)}\]

The additive model:

\[\large y_{ij} = \color{#446E9B}{\mu} + \color{#D47500}{\alpha_i} + \color{#3CB521}{\epsilon_{ij}}\]

where

\[\large \epsilon_{ij} \sim normal(0, \sigma^2)\]

The additive model

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

\[\large \epsilon_{ij} \sim normal(0, \sigma^2)\]

Notes:

• \(\mu\) is the grand mean of the population, estimated by \(\bar{y}.\)

• \(\alpha_i\) is the effect of treatment i, estimated by \(\bar{y}_i - \bar{y}.\)

  • It is the deviation of the group mean from the grand mean

  • If all \(\alpha_i = 0\), there is no treatment effect

  • Thus, we can re-write the null hypothesis \(H_0 : \alpha_1 = \alpha_2=... =\alpha_a = 0\)

• \(\epsilon_{ij}\) is the residual error, estimated by \(y_{ij} - \bar{y}_i\)

  • It is the unexplained (random) deviation of the observation from the group mean

• Note that this is just a re-arrangement of the linear model

Sums of squares

Back to the F-test. Where to the variances for the test come from and how do they related to the null hypothesis?

• To answer those questions, we need to quantify where the variation in our data comes from

• To do that, we calculate sums of squares

First, we calculate the sums of squares among, which quantifies variation among treatment groups⁵

\[\large SS_A = n \sum_i (\bar{y}_i - \bar{y}.)^2\]

Motivating example

Question: Is there a difference in growth among the four treatment groups?

Sums of squares

Variation among groups

\[\large SS_A = n \sum_i (\bar{y}_i - \bar{y}.)^2\]

Next, we calculate the sums of squares within, which quantifies variation within treatment groups

\[\large SS_W = \sum_i \sum_j (y_{ij} - \bar{y}_i)^2\]

Motivating example

Question: Is there a difference in growth among the four treatment groups?

Sums of squares

Variation among groups

\[\large SS_A = n \sum_i (\bar{y}_i - \bar{y}.)^2\]

Variation within groups

\[\large SS_W = \sum_i \sum_j (y_{ij} - \bar{y}_i)^2\]

Note that we can also calculate the total variation in our observations⁶

\[\large SS_T = \sum_i \sum_j (y_{ij} - \bar{y}.)^2= SS_A + SS_W\]

Motivating example

Question: Is there a difference in growth among the four treatment groups?

Mean squares

The F-test is based on variances. Sums of squares quantify variation but they are not technically variances.

To covert the sums of squares to variances, divide by the degrees of freedom:⁷

Mean squares among

\[\large MS_A = \frac{SS_A}{a-1}\]

Mean squares within

\[\large MS_W = \frac{SS_W}{a(n-1)}\]

f-statistic

From the mean squares, we can finally calculate the F-statistic:

\[\large F = \frac{MS_A}{MS_W}\]

Note:

• F is a ratio measuring the variance among groups to the variance within groups

• If there is a large treatment effect, what happens to \(MS_A\)?

• If there is little residual variation, what happens to \(MS_W\)?

• Large values of \(F\) indicate treatment effects are large relative to residual variation, but can we conclude there is a treatment effect?

The f-distribution

If the null hypothesis is true (\(\sigma^2_A = \sigma^2_B\)), the ratio of sample variances follows an F-distribution

Properties:

• \(F > 0\)

• \(F\)-distribution is not symmetrical

• Shape of distribution depends on an ordered pair of degrees of freedom, \(df_A\) and \(df_B\)

The f-distribution

Just like the t-distribution, we can use the F-distribution to calculate p-values and test the null hypothesis that the variances are equal

Analysis of variance

Using the F-statistic to test whether there is a treatment effect is commonly referred to as Analysis of Variance (ANOVA)

• Especially in experimental settings, ANOVA is very commonly used

• When treated as an ANOVA, the model is usually presented in a table providing the degrees of freedom, sums of squares, mean squares, F-statistic, and p-value

• The linear model approach is increasingly common, especially in observational settings

• From a practical standpoint, the linear model and ANOVA approaches provide exactly the same information, just presented in different ways (just like the linear model vs t-test we saw previously)

Worked example

Suppose we are interested in the effect of elevation on the abundance of Canada Warblers

	Elevation
Replicate	Low	Medium	High
1	1	2	4
2	3	0	7
3	0	4	5
4	2	3	5

Worked example

data(warblerdata)

fit.lm <- lm(Count ~ Elevation, data = warblerdata)

summary(fit.lm)

Call:
lm(formula = Count ~ Elevation, data = warblerdata)

Residuals:
   Min     1Q Median     3Q    Max 
-2.250 -0.688 -0.250  0.938  1.750 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)        5.250      0.717    7.32  4.5e-05 ***
ElevationLow      -3.750      1.014   -3.70   0.0049 ** 
ElevationMedium   -3.000      1.014   -2.96   0.0160 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.43 on 9 degrees of freedom
Multiple R-squared:  0.63,  Adjusted R-squared:  0.548 
F-statistic: 7.66 on 2 and 9 DF,  p-value: 0.0114

Worked example

data(warblerdata)

fit.aov <- aov(Count ~ Elevation, data = warblerdata)

summary(fit.aov)

            Df Sum Sq Mean Sq F value Pr(>F)  
Elevation    2   31.5   15.75    7.66  0.011 *
Residuals    9   18.5    2.06                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note:

• The row named for the treatment variable (e.g., Elevation) is the among group variation

• The row named Residuals is the within group variation

• The degrees of freedom, F-statistic, and p-value are exactly the same as returned by lm()⁸

Anova table from lm

In case you want more detail about the relationship between lm and aov output:

• lm() function also returns residuals (e.g., \(y_i - E[y_i]\))

fit.lm$residual

    1     2     3     4     5     6     7     8     9    10    11    12 
-0.50  1.50 -1.50  0.50 -0.25 -2.25  1.75  0.75 -1.25  1.75 -0.25 -0.25 

• Residual sum of squares

sum(fit.lm$residuals^2)

[1] 18.5

• Residual mean square

sum(fit.lm$residuals^2)/9

[1] 2.056

Anova table from lm

In case you want more detail about the relationship between lm and aov output:

What about among group variation?

fit.lm$fitted.values

   1    2    3    4    5    6    7    8    9   10   11   12 
1.50 1.50 1.50 1.50 2.25 2.25 2.25 2.25 5.25 5.25 5.25 5.25 

• Treatment sum of squares

sum((fit.lm$fitted.values - mean(fit.lm$fitted.values))^2)

[1] 31.5

• Treatment mean square

sum((fit.lm$fitted.values - mean(fit.lm$fitted.values))^2)/2

[1] 15.75

Anova table from lm

In case you want more detail about the relationship between lm and aov output:

F-statistic and p-value⁹

MSa <- sum((fit.lm$fitted.values - mean(fit.lm$fitted.values))^2)/2
MSe <- sum(fit.lm$residuals^2)/9

(F <- MSa/MSe)

[1] 7.662

(p <- pf(F, 2, 9, lower.tail = FALSE))

[1] 0.0114

Summary

1. Models with categorical variables can be fit using either lm() or aov()

2. The test traditionally known as ANOVA is just a linear model, slightly modified from the one used analyze binary predictor variables

3. R codes categorical predictors using dummy variables and lm coefficients correspond to the difference between the reference level and the other treatment levels

4. lm() or aov() provide the exact same information, just presented in different ways

5. ANOVA output is commonly used for manipulative experiments, linear model output for observational studies¹⁰

6. But how do we tell which groups differ from each other (aside from just the reference level)…

Looking ahead

Next time: Multiple comparisons

Reading: Fieberg chp. 3.9 and 3.12

Footnotes

1. We will discuss this issue in detail in the next lecture

2. Sometimes referred to as the design matrix

3. For a model with \(b\) parameters and \(n\) observations, the model matrix \(X\) is a \(n \times b\) matrix. The parameters themselves form a \(b \times 1\) matrix denoted \(B\). Using matrix multiplication, \(X * B\) returns a \(n \times 1\) matrix with each value corresponding to \(E[y_i]\) for each observation

4. Notice that this is just an algebraic manipulation of the linear model. We’ve re-written the equation but it returns the same value as the previous equation we used

5. Note what this equation is doing: We take the difference between each group mean and the grand mean and square them so there are no negative values (essentially, the length of the orange bars on the previous graph regardless of direction). Then we add them together to get the total length of those lines. Multiplying by \(n\) is basically a way to account for how many observations this variation in based on. The equation here assumes a balanced design (same \(n\) for each group), though \(SS_A\) can be calculated for unbalanced designs

6. And, importantly, that \(SS_A\) and \(SS_W\) partition all of the variation into two categories: variation due to differences among treatments (\(SS_A\)) and variation not due to treatments \(SS_W\)

7. Degrees of freedom are measures of how much information we have, specially how many independent pieces of information we have to estimate parameters. Generally, degrees of freedom are equal to the number of observations minus the number of parameters we need to estimate from those observations. Another way to view degrees of freedom. Imagine you know the mean of 10 observations is \(\bar{y}\). Based on this information, you have no way of knowing the value of the fist 9 observations - they are free to be any value. But once you know the value of observations 1-9, the value of observation 10 is also known - it is not free to vary. So to estimate the mean from \(n=10\), you have \(df = 9\)

8. Which makes sense, because they are from the exact same model, just presented in two different ways

9. The pf() function calculates the area under an F-distribution with \(df1\) and \(df2\) degrees of freedom. lower.tail = FALSE returns only the area under the upper tail, which is what we want since we testing whether \(MS_A\) is greater than \(MS_W\)

10. This is just convention. Since they are the same model, you could use either one regardless of the type of experiment