effect {effects}R Documentation

Functions For Constructing Effect Plots

Description

effect and Effect construct an "eff" object for a term (usually a high-order term) in a linear (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 (mlm), 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 requires specifying a the components of the term, for example c("a", "b") rather than the term itself. Effect is consequently more flexible and robust than effect, and will work with some models for which effect fails, such as models with nested terms (see the examples). The effect function works by constructing a call to Effect.

The function can also be used with some mixed-effects models, including linear and generalized linear model fits from lmer and glmer from the lme4 package and lme from the nlme package; ordinal logistic mixed effects fit with clmm2 from the ordinal package, and polytomous latent-class models fit by the poLCA function in the poLCA package. The analysis in mixed effects is for the fixed effects only, not for random effects.

There is a default method for Effect that should work with any model object that has a linear predictor and responds to the coef, model.frame, formula, and vcov functions.

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").

Usage

effect(term, mod, vcov.=vcov, ...)

Effect(focal.predictors, mod, ...)

## S3 method for class 'lm'
Effect(focal.predictors, mod, xlevels = list(), default.levels = NULL, given.values,
    vcov.=vcov, se = TRUE, confidence.level = 0.95, 
    transformation = list(link = family(mod)$linkfun, inverse = family(mod)$linkinv), 
    typical = mean, offset = mean, 
    partial.residuals=FALSE, quantiles=seq(0.2, 0.8, by=0.2),
    x.var=NULL,  ...)
        
## S3 method for class 'gls'
Effect(focal.predictors, mod, xlevels = list(), default.levels=NULL, 
    given.values, vcov.=vcov, se = TRUE, confidence.level = 0.95, transformation = NULL, 
    typical = mean, ...)
    
## S3 method for class 'multinom'
Effect(focal.predictors, mod, 
    confidence.level=.95, xlevels=list(), default.levels=NULL, 
    given.values, vcov.=vcov, se=TRUE, typical=mean, ...)
                            
## S3 method for class 'polr'
Effect(focal.predictors, mod, 
    confidence.level=.95, xlevels=list(), default.levels=NULL, 
    given.values, vcov.=vcov, se=TRUE, typical=mean, latent=FALSE, ...)
    
## S3 method for class 'clm2'
Effect(focal.predictors, mod, ...)
    
## S3 method for class 'mer'
Effect(focal.predictors, mod, 
    KR=FALSE, ...)

## S3 method for class 'merMod'
Effect(focal.predictors, mod, 
    KR=FALSE, ...)

## S3 method for class 'lme'
Effect(focal.predictors, mod, ...)

## S3 method for class 'poLCA'
Effect(focal.predictors, mod, ...)

## S3 method for class 'mlm'
Effect(focal.predictors, mod, response, ...)

## Default S3 method:
Effect(focal.predictors, mod, xlevels = list(), 
    default.levels = NULL, given.values,
    vcov. = vcov, se = TRUE, confidence.level = 0.95, 
    transformation = list(link = I, inverse = I), 
    typical = mean, offset = mean, ...)
    
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, ...)

Arguments

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 (:) or asterisks (*) may be used for interactions.

focal.predictors

a character vector of one or more predictors in the model.

mod

an object of class "lm", "gls", "glm", "multinom", "polr", "mer" (or "merMod"), "lme" or "poLCA".

xlevels

this argument is used to set the number of levels for any focal predictor that is not a factor. If xlevels=NULL, the default, then the number and values of levels for any numeric predictor is determined by grid.pretty. If xlevels=n is an integer, then each numeric predictor is represented by n equally spaced levels. More generally, xlevels can be a named list of values at which to set each numeric predictor. For example, xlevels=list(x1=c(2, 4, 7), x2=5) would use the values 2, 4 and 7 for the levels of x1, 5 equally spaced levels for the levels of x2, and use the default for any other numeric predictors. If partial residuals are computed, then the focal predictor that is to appear on the horizontal axis of an effect plot is evaluated at 100 equally spaced values along its full range, and, by default, other numeric predictors are evaluated at the quantiles specified in the quantiles argument, unless their values are given explicitly in xlevels.

default.levels

ignored, but included for compatibility with pre-July 2013 versions of this package. Use xlevels instead.

given.values

a numeric vector of named elements, setting particular columns of the model matrix to specific values for predictors that are not focal predictors; if specified, this argument takes precedence over the application of the function given in the typical argument (below). Care must be taken in specifying these values — e.g., for a factor, the values of all contrasts should be given and these should be consistent with each other.

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, vcov, which will result in the function call vcov(mod) to get the variance covariance matrix. You can use the name of any function with mod as the value of its first argument that returns an estimated sample covariance matrix, such as the hccm function in the car package that returns a heteroscedasticity corrected estimate with linear models.

se

if TRUE, the default, calculate standard errors and confidence limits for the effects. For mer, merMod, and lme objects, the normal distribution is used to get confidence limits.

confidence.level

level at which to compute confidence limits based on the standard-normal distribution; the default is 0.95.

KR

if TRUE and the pbkrtest package is installed, use the Kenward-Roger coefficient covariance matrix to compute effect standard errors for linear mixed models fit with lmer in the lme4 package. The default is FALSE because the computation can be very slow.

transformation

a two-element list with elements link and inverse. For a generalized linear model, these are by default the link function and inverse-link (mean) function. For a linear model, these default to NULL. If NULL, the identify function, I, is used; this effect can also be achieved by setting the argument to NULL. The inverse-link may be used to transform effects when they are printed or plotted; the link may be used in positioning axis labels (see below). If the link is not given, an attempt will be made to approximate it from the inverse-link.

typical

a function to be applied to the columns of the model matrix over which the effect is "averaged"; the default is mean.

offset

a function to be applied to the offset values (if there is an offset) in a linear or generalized linear model, or a mixed-effects model fit by lmer or glmer; or a numeric value, to which the offset will be set. The default is the mean function, and thus the offset will be set to its mean. Note: Only offsets defined by the offset argument to lm, glm, lmer, or glmer will be handled correctly; use of the offset function in the model formula is not supported.

partial.residuals

if TRUE, partial residuals will be computed for an effect in a linear or generalized linear model; if FALSE (the default), partial residuals are suppressed.

quantiles

quantiles at which to evaluate numeric focal predictors not on the horizontal axis, used only when partial residuals are displayed; superceded if the xlevels argument gives specific values for a predictor.

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 NULL, in which case the first numeric predictor in the effect will be used if partial residuals are to be computed. This argument is intended to be used when partial.residuals is TRUE; otherwise, the variable on the horizontal axis can be chosen when the effect object is plotted: see plot.eff.

latent

if TRUE, effects in a proportional-odds logit model are computed on the scale of the latent response; if FALSE (the default) effects are computed as individual-level probabilities and logits.

x

an object of class "eff", "effpoly", or "efflatent".

transform

a transformation to be applied to the effects and confidence limits, by default taken from the inverse link function saved in the "eff" object.

row.names, optional

not used.

response

for an mlm, a vector containing the name(s) or indices of one or more response variable(s). The default is to use all responses in the model.

object

an object of class "eff" for which the covariance matrix of the effects is desired.

...

arguments to be passed down.

Details

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.

Value

For lm, glm, mer and lme, effect and Effect return an "eff" object, and for multinom, polr and clmm2, an "effpoly" object, with the components listed below. For mlm with one response specified, an "eff" object, otherwise a class "efflist" object, 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 "effpoly" objects) levels of the polytomous response variable.

variables

a list with information about each predictor, including its name, whether it is a factor, and its levels or values.

fit

(for "eff" objects) a one-column matrix of fitted values, representing the effect on the scale of the linear predictor; this is a ravelled table, representing all combinations of predictor values.

prob

(for "effpoly" objects) a matrix giving fitted probabilities for the effect for the various levels of the the response (columns) and combinations of the focal predictors (rows).

logit

(for "effpoly" objects) a matrix giving fitted logits for the effect for the various levels of the the response (columns) and combinations of the focal predictors (rows).

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 gls models, this is now necessarily 0.

offset

value to which the offset is fixed; 0 if there is no offset.

model

(for "effpoly" objects) "multinom" or "polr", as appropriate.

vcov

(for "eff" objects) a covariance matrix for the effect, on the scale of the linear predictor.

se

(for "eff" objects) a vector of standard errors for the effect, on the scale of the linear predictor.

se.prob, se.logit

(for "effpoly" objects) matrices of standard errors for the effect, on the probability and logit scales.

lower, upper

(for "eff" objects) one-column matrices of confidence limits, on the scale of the linear predictor.

lower.prob, upper.prob, lower.logit, upper.logit

(for "effpoly" objects) matrices of confidence limits for the fitted logits and probabilities; the latter are computed by transforming the former.

confidence.level

for the confidence limits.

transformation

(for "eff" objects) a two-element list, with element link giving the link function, and element inverse giving the inverse-link (mean) function.

fitted.rounded

partial fitted values at the observed values of the predictor x.var and the values of the other predictors that appear in the effect display; predictors not in the effect are held constant to typical values. This and the following two elements are NULL if partial residuals aren't computed and pertain only to linear or generalized linear models.

fitted

partial fitted values for the observed values of all predictors that appear in the effect display; predictors not in the effect are held constant to typical values.

partial.residuals.raw

partial residuals for the effect computed at the actual values of all focal predictors.

partial.residuals.adjusted

partial residuals for the effect computed at the panel-rounded values of focal predictors, except for the predictor corresponding to xvar, which will appear on the horizontal axis of the plotted effect.

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 NULL if partial residuals aren't computed.

family

for a "glm" model, the name of the distributional family of the model; for an "lm" model, this is "gaussian"; otherwise NULL. The family controls how partial residuals are smoothed in plots.

effectList returns 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) to fit 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.

Warnings and Limitations

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. Similarly variables of class date or "times", which are usually differences between "dates" variables, should be converted to numeric variables outside 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"). In addition, factors cannont be declared on the fly (e.g., using y ~ a + factor(b).

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.

With multinom models with two or more crossed factors with an empty cell, the 'plot' command with style="stacked" apparently does not work because of a bug in the barchart function in lattice. However, the default style="lines" does 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 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 logistic link. The effect methods can also be used with objects created using clm2 fitting ordinal response logistic models 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.

Author(s)

John Fox jfox@mcmaster.ca, Sanford Weisberg sandy@umn.edu and Jangman Hong.

References

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/>.

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.

See Also

print.eff, summary.eff, plot.eff, print.summary.eff, print.effpoly, summary.effpoly, plot.effpoly, print.efflist, summary.efflist, plot.efflist, xyplot, densityplot

Examples

# Note: Some of these examples are marked as "don't test"
#       to reduce the execution times of the examples
#       when the package is checked.

mod.cowles <- glm(volunteer ~ sex + neuroticism*extraversion,
    data=Cowles, family=binomial)
eff.cowles <- allEffects(mod.cowles, xlevels=list(extraversion=seq(0, 24, 6)),
    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', ylab="Prob(Volunteer)")

plot(eff.cowles, 'neuroticism:extraversion', ylab="Prob(Volunteer)",
    ticks=list(at=c(.1,.25,.5,.75,.9)))

plot(eff.cowles, 'neuroticism:extraversion', multiline=TRUE,
    ylab="Prob(Volunteer)")

plot(effect('sex:neuroticism:extraversion', mod.cowles,
    xlevels=list(extraversion=seq(0, 24, 6))), multiline=TRUE)


# a nested model:

mod <- lm(log(prestige) ~ income:type + education, data=Prestige)

# does not work: effect("income:type", mod, transformation=list(link=log, inverse=exp))

plot(Effect(c("income", "type"), mod, transformation=list(link=log, inverse=exp),
	ylab="prestige")) # works


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),
		given.values=c(gendermale=0.5)),
	style="stacked", colors=c("blue", "red", "orange"), rug=FALSE)
    
 
plot(effect("Europe*political.knowledge", mod.beps, # equivalent
    	xlevels=list(political.knowledge=0:3),
		given.values=c(gendermale=0.5)),
	style="stacked", colors=c("blue", "red", "orange"), rug=FALSE)

}

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), style="stacked") 


plot(effect("country*poly(age, 3)", mod.wvs), style="stacked") # equivalent

plot(effect("country*poly(age, 3)", latent=TRUE, mod.wvs))

}


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],
    transform.x=list(income=list(trans=log10, inverse=function(x) 10^x)),
    ticks.x=list(income=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, transform.x=list(income=c(trans=log, inverse=exp)),
     ticks.x=list(income=list(at=c(1000, 2000, 5000, 10000, 20000))),
     xlab="income, log-scale")

# y-axis:  scale is log, tick labels are arithmetic
# x-axis:  scale is arithmetic, tick labels are arithmetic
plot(eff.trans, ylab="prestige")

# y-axis:  scale is arithmetic, tick labels are arithmetic
# x-axis:  scale is arithmetic, tick labels are arithmetic
plot(eff.trans, type="response", ylab="prestige")

# y-axis:  scale is log, tick labels are arithmetic
# x-axis:  scale is log, tick labels are arithmetic
plot(eff.trans, transform.x=list(income=c(trans=log, inverse=exp)),
     ticks.x=list(income=list(at=c(1000, 2000, 5000, 10000, 20000))),
     xlab="income, log-scale", ylab="prestige, log-scale",
     main="Both effect and X in log-scale")

# y-axis:  scale is arithmetic, tick labels are airthmetic
# x-axis:  scale is log, tick labels are arithmetic
plot(eff.trans, transform.x=list(income=c(trans=log, inverse=exp)),
     ticks.x=list(income=list(at=c(1000, 2000, 5000, 10000, 20000))),
     type="link",
     xlab="income, log-scale", ylab="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), grid=TRUE) # equivalent

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), grid=TRUE)  # equivalent
}
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), 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, rug=FALSE, grid=TRUE)
    
	plot(Effect("educ", mod, response="read"))
    
  detach(package:heplots)
}

# 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, partial.residuals=TRUE)) # standard C+R plots

mod.prestige.2 <- lm(prestige ~ type*(income + education), data=Prestige)
plot(allEffects(mod.prestige.2, partial.residuals=TRUE))

mod.prestige.3 <- lm(prestige ~ type + income*education, data=Prestige)
plot(Effect(c("income", "education"), mod.prestige.3, partial.residuals=TRUE),
  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)

.save.strip <- setStrip() # change color of lattice strips
plot(Effect(c("x1", "x2"), mod.1, partial.residuals=TRUE)) # correct model
plot(Effect(c("x1", "x2"), mod.2, partial.residuals=TRUE)) # wrong model
plot(Effect(c("x1", "x2"), mod.3, partial.residuals=TRUE)) # wrong model
restoreStrip(.save.strip)
remove(.save.strip)
[Package effects version 3.0-4 Index]