I've been playing around with contrasts in lm by specifying the contrasts argument. So, I want to specify a specific contrast to be tested Say: > y _ rnorm(100) > x _ cut(rnorm(100, mean=y, sd=0.25),c(-3,-1.5,0,1.5,3)) > reg _ lm(y ~ x, contrasts=list(x=c(1,0,0,-1))) > coef(reg)[2] x1 -1.814101 I was surprised to see that I get a different estimate for the contrast specified this way than if I do > reg2 _ lm(y ~ x, contrasts=list(x=contr.sum)) > coef(reg2)[2] x1 -1.816682 even though the first contrast in the same in both cases: > reg$contrasts $x [,1] [,2] [,3] (-3,-1.5] 1 -0.2697521 0.42099143 (-1.5,0] 0 -0.3256196 -0.80247857 (0,1.5] 0 0.8651239 -0.03950429 (1.5,3] -1 -0.2697521 0.42099143 > reg2$contrasts $x [,1] [,2] [,3] (-3,-1.5] 1 0 0 (-1.5,0] 0 1 0 (0,1.5] 0 0 1 (1.5,3] -1 -1 -1 Questions: 1) How are columns 2 and 3 of the first contrast matrix being created and why? 2) What is the proper way of doing custom contrasts? I'm interested in a 'general' solution. Thanks, Greg LEGAL NOTICE Unless expressly stated otherwise, this message is confidential and may be privileged. It is intended for the addressee(s) only. Access to this E-mail by anyone else is unauthorized. If you are not an addressee, any disclosure or copying of the contents of this E-mail or any action taken (or not taken) in reliance on it is unauthorized and may be unlawful. If you are not an addressee, please inform the sender immediately. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Thu, 30 Aug 2001, Warnes, Gregory R wrote:> I've been playing around with contrasts in lm by specifying the contrasts > argument. So, I want to specify a specific contrast to be tested > > Say: > > > y _ rnorm(100) > > x _ cut(rnorm(100, mean=y, sd=0.25),c(-3,-1.5,0,1.5,3)) > > reg _ lm(y ~ x, contrasts=list(x=c(1,0,0,-1))) > > coef(reg)[2] > x1 > -1.814101 > > I was surprised to see that I get a different estimate for the contrast > specified this way than if I doThat's not an estimate of the contrast, just of one element of the contrast. The full answer is something like> coef(reg)(Intercept) x1 x2 x3 0.09832456 -1.84557466 0.74762184 0.34896487> > > reg2 _ lm(y ~ x, contrasts=list(x=contr.sum)) > > coef(reg2)[2] > x1 > -1.816682 > > even though the first contrast in the same in both cases: > > > reg$contrasts > $x > [,1] [,2] [,3] > (-3,-1.5] 1 -0.2697521 0.42099143 > (-1.5,0] 0 -0.3256196 -0.80247857 > (0,1.5] 0 0.8651239 -0.03950429 > (1.5,3] -1 -0.2697521 0.42099143 > > reg2$contrasts > $x > [,1] [,2] [,3] > (-3,-1.5] 1 0 0 > (-1.5,0] 0 1 0 > (0,1.5] 0 0 1 > (1.5,3] -1 -1 -1Why were you surprised? If I do lm(y ~ x + z) lm(y ~ x + w) for continuous x and different continuous z and w I should not expect the same regression coefficient for x. That's what you've done. The exception is where z and w are orthogonal to x, but columns of the model matrix are not in general orthogonal, especially with unbalanced designs (as here) and several of the standard contrasts functions (as in your second example).> Questions: > > 1) How are columns 2 and 3 of the first contrast matrix being created and > why?by contrasts<-(), to complete a span of the space. It uses orthogonal terms,> 2) What is the proper way of doing custom contrasts? I'm interested in a > 'general' solution.function C(). I suspect you actually intended lm(y ~ C(x, contr.sum, 1)) and confused yourself by not looking at the whole answer. I do recommend Bill Venables' account of coding in chapter 6 of MASS: there seem to be several misconceptions to clear up here. One thing to watch is that R's default, contr.treatment, are *not contrasts* in the conventional meaning. -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Thanks Brian,
Reading MASS 6.2 really helped. I've now written a function, contrast.lm,
that does what I want:
> y _ rnorm(100)
> x _ cut(rnorm(100, mean=y, sd=0.25),c(-3,-1.5,0,1.5,3))
> reg _ lm(y ~ x, contrasts=list(x=contr.sum))
>
> # contrast mean of 1st group vs mean of 4th group
> contrast.lm(reg, x, c(1,0,0,-1) )
Estimate Std. Error t value Pr(>|t|)
x c=( 1 0 0 -1 ) -3.615228 0.2841731 -12.72192 0
> # estimate should be equal to:
> gm[1] - gm[4]
(-3,-1.5]
-3.615228
> contrast.lm(reg,x,cmat)
Estimate Std. Error t value Pr(>|t|)
x1 vs 4 3.615228 0.2841731 12.72192 0
x1+2 vs 3+4 2.500153 0.1504516 16.61766 0
x1 vs 2+3+4 2.563423 0.2406548 10.65187 0
>
A primary purpose for writing this function is to let my colleagues are used
to SAS specify contrasts in a familiar way. I'm attaching the code to
contrast.lm() below in case it might be useful to someone else. I'll also
include it in the next version of the gregmisc library. Comments or
suggestions welcome.
-Greg
-----R code below this line--------
# Function to compute (and test) arbitrary contrasts for regression objects
# Usage:
#
# y _ rnorm(100)
# x _ cut(rnorm(100, mean=y, sd=0.25),c(-3,-1.5,0,1.5,3))
# reg _ lm(y ~ x, contrasts=list(x=contr.sum))#
#
# contrast.lm(reg, x, c(1,0,0,-1) ) # one contrast
#
# cmat <- rbind( "1 vs 4" =c(-1 , 0 , 0 , 1 ),
# "1+2 vs 3+4"=c(-1/2,-1/2, 1/2, 1/2),
# "1 vs 2+3+4"=c(-3/3, 1/3, 1/3, 1/3))
# contrast.lm(reg,x,cmat) # set of contrasts
#
contrast.lm _ function(reg, varname, coeff)
{
# make sure we have the NAME of the variable
if(!is.character(varname))
varname <- deparse(substitute(varname))
# make coeff into a matrix
if(!is.matrix(coeff))
coeff <- matrix(coeff, nrow=1)
# make sure columns are labeled
if (is.null(rownames(coeff)))
{
rn <- vector(length=nrow(coeff))
for(i in 1:nrow(coeff))
rn[i] <- paste(" c=(",paste(coeff[i,],collapse="
"), ")")
rownames(coeff) <- rn
}
# now convert into the proper form for the contrast matrix
cmat <- ginv(coeff) # requires library(MASS)
colnames(cmat) <- rownames(coeff)
# call lm with the specified contrast
m <- reg$call
if(is.null(m$contrasts))
m$contrasts <- list()
m$contrasts[[varname]] <- cmat
r <- eval(m, parent.frame())
# now return the correct elements ....
retval <- coef(summary(r))
rn <- paste(varname,rownames(coeff),sep="")
ind <- match(rn,rownames(retval))
retval[ind,,drop=F]
}
LEGAL NOTICE
Unless expressly stated otherwise, this message is confidential and may be
privileged. It is intended for the addressee(s) only. Access to this E-mail by
anyone else is unauthorized. If you are not an addressee, any disclosure or
copying of the contents of this E-mail or any action taken (or not taken) in
reliance on it is unauthorized and may be unlawful. If you are not an addressee,
please inform the sender immediately.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._