Re: [R] testing for significantly different slopes

From: Gregory Warnes <gregory.warnes_at_mac.com>
Date: Wed, 05 Mar 2008 22:51:02 -0500

>> 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



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Thu 06 Mar 2008 - 09:29:53 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 06 Mar 2008 - 09:30:19 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive