(If you use this information, please cite it.)

Introduction

This file shows how to fit and use linear models combining categorical and continuous predictors. Of focal importance is the presence of slopes that may differ among levels of the grouping factor.

We use primarily two packages in addition to base and stats: car and multcomp

Install those now if they are not available. Click on the Packages tab and then on Install. You can also run install.packages(c("car", "multcomp"))

Filaree plants were grown in three different temperature ranges (categorical groups) over a number of days (continuous predictor), and their weight was measured periodically and log transformed (response or dependent variable). The study had three goals:

  1. Determine if temperature affects the relative growth rate of filaree,
  2. estimate the mean mass of plants at different ages, and
  3. compare the mass of the plants at different ages within and across groups.

Under the assumption of constant relative growth rate in early growth stages,

SLOPE = RELATIVE GROWTH RATE

In a later post I will show that the data have enough information to reject the linear model, but that is another story.

filaree <- read.csv("xmpl_filaree.csv", header=TRUE)

filaree
##    Tgroup day    lnwt day.22.4
## 1       1   9 2.76312    -13.4
## 2       1   9 2.83220    -13.4
## 3       1   9 2.80917    -13.4
## 4       1  15 3.40785     -7.4
## 5       1  15 3.66113     -7.4
## 6       1  15 3.40785     -7.4
## 7       1  22 4.21376     -0.4
## 8       1  22 4.02955     -0.4
## 9       1  22 4.09863     -0.4
## 10      1  30 4.95059      7.6
## 11      1  30 4.69730      7.6
## 12      1  30 4.74336      7.6
## 13      1  37 5.31901     14.6
## 14      1  37 5.22690     14.6
## 15      1  37 5.15782     14.6
## 16      2   8 3.10851    -14.4
## 17      2   8 3.01641    -14.4
## 18      2   8 3.01641    -14.4
## 19      2  18 4.21376     -4.4
## 20      2  18 4.32889     -4.4
## 21      2  18 4.05258     -4.4
## 22      2  23 4.88181      0.6
## 23      2  23 4.65125      0.6
## 24      2  23 4.81243      0.6
## 25      2  30 5.15782      7.6
## 26      2  30 5.29598      7.6
## 27      2  30 5.29598      7.6
## 28      2  36 5.52624     13.6
## 29      2  36 5.48019     13.6
## 30      2  36 5.48019     13.6
## 31      3   7 2.99338    -15.4
## 32      3   7 3.01641    -15.4
## 33      3   7 3.13154    -15.4
## 34      3  14 4.23678     -8.4
## 35      3  14 4.25981     -8.4
## 36      3  14 4.09863     -8.4
## 37      3  22 4.74336     -0.4
## 38      3  22 4.83546     -0.4
## 39      3  22 4.74336     -0.4
## 40      3  29 5.24993      6.6
## 41      3  29 5.27295      6.6
## 42      3  29 5.36506      6.6
## 43      3  36 5.89466     13.6
## 44      3  36 5.81558     13.6
## 45      3  36 5.82558     13.6
str(filaree)
## 'data.frame':    45 obs. of  4 variables:
##  $ Tgroup  : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ day     : int  9 9 9 15 15 15 22 22 22 30 ...
##  $ lnwt    : num  2.76 2.83 2.81 3.41 3.66 ...
##  $ day.22.4: num  -13.4 -13.4 -13.4 -7.4 -7.4 -7.4 -0.4 -0.4 -0.4 7.6 ...
filaree <- filaree %>%
  mutate(group = factor(Tgroup)) %>%
  select(group, day, lnwt)

str(filaree)
## 'data.frame':    45 obs. of  3 variables:
##  $ group: Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ day  : int  9 9 9 15 15 15 22 22 22 30 ...
##  $ lnwt : num  2.76 2.83 2.81 3.41 3.66 ...
summary(filaree)
##  group       day            lnwt      
##  1:15   Min.   : 7.0   Min.   :2.763  
##  2:15   1st Qu.:14.0   1st Qu.:3.661  
##  3:15   Median :22.0   Median :4.697  
##         Mean   :22.4   Mean   :4.425  
##         3rd Qu.:30.0   3rd Qu.:5.250  
##         Max.   :37.0   Max.   :5.895
scatterplot(lnwt ~ day | group, filaree, smooth = FALSE)

options(contrasts =c(“contr.treatment”, “contr.poly”))

options(contrasts =c("contr.treatment", "contr.poly"))

fullm <- lm(lnwt ~ group * day, filaree)

summary(fullm)  # meaning of parameters?
## 
## Call:
## lm(formula = lnwt ~ group * day, data = filaree)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30047 -0.12657 -0.02163  0.12080  0.32725 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.130009   0.112944  18.859  < 2e-16 ***
## group2      0.385378   0.163615   2.355  0.02363 *  
## group3      0.522179   0.155178   3.365  0.00173 ** 
## day         0.086632   0.004566  18.972  < 2e-16 ***
## group2:day  0.002028   0.006585   0.308  0.75970    
## group3:day  0.005034   0.006372   0.790  0.43431    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1778 on 39 degrees of freedom
## Multiple R-squared:  0.9689, Adjusted R-squared:  0.965 
## F-statistic: 243.4 on 5 and 39 DF,  p-value: < 2.2e-16
anova(fullm)
## Analysis of Variance Table
## 
## Response: lnwt
##           Df Sum Sq Mean Sq   F value    Pr(>F)    
## group      2  2.600   1.300   41.1424 2.463e-10 ***
## day        1 35.832  35.832 1133.8739 < 2.2e-16 ***
## group:day  2  0.020   0.010    0.3168    0.7303    
## Residuals 39  1.232   0.032                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow = c(2,2)) # make 2 rows and rwo columns of graphs in one page

plot(fullm, col = filaree$group)  # LATER! Ignore lack fo fit for now.

Equality of slopes is not rejected, so we should stop the test of slope here. We continue to show what to do if the equality of all slopes had been rejected

Compare slopes

We make all possible comparisons, so we need to adjust \(\alpha\) to control family-wise error rate. We use the Bonferroni correction, which is very conservative: real \(\alpha\) is at most nominal \(\alpha\).

We use a “trick” to get the coefficients to extract the comparisons we want.

First we get the coefficients (L-vectors) that give use the intercepts for each group and the predicted lnwt for each group when day = 1, which is the sum of intercept and slope for each group. Then we remove the intercepts from the estimated lnwt for day = 1 and the result is a set of coefficients to extract each slope.

(data4mmat <- expand.grid(group = unique(filaree$group), day = c(0, 1)))
##   group day
## 1     1   0
## 2     2   0
## 3     3   0
## 4     1   1
## 5     2   1
## 6     3   1
(mym <- model.matrix(~ group * day, data4mmat, contrasts = fullm$contrasts))
##   (Intercept) group2 group3 day group2:day group3:day
## 1           1      0      0   0          0          0
## 2           1      1      0   0          0          0
## 3           1      0      1   0          0          0
## 4           1      0      0   1          0          0
## 5           1      1      0   1          1          0
## 6           1      0      1   1          0          1
## attr(,"assign")
## [1] 0 1 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
mym[4:6, 1:3] <- 0

mym # see the final matrix of L-vectors for all parameters
##   (Intercept) group2 group3 day group2:day group3:day
## 1           1      0      0   0          0          0
## 2           1      1      0   0          0          0
## 3           1      0      1   0          0          0
## 4           0      0      0   1          0          0
## 5           0      0      0   1          1          0
## 6           0      0      0   1          0          1
## attr(,"assign")
## [1] 0 1 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"

Now, we use the rows of the matrix to create L-vectors that compare slopes.

L4B11 <- mym[4, ]
L4B12 <- mym[5, ]
L4B13 <- mym[6, ]

L4.B11.B12 <- L4B11 - L4B12 # compares Bhat11 vs Bhat12
L4.B11.B13 <- L4B11 - L4B13
L4.B12.B13 <- L4B12 - L4B13

Put together the L-vectors for comparisons in a hypotheses matrix and submit to multcomp.

H0mat <- rbind(B11hat_B12hat =  L4.B11.B12,
              B11hat_B13hat =   L4.B11.B13,
              B12hat_B13hat =   L4.B12.B13)

tH0 <- glht(fullm, linfct = H0mat)
summary(tH0, test = adjusted("bonferroni")) # No RGR's are different from each other.
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = lnwt ~ group * day, data = filaree)
## 
## Linear Hypotheses:
##                     Estimate Std. Error t value Pr(>|t|)
## B11hat_B12hat == 0 -0.002028   0.006585  -0.308        1
## B11hat_B13hat == 0 -0.005034   0.006372  -0.790        1
## B12hat_B13hat == 0 -0.003006   0.006501  -0.462        1
## (Adjusted p values reported -- bonferroni method)

Other tests

Most of the questions that we are usually interested in are implemented as tests of the hypothesis that certain linear combinations of parameters are zero, greater than or equal to zero or less than or equal to zero. The last two options are more useful, because they allow us to implement the idea of Practically Significant Difference.

We are also interested in making confidence intervals for estimated quantities that usually are linear combinations of estimated parameters.

The column vector that R prints under the “Estimate” column of a summary is the vector of estimated parameters that we have to operate with and implement tests of hypotheses.

For now, we will ignore some hypotheses or questions that involve nonlinear functions of estimated parameters, like finding an X for which the expected values of Y is a given number. For example a person trying to decide how long it will take plants to reach market size may want to estimate the age at which plants reach a certain mass.)

Suppose that we want to determine if mean mass of plants in group 3 on day 25 reaches a required lnwt = 5.

We hypothesize \(H_0: \mu_{3, Y|X = 25} \leq 5\) v. \(H_A: \mu_{3, Y|X = 25} \gt 5\)

A null hypothesis that an estimated value is equal to a fixed number is tested with a confidence interval. Because the hypothesis is one-sided, we put all p-value in one side of the distribution. In this case, if the 90% interval is all above 5 we reject the null.

We estimate \(\mu_{3, Y|X = 25}\) with \(\hat{\beta}_{03} + 25 \ \hat{\beta}_{13}\)


\[\hat{\mu}_{3, Y|X = 25} = \hat{\beta}_{03} + 25 \ \hat{\beta}_{13} = L_{\beta_{03}} \ \mathbf{\hat{\beta}}+ L_{\beta_{13}} \ \hat{\mathbf{\beta}}\]
where \(\hat{\mathbf{\beta}}\) is the vector of estimated parameters that R gives us (coef(modelname)).

# First, make L vectors for the intercepts

L4B01 <- mym[1, ]
L4B02 <- mym[2, ]
L4B03 <- mym[3, ]

# Estimate the value on the line for day 25

mu.hat.3.25 <- L4B03 %*% coef(fullm) + 25 * L4B13 %*% coef(fullm)

# factor out beta.hat to get the L-vector for the estimate

L4.3.25 <- L4B03 + 25 * L4B13

# 90% CI

dfe <- anova(fullm)$Df[4]
tval <- qt(p = 0.95, df = dfe)
s.e. <- sqrt(L4.3.25 %*% vcov(fullm) %*% L4.3.25)
mu.hat.3.25 + s.e. * tval
##          [,1]
## [1,] 5.025247
mu.hat.3.25 - s.e. * tval
##          [,1]
## [1,] 4.862411
# Use glht

H0mat2 <- rbind(lnwtg3day25 =  L4.3.25)

tH02 <- glht(fullm, linfct = H0mat2, alternative = "greater", rhs = 5)
summary(tH02, test = adjusted("bonferroni")) # No RGR's are different from each other.
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = lnwt ~ group * day, data = filaree)
## 
## Linear Hypotheses:
##                  Estimate Std. Error t value Pr(>t)
## lnwtg3day25 <= 5  4.94383    0.04832  -1.162  0.874
## (Adjusted p values reported -- bonferroni method)

options(contrasts =c(“contr.sum”, “contr.poly”))

Change the parameterization to contr.sum. This means that instead of calculating the parameters as deviations from the first level, R calculates them as deviations from the overall mean of parameter values. The constraint is that the sum of the deviations or “effects” have to be zero.

Change contrasts option and rerun model.

options(contrasts =c("contr.sum", "contr.poly"))

fullmSum <- lm(lnwt ~ group * day, filaree)

summary(fullmSum)  # meaning of numbers printed is different from before
## 
## Call:
## lm(formula = lnwt ~ group * day, data = filaree)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30047 -0.12657 -0.02163  0.12080  0.32725 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.4325281  0.0650587  37.390  < 2e-16 ***
## group1      -0.3025192  0.0921127  -3.284  0.00217 ** 
## group2       0.0828589  0.0943598   0.878  0.38526    
## day          0.0889857  0.0026481  33.603  < 2e-16 ***
## group1:day  -0.0023541  0.0037367  -0.630  0.53237    
## group2:day  -0.0003258  0.0038099  -0.086  0.93228    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1778 on 39 degrees of freedom
## Multiple R-squared:  0.9689, Adjusted R-squared:  0.965 
## F-statistic: 243.4 on 5 and 39 DF,  p-value: < 2.2e-16
anova(fullmSum) # Note that anova is identical, regardless of contrast type.
## Analysis of Variance Table
## 
## Response: lnwt
##           Df Sum Sq Mean Sq   F value    Pr(>F)    
## group      2  2.600   1.300   41.1424 2.463e-10 ***
## day        1 35.832  35.832 1133.8739 < 2.2e-16 ***
## group:day  2  0.020   0.010    0.3168    0.7303    
## Residuals 39  1.232   0.032                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow = c(2,2)) # make 2 rows and rwo columns of graphs in one page

plot(fullmSum, col = filaree$group)  # LATER!

Note that instead of showing groups 2 and 3 it now shows groups 1 and 2. Estimated parameters are in the order they appear in the Estimate column:

\[\hat{\beta_0} \\ \Delta \hat{\beta}_{01} \\ \Delta \hat{\beta}_{02} \\ \hat{\beta}_1 \\ \Delta \hat{\beta}_{11} \\ \Delta \hat{\beta}_{12}\]

where, using i to represent the group number, \(\hat{\beta}_{0i}\) is the intercept of group i, \(\hat{\beta}_{1i}\) is the slope of group i and

\[\hat{\beta}_0 = \sum_{i = 1}^{i = 3} \hat{\beta}_{0i} / 3 \\[20pt] \hat{\beta}_1 = \sum_{i = 1}^{i = 3} \hat{\beta}_{1i} / 3 \\[20pt] \hat{\beta}_{0i} = \hat{\beta}_0 + \Delta \hat{\beta}_{0i} \\[20pt] \hat{\beta}_{1i} = \hat{\beta}_1 + \Delta \hat{\beta}_{1i}\]

A result of these definitions to keep in mind is that because of the first and second equations in the last set above, the following is true:

\[\Delta \hat{\beta}_{03} = -(\Delta \hat{\beta}_{01} + \Delta \hat{\beta}_{02}) \\[20pt] \Delta \hat{\beta}_{13} = -(\Delta \hat{\beta}_{11} + \Delta \hat{\beta}_{12}) \\[20pt] \text{thus, } \qquad \quad \hat{\beta}_{13} = \hat{\beta}_1 - \Delta \hat{\beta}_{11} - \Delta \hat{\beta}_{12}\]

(data4mmatSum <- expand.grid(group = unique(filaree$group), day = c(0, 1)))
##   group day
## 1     1   0
## 2     2   0
## 3     3   0
## 4     1   1
## 5     2   1
## 6     3   1
(mymSum <- model.matrix(~ group * day, data4mmatSum, contrasts = fullmSum$contrasts))
##   (Intercept) group1 group2 day group1:day group2:day
## 1           1      1      0   0          0          0
## 2           1      0      1   0          0          0
## 3           1     -1     -1   0          0          0
## 4           1      1      0   1          1          0
## 5           1      0      1   1          0          1
## 6           1     -1     -1   1         -1         -1
## attr(,"assign")
## [1] 0 1 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.sum"
mymSum[4:6, 1:3] <- 0

mymSum # see the final matrix of L-vectors for all parameters
##   (Intercept) group1 group2 day group1:day group2:day
## 1           1      1      0   0          0          0
## 2           1      0      1   0          0          0
## 3           1     -1     -1   0          0          0
## 4           0      0      0   1          1          0
## 5           0      0      0   1          0          1
## 6           0      0      0   1         -1         -1
## attr(,"assign")
## [1] 0 1 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.sum"

Now, we use the rows of the matrix to create L-vectors that compare slopes.

L4B11s <- mymSum[4, ]
L4B12s <- mymSum[5, ]
L4B13s <- mymSum[6, ]

L4.B11.B12s <- L4B11s - L4B12s # compares Bhat11 vs Bhat12
L4.B11.B13s <- L4B11s - L4B13s
L4.B12.B13s <- L4B12s - L4B13s

Put together the L-vectors for comparisons in a hypotheses matrix and submit to multcomp.

H0matSum <- rbind(B11hat_B12hat =  L4.B11.B12s,
              B11hat_B13hat =   L4.B11.B13s,
              B12hat_B13hat =   L4.B12.B13s)

tH0Sum <- glht(fullmSum, linfct = H0matSum)
summary(tH0Sum, test = adjusted("bonferroni")) # No RGR's are different from each other.
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = lnwt ~ group * day, data = filaree)
## 
## Linear Hypotheses:
##                     Estimate Std. Error t value Pr(>|t|)
## B11hat_B12hat == 0 -0.002028   0.006585  -0.308        1
## B11hat_B13hat == 0 -0.005034   0.006372  -0.790        1
## B12hat_B13hat == 0 -0.003006   0.006501  -0.462        1
## (Adjusted p values reported -- bonferroni method)

We can see that the results are the same as before. All the time, regardless of the options(), the estimated betas are the same, what differs is what R shows and uses for the calculations.

Differences in expected plant sizes

The last part type of test that we may be interested in is a comparison of expected plant sizes for different groups at different times. Imagine that the greenhouse operator calculated that the cost of growing plants in the high temperature for 15 days is the same as growing them in the low temperature for 23 days. From all other points of view, it is preferable to have many short cycles. Market research indicates that plants that are within 0.25 units of log weight do not differ in price. The question we need to answer is: are plants grown 15 days in high temperature more than 0.2 units smaller than those grown for 23 days in high temperature?

\[H_0: \qquad \quad \mu_{1, Y|X = 23} - \mu_{3, Y|X = 15} \ge 0.20 \\[15pt] H_A: \qquad \quad \mu_{1, Y|X = 23} - \mu_{3, Y|X = 15} \lt 0.20\]

First we get the L-vector to estimate the difference and then we submit it to glht to test if the estimate is significantly greater than 0.2. For clarity, we use the symbol \(L_{01}\) to represent the L-vector to get the estimated intercept for group 1, etc., \(L_{H_0}\) to represent the L-vector that estimates the difference we are testing, and \(\mathbf{\beta}\).

\[\begin{aligned} \hat{\mu}_{1, Y|X = 23} - \hat{\mu}_{3, Y|X = 15} &= L_{01} \ \mathbf{\hat{\beta}} + 23 \ L_{11} \ \mathbf{\hat{\beta}} - L_{03} \ \mathbf{\hat{\beta}} + 15 \ L_{13} \ \mathbf{\hat{\beta}} \\[10pt] &=(L_{01} + 23 \ L_{11} - L_{03} + 15 \ L_{13}) \ \mathbf{\hat{\beta}} \\[20pt] \text{thus, } \qquad L_{H_0} &= L_{01} + 23 \ L_{11} - L_{03} + 15 \ L_{13} \end{aligned}\]

L4B01s <- mymSum[1, ] # define the L-vectors for the intercepts
L4B02s <- mymSum[2, ]
L4B03s <- mymSum[3, ]


L4H0 <- L4B01s + 23 * L4B11s - L4B03s - 15 * L4B13s

H0matSizeDiff <- rbind(mu1.23.mu3.15 =  L4H0)

tH0SizeDiff <- glht(fullmSum, linfct = H0matSizeDiff, alternative = "less", rhs = 0.25)
summary(tH0SizeDiff, test = adjusted("bonferroni")) # Reject H0
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = lnwt ~ group * day, data = filaree)
## 
## Linear Hypotheses:
##                       Estimate Std. Error t value Pr(<t)  
## mu1.23.mu3.15 >= 0.25  0.09536    0.07126   -2.17 0.0181 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)

The null hypothesis is rejected, which means that there is enough information in the data to make a good guess that plants grown for 15 days in high temperature as at least as heavy as those grown in low temperature for 23 days. The grower concludes that it is better to have shorter cycles with hotter greenhouses.