effect {effects} | R Documentation |
Effect
and effect
construct an "eff"
object for a term (usually a high-order term)
in a linear model (fit by lm
or gls
) or generalized linear model (fit by glm
), or an
"effpoly"
object for a term in a
multinomial or proportional-odds logit model (fit respectively by multinom
or polr
),
absorbing the lower-order terms marginal
to the term in question, and averaging over other terms in the model. For multivariate linear models
(of class "mlm"
, fit by lm
), the function constructs a list of "eff"
objects separately for the various
response variables.
effect
builds the required object by specifying explicity a focal term like "a:b"
for an a
by b
interaction. Effect
specifies the predictors in
the term, for example c("a", "b")
, rather than the term itself. Effect
is consequently
more flexible and robust than
effect
, and will succeed with some models for which effect
fails. The effect
function
works by constructing a call to Effect
.
The Effect
and effect
functions can also be used with many other models; see Effect.default
and a vignette for this package.
allEffects
identifies all of the high-order terms in a model and returns
a list of "eff"
or "effpoly"
objects (i.e., an object of type "efflist"
).
For information on computing and displaying predictor effects, see predictorEffect
and plot.predictoreff
.
For further information about plotting effects, see plot.eff
.
effect(term, mod, vcov.=vcov, ...) ## Default S3 method: effect(term, mod, vcov.=vcov, ...) Effect(focal.predictors, mod, ...) ## S3 method for class 'lm' Effect(focal.predictors, mod, xlevels=list(), fixed.predictors, vcov. = vcov, se=TRUE, transformation = list(link = family(mod)$linkfun, inverse = family(mod)$linkinv), residuals=FALSE, quantiles=seq(0.2, 0.8, by=0.2), x.var=NULL, ..., #legacy arguments: given.values, typical, offset, confint, confidence.level, partial.residuals) ## S3 method for class 'multinom' Effect(focal.predictors, mod, xlevels=list(), fixed.predictors, vcov. = vcov, se=TRUE, ..., #legacy arguments: confint, confidence.level, given.values, typical) ## S3 method for class 'polr' Effect(focal.predictors, mod, xlevels=list(), fixed.predictors, vcov.=vcov, se=TRUE, latent=FALSE, ..., #legacy arguments: confint, confidence.level, given.values, typical) ## S3 method for class 'svyglm' Effect(focal.predictors, mod, fixed.predictors, ...) allEffects(mod, ...) ## Default S3 method: allEffects(mod, ...) ## S3 method for class 'eff' as.data.frame(x, row.names=NULL, optional=TRUE, transform=x$transformation$inverse, ...) ## S3 method for class 'effpoly' as.data.frame(x, row.names=NULL, optional=TRUE, ...) ## S3 method for class 'efflatent' as.data.frame(x, row.names=NULL, optional=TRUE, ...) ## S3 method for class 'eff' vcov(object, ...)
term |
the quoted name of a term, usually, but not necessarily, a high-order
term in the model. The term must be given exactly as it appears in the printed
model, although either colons ( |
focal.predictors |
a character vector of one or more predictors in the model in any order. |
mod |
an object of the appropriate class. If no method exists for that class, |
xlevels |
this argument is used to set the number of levels for any
focal predictor that is not a factor. If |
fixed.predictors |
an optional list of specifications affecting the values at which fixed predictors for an effect are set, potentially including:
|
vcov. |
A function or the name of a function that will be used to get the estimated variance-covariance
matrix of the estimated coefficients. This will ordinarily be the default,
|
se |
|
transformation |
a two-element list with elements |
residuals |
if |
quantiles |
quantiles at which to evaluate numeric focal predictors not on the
horizontal axis, used only when partial residuals are displayed; superceded if the |
x.var |
the name or index of the numeric predictor to define the horizontal axis of an effect
plot for a linear or generalized linear model; the default is |
latent |
if |
x |
an object of class |
transform |
a transformation to be applied to the effects and confidence limits,
by default taken from the inverse link function saved in the |
row.names, optional |
not used. |
object |
an object of class |
... |
arguments to be passed down. |
confint, confidence.level, given.values, typical, offset, partial.residuals |
legacy arguments retained for backwards compatability; if present, these arguments take precedence
over |
Normally, the functions to be used directly are allEffects
, to return
a list of high-order effects, and the generic plot
function to plot the effects.
(see plot.efflist
, plot.eff
, and plot.effpoly
).
Alternatively, Effect
can be used to vary a subset of predictors over their ranges,
while other predictors are held to typical values.
Plots are drawn using the xyplot
(or in some cases,
the densityplot
) function in the
lattice package. Effects may also be printed (implicitly or explicitly via
print
) or summarized (using summary
)
(see print.efflist
, summary.efflist
,
print.eff
, summary.eff
, print.effpoly
, and summary.effpoly
).
If asked, the effect
function will compute effects for terms that have
higher-order relatives in the model, averaging over those terms (which rarely makes sense), or for terms that
do not appear in the model but are higher-order relatives of terms that do.
For example, for the model Y ~ A*B + A*C + B*C
, one could
compute the effect corresponding to the absent term A:B:C
, which absorbs the constant, the
A
, B
, and C
main effects, and the three two-way interactions. In either of these
cases, a warning is printed.
The as.data.frame
methods convert effect objects to data frames to facilitate the construction
of custom displays. In the case of "eff"
objects, the se
element in the data frame is always
on the scale of the linear predictor, and the transformation used for the fit and confidence limits is saved in
a "transformation"
attribute.
See predictorEffects
for an alternative paradigm for getting effects.
For lm
, glm
,svyglm
, mer
and lme
, effect
and Effect
return
an "eff"
object, and for multinom
,
polr
, clm
, clmm
and clm2
, an "effpoly"
object, with the components listed below.
For an "mlm"
object with one response specified, an "eff"
object is returned, otherwise an "efflist"
object
is returned, containing one "eff"
object for each response
.
term |
the term to which the effect pertains. |
formula |
the complete model formula. |
response |
a character string giving the name of the response variable. |
y.levels |
(for |
variables |
a list with information about each predictor, including its name, whether it is a factor, and its levels or values. |
fit |
(for |
prob |
(for |
logit |
(for |
x |
a data frame, the columns of which are the predictors in the effect, and the rows of which give all combinations of values of these predictors. |
model.matrix |
the model matrix from which the effect was calculated. |
data |
a data frame with the data on which the fitted model was based. |
discrepancy |
the percentage discrepancy for the ‘safe’ predictions of the original fit; should be very close to 0.
Note: except for |
offset |
value to which the offset is fixed; |
model |
(for |
vcov |
(for |
se |
(for |
se.prob, se.logit |
(for |
lower, upper |
(for |
lower.prob, upper.prob, lower.logit, upper.logit |
(for |
confidence.level |
for the confidence limits. |
transformation |
(for |
residuals |
(working) residuals for linear or generalized linear models, to be used by
|
x.var |
the name of the predictor to appear on the horizontal axis of an effect plot made from the
returned object; will usually be |
family |
for a |
allEffects
returns an "efflist"
object, a list of "eff"
or "effpoly"
objects
corresponding to the high-order terms of the model.
If mod
is of class "poLCA"
(from the poLCA
package), representing a
polytomous latent class model, effects are computed for the predictors given the
estimated latent classes.
The result is of class "eff"
if the latent class model has 2 categories
and of class "effpoly"
with more than 2 categories.
The Effect
function handles factors and covariates differently, and is likely to be confused if one is changed to the other
in a model formula. Consequently, formulas that include calls to as.factor
, factor
, or numeric
(as, e.g., in y ~ as.factor(income)
) will cause errors. Instead, create the modified variables outside of the model
formula (e.g., fincome <- as.factor(income)
) and use these in the model formula.
Factors cannot have colons in level names (e.g., "level:A"
); the effect
function will confuse the
colons with interactions; rename levels to remove or replace the colons (e.g., "level.A"
).
The functions in the effects package work properly with predictors that are numeric or factors; consequently, e.g., convert logical predictors to factors, and dates to numeric.
Empty cells in crossed-factors are now permitted for "lm"
, "glm"
, and "multinom"
models.
For "multinom"
models with two or more crossed factors with an empty cell, stacked area plots
apparently do not work because of a bug in the barchart
function in the lattice package. However, the default
line plots do work.
Offsets in linear and generalized linear models are supported, as are offsets in mixed models fit by
lmer
or glmer
, but must be supplied through the offset
argument to lm
, glm
, lmer
or glmer
;
offsets supplied via calls to the offset
function on the right-hand side
of the model formula are not supported.
Fitting ordinal mixed-models using clmm
or clmm2
permits many options, including a variety of link functions,
scale functions, nominal regressors, and various methods for setting thresholds. Effects are currently generated
only for the default values of the arguments scale
, nominal
, link
and threshold
,
which is equivalent to fitting an ordinal response mixed effects model with a logit link.
The effect methods can also be
used with objects created using clm
or clm2
fitting ordinal response models with the same links permitted by polr with no random effects, with
results similar to those from polr
in the MASS package.
Calling any of these functions from within a user-written function may result in errors due
to R's scoping rules. See the vignette embedding.pdf
for the car package
for a solution to this problem.
John Fox jfox@mcmaster.ca, Sanford Weisberg sandy@umn.edu and Jangman Hong.
Fox, J. (1987). Effect displays for generalized linear models. Sociological Methodology 17, 347–361.
Fox, J. (2003) Effect displays in R for generalised linear models. Journal of Statistical Software 8:15, 1–27, <http://www.jstatsoft.org/v08/i15/>.
Fox, J. and R. Andersen (2006). Effect displays for multinomial and proportional-odds logit models. Sociological Methodology 36, 225–255.
Fox, J. and J. Hong (2009). Effect displays in R for multinomial and proportional-odds logit models: Extensions to the effects package. Journal of Statistical Software 32:1, 1–24, <http://www.jstatsoft.org/v32/i01/>.
Fox, J. and S. Weisberg (forthcoming). Visualizing Fit and Lack of Fit in Complex Regression Models: Effect Plots with Partial Residuals. Journal of Statistical Software.
Hastie, T. J. (1992). Generalized additive models. In Chambers, J. M., and Hastie, T. J. (eds.) Statistical Models in S, Wadsworth.
Weisberg, S. (2014). Applied Linear Regression, 4th edition, Wiley, http://z.umn.edu/alr4ed.
LegacyArguments
. For information on printing, summarizing, and plotting effects:
print.eff
, summary.eff
, plot.eff
,
print.summary.eff
,
print.effpoly
, summary.effpoly
, plot.effpoly
,
print.efflist
, summary.efflist
,
plot.efflist
, xyplot
,
densityplot
.
mod.cowles <- glm(volunteer ~ sex + neuroticism*extraversion, data=Cowles, family=binomial) eff.cowles <- allEffects(mod.cowles, xlevels=list(extraversion=seq(0, 24, 6)), fixed.predictors=list(given.values=c(sexmale=0.5))) eff.cowles as.data.frame(eff.cowles[[2]]) # the following are equivalent: eff.ne <- effect("neuroticism*extraversion", mod.cowles) Eff.ne <- Effect(c("neuroticism", "extraversion"), mod.cowles) all.equal(eff.ne$fit, Eff.ne$fit) plot(eff.cowles, 'sex', axes=list(y=list(lab="Prob(Volunteer)"))) plot(eff.cowles, 'neuroticism:extraversion', axes=list(y=list(lab="Prob(Volunteer)", ticks=list(at=c(.1,.25,.5,.75,.9))))) plot(Effect(c("neuroticism", "extraversion"), mod.cowles, se=list(type="Scheffe"), xlevels=list(extraversion=seq(0, 24, 6)), fixed.predictors=list(given.values=c(sexmale=0.5))), axes=list(y=list(lab="Prob(Volunteer)", ticks=list(at=c(.1,.25,.5,.75,.9))))) plot(eff.cowles, 'neuroticism:extraversion', lines=list(multiline=TRUE), axes=list(y=list(lab="Prob(Volunteer)"))) plot(effect('sex:neuroticism:extraversion', mod.cowles, xlevels=list(extraversion=seq(0, 24, 6))), lines=list(multiline=TRUE)) # a nested model: mod <- lm(log(prestige) ~ income:type + education, data=Prestige) plot(Effect(c("income", "type"), mod, transformation=list(link=log, inverse=exp)), axes=list(y=list(lab="prestige"))) if (require(nnet)){ mod.beps <- multinom(vote ~ age + gender + economic.cond.national + economic.cond.household + Blair + Hague + Kennedy + Europe*political.knowledge, data=BEPS) plot(effect("Europe*political.knowledge", mod.beps, xlevels=list(political.knowledge=0:3))) plot(Effect(c("Europe", "political.knowledge"), mod.beps, xlevels=list(Europe=1:11, political.knowledge=0:3), fixed.predictors=list(given.values=c(gendermale=0.5))), lines=list(col=c("blue", "red", "orange")), axes=list(x=list(rug=FALSE), y=list(style="stacked"))) plot(effect("Europe*political.knowledge", mod.beps, # equivalent xlevels=list(Europe=1:11, political.knowledge=0:3), fixed.predictors=list(given.values=c(gendermale=0.5))), lines=list(col=c("blue", "red", "orange")), axes=list(x=list(rug=FALSE), y=list(style="stacked"))) } if (require(MASS)){ mod.wvs <- polr(poverty ~ gender + religion + degree + country*poly(age,3), data=WVS) plot(effect("country*poly(age, 3)", mod.wvs)) plot(Effect(c("country", "age"), mod.wvs), axes=list(y=list(style="stacked"))) plot(effect("country*poly(age, 3)", mod.wvs), axes=list(y=list(style="stacked"))) # equivalent plot(effect("country*poly(age, 3)", latent=TRUE, mod.wvs)) plot(effect("country*poly(age, 3)", latent=TRUE, mod.wvs, se=list(type="scheffe"))) # Scheffe-type confidence envelopes } mod.pres <- lm(prestige ~ log(income, 10) + poly(education, 3) + poly(women, 2), data=Prestige) eff.pres <- allEffects(mod.pres, xlevels=50) plot(eff.pres) plot(eff.pres[1], axes=list(x=list(income=list( transform=list(trans=log10, inverse=function(x) 10^x), ticks=list(at=c(1000, 2000, 5000, 10000, 20000)) )))) # linear model with log-response and log-predictor # to illustrate transforming axes and setting tick labels mod.pres1 <- lm(log(prestige) ~ log(income) + poly(education, 3) + poly(women, 2), data=Prestige) # effect of the log-predictor eff.log <- Effect("income", mod.pres1) # effect of the log-predictor transformed to the arithmetic scale eff.trans <- Effect("income", mod.pres1, transformation=list(link=log, inverse=exp)) #variations: # y-axis: scale is log, tick labels are log # x-axis: scale is arithmetic, tick labels are arithmetic plot(eff.log) # y-axis: scale is log, tick labels are log # x-axis: scale is log, tick labels are arithmetic plot(eff.log, axes=list(x=list(income=list( transform=list(trans=log, inverse=exp), ticks=list(at=c(5000, 10000, 20000)), lab="income, log-scale")))) # y-axis: scale is log, tick labels are arithmetic # x-axis: scale is arithmetic, tick labels are arithmetic plot(eff.trans, axes=list(y=list(lab="prestige"))) # y-axis: scale is arithmetic, tick labels are arithmetic # x-axis: scale is arithmetic, tick labels are arithmetic plot(eff.trans, axes=list(y=list(type="response", lab="prestige"))) # y-axis: scale is log, tick labels are arithmetic # x-axis: scale is log, tick labels are arithmetic plot(eff.trans, axes=list( x=list(income=list( transform=list(trans=log, inverse=exp), ticks=list(at=c(1000, 2000, 5000, 10000, 20000)), lab="income, log-scale")), y=list(lab="prestige, log-scale")), main="Both response and X in log-scale") # y-axis: scale is arithmetic, tick labels are arithmetic # x-axis: scale is log, tick labels are arithmetic plot(eff.trans, axes=list( x=list( income=list(transform=list(trans=log, inverse=exp), ticks=list(at=c(1000, 2000, 5000, 10000, 20000)), lab="income, log-scale")), y=list(type="response", lab="prestige"))) if (require(nlme)){ # for gls() mod.hart <- gls(fconvict ~ mconvict + tfr + partic + degrees, data=Hartnagel, correlation=corARMA(p=2, q=0), method="ML") plot(allEffects(mod.hart)) detach(package:nlme) } if (require(lme4)){ data(cake, package="lme4") fm1 <- lmer(angle ~ recipe * temperature + (1|recipe:replicate), cake, REML = FALSE) plot(Effect(c("recipe", "temperature"), fm1)) plot(effect("recipe:temperature", fm1), axes=list(grid=TRUE)) # equivalent (plus grid) if (any(grepl("pbkrtest", search()))) detach(package:pbkrtest) detach(package:lme4) } if (require(nlme) && length(find.package("lme4", quiet=TRUE)) > 0){ data(cake, package="lme4") cake$rep <- with(cake, paste( as.character(recipe), as.character(replicate), sep="")) fm2 <- lme(angle ~ recipe * temperature, data=cake, random = ~ 1 | rep, method="ML") plot(Effect(c("recipe", "temperature"), fm2)) plot(effect("recipe:temperature", fm2), axes=list(grid=TRUE)) # equivalent (plus grid) } detach(package:nlme) if (require(poLCA)){ data(election) f2a <- cbind(MORALG,CARESG,KNOWG,LEADG,DISHONG,INTELG, MORALB,CARESB,KNOWB,LEADB,DISHONB,INTELB)~PARTY*AGE nes2a <- poLCA(f2a,election,nclass=3,nrep=5) plot(Effect(c("PARTY", "AGE"), nes2a), axes=list(y=list(style="stacked"))) } # mlm example if (require(heplots)) { data(NLSY, package="heplots") mod <- lm(cbind(read,math) ~ income+educ, data=NLSY) eff.inc <- Effect("income", mod) plot(eff.inc) eff.edu <- Effect("educ", mod) plot(eff.edu, axes=list(x=list(rug=FALSE), grid=TRUE)) plot(Effect("educ", mod, response="read")) detach(package:heplots) } # svyglm() example (adapting an example from the survey package) if (require(survey)){ data(api) dstrat<-svydesign(id=~1, strata=~stype, weights=~pw, data=apistrat, fpc=~fpc) mod <- svyglm(sch.wide ~ ell + meals + mobility, design=dstrat, family=quasibinomial()) plot(allEffects(mod), axes=list(y=list(lim=log(c(0.4, 0.99)/c(0.6, 0.01)), ticks=list(at=c(0.4, 0.75, 0.9, 0.95, 0.99))))) } # component + residual plot examples Prestige$type <- factor(Prestige$type, levels=c("bc", "wc", "prof")) mod.prestige.1 <- lm(prestige ~ income + education, data=Prestige) plot(allEffects(mod.prestige.1, residuals=TRUE)) # standard C+R plots plot(allEffects(mod.prestige.1, residuals=TRUE, se=list(type="scheffe"))) # with Scheffe-type confidence bands mod.prestige.2 <- lm(prestige ~ type*(income + education), data=Prestige) plot(allEffects(mod.prestige.2, residuals=TRUE)) mod.prestige.3 <- lm(prestige ~ type + income*education, data=Prestige) plot(Effect(c("income", "education"), mod.prestige.3, residuals=TRUE), partial.residuals=list(span=1)) # artificial data set.seed(12345) x1 <- runif(500, -75, 100) x2 <- runif(500, -75, 100) y <- 10 + 5*x1 + 5*x2 + x1^2 + x2^2 + x1*x2 + rnorm(500, 0, 1e3) Data <- data.frame(y, x1, x2) mod.1 <- lm(y ~ poly(x1, x2, degree=2, raw=TRUE), data=Data) # raw=TRUE necessary for safe prediction mod.2 <- lm(y ~ x1*x2, data=Data) mod.3 <- lm(y ~ x1 + x2, data=Data) plot(Effect(c("x1", "x2"), mod.1, residuals=TRUE)) # correct model plot(Effect(c("x1", "x2"), mod.2, residuals=TRUE)) # wrong model plot(Effect(c("x1", "x2"), mod.3, residuals=TRUE)) # wrong model