gmodels (version 2.18.1.1)

fit.contrast: Compute and test arbitrary contrasts for regression objects

Description

Compute and test arbitrary contrasts for regression objects.

Usage

fit.contrast(model, varname, coeff, ... )
# S3 method for lm
fit.contrast(model, varname, coeff, showall=FALSE,
             conf.int=NULL, df=FALSE, ...)
# S3 method for lme
fit.contrast(model, varname, coeff, showall=FALSE,
             conf.int=NULL, df=FALSE, ...)

Value

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.

Arguments

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 TRUE, all model cofficients will be returned. If FALSE (the default), only the coefficients corresponding to the specified contrast will be returned.

conf.int

numeric value on (0,1) or NULL. If a numeric value is specified, confidence intervals with nominal coverage probability conf.int will be computed. If NULL, confidence intervals will not be computed.

df

boolean indicating whether to return a column containing the degrees of freedom.

...

optional arguments provided by methods.

Author

Gregory R. Warnes

Details

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.

References

Venables & Ripley, Section 6.2

See Also

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

Examples

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)