[R] testing for significantly different slopes

[Gregory Warnes]

>> What's the problem?
>
> The problem is that I would like to do a pair-wise comparison  
> between the
> multiple slopes. For example with this model:
>
> lm1 <- lm (Sepal.Length ~ Species/Sepal.Width -1, data=iris)
>
> # truncated output from summary(lm1)
> # just the slope terms
> Speciessetosa:Sepal.Width       0.6905     0.1657   4.166 5.31e-05 ***
> Speciesversicolor:Sepal.Width   0.8651     0.2002   4.321 2.88e-05 ***
> Speciesvirginica:Sepal.Width    0.9015     0.1948   4.628 8.16e-06 ***
>
> If I wanted to test the hypothesis that Speciessetosa:Sepal.Width was
> different than Speciesversicolor:Sepal.Width, what course of action  
> should I
> take?
>
> I have found an example in the gmodels package, using the estimable()
> function, but the documentation is not clear to me.

The gmodels::estimable function will allow you to test any hypothesis  
constructed from the model coefficients.  Each row of the contrast  
matrix 'cm' is applied to the vector of coefficients 'beta'  
individually and tested (by default) against the null:

	cm[i] %*% beta == 0

This is accomplished using t-tests using the appropriate  
transformation of the model parameter covariance matrix to determine  
the correct standard deviation.

> Here is the example:
>
> library(gmodels)
>
> # example from manual
> lm1 <- lm (Sepal.Length ~ Sepal.Width + Species + Sepal.Width:Species,
> data=iris)
>
> cm <- rbind(
> 'Setosa vs. Versicolor'   = c(0, 0, 1, 0, 1, 0),
> 'Setosa vs. Virginica'    = c(0, 0, 0, 1, 0, 1),
> 'Versicolor vs. Virginica'= c(0, 0, 1,-1, 1,-1)
> )
>
> estimable(lm1, cm)
>
> This *appears* to test what I am after, but I am not clear on how  
> the 'cm'
> argument works.


This example code tests whether the intercept *and* slope are equal  
across species, with a little extra  complexity because the  
species=Setosa is represeted as the 'baseline' group that is included  
in the intercept term.

To test whether just the slope  Speciesversicolor:Sepal.Width is  
equal to the slope Speciesvirginica:Sepal.Width, you can do a single  
contrast row:

 > estimable(lm1, c(0, 0, 0, 0, 1, -1) )
                   Estimate Std. Error    t value  DF  Pr(>|t|)
(0 0 0 0 1 -1) -0.03645676  0.2793277 -0.1305161 144 0.8963403

For the one-contrast case, one can use the shortcut:

 > estimable(lm1, c("Speciesversicolor:Sepal.Width"=1,  
"Speciesvirginica:Sepal.Width"=-1) )
                   Estimate Std. Error    t value  DF  Pr(>|t|)
(0 0 0 0 1 -1) -0.03645676  0.2793277 -0.1305161 144 0.8963403

Now, to each of the pairwise tests in a single call, construct a  
matrix with one row per contrast, and one column per model parameter:

 > cm <- rbind(
+ 'Setosa vs. Versicolor Slope'   = c(0, 0, 0,  1, -1, 0),
+ 'Setosa vs. Virginica Slope'    = c(0, 0, 0,  1,  0, -1),
+ 'Versicolor vs. Virginica Slope'= c(0, 0, 0,  0,  1, -1)
+ )
 >
 > colnames(cm) <- names(coef(lm1))
 >
 > #print(cm) # omitted here for brevity, but instructive
 >
 > estimable(lm1, cm)
                                   Estimate Std. Error    t value   
DF  Pr(>|t|)
Setosa vs. Versicolor Slope    -0.17458800  0.2598919 -0.6717717 144  
0.5028054
Setosa vs. Virginica Slope     -0.21104476  0.2557557 -0.8251811 144  
0.4106337
Versicolor vs. Virginica Slope -0.03645676  0.2793277 -0.1305161 144  
0.8963403

I hope this helps.

-Greg