contrast.Design {Design}  R Documentation 
This function computes one or more contrasts of the estimated
regression coefficients in a fit from one of the functions in Design,
along with standard errors, confidence limits, t or Z statistics, Pvalues.
General contrasts are handled by obtaining the design matrix for two
sets of predictor settings (a
, b
) and subtracting the
corresponding rows of the two design matrics to obtain a new contrast
design matrix for testing the a
 b
differences. This allows for
quite general contrasts (e.g., estimated differences in means between
a 30 year old female and a 40 year old male).
This can also be used
to obtain a series of contrasts in the presence of interactions (e.g.,
female:male log odds ratios for several ages when the model contains
age by sex interaction). Another use of contrast
is to obtain
centerweighted (Type III test) and subjectweighted (Type II test)
estimates in a model containing treatment by center interactions. For
the latter case, you can specify type="average"
and an optional
weights
vector to average the withincenter treatment contrasts.
The design contrast matrix computed by contrast.Design
can be used
by the bootplot
and confplot
functions to obtain bootstrap
nonparametric confidence intervals for contrasts.
By omitting the b
argument, contrast
can be used to obtain
an average or weighted average of a series of predicted values, along
with a confidence interval for this average. This can be useful for
"unconditioning" on one of the predictors (see the next to last
example).
Specifying type="joint"
, and specifying at least as many contrasts
as needed to span the space of a complex test, one can make
multiple degree of freedom tests flexibly and simply. Redundant
contrasts will be ignored in the joint test. See the examples below.
These include an example of an "incomplete interaction test" involving
only two of three levels of a categorical variable (the test also tests
the main effect).
When more than one contrast is computed, the list created by
contrast.Design
is suitable for plotting (with error bars or bands)
with xYplot
or Dotplot
(see the last example before the
type="joint"
examples).
contrast(fit, ...) ## S3 method for class 'Design': contrast(fit, a, b, cnames=NULL, type=c("individual", "average", "joint"), weights="equal", conf.int=0.95, tol=1e7, ...) ## S3 method for class 'contrast.Design': print(x, X=FALSE, fun=function(u)u, jointonly=FALSE, ...)
fit 
a fit of class "Design"

a 
a list containing settings for all predictors that you do not wish to
set to default (adjustto) values. Usually you will specify two
variables in this list, one set to a constant and one to a sequence of
values, to obtain contrasts for the sequence of values of an
interacting factor. The gendata function will generate the
necessary combinations and default values for unspecified predictors.

b 
another list that generates the same number of observations as a ,
unless one of the two lists generates only one observation. In that
case, the design matrix generated from the shorter list will have its
rows replicated so that the contrasts assess several differences
against the one set of predictor values. This is useful for comparing
multiple treatments with control, for example. If b is missing, the
design matrix generated from a is analyzed alone.

cnames 
vector of character strings naming the contrasts when
type!="average" . Usually cnames is not necessary as
contrast.Design tries to name the contrasts by examining which
predictors are varying consistently in the two lists. cnames will
be needed when you contrast "noncomparable" settings, e.g., you
compare list(treat="drug", age=c(20,30)) with
list(treat="placebo"), age=c(40,50))

type 
set type="average" to average the individual contrasts (e.g., to
obtain a Type II or III contrast). Set type="joint" to jointly
test all nonredundant contrasts with a multiple degree of freedom test
and no averaging.

weights 
a numeric vector, used when type="average" , to obtain weighted contrasts

conf.int 
confidence level for confidence intervals for the contrasts 
tol 
tolerance for qr function for determining which
contrasts are redundant, and for inverting the covariance matrix
involved in a joint test 
... 
unused 
x 
result of contrast 
X 
set X=TRUE to print design matrix used in computing the contrasts (or
the average contrast)

fun 
a function to transform the contrast, SE, and lower and upper
confidence limits before printing. For example, specify fun=exp to
antilog them for logistic models.

jointonly 
set to FALSE to omit printing of individual
contrasts 
a list of class "contrast.Design"
containing the elements
Contrast
, SE
, Z
, var
, df.residual
Lower
, Upper
, Pvalue
, X
, cnames
, redundant
, which denote the contrast
estimates, standard errors, Z or tstatistics, variance matrix,
residual degrees of freedom (this is NULL
if the model was not
ols
), lower and upper confidence limits, 2sided Pvalue, design
matrix, contrast names (or NULL
), and a logical vector denoting
which contrasts are redundant with the other contrasts. If there are
any redundant contrasts, when the results of contrast
are
printed, and asterisk is printed at the start of the corresponding lines.
Frank Harrell
Department of Biostatistics
Vanderbilt University School of Medicine
f.harrell@vanderbilt.edu
predict.Design
, gendata
, bootcov
, summary.Design
,
anova.Design
, plot.Design
, qr
set.seed(1) age < rnorm(200,40,12) sex < factor(sample(c('female','male'),200,TRUE)) logit < (sex=='male') + (age40)/5 y < ifelse(runif(200) <= plogis(logit), 1, 0) f < lrm(y ~ pol(age,2)*sex) # Compare a 30 year old female to a 40 year old male # (with or without age x sex interaction in the model) contrast(f, list(sex='female', age=30), list(sex='male', age=40)) # For a model containing two treatments, centers, and treatment # x center interaction, get 0.95 confidence intervals separately # by cente center < factor(sample(letters[1:8],500,TRUE)) treat < factor(sample(c('a','b'), 500,TRUE)) y < 8*(treat=='b') + rnorm(500,100,20) f < ols(y ~ treat*center) lc < levels(center) contrast(f, list(treat='b', center=lc), list(treat='a', center=lc)) # Get 'Type III' contrast: average b  a treatment effect over # centers, weighting centers equally (which is almost always # an unreasonable thing to do) contrast(f, list(treat='b', center=lc), list(treat='a', center=lc), type='average') # Get 'Type II' contrast, weighting centers by the number of # subjects per center. Print the design contrast matrix used. k < contrast(f, list(treat='b', center=lc), list(treat='a', center=lc), type='average', weights=table(center)) print(k, X=TRUE) # Note: If other variables had interacted with either treat # or center, we may want to list settings for these variables # inside the list()'s, so as to not use default settings # For a 4treatment study, get all comparisons with treatment 'a' treat < factor(sample(c('a','b','c','d'), 500,TRUE)) y < 8*(treat=='b') + rnorm(500,100,20) dd < datadist(treat,center); options(datadist='dd') f < ols(y ~ treat*center) lt < levels(treat) contrast(f, list(treat=lt[1]), list(treat=lt[ 1]), cnames=paste(lt[1],lt[1],sep=':'), conf.int=1.05/3) # Compare each treatment with average of all others for(i in 1:length(lt)) { cat('Comparing with',lt[i],'\n\n') print(contrast(f, list(treat=lt[i]), list(treat=lt[ i]), type='average')) } options(datadist=NULL) # Six ways to get the same thing, for a variable that # appears linearly in a model and does not interact with # any other variables. We estimate the change in y per # unit change in a predictor x1. Methods 4, 5 also # provide confidence limits. Method 6 computes nonparametric # bootstrap confidence limits. Methods 26 can work # for models that are nonlinear or nonadditive in x1. # For that case more care is needed in choice of settings # for x1 and the variables that interact with x1. ## Not run: coef(fit)['x1'] # method 1 diff(predict(fit, gendata(x1=c(0,1)))) # method 2 g < Function(fit) # method 3 g(x1=1)  g(x1=0) summary(fit, x1=c(0,1)) # method 4 k < contrast(fit, list(x1=1), list(x1=0)) # method 5 print(k, X=TRUE) fit < update(fit, x=TRUE, y=TRUE) # method 6 b < bootcov(fit, B=500, coef.reps=TRUE) bootplot(b, X=k$X) # bootstrap distribution and CL # In a model containing age, race, and sex, # compute an estimate of the mean response for a # 50 year old male, averaged over the races using # observed frequencies for the races as weights f < ols(y ~ age + race + sex) contrast(f, list(age=50, sex='male', race=levels(race)), type='average', weights=table(race)) ## End(Not run) # Plot the treatment effect (drug  placebo) as a function of age # and sex in a model in which age nonlinearly interacts with treatment # for females only set.seed(1) n < 800 treat < factor(sample(c('drug','placebo'), n,TRUE)) sex < factor(sample(c('female','male'), n,TRUE)) age < rnorm(n, 50, 10) y < .05*age + (sex=='female')*(treat=='drug')*.05*abs(age50) + rnorm(n) f < ols(y ~ rcs(age,4)*treat*sex) d < datadist(age, treat, sex); options(datadist='d') # show separate estimates by treatment and sex plot(f, age=NA, treat=NA, sex='female') plot(f, age=NA, treat=NA, sex='male') ages < seq(35,65,by=5); sexes < c('female','male') w < contrast(f, list(treat='drug', age=ages, sex=sexes), list(treat='placebo', age=ages, sex=sexes)) xYplot(Cbind(Contrast, Lower, Upper) ~ age  sex, data=w, ylab='Drug  Placebo') xYplot(Cbind(Contrast, Lower, Upper) ~ age, groups=sex, data=w, ylab='Drug  Placebo', method='alt bars') options(datadist=NULL) # Examples of type='joint' contrast tests set.seed(1) x1 < rnorm(100) x2 < factor(sample(c('a','b','c'), 100, TRUE)) dd < datadist(x1, x2); options(datadist='dd') y < x1 + (x2=='b') + rnorm(100) # First replicate a test statistic from anova() f < ols(y ~ x2) anova(f) contrast(f, list(x2=c('b','c')), list(x2='a'), type='joint') # Repeat with a redundancy; compare a vs b, a vs c, b vs c contrast(f, list(x2=c('a','a','b')), list(x2=c('b','c','c')), type='joint') # Get a test of association of a continuous predictor with y # First assume linearity, then cubic f < lrm(y>0 ~ x1 + x2) anova(f) contrast(f, list(x1=1), list(x1=0), type='joint') # a minimum set of contrasts xs < seq(2, 2, length=20) contrast(f, list(x1=0), list(x1=xs), type='joint') # All contrasts were redundant except for the first, because of # linearity assumption f < lrm(y>0 ~ pol(x1,3) + x2) anova(f) contrast(f, list(x1=0), list(x1=xs), type='joint') print(contrast(f, list(x1=0), list(x1=xs), type='joint'), jointonly=TRUE) # All contrasts were redundant except for the first 3, because of # cubic regression assumption # Now do something that is difficult to do without cryptic contrast # matrix operations: Allow each of the three x2 groups to have a different # shape for the x1 effect where x1 is quadratic. Test whether there is # a difference in mean levels of y for x2='b' vs. 'c' or whether # the shape or slope of x1 is different between x2='b' and x2='c' regardless # of how they differ when x2='a'. In other words, test whether the mean # response differs between group b and c at any value of x1. # This is a 3 d.f. test (intercept, linear, quadratic effects) and is # a better approach than subsetting the data to remove x2='a' then # fitting a simpler model, as it uses a better estimate of sigma from # all the data. f < ols(y ~ pol(x1,2) * x2) anova(f) contrast(f, list(x1=xs, x2='b'), list(x1=xs, x2='c'), type='joint') # Note: If using a spline fit, there should be at least one value of # x1 between any two knots and beyond the outer knots.