1% $Id: fit.contrast.Rd 2184 2018-06-20 20:46:34Z warnes $
2%
3\name{fit.contrast}
4\alias{fit.contrast}
5\alias{fit.contrast.lm}
6\alias{fit.contrast.lme}
7%%\alias{fit.contrast.mer}
8\title{Compute and test arbitrary contrasts for regression objects}
9\description{
10 Compute and test arbitrary contrasts for regression objects.
11}
12\usage{
13fit.contrast(model, varname, coeff, ... )
14\method{fit.contrast}{lm}(model, varname, coeff, showall=FALSE,
15             conf.int=NULL, df=FALSE, ...)
16\method{fit.contrast}{lme}(model, varname, coeff, showall=FALSE,
17             conf.int=NULL, df=FALSE, ...)
18%% \method{fit.contrast}{mer}(model, varname, coeff, showall=FALSE,
19%%             conf.int=NULL, sim.mer = TRUE, n.sim = 1000, ...)
20}
21\arguments{
22  \item{model}{regression (lm,glm,aov,lme) object for which the
23    contrast(s) will be computed.}
24  \item{varname}{variable name}
25  \item{coeff}{vector or matrix specifying contrasts (one per row).}
26  \item{showall}{return all regression coefficients. If \code{TRUE}, all
27    model cofficients will be returned.  If \code{FALSE} (the default),
28    only the coefficients corresponding to the specified contrast will
29    be returned.}
30  \item{conf.int}{numeric value on (0,1) or NULL.  If a numeric value is
31    specified, confidence intervals with nominal coverage probability
32    \code{conf.int} will be computed.  If \code{NULL}, confidence
33    intervals will not be computed.}
34  \item{df}{boolean indicating whether to return a column containing the
35    degrees of freedom.}
36  \item{\dots}{optional arguments provided by methods.}
37%%  \item{sim.mer}{Logical value. If TRUE p-values and confidence
38%%    intervals will be estimated using \code{mcmcsamp}. This option only takes effect for mer
39%%    objects.}
40%%  \item{n.sim}{Number of samples to use in \code{mcmcsamp}.}
41  }
42
43\details{
44  Computes the specified contrast(s) by re-fitting the model with the
45  appropriate arguments.  A contrast of the form \code{c(1,0,0,-1)}
46  would compare the mean of the first group with the mean of the fourth group.
47}
48\value{
49  Returns a matrix containing estimated coefficients, standard errors, t
50  values, two-sided p-values. If \code{df} is TRUE, an additional column
51  containing the degrees of freedom is included.  If \code{conf.int} is
52  specified lower and upper confidence limits are also returned.}
53\references{Venables & Ripley, Section 6.2}
54
55\author{ Gregory R. Warnes \email{greg@warnes.net}}
56
57\seealso{ \code{\link{lm}}, \code{\link{contrasts}},
58  \code{\link{contr.treatment}},  \code{\link{contr.poly}},
59  Computation and testing of General Linear Hypothesis:
60  \code{\link{glh.test}}, Computation and testing of estimable functions
61  of model coefficients: \code{\link{estimable}}, \code{\link{make.contrasts}}
62  }
63
64\examples{
65y <- rnorm(100)
66x <-  cut(rnorm(100, mean=y, sd=0.25),c(-4,-1.5,0,1.5,4))
67reg <- lm(y ~ x)
68summary(reg)
69
70# look at the group means
71gm <- sapply(split(y,x),mean)
72gm
73
74
75# mean of 1st group vs mean of 4th group
76fit.contrast(reg, x, c(    1,    0,    0,   -1) )
77# estimate should be equal to:
78gm[1] - gm[4]
79
80# mean of 1st and 2nd groups vs mean of 3rd and 4th groups
81fit.contrast(reg, x, c( -1/2, -1/2,  1/2,  1/2) )
82# estimate should be equal to:
83sum(-1/2*gm[1], -1/2*gm[2], 1/2*gm[3], 1/2*gm[4])
84
85# mean of 1st group vs mean of 2nd, 3rd and 4th groups
86fit.contrast(reg, x, c( -3/3,  1/3,  1/3,  1/3) )
87# estimate should be equal to:
88sum(-3/3*gm[1], 1/3*gm[2], 1/3*gm[3], 1/3*gm[4])
89
90# all at once
91cmat <- rbind( "1 vs 4"    =c(-1, 0, 0, 1),
92               "1+2 vs 3+4"=c(-1/2,-1/2, 1/2, 1/2),
93               "1 vs 2+3+4"=c(-3/3, 1/3, 1/3, 1/3))
94fit.contrast(reg,x,cmat)
95
96#
97x2 <- rnorm(100,mean=y,sd=0.5)
98reg2 <- lm(y ~ x + x2 )
99fit.contrast(reg2,x,c(-1,0,0,1))
100
101#
102# Example for Analysis of Variance
103#
104
105set.seed(03215)
106Genotype <- sample(c("WT","KO"), 1000, replace=TRUE)
107Time <- factor(sample(1:3, 1000, replace=TRUE))
108y <- rnorm(1000)
109data <- data.frame(y, Genotype, Time)
110
111
112# Compute Contrasts & obtain 95\% confidence intervals
113
114model <- aov( y ~ Genotype + Time + Genotype:Time, data=data )
115
116fit.contrast( model, "Genotype", rbind("KO vs WT"=c(-1,1) ), conf=0.95 )
117
118fit.contrast( model, "Time",
119            rbind("1 vs 2"=c(-1,1,0),
120                  "2 vs 3"=c(0,-1,1)
121                  ),
122            conf=0.95 )
123
124
125cm.G <- rbind("KO vs WT"=c(-1,1) )
126cm.T <- rbind("1 vs 2"=c(-1,1,0),
127              "2 vs 3"=c(0,-1,1) )
128
129# Compute contrasts and show SSQ decompositions
130
131model <- aov( y ~ Genotype + Time + Genotype:Time, data=data,
132              contrasts=list(Genotype=make.contrasts(cm.G),
133                             Time=make.contrasts(cm.T) )
134            )
135
136summary(model, split=list( Genotype=list( "KO vs WT"=1 ),
137                           Time = list( "1 vs 2" = 1,
138                                        "2 vs 3" = 2 ) ) )
139
140
141# example for lme
142library(nlme)
143data(Orthodont)
144fm1 <- lme(distance ~ Sex, data = Orthodont,random=~1|Subject)
145
146# Contrast for sex.  This example is equivalent to standard treatment
147# contrast.
148#
149fit.contrast(fm1, "Sex", c(-1,1), conf.int=0.95 )
150#
151# and identical results can be obtained using lme built-in 'intervals'
152#
153intervals(fm1)
154
155# Cut age into quantile groups & compute some contrasts
156Orthodont$AgeGroup <- gtools::quantcut(Orthodont$age)
157fm2 <- lme(distance ~ Sex + AgeGroup, data = Orthodont,random=~1|Subject)
158#
159fit.contrast(fm2, "AgeGroup", rbind("Linear"=c(-2,-1,1,2),
160                                    "U-Shaped"=c(-1,1,1,-1),
161                                    "Change-Point at 11"=c(-1,-1,1,1)),
162                              conf.int=0.95)
163
164
165}
166\keyword{ models }
167\keyword{ regression }
168
169