effect: Functions For Constructing Effect Plots

Description

effect constructs 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. The function can also be used with mixed-effects models fit with lmer from the lme4 package, or fit with lme from the nlme package. In mixed effects models the analysis is for the fixed effects only, not for random effects. Effect also constructs an an "eff" object, but rather than focusing on a term in the model, it focuses on a subset of the predictors. 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). 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, ...)

## S3 method for class 'lm':
effect(term, mod, xlevels=list(), default.levels=10, given.values, 
	se=TRUE, confidence.level=.95, 
    transformation=list(link=family(mod)$linkfun, inverse=family(mod)$linkinv), 
    typical=mean, offset=mean, ...)
    
## S3 method for class 'gls':
effect(term, mod, xlevels=list(), default.levels=10, given.values,
	se=TRUE, confidence.level=.95, transformation=NULL, typical=mean, ...)
    
## S3 method for class 'multinom':
effect(term, mod, confidence.level=.95, xlevels=list(), default.levels=10, 
	given.values, se=TRUE, typical=mean, ...)
		
## S3 method for class 'polr':
effect(term, mod, confidence.level=.95, xlevels=list(), default.levels=10, 
	given.values, se=TRUE, typical=mean, latent=FALSE, ...)
	
## S3 method for class 'mer':
effect(term, mod, ...)

## S3 method for class 'lme':
effect(term, mod, ...)

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

## S3 method for class 'lm':
Effect(focal.predictors, mod, xlevels=list(), default.levels=10, given.values, 
    se=TRUE, confidence.level=0.95, 
    transformation=list(link=family(mod)$linkfun, inverse=family(mod)$linkinv), 
    typical=mean, offset=mean, ...)
        
## S3 method for class 'gls':
Effect(focal.predictors, mod, xlevels = list(), default.levels = 10, 
    given.values, 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=10, 
    given.values, se=TRUE, typical=mean, ...)
                            
## S3 method for class 'polr':
Effect(focal.predictors, mod, 
    confidence.level=.95, xlevels=list(), default.levels=10, 
    given.values, se=TRUE, typical=mean, latent=FALSE, ...)
    
## S3 method for class 'mer':
Effect(focal.predictors, mod, ...)

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

## S3 method for class 'mer':
allEffects(mod, ...)

## S3 method for class 'lme':
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 f

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

xlevels

an optional list of values at which to set covariates, with components of the form covariate.name = vector.of.values.

default.levels

number of values for covariates that are not specified explicitly via xlevels; covariate values set by default are evenly spaced between the minimum and maximum values in the data.

given.values

a numeric vector of named elements, setting particular columns of the model matrix to specific values for terms not appearing in an effect; if specified, this argument takes precedence over the application of the function given in the

if TRUE, the default, calculate standard errors and confidence limits for the effects. For mer 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.

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

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

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.

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.

object

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

...

arguments to be passed down.

Value

For lm, glm, mer and lme, effect and Effect return an "eff" object, and for multinom and polr, an "effpoly" object, with the following components:
termthe term to which the effect pertains.
formulathe complete model formula.
responsea character string giving the name of the response variable.
y.levels(for "effpoly" objects) levels of the polytomous response variable.
variablesa 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).
xa 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.matrixthe model matrix from which the effect was calculated.
dataa data frame with the data on which the fitted model was based.
discrepancythe 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.
offsetvalue 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.levelfor 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.
effectList returns a list of "eff" or "effpoly" objects corresponding to the high-order terms of the model.

Warnings and Limitations

The effect function handles factors and covariates differently, and becomes 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 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 do not handle logical predictors; convert logical predictors to factors. 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, the 'plot' command with style="stacked" apparently does not work because of a bug in the barchart function in lattice. However, the stype="lines", the default, 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.

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.

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.

Examples

Run this code

mod.cowles <- glm(volunteer ~ sex + neuroticism*extraversion, 
    data=Cowles, family=binomial)
eff.cowles <- allEffects(mod.cowles, xlevels=list(neuroticism=0:24, 
    extraversion=seq(0, 24, 6)), given.values=c(sexmale=0.5))
eff.cowles

	# 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(neuroticism=0:24, extraversion=seq(0, 24, 6))), multiline=TRUE)
    
as.data.frame(eff.cowles[[2]])
    
# 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
    
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(Europe=1:11, political.knowledge=0:3)))  
plot(effect("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(c("Europe","political.knowledge"), mod.beps,  # equivalent
    	xlevels=list(Europe=1:11, political.knowledge=0:3),
		given.values=c(gendermale=0.5)),
	style="stacked", colors=c("blue", "red", "orange"), rug=FALSE)
    
mod.wvs <- polr(poverty ~ gender + religion + degree + country*poly(age,3),
	data=WVS)
plot(effect("country*poly(age, 3)", mod.wvs))
plot(effect("country*poly(age, 3)", mod.wvs), style="stacked")
plot(Effect(c("country", "age"), 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, default.levels=50)
plot(eff.pres)
plot(eff.pres[1], 
    transform.x=list(income=list(trans=log10, inverse=function(x) 10^x)),
    ticks.x=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(c("income"), mod.pres1, default.levels=50)
# effect of the log-predictor transformed to the arithmetic scale
eff.trans <- Effect(c("income"), mod.pres1, default.levels=50,
       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(at=c(1000, 2000, 5000, 10000, 20000)),
       xlab="income, log-scale")

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

# y-axis:  scale is arithmetic, tick labels are airthmetic
# x-axis:  scale is arithmetic, tick labels are arithmetic
plot(eff.trans, rescale.axis=FALSE, 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(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(at=c(1000, 2000, 5000, 10000, 20000)),
       rescale.axis=FALSE,
       xlab="income, log-scale", ylab="prestige")

library(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))

	library(lme4)
	data(cake, package="lme4")
	fm1 <- lmer(angle ~ recipe * temperature + (1|recipe:replicate), cake,
                   REML = FALSE)
	plot(effect("recipe:temperature", fm1), grid=TRUE)
    plot(Effect(c("recipe", "temperature"), fm1)) # equivalent

	detach(package:lme4) # if previously attached
	library(nlme)
	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("recipe:temperature", fm2), grid=TRUE)
    plot(Effect(c("recipe", "temperature"), fm2)) # equivalent

Run the code above in your browser using DataLab