[R] Comparison of age categories using contrasts
Dear listers, I would like to compare the levels of a factor with 8 age categories (0,10] (10,20] (20,30] (30,40] (40,50] (50,60] (60,70] (70,90] (however, the factor has not been ordered yet). The default in glm is cont.treatment (for unordered factors) and that leads to compare each level to the first one. I would rather prefer to compare the 2nd to the 1st, the 3rd to the 2nd, the 4th to the 3rd, etc... My understanding is that cont.poly may make the trick, eg specified like this: mod3<-glm(AE~agecat, family=binomial,data=qinghai2, contrasts=list(agecat="contr.poly")) but I am not sure to be right. Would be grateful if a true statistician can confirm or fire me... and before definitive fire tell me how to manage with this... __ 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.
[R] Comparison of age categories using contrasts
Dear listers, I would like to compare the levels of a factor with 8 age categories (0,10] (10,20] (20,30] (30,40] (40,50] (50,60] (60,70] (70,90] (however, the factor has not been ordered yet). The default in glm is cont.treatment (for unordered factors) and that leads to compare each level to the first one. I would rather prefer to compare the 2nd to the 1st, the 3rd to the 2nd, the 4th to the 3rd, etc... My understanding is that cont.poly may make the trick, eg specified like this: mod3<-glm(AE~agecat, family=binomial,data=qinghai2, contrasts=list(agecat="contr.poly")) but I am not sure to be right. Would be grateful if a true statistician can confirm or fire me... and before definitive fire tell me how to manage with this... __ 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.
Re: [R] Comparison of age categories using contrasts
Hi Patrick, >> The default in glm is cont.treatment (for unordered factors) and that >> leads to compare each level to the first one. I would rather prefer to >> compare the 2nd to the 1st, the 3rd to the 2nd, the 4th to the 3rd, >> etc... The functions ?C and ?contrasts allow you to set up your own contrast matrix. A function to carry out the type of comparisons you are interested in has been written by Venables & Ripley. See contr.sdif in package MASS. Regards, Mark. Patrick Giraudoux wrote: > > Dear listers, > > I would like to compare the levels of a factor with 8 age categories > (0,10] (10,20] (20,30] (30,40] (40,50] (50,60] (60,70] (70,90] (however, > the factor has not been ordered yet). The default in glm is > cont.treatment (for unordered factors) and that leads to compare each > level to the first one. I would rather prefer to compare the 2nd to the > 1st, the 3rd to the 2nd, the 4th to the 3rd, etc... My understanding is > that cont.poly may make the trick, eg specified like this: > > mod3<-glm(AE~agecat, family=binomial,data=qinghai2, > contrasts=list(agecat="contr.poly")) > > but I am not sure to be right. > > Would be grateful if a true statistician can confirm or fire me... and > before definitive fire tell me how to manage with this... > > __ > 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-tp22031426p22032444.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.
Re: [R] Comparison of age categories using contrasts
Mark Difford wrote: Hi Patrick, The default in glm is cont.treatment (for unordered factors) and that leads to compare each level to the first one. I would rather prefer to compare the 2nd to the 1st, the 3rd to the 2nd, the 4th to the 3rd, etc... The functions ?C and ?contrasts allow you to set up your own contrast matrix. A function to carry out the type of comparisons you are interested in has been written by Venables & Ripley. See contr.sdif in package MASS. Regards, Mark. Patrick Giraudoux wrote: Dear listers, I would like to compare the levels of a factor with 8 age categories (0,10] (10,20] (20,30] (30,40] (40,50] (50,60] (60,70] (70,90] (however, the factor has not been ordered yet). The default in glm is cont.treatment (for unordered factors) and that leads to compare each level to the first one. I would rather prefer to compare the 2nd to the 1st, the 3rd to the 2nd, the 4th to the 3rd, etc... My understanding is that cont.poly may make the trick, eg specified like this: mod3<-glm(AE~agecat, family=binomial,data=qinghai2, contrasts=list(agecat="contr.poly")) but I am not sure to be right. Would be grateful if a true statistician can confirm or fire me... and before definitive fire tell me how to manage with this... My preference is to not reformulate the model to get desired contrasts but to think of the model as something that uses a convenient coding, then do after-fitting contrasts. In the Design package you can fit a model (say, using glmD) then use contrast(fit, list(predictor settings 1), list(predictor settings 2)) to get what you want. ?contrast.Design provides the details. I'll bet that one of John Fox's packages also provides a good way to do this. Frank -- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University __ 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.
Re: [R] Comparison of age categories using contrasts
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, -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 > -Original Message- > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- > project.org] On Behalf Of Patrick Giraudoux > Sent: Sunday, February 15, 2009 9:23 PM > To: r-h...@stat.math.ethz.ch > Subject: [R] Comparison of age categories using contrasts > > Dear listers, > > I would like to compare the levels of a factor with 8 age categories > (0,10] (10,20] (20,30] (30,40] (40,50] (50,60] (60,70] (70,90] > (however, > the factor has not been ordered yet). The default in glm is > cont.treatment (for unordered factors) and that leads to compare each > level to the first one. I would rather prefer to compare the 2nd to the > 1st, the 3rd to the 2nd, the 4th to the 3rd, etc... My understanding is > that cont.poly may make the trick, eg specified like this: > > mod3<-glm(AE~agecat, family=binomial,data=qinghai2, > contrasts=list(agecat="contr.poly")) > > but I am not sure to be right. > > Would be grateful if a true statistician can confirm or fire me... and > before definitive fire tell me how to manage with this... > > __ > 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. __ 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.
Re: [R] Comparison of age categories using contrasts
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 [[alternative HTML version deleted]] __ 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.
Re: [R] Comparison of age categories using contrasts
On Mon, Feb 16, 2009 at 5:28 PM, Patrick Giraudoux
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.
-0.3449623 d 4.3510407 0.120 3.975904726 4.7261766 23.04 0.
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.
Re: [R] 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
> 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.
> -0.3449623 d 4.3510407 0.120 3.975904726 4.7261766 23.04 0.
>
>
> 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-proje
Re: [R] Comparison of age categories using contrasts
On 2/16/2009 10:18 PM, Dylan Beaudette wrote:
> On Mon, Feb 16, 2009 at 5:28 PM, Patrick Giraudoux
> 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.
> -0.3449623 d 4.3510407 0.120 3.975904726 4.7261766 23.04 0.
>
> 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.
d 4.4308653 0.2731223 3.8884207 4.9733098 16.22 0.
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.
b -4.0635198 0.2782704 -4.6161888 -3.5108507 -14.60 0.
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.
Re: [R] Comparison of age categories using contrasts
Hi Dylan, Chuck,
>> contrast(l, a=list(f=levels(d$f)[1:3], x=0), b=list(f=levels(d$f)[4],
>> x=0))
There is a subtlety here that needs to be emphasized. Setting the
interacting variable (x) to zero is reasonable in this case, because the
mean value of rnorm(n) is zero. However, in the real world of real data it
is somewhat unusual for zero to be a reasonable value for an interacting
variable. In fact, it is this setting-to-zero that causes multcomp's "bells
to ring" with a word of warning.
##
glht(l, linfct=mcp(f="Seq"))$linfct
Regards, Mark.
Chuck Cleland wrote:
>
> On 2/16/2009 10:18 PM, Dylan Beaudette wrote:
>> On Mon, Feb 16, 2009 at 5:28 PM, Patrick Giraudoux
>> 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.
>> -0.3449623 d 4.3510407 0.120 3.975904726 4.7261766 23.04 0.
>>
>> 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.
> d 4.4308653 0.2731223 3.8884207 4.9733098 16.22 0.
>
> 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.
> b -4.0635198 0.2782704 -4.6161888 -3.5108507 -14.60 0.
> 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/mailm
Re: [R] Comparison of age categories using contrasts
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
>> 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)
>> y
