On 3/5/2008 10:09 AM, jebyrnes wrote:
> Huh. Very interesting. I haven't really worked with manipulating contrast
> matrices before, save to do a prior contrasts. Could you explain the matrix
> you laid out just a bit more so that I can generalize it to my case?
Each column corresponds to one of the coefficients in the model, and
each row specifies a particular contrast. The numbers in the matrix
indicate how the model coefficients are combined to indicate a
particular difference in means.
For example, the first row indicates that the third coefficient
(woolB) is multiplied by -1. The baseline categories are A and L for
the wool and tension factors, so the woolB effect in fm is the simple
effect of B vs. A in the baseline category of the tension factor.
Multiplying this coefficient by -1 produces an A vs. B comparison in the
baseline category of the tension factor.
When I want to check that contrasts are as intended, I use contrast()
in the contrast package (by Steve Weston, Jed Wing, & Max Kuhn). That
function allows you to specify factor levels by name to construct the
contrast. For example:
library(contrast)
# M vs. H at B
contrast(fm, a=list(tension = "M", wool = "B"),
b=list(tension = "H", wool = "B"))
lm model parameter contrast
Contrast S.E. Lower Upper t df Pr(>|t|)
10 5.157299 -0.3694453 20.36945 1.94 48 0.0584
It also allows you to print the design matrix for a contrast:
contrast(fm, a=list(tension = "M", wool = "B"),
b=list(tension = "H", wool = "B"))$X
(Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB
1 0 1 -1 0 1 -1
> Chuck Cleland wrote:
>>
>> One approach would be to use glht() in the multcomp package. You
>> need to work out how to formulate the matrix of coefficients that give
>> the desired contrasts. Here is an example using the warpbreaks data
>> frame:
>>
>> fm <- lm(breaks ~ tension*wool, data=warpbreaks)
>>
>> # names(coef(fm))
>> # (Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB
>>
>> cm <- rbind(
>> "A vs. B at L" = c(0, 0, 0,-1, 0, 0),
>> "A vs. B at M" = c(0, 0, 0,-1,-1, 0),
>> "A vs. B at H" = c(0, 0, 0,-1, 0,-1),
>> "M vs. L at A" = c(0, 1, 0, 0, 0, 0),
>> "M vs. H at A" = c(0, 1,-1, 0, 0, 0),
>> "L vs. H at A" = c(0, 0,-1, 0, 0, 0),
>> "M vs. L at B" = c(0, 1, 0, 0, 1, 0),
>> "M vs. H at B" = c(0, 1,-1, 0, 1,-1),
>> "L vs. H at B" = c(0, 0,-1, 0, 0,-1))
>>
>> library(multcomp)
>>
>> summary(glht(fm, linfct = cm), test = adjusted(type="none"))
>>
>> Simultaneous Tests for General Linear Hypotheses
>>
>> Fit: lm(formula = breaks ~ tension * wool, data = warpbreaks)
>>
>> Linear Hypotheses:
>> Estimate Std. Error t value p value
>> A vs. B at L == 0 16.3333 5.1573 3.167 0.002677 **
>> A vs. B at M == 0 -4.7778 5.1573 -0.926 0.358867
>> A vs. B at H == 0 5.7778 5.1573 1.120 0.268156
>> M vs. L at A == 0 -20.5556 5.1573 -3.986 0.000228 ***
>> M vs. H at A == 0 -0.5556 5.1573 -0.108 0.914665
>> L vs. H at A == 0 20.0000 5.1573 3.878 0.000320 ***
>> M vs. L at B == 0 0.5556 5.1573 0.108 0.914665
>> M vs. H at B == 0 10.0000 5.1573 1.939 0.058392 .
>> L vs. H at B == 0 9.4444 5.1573 1.831 0.073270 .
>> ---
>> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>> (Adjusted p values reported -- none method)
--
Chuck Cleland, Ph.D.
NDRI, Inc.
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
______________________________________________
R-help@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.
