Eﬀect Displays in R for Generalised Linear Models

This paper describes the implementation in R of a method for tabular or graphical display of terms in a complex generalised linear model. By complex, I mean a model that contains terms related by marginality or hierarchy, such as polynomial terms, or main eﬀects and interactions. I call these tables or graphs eﬀect displays . Eﬀect displays are constructed by identifying high-order terms in a generalised linear model. Fitted values under the model are computed for each such term. The lower-order ‘relatives’ of a high-order term (e


Background and Motivation
This paper describes the implementation in R of a method for tabular or graphical display of terms in a complex generalised linear model.By complex, I mean a model that contains terms related by marginality or hierarchy (in the sense of Nelder, 1977), such as polynomial terms, or main effects and interactions.I call these tables or graphs effect displays.
I assume that readers are familiar with R, which is a free, open-source implementation of the S statistical programming language and computing environment.Extensive information on R is available at <http://www.r-project.org/>.The software described in the paper is distributed as the effects package for R, available on the Comprehensive R Archive Network (CRAN) at <http://cran.r-project.org/>.
The general approach implemented in the effects package is a modest extension of methods introduced by Fox (1987).Similar, if less general, ideas have a long history, dating at least to Fisher's (1936) "adjusted means" in the analysis of covariance.Goodnight and Harvey's (1978) "least-squares means" in analysis of variance and covariance, and Searle, Speed, and Milliken's (1980) "estimated population marginal means" are other examples, though all of these approaches are restricted to interactions among factors (i.e., categorical predictors) in a linear model.Recently King, Tomz, and Wittenberg (2000) and Tomz, Wittenberg, and King (2003) have presented similar ideas, but their approach, based on Monte-Carlo simulation of a model, is much more complicated than the one discussed here.
To motivate effect displays, consider the following data on police treatment of individuals arrested in Toronto for simple possession of small quantities of marijuana.(The data discussed here are part of a larger data set featured in a series of articles in the Toronto Star newspaper.)Under these circumstances police have the option of releasing an arrestee with a summons to appear in court -similar to a traffic ticket; alternatively, the individual may be brought to the police station for questioning and possible indictment.The principal question of interest is whether and how the probability of release is influenced by the subject's race, age, and other characteristics.The data are in the data frame Arrests in the effects package: > library(effects) > data(Arrests) > dim(Arrests) [1]  Most of the variables in the data set are self explanatory, with the following exceptions: • colour is either Black or White.The original data included the additional categories "Brown" and "Other," but their meaning is ambiguous and their use relatively infrequent.
• The observations span the years 1997 through (part of) 2002.A few arrests in 1996 were eliminated.
In the analysis reported below, year is treated as a factor (i.e., as a categorical predictor).
• When suspects are stopped by the police, their names are checked in six data bases -of previous arrests, previous convictions, parole status, and so on.The variable checks records the number of data bases on which an individual's name appeared.
Preliminary analysis of the data suggested a logit model including interactions between colour and year and between colour and age, and main effects of employed, citizen, and checks.The effects of age and checks appear to be reasonably linear on the logit scale and are modelled as such: > Arrests$year <-as.factor(Arrests$year)> arrests.It is relatively difficult (but not impossible) to discern from inspection of the coefficients how colour, year, and age combine to influence the probability of release.(In contrast, the effects of employed, citizen, and checks are quite clear from the coefficients, although some mental arithmetic is required unless one is comfortable with interpreting effects on the logit scale.)The model employs the default 0/1 dummy coding for factors that R terms "treatment" contrasts (but the coding employed is immaterial to the effect displays described below).

Effect Displays
A general principle of interpretation for statistical models containing terms that are marginal to others is that high-order terms should be combined with their lower-order relatives -for example, an interaction between two factors should be combined with the main effects marginal to the interaction.In conformity with this principle, Fox (1987) suggests identifying the high-order terms in a generalised linear model.Fitted values under the model are computed for each such term.The lower-order 'relatives' of a high-order term (e.g., main effects marginal to an interaction) are absorbed into the term, allowing the predictors appearing in the high-order term to range over their values.The values of other predictors are fixed at typical values: for example, a covariate could be fixed at its mean or median, a factor at its proportional distribution in the data, or to equal proportions in its several levels.This procedure is illustrated below.
Consider the generalised linear model with linear predictor η = Xβ and link function g(µ) = η, where µ is the expectation of the response vector y.Suppose that we have an estimate b of β, along with the estimated covariance matrix V (b) of b.Let the rows of X * include all combinations of values of predictors appearing in a high-order term, along with typical values of the remaining predictors.The structure of X * with respect, for example, to interactions, is the same as that of the model matrix X.Then the fitted values η * = X * b represent the effect in question, and a table or graph of these values -or, alternatively, of the fitted values transformed to the scale of the response, g −1 ( η * ) -is an effect display.The standard errors of η * , available as the square-root diagonal entries of X * V (b)X * , may be used to compute pointwise confidence intervals for the effects, the end-points of which may then also be transformed to the scale of the response.
For example, colour:age is a high-order term of the logit model fit to the Arrests data in the preceding section.The lower-order relatives of this term are the constant, and the terms for colour and age.The arrestees in the data set range in age from 12 to 66, with the .025and .975quantiles at 15 and 45, respectively.As mentioned, there are two levels of colour, so if we set age to all integer values between 15 and 45, there are 62 combinations of values of these two predictors.Other predictors in the model (employed, citizen, checks, and year) will be set to typical values, as previously explained.
In particular: • Column 1 of X * represents the constant.
• Column 2 reflects the 79 percent of arrestees who were at level Yes of employed, and hence had values of 1 on the treatment-coded contrast for this factor; 0.79 is therefore also the mean of the contrast.Note that this column, along with other constant columns in X * , is in effect absorbed in the constant term, and therefore influences only the average level of the computed effects (on the scale of the linear predictor).
• Column 3 reflects the 85 percent of arrestees who were in level Yes of citizen.The vertical axis is labelled on the probability scale, and a 95-percent pointwise confidence interval is drawn around the estimated effect.
• Column 5 repeats the two values 0 and 1 for the contrast for colour (to be taken in combination with the values of age in column 11).
• Columns 6 through 10 represent the contrasts for year, and contain the proportions of arrestees in years 1998 through 2002; this reflects the use of the first level of year, 1997, as the baseline level.
• Column 11 contains the twice-repeated integer values of age, from 15 through 65.
• Columns 12 through 16 are for the interaction of colour with year (which is absorbed in the colour term).
• Column 17 is for the colour by age interaction.
Although it is not difficult to construct X * directly in this manner, it is tedious to do so.A graphical effect display for the colour by age interaction, computed by the software described in Section 3, appears in Figure 1: Apparently age has quite a different relationship to the probability of release for blacks and whites: older blacks are more likely to be released than younger blacks, while older whites are less likely to be released than younger whites; the relationship between age and the probability of release is also steeper for blacks than for whites.
Notice from the unequal spacing of the tick marks on the vertical axis of Figure 1 that although the axis is labelled on the scale of the response (i.e., the probability scale), the effects are plotted on the scale of the linear predictor (the logit scale); consequently, the lines plotted on the display are straight.I return to this point below.

Effect-Display Software
The software in the effects package consists of several related functions, described in more specific detail in the help pages for the package given as an appendix to this paper: • The function effect returns an object of class effect, containing information for constructing an effect display.The essential input to effect includes a linear (lm) or generalised-linear (glm) model object, and a term for which the effect is to be computed.This term will usually be a high-order term of the model, but it is possible to compute effects for lower-order terms (averaging over their higher-order relatives), and for terms that are not in the model but that have relatives that are in the model.For example, in the model fit to the Arrests data, the two way interactions colour:year and colour:age appear in the model, but not the three-way interaction colour:year:age (nor the two-way interaction year:age).Computing the colour:year:age effect combines the colour:year and colour:age interactions.In my experience, the result of combining lower-order interactions for overlapping predictors (here colour in colour:year and colour:age) is not always obvious, and such a display can therefore prove informative.For example, in Figure 2 the lines for each of blacks and whites are parallel across panels, reflecting the absence of a year:age interaction, but the relative heights of the lines vary, reflecting the colour:year interaction: In the earlier years, the line for blacks is mostly below that for whites, suggesting that (holding the other factors constant) at most ages whites were more likely than blacks to be released (particularly when we take into account the fact that most arrestees were young); in the last two years, however, the line for blacks is mostly above that for whites.
• Figure 2 is computed by the following command (note the warning message): > plot(effect("colour:year:age", arrests.mod,xlevels=list(age=15:45)), + multiline=TRUE, ylab="Probability(released)", rug=FALSE) Warning message: colour:year:age does not appear in the model in: effect("colour:year:age", arrests.mod,xlevels = list(age = 15:45)) • In this command, I have supplied an optional argument to the effect function: xlevels= list(age=15:45)) specifies values for age, which otherwise would default to 10 equally spaced values across the full range of the predictor.As mentioned, the ages 15 through 45 encompass 95 percent of the observations, excluding therefore a small number of very young and very old individuals.
• There are also optional arguments supplied to plot: -multiline=TRUE specifies that separate graphs should not be drawn for each level of colour.By default, confidence envelopes are suppressed in a multi-line effect display.
-ylab="Probability(released)" gives a non-default label for the vertical axis (the default would simply be the name of the response, released).
-rug=FALSE suppresses a "rug plot" (one-dimensional scatterplot) for age, which would otherwise appear on the horizontal axis of the graph.In some instances -but not here, where the data set is quite large -displaying the distribution of the variable on the horizontal axis in this manner can give the viewer a rough sense of where the data are located.
• The function all.effectstakes a linear or generalised linear model object as its required argument, finds all high-order terms in the model, and returns a list of effects corresponding to these terms.
• There are print, summary, and plot methods both for effect objects and for effect.listobjects (returned by effect and all.effects,respectively).Effect plots are created using Trellis graphics (via the lattice package in R).The plot method for effect.listobjects presents a text menu from which the user can select effects to graph.For example, to produce Figure 1 The plot method for effect objects allows one to graph effects on the response scale (rather than the scale of the linear predictor) or to transform the vertical axis of the plot in an arbitrary manner.Plotting effects on the response scale makes them more complex, however: Except in the case of the identity link, effects that are linear on the scale of the linear predictor are not linear on the scale of the response.Consequently, for example, profiles of effects are not parallel on the response scale even in the absence of interactions.As well, plotting on the scale of the linear predictor makes the choice of "typical" values for excluded effects less critical, in that these values only affect the labeling of the vertical axis (i.e., the level of the effects) and not the configuration of the display.Applied to the preceding example, plotting on the response (probability) scale produces Figure 3

Some Details
Many aspects of effect displays can be controlled through optional arguments to the effect function and to the plot, print, and summary methods for effect objects.Defaults are set in a manner intended to produce pleasing displays.For example, the predictor with the largest number of values (or levels) is plotted by default on the horizontal axis; similarly, in a multi-line display, the predictor with the smallest number of values is used to define the lines, different colours and line types are selected to represent the values of the predictor, and a legend is drawn at the top of the plot (as in Figures 2 and 3).All this can be modified via optional arguments.Some thought should be given to the scale on which effects are displayed.The default behaviour is to compute and graph effects on the scale of the linear predictor but to label the vertical axis on the scale of the response.Arbitrary transformations are supported, however.For example, if a linear model is fit to the log of the response (e.g., log dollar income), then the exponential function could be used to transform ticks mark labels to the original scale (dollars).See the transformation argument to effect, and the type and rescale.axisarguments to the plot method for effect objects.
The effect objects contain a variety of information that can be used to construct custom displays, and print and summary methods for effect objects present tables of the results (note that an interaction may be specified either with asterisks or colons -they are treated as equivalent since lower-order terms are in any event absorbed):  By default, the print and summary methods express effects on the scale of the response, but this behaviour may be modified via the argument type.In general, I find tabular displays of effects less satisfactory than graphical displays, but when the predictors are factors with relatively few levels, tabular displays may suffice.
To calculate fitted values for effect displays, the effect function implements the strategy for "safe prediction" described by Hastie (1992: Sec. 7.3.3;see also Venables and Ripley, 2002: Sec. 6.4).This is an issue for terms in a linear or generalised linear model with bases that depend upon the data, such as orthogonal polynomial regressors or regression splines.In such cases, naive generation of a new model matrix to be used along with the coefficients for the original fit produces incorrect results.
A simple example of effect plots incorporating regression splines, orthogonal-polynomial regressors, and a transformed predictor appears in Figure 4.The data for this example pertain to the rated prestige of 102 Canadian occupations.The prestige of the occupations is regressed on three predictors, all derived from the 1971 Census of Canada: the average income of occupational incumbents, in dollars (represented in the model as the log of income); the average education of occupational incumbents, in years (represented by a B-spline with three degrees of freedom); and the percentage of occupational incumbents who were women (represented by an orthogonal polynomial of degree two): > data(Prestige) > library(splines) # for bs > prestige.mod<-lm(prestige ~log(income) + bs(education, df=3) + poly(women, 2), + data=Prestige) > summary(prestige.mod)Call: lm(formula = prestige ~log(income) + bs(education, df = 3) + poly(women, 2), data = Prestige) > plot(all.effects(prestige.mod,default.levels=50))1:log(income) 2:bs(education, df = 3) 3:poly(women, 2) Selection: 1 1:log(income) 2:bs(education, df = 3) 3:poly(women, 2) Selection: 2 1:log(income) 2:bs(education, df = 3) 3:poly(women, 2) Selection: 3 1:log(income) 2:bs(education, df = 3) 3:poly(women, 2) Selection: 0 > When a model includes interactions between or among covariates, consideration should be given to the values at which the covariates are set in effect displays and to the form of the display.Consider the following example, taken from research by Cowles and Davis (1987) on volunteering for psychological experiments (discussed in greater detail in Fox, 1987).The response variable in the study is dichotomous: whether or not each of 1421 subjects volunteered to participate in an experiment.The authors modelled the data in a logistic regression of volunteering on the factor sex, the personality dimensions neuroticism and extraversion, and the product of neuroticism and extraversion.The two covariates each can take on integer values between 0 and 24.

Format
This data frame contains the following columns: neuroticism scale from Eysenck personality inventory.
extraversion scale from Eysenck personality inventory.sex a factor with levels: female; male.
volunteer volunteeing, a factor with levels: no; yes.

Source
Cowles, M. and C. Davis (1987) The subject matter of psychology: Volunteers.British Journal of Social Psychology 26, 97-102.

effect
Functions For Constructing Effect Plots Description effect constructs an "effect" object for a term (usually a high-order term) in a linear or generalized linear model, absorbing the lower-order terms marginal to the term in question, and averaging over other terms in the model.
all.effects identifies all of the high-order terms in a model and returns a list of "effect" objects (i.e., an object of type "effect.list").se if TRUE, the default, calculate standard errors and confidence limits for the effects.confidence.levellevel 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 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.

Usage
typical a function to be applied to the columns of the model matrix over which the effect is "averaged"; the default is mean.
... arguments to be passed down.
x an object of type "effect".row.names,optional not used.

Details
Normally, the functions to be used directly are all.effects,to return a list of high-order effects, and the generic plot function to plot the effects.(see plot.effect.listand plot.effect).Plots are drawn using the xyplot function in the lattice package.Effects may also be printed (implicitly or explicitly via print) or summarized (using summary) (see print.effect.list,summary.effect.list,print.effect,and summary.effect).
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.
In calculating effects, the strategy for 'safe' prediction described in Hastie (1992: Sec. 7.3.3) is employed.
Value effect returns an "effect" object with the following components: term the term to which the effect pertains.formula the complete model formula.
response a character string giving the response variable.
variables a list with information about each predictor, including its name, whether it is a factor, and its levels or values.fit 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.
x a data frame, the columns of which are the predictors, and the rows of which give all combinations of values of the 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.
se a vector of standard errors for the effect, on the scale of the linear predictor.
lower, upper one-column matrices of confidence limits, on the scale of the linear predictor.confidence.levelcorresponding to the confidence limits.transformation a two-element list, with element link giving the link function, and element inverse giving the inverse-link (mean) function.

Arguments
x an object of type "effect", "summary.effect",or "effect.list",as appropriate.
object an object of type "effect" or "effect.list",as appropriate.
type if "response" (the default), effects are printed or the vertical axis is labelled on the scale of the response variable; if "link", effects are printed or the vertical axis labelled on the scale of the linear predictor.
x.var the index (number) of the covariate or factor to place on the horizontal axis of each panel of the effect plot.rescale.axisif TRUE (the default), the tick marks on the vertical axis are labelled on the response scale (e.g., the probability scale for effects computed on the logit scale for a binomial GLM).
selection the optional index (number) or quoted name of the effect in an effect list to be plotted; if not supplied, a menu of high-order terms is presented.
... arguments to be passed down.

Details
In a generalized linear model, by default, the print and summary methods for effect objects print the computed effects on the scale of the response variable using the inverse of the link function.In a logit model, for example, this means that the effects are expressed on the probability scale.
By default, effects in a GLM are plotted on the scale of the linear predictor, but the vertical axis is labelled on the response scale.This preserves the linear structure of the model while permitting interpretation on what is usually a more familiar scale.This approach may also be used with linear models, for example to display effects on the scale of the response even if the data are analyzed on a transformed scale, such as log or square-root.

Value
The summary method for "effect" objects returns a "summary.effect"object with the following components (those pertaining to confidence limits need not be present): header a character string to label the effect.
effect an array containing the estimated effect.
lower.header a character string to label the lower confidence limits.
lower an array containing the lower confidence limits.
upper.header a character string to label the upper confidence limits.
upper an array containing the upper confidence limits.

Figure 1 :
Figure 1: Effect display for the interaction of colour and age in the logit model fit to the Arrests data.The vertical axis is labelled on the probability scale, and a 95-percent pointwise confidence interval is drawn around the estimated effect.

Figure 3 :
Figure3: Plotting on the scale of the response (the probability scale) -compare to Figure2.Note that the tick marks on the vertical axis are equally spaced; that the lines in the plot are not quite straight; and that corresponding lines are not quite parallel across panels.

Figure 4 :
Figure 4: Effect plots for the predictors of prestige in the Canadian occupational prestige data.The model includes the log of income, a B-spline in education, and a quadratic in women.

Figure 5 :Figure 6 :
Figure 5: Default effect display for the interaction between two covariatesneuroticism and extraversion in Cowles and Davis's volunteering data.
Figure 2: A graph of the colour by year by age effect, which corresponds to a term not in the model.The model includes colour:year and colour:age interactions, but not the colour:year:age interaction. : effect("colour:year:age", arrests.mod,xlevels = list(age = 15:45)) Type of occupation.A factor with levels (note: out of order): bc, Blue Collar; prof, Professional, Managerial, and Technical; wc, White Collar.Census of Canada.Vol. 3, Part 6. Statistics Canada [pp. 19-1-19-21].Personal communication from B. Blishen, W. Carroll, and C. Moore, Departments of Sociology, York University and University of Victoria. type The default is the predictor with the largest number of levels or values.z.varthe index (number) of the covariate or factor for which individual lines are to be drawn in each panel of the effect plot.The default is the predictor with the smallest number of levels or values.This argument is only used if multiline = TRUE.multilineifTRUE,eachpanel of the display represents combinations of values of two predictors, with one predictor (corresponding to x.var) on the horzontal axis, and the other (corresponding to z.var) used to define lines in the graph; defaults to TRUE if there are no standard errors in the object being plotted, and FALSE otherwise.rugifTRUE,thedefault,a rug plot is shown giving the marginal distribution of the predictor on the horizontal axis, if this predictor is a covariate.xlabthelabelforthehorizontal axis of the effect plot; if missing, the function will use the name of the predictor on the horizontal axis.ylab the label for the vertical axis of the effect plot; the default is the response variable for the model from which the effect was computed.isusedtoplot effects, colors[2] to plot confidence bands.In a mulitline plot, the successive colors correspond to the levels of the z.var covariate or factor.symbols,linescorresponding to the levels of the z.var covariate or factor on a multiline plot.These arguments are used only if multiline = TRUE; in this case a legend is drawn at the top of the display.NULL, the default, then the vertical axes are scaled from the data.factor.namesalogicalvalue, default TRUE, that controls the inclusion of factor names in conditioning-variable labels.ticksatwo-item list controlling the placement of tick marks on the vertical axis, with elements at and n.If at=NULL (the default), the program attempts to find 'nice' locations for the ticks, and the value of n (default, 5) gives the approximate number of tick marks desired; if at is non-NULL, then the value of n is ignored.alternatingif TRUE (the default), the tick labels alternate by panels in multi-panel displays from left to right and top to bottom; if FALSE, tick labels appear at the bottom and on the left.