[R] Contrasts in lm

[BXC (Bendix Carstensen)]

Just yesterday I wrote the following function. But you have to set up the 
contrast matrix for yourself and make sure it is sensible.

ci.mat <- 
function( alpha = 0.05 ) 
{
# Gives 1 2 x 3 matrix to multiply onto the first two
# colums of coefficients to give c.i.s
# BxC, 10/00
rbind( c(1,1,1), qnorm(1-alpha/2)*c(0,-1,1) )
}  

contr <- 
function( obj, cm, alpha=0.05 )
{
# Function to compute arbitrary contrasts with c.i.
# from a linear model ( lm, glm or nlme )
# BxC, 12/01.
if ( "lme" %in% class( obj ) ) 
  {
  cf  <- summary(obj)$tTable
  rho <- summary(obj)$cor  
  vcv <- rho * outer(cf[,2],cf[,2])
  }
if ( "lm" %in% class( obj ) )
  {
  cf  <- summary(obj)$coefficients
  vcv <- summary(obj)$cov.unscaled  
  }
if ( !dim( cm )[2]==dim( cf )[1] ) stop(  
   paste( "\n Dimension of ",
          deparse( substitute( cm ) ), ": ", paste( dim(cm), collapse="x" ),
          ", not compatible with no of parameters in ",
          deparse( substitute( obj ) ), ": ", dim(cf)[1], sep="" ) )
ct <- cm %*% cf[,1]
vc <- sqrt( diag( cm %*% vcv %*% t(cm) ) )
cbind( ct, vc ) %*% ci.mat( alpha=alpha )
}

----------------------
Bendix Carstensen
Senior Statistician
Steno Diabetes Centre
Niels Steensens Vej 2
DK-2820 Gentofte
Denmark
tel: +45 44 43 87 38
mob: +45 28 25 87 38
fax: +45 44 43 73 13
bxc at novonordisk.com
www.biostat.ku.dk/~bxc
----------------------



> -----Original Message-----
> From: Søren Højsgaard [mailto:Soren.Hojsgaard at agrsci.dk]
> Sent: 6. december 2001 10:13
> To: r-help at stat.math.ethz.ch
> Subject: [R] Contrasts in lm
> 
> 
> Dear all,
> 
> In SAS (GLM and MIXED) estimable functions (linear functions of the
> parameters) can be specified in the ESTIMATE and CONTRAST statements.
> 
> Has anyone written a similar "utility" for use in connection with lm?
> 
> Thanks in advance
> 
> Søren Højsgaard 
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
> -.-.-.-.-.-.-.-.-
> 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
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
> _._._._._._._._._
> 
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._