On 2/16/2009 10:18 PM, Dylan Beaudette wrote:
> On Mon, Feb 16, 2009 at 5:28 PM, Patrick Giraudoux
> <patrick.giraud...@univ-fcomte.fr> wrote:
>> Greg Snow a écrit :
>>> One approach is to create your own contrasts matrix:
>>>
>>>
>>>> mycmat <- diag(8)
>>>> mycmat[ row(mycmat) == col(mycmat) + 1 ] <- -1
>>>> mycmati <- solve(mycmat)
>>>> contrasts(agefactor) <- mycmati[,-1]
>>>>
>>> Now when you use agefactor, the intercept will be the first age group and
>>> the slopes will be the differences between the pairs of groups (make sure
>>> that the order of the levels of agefactor is correct).
>>>
>>> The difference between this method and the contr.sdif function in MASS is
>>> how the intercept will end up being interpreted (and the dimnames).
>>>
>>> Hope this helps,
>>>
>>>
>> Actually, exactly what I needed including the reference to contr.sdif in
>> MASS I did not spot before (although I am a faithful reader of the
>> yellow book... but so many things still escape to me). Again thanks a lot.
>>
>> Patrick
>>
>
> Glad you were able to solve your problem. Frank replied earlier with
> the suggestion to use contrast.Design() to perform the tests after the
> initial model has been fit. Perhaps I am a little denser than normal,
> but I cannot figure out how to apply contrast.Design() to a similar
> model with several levels of a factor. Example data:
>
> # need these
> library(lattice)
> library(Design)
>
> # replicate an important experimental dataset
> set.seed(10101010)
> x <- rnorm(100)
> y1 <- x[1:25] * 2 + rnorm(25, mean=1)
> y2 <- x[26:50] * 2.6 + rnorm(25, mean=1.5)
> y3 <- x[51:75] * 2.9 + rnorm(25, mean=5)
> y4 <- x[76:100] * 3.5 + rnorm(25, mean=5.5)
>
> d <- data.frame(x=x, y=c(y1,y2,y3,y4), f=factor(rep(letters[1:4], each=25)))
>
> # plot
> xyplot(y ~ x, groups=f, data=d, auto.key=list(columns=4, title='Beard
> Type', lines=TRUE, points=FALSE, cex=0.75), type=c('p','r'),
> ylab='Number of Pirates', xlab='Distance from Land')
>
> # standard comparison to base level of f
> summary(lm(y ~ x * f, data=d))
>
> # compare to level 4 of f
> summary(lm(y ~ x * C(f, base=4), data=d))
>
> # now try with contrast.Design():
> dd <- datadist(d)
> options(datadist='dd')
> l <- ols(y ~ x * f, data=d)
>
> # probably the wrong approach, as the results do not look too familiar:
> contrast(l, list(f=levels(d$f)))
>
> x f Contrast S.E. Lower Upper t Pr(>|t|)
> -0.3449623 a 0.3966682 0.1961059 0.007184856 0.7861515 2.02 0.0460
> -0.3449623 b 0.5587395 0.1889173 0.183533422 0.9339456 2.96 0.0039
> -0.3449623 c 4.1511677 0.1889170 3.775962254 4.5263730 21.97 0.0000
> -0.3449623 d 4.3510407 0.1888820 3.975904726 4.7261766 23.04 0.0000
>
> This is adjusting the output to a given value of 'x'... Am I trying to
> use contrast.Design() for something that it was not intended for? Any
> tips Frank?
Does this help?
# Compare with summary(lm(y ~ x * f, data=d))
contrast(l, a=list(f=levels(d$f)[2:4], x=0),
b=list(f=levels(d$f)[1], x=0))
f Contrast S.E. Lower Upper t Pr(>|t|)
b 0.3673455 0.2724247 -0.1737135 0.9084046 1.35 0.1808
c 4.1310015 0.2714011 3.5919754 4.6700275 15.22 0.0000
d 4.4308653 0.2731223 3.8884207 4.9733098 16.22 0.0000
Error d.f.= 92
# Compare with summary(lm(y ~ x * C(f, base=4), data=d))
contrast(l, a=list(f=levels(d$f)[1:3], x=0),
b=list(f=levels(d$f)[4], x=0))
f Contrast S.E. Lower Upper t Pr(>|t|)
a -4.4308653 0.2731223 -4.9733098 -3.8884207 -16.22 0.0000
b -4.0635198 0.2782704 -4.6161888 -3.5108507 -14.60 0.0000
c -0.2998638 0.2772684 -0.8505427 0.2508151 -1.08 0.2823
Error d.f.= 92
> Cheers,
> Dylan
>
> ______________________________________________
> 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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
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.
