Comparison of age categories using contrasts
Hi Dylan,
Am I trying to use contrast.Design() for something that it was not intended for? ...
I think Prof. Harrell's main point had to do with how interactions are
handled. You can also get the kind of contrasts that Patrick was interested
in via multcomp. If we do this using your artificial data set we see that
the contrasts are the same as those got by fitting the model using
contr.sdif, but a warning is generated about interactions in the model &c.
[see example code]
Part of Prof. Harrell's "system" is that in generating contrasts via
contr.Design an appropriate value is automatically chosen for the
interacting variable (in this case the median value of x) so that a
reasonable default set of contrasts is calculated. This can of course be
changed.
Coming to your question [?] about how to generate the kind of contrasts that
Patrick wanted using contrast.Design. Well, it is not that straightforward,
though I may have missed something in the documentation to the function. In
the past I have cobbled them together using the kind of hack given below:
## Exampe code
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)))
# now try with contrast.Design():
library(multcomp)
dd <- datadist(d)
options(datadist='dd')
l <- ols(y ~ x * f, data=d)
## model fitted using successive difference contrasts
library(MASS)
T.lm <- lm(y ~ x * f, contrasts=list(f="contr.sdif"), data=d)
summary(T.lm)
## show the warning: model fitted using ols() and treatment contrasts
summary(glht(l, linfct=mcp(f="Seq")))
## "custom" successive difference contrasts using contrast.Design
TFin <- {}
for (i in 1:(length(levels(d$f))-1)) {
tcont <- contrast(l, a=list(f=levels(d$f)[i]),
b=list(f=levels(d$f)[i+1]))
TFin <- rbind(TFin, tcont)
rownames(TFin)[i] <- paste(levels(d$f)[i], levels(d$f)[i+1], sep="-")
}
TFin[,1:9]
Regards, Mark.
Dylan Beaudette-2 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?
Cheers,
Dylan
______________________________________________ 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.
View this message in context: http://www.nabble.com/Comparison-of-age-categories-using-contrasts-tp22031426p22056520.html Sent from the R help mailing list archive at Nabble.com.