Comparison of age categories using contrasts

contrast(l, a=list(f=levels(d$f)[1:3], x=0), b=list(f=levels(d$f)[4],  
x=0))

On 2/16/2009 10:18 PM, Dylan Beaudette wrote:

On Mon, Feb 16, 2009 at 5:28 PM, Patrick Giraudoux
<patrick.giraudoux at 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

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