Hi Dylan, Chuck, Mark Difford wrote: >> 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:
Here is a much simpler route, and a correction to my earlier posting (helped by a little off-list prompting from Prof. Harrell). All that is required to get successive difference (i.e. sdif) contrasts using contrast.Design is the following: contrast(l, a=list(f=c("a","b","c")), b=list(f=c("b","c","d"))) ## Run a full example 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))) dd <- datadist(d) options(datadist='dd') l <- ols(y ~ x * f, data=d) ## compare with multcomp by setting x=0 summary(glht(l, linfct=mcp(f="Seq"))) contrast(l, b=list(x=0, f=c("a","b","c")), a=list(x=0, f=c("b","c","d"))) Regards, and apologies for any confusion, Mark. Mark Difford wrote: > > 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.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? >> >> 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. >> >> > > -- View this message in context: http://www.nabble.com/Comparison-of-age-categories-using-contrasts-tp22031426p22060776.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ 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.