[R] contrasts in lm

[Warnes, Gregory R]

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