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