Compute and test arbitrary contrasts for regression objects
Compute and test arbitrary contrasts for regression objects.
fit.contrast(model, varname, coeff, ... )
## S3 method for class 'lm'
fit.contrast(model, varname, coeff, showall=FALSE,
conf.int=NULL, df=FALSE, ...)
## S3 method for class 'lme'
fit.contrast(model, varname, coeff, showall=FALSE,
conf.int=NULL, df=FALSE, ...)model |
regression (lm,glm,aov,lme) object for which the contrast(s) will be computed. |
varname |
variable name |
coeff |
vector or matrix specifying contrasts (one per row). |
showall |
return all regression coefficients. If |
conf.int |
numeric value on (0,1) or NULL. If a numeric value is
specified, confidence intervals with nominal coverage probability
|
df |
boolean indicating whether to return a column containing the degrees of freedom. |
... |
optional arguments provided by methods. |
Computes the specified contrast(s) by re-fitting the model with the
appropriate arguments. A contrast of the form c(1,0,0,-1)
would compare the mean of the first group with the mean of the fourth group.
Returns a matrix containing estimated coefficients, standard errors, t
values, two-sided p-values. If df is TRUE, an additional column
containing the degrees of freedom is included. If conf.int is
specified lower and upper confidence limits are also returned.
Gregory R. Warnes greg@warnes.net
Venables & Ripley, Section 6.2
lm, contrasts,
contr.treatment, contr.poly,
Computation and testing of General Linear Hypothesis:
glh.test, Computation and testing of estimable functions
of model coefficients: estimable, make.contrasts
y <- rnorm(100)
x <- cut(rnorm(100, mean=y, sd=0.25),c(-4,-1.5,0,1.5,4))
reg <- lm(y ~ x)
summary(reg)
# look at the group means
gm <- sapply(split(y,x),mean)
gm
# mean of 1st group vs mean of 4th group
fit.contrast(reg, x, c( 1, 0, 0, -1) )
# estimate should be equal to:
gm[1] - gm[4]
# mean of 1st and 2nd groups vs mean of 3rd and 4th groups
fit.contrast(reg, x, c( -1/2, -1/2, 1/2, 1/2) )
# estimate should be equal to:
sum(-1/2*gm[1], -1/2*gm[2], 1/2*gm[3], 1/2*gm[4])
# mean of 1st group vs mean of 2nd, 3rd and 4th groups
fit.contrast(reg, x, c( -3/3, 1/3, 1/3, 1/3) )
# estimate should be equal to:
sum(-3/3*gm[1], 1/3*gm[2], 1/3*gm[3], 1/3*gm[4])
# all at once
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))
fit.contrast(reg,x,cmat)
#
x2 <- rnorm(100,mean=y,sd=0.5)
reg2 <- lm(y ~ x + x2 )
fit.contrast(reg2,x,c(-1,0,0,1))
#
# Example for Analysis of Variance
#
set.seed(03215)
Genotype <- sample(c("WT","KO"), 1000, replace=TRUE)
Time <- factor(sample(1:3, 1000, replace=TRUE))
y <- rnorm(1000)
data <- data.frame(y, Genotype, Time)
# Compute Contrasts & obtain 95% confidence intervals
model <- aov( y ~ Genotype + Time + Genotype:Time, data=data )
fit.contrast( model, "Genotype", rbind("KO vs WT"=c(-1,1) ), conf=0.95 )
fit.contrast( model, "Time",
rbind("1 vs 2"=c(-1,1,0),
"2 vs 3"=c(0,-1,1)
),
conf=0.95 )
cm.G <- rbind("KO vs WT"=c(-1,1) )
cm.T <- rbind("1 vs 2"=c(-1,1,0),
"2 vs 3"=c(0,-1,1) )
# Compute contrasts and show SSQ decompositions
model <- aov( y ~ Genotype + Time + Genotype:Time, data=data,
contrasts=list(Genotype=make.contrasts(cm.G),
Time=make.contrasts(cm.T) )
)
summary(model, split=list( Genotype=list( "KO vs WT"=1 ),
Time = list( "1 vs 2" = 1,
"2 vs 3" = 2 ) ) )
# example for lme
library(nlme)
data(Orthodont)
fm1 <- lme(distance ~ Sex, data = Orthodont,random=~1|Subject)
# Contrast for sex. This example is equivalent to standard treatment
# contrast.
#
fit.contrast(fm1, "Sex", c(-1,1), conf.int=0.95 )
#
# and identical results can be obtained using lme built-in 'intervals'
#
intervals(fm1)
# Cut age into quantile groups & compute some contrasts
Orthodont$AgeGroup <- gtools::quantcut(Orthodont$age)
fm2 <- lme(distance ~ Sex + AgeGroup, data = Orthodont,random=~1|Subject)
#
fit.contrast(fm2, "AgeGroup", rbind("Linear"=c(-2,-1,1,2),
"U-Shaped"=c(-1,1,1,-1),
"Change-Point at 11"=c(-1,-1,1,1)),
conf.int=0.95)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.