Unit VIII Discussion Board RCH

7.1 Linear Regression Models

As mentioned, linear least-squares regression is typically taken up in a basic statistics course. The normal linear regression model is written

yi=β0+β1x1i+β2x2i+⋯+βkxki+εi=E(yi)+εi

(7.1)

where yi is the value of the response variable for the ith of n independently sampled observations; x1i, x2i,…, xki are the values of k explanatory variables; and the errors εi are normally and independently distributed with 0 means and constant variance, εi ∼ NID(0,σε2). Both y and the xs are numeric variables, and the model assumes that the average value E(y) of yis a linear function—that is, a simple weighted sum—of the xs. 1  If there is just one x (i.e., if k = 1), then  Equation 7.1  is called the linear simple regression model; if there are more than one x (k ≥ 2), then it is called the linear multiple regression model.

The normal linear model is optimally estimated by the method of least squares, producing the fitted model

yi=b0+b1x1i+b2x2i+⋯+bkxki+ei=y^i+ei

where y^i is the fitted value and ei the residual for observation i. The least-squares criterion

FIGURE 7.1: The Linear Regression dialog for Duncan’s occupational prestige data.

selects the values of the bs that minimize the sum of squared residuals, ∑ei2. The least-squares regression coefficients are easily computed, and, in addition to having desirable statistical properties under the model (such as efficiency and unbias), statistical inference based on the least-squares estimates is very simple (see, e.g., the references given at the beginning of the chapter).

The simplest way to fit a linear regression model in the R Commander is by the Linear Regression dialog. To illustrate, I’ll use Duncan’s occupational prestige data (introduced in  Chapter 4 ). Duncan’s data set resides in the car package, and so I can read the data into the R Commander via Data > Data in packages > Read data from an attached package (see  Section 4.2.4 ). Then selecting Statistics > Fit models > Linear regression produces the dialog in  Figure 7.1 . To complete the dialog, I click on prestige in the Response variable list, and Ctrl-click on education and income in the Explanatory variables list. Finally, pressing the OK button produces the output shown in  Figure 7.2 .

The commands generated by the Linear Regression dialog use the lm (linear model) function in R to fit the model, creating RegModel.1, and then summarize the model to produce printed output. The summary output includes information about the distribution of the residuals; coefficient estimates, their standard errors, t statistics for testing the null hypothesis that each population regression coefficient is 0, and the two-sided p-values for these tests; the standard deviation of the residuals (“residual standard error”) and residual degrees of freedom; the squared multiple correlation, R2, for the model and R2 adjusted for degrees of freedom; and the omnibus F test for the hypothesis that all population slope coefficients (here the coefficients of education and income) are 0 (H0: β1 = β2 = 0, for the example).

This is more or less standard least-squares regression output, similar to printed output produced by almost all statistical packages. What is unusual is that in addition to the printout in  Figure 7.2 , the R Commander creates and retains a linear model object on which I can perform further computations, as illustrated later in this chapter.

The Model button in the R Commander toolbar now reads RegModel.1, rather than <No active model>, as it did at the beginning the session. Just as you can choose among data sets residing in memory (if there are more than one) by pressing the Data set button in the toolbar, you can similarly choose among statistical models (if there are more than one) by pressing the Model button. Equivalently, you can pick Models > Select active model from the R Commander menus. Moreover, the R Commander takes care of coordinating data sets and models, by associating each statistical model with the data set to which it is fit. Consequently, selecting a statistical model makes the data set to which it was fit the active data set, if that isn’t already the case.

FIGURE 7.2: Output from Duncan’s regression of occupational prestige on income and education, produced by the Linear Regression dialog.

The variable lists in the Linear Regression dialog in  Figure 7.1  include only numeric variables. For example, the factor type (type of occupation) in Duncan’s data set, with levels “bc” (blue-collar), “wc” (white-collar), and “prof”(professional, technical, or managerial), doesn’t appear in either variable list. Moreover, the explanatory variables that are selected enter the model linearly and additively. The Linear Model dialog, described in the next section, is capable of fitting a much wider variety of regression models.

In completing the Linear Regression dialog in  Figure 7.1 , I left the name of the model at its default, RegModel.1. The R Commander generates unique model names automatically during a session, each time incrementing the model number (here 1).

I also left the Subset expression at its default, <all valid cases>. Had I instead entered type == “bc”, 2  for example, the regression model would have been fit only to blue-collar occupations. As in this example, the subset expression can be a logical expression, returning the value TRUE or FALSE for each case (see  Section 4.4.2 ), a vector of case indices to include, 3  or a negative vector of case indices to exclude. For example, 1:25 would include the first 25 occupations, while -c(6, 16) would exclude occupations 6 and 16. 4  All of the statistical modeling dialogs in the R Commander allow subsets of cases to be specified in this manner.

7.2 Linear Models with Factors*

Like the Linear Regression dialog described in the preceding section, the Linear Model dialog can fit additive linear regression models, but it is much more flexible: The Linear Model dialog accommodates transformations of the response and explanatory variables, factors as well as numeric explanatory variables on the right-hand-side of the regression model, nonlinear functions of explanatory variables expressed as polynomials and regression splines, and interactions among explanatory variables. All this is accomplished by allowing the user to specify the model as an R linear-model formula. Linear-model formulas in R are inherited from the S programming language (Chambers and Hastie, 1992), and are a version of notation for expressing linear models originally introduced by Wilkinson and Rogers (1973).

7.2.1 Linear-Model Formulas

An R linear-model formula is of the general form response-variable ∼ linear-predictor . The tilde (~) in a linear-model formula can be read as “is regressed on.” Thus, in this general form, the response variable is regressed on a linear predictor comprising the terms in the right-hand side of the model.

The left-hand side of the model formula, response-variable , is an R expression that evaluates to the numeric response variable in the model, and is usually simply the name of the response variable—for example, prestige in Duncan’s regression. You can, however, transform the response variable directly in the model formula (e.g., log10(income)) or compute the response as a more complex arithmetic expression (e.g., log(investment.income + hourly.wage.rate*hours.worked). 5

The formulation of the linear predictor on the right-hand side of a model formula is more complex. What are normally arithmetic operators (+, -, *, /, and ^) in R expressions have special meanings in a model formula, as do the operators : (colon) and %in%. The numeral 1 (one) may be used to represent the regression constant (i.e., the intercept) in a model formula; this is usually unnecessary, however, because an intercept is included by default. A period (.) represents all of the variables in the data set with the exception of the response. Parentheses may be used for grouping, much as in an arithmetic expression.

In the large majority of cases, you’ll be able to formulate a model using only the operators + (interpreted as “and”) and * (interpreted as “crossed with”), and so I’ll emphasize these operators here. The meaning of these and the other model-formula operators are summarized and illustrated in  Table 7.1 . Especially on first reading, feel free to ignore everything in the table except +, :, and * (and : is rarely used directly).

A final formula subtlety: As I’ve explained, the arithmetic operators take on special meanings on the right-hand side of a linear-model formula. A consequence is that you can’t use these operators directly for arithmetic. For example, fitting the model savings ~ wages + interest + dividends estimates a separate regression coefficient for each of wages, interest, and dividends. Suppose, however, that you want to estimate a single coefficient for the sum of these variables—in effect, setting the three coefficients equal to each other. The solution is to “protect” the + operator inside a call to the I (identity or inhibit) function, which simply returns its argument unchanged: 6  savings ∼ I(wages + interest + dividends). This formula works as desired because arithmetic operators like + have their usual meaning within a function call on the right-hand side of the formula—implying, incidentally, that savings ∼ log10(wages + interest + dividends) also works as intended, estimating a single coefficient for the log base 10 of the sum of wages, interest, and dividends.

TABLE 7.1: Operators and other symbols used on the right-hand side of R linear-model formulas.

Operator	Meaning	Example	Interpretation
+	and	x1 + x2	x1 and x2
:	interaction	x1:x2	interaction of x1 and x2
*	crossing	x1*x2	x1 crossed with x2 (i.e., x1 + x2 + x1:x2)
-	remove	x1 - 1	regression through the origin (for numeric x1)
^k	cross to order k	(x1 + x2 + x3)^2	same as x1*x2 + x1*x3 + x2*x3
%.in%.	nesting	province %in% country	province nested in country
/	nesting	country/province	same as country + province %in% country
Symbol	Meaning	Example	Interpretation
1	intercept	x1 - 1	suppress the intercept
.	everything but the response	y ~ .	regress y on everything else
( )	grouping	x1*(x2 + x3)	same as x1*x2 + x1*x3

The symbols x1, x2, and x3 represent explanatory variables and could be either numeric or factors.

7.2.2 The Principle of Marginality

Introduced by Nelder (1977), the principle of marginality is a rule for formulating and interpreting linear (and similar) statistical models. According to the principle of marginality, if an interaction, say x1:x2, is included in a linear model, then so should the main effects, x1 and x2, that are marginal to—that is lower-order relatives of—the interaction. Similarly, the lower-order interactions x1:x2, x1:x3, and x2:x3 are marginal to the three-way interaction x1:x2:x3. The regression constant (1 in an R model formula) is marginal to every other term in the model. 7

It is in most circumstances difficult in R to formulate models that violate the principle of marginality, and trying to do so can produce unintended results. For example, although it may appear that the model y ∼ f*x - x - 1, where f is a factor and x is a numeric explanatory variable, 8  violates the principle of marginality by removing the regression constant and x slope, the model that R actually fits includes a separate intercept and slope for each level of the factor f. Thus, the model y ∼ f*x - x - 1 is equivalent to (i.e., an alternative parametrization of) y ∼ f*x. It is almost always best to stay away from such unusual model formulas.

7.2.3 Examples Using the Canadian Occupational Prestige Data

For concreteness, I’ll formulate several linear models for the Canadian occupational prestige data (introduced in  Section 4.2.3  and described in  Table 4.2  on  page 61 ), regressing prestige on income, education, women (gender composition), and type (type of occupation). The last variable is a factor (categorical variable) and so it cannot enter into the linear model directly. When a factor is included in a linear-model formula, R generates contrasts to represent the factor—one fewer than the number of levels of the factor. I’ll explain how this works in greater detail in Section 7.2.4, but the default in the R Commander (and R more generally) is to use 0/1 dummy-variable regressors, also called indicator variables.

A version of the Canadian occupational prestige data resides in the data frame Prestige in the car package,9 and it’s convenient to read the data into the R Commander from this source via Data > Data in packages > Read data from an attached package. Prestige replaces Duncan as the active data set.

Recall that 4 of the 102 occupations in the Prestige data set have missing values (NA) for occupational type. Because I will fit several regression models to the Prestige data, not all of which include type, I begin by filtering the data set for missing values, selecting Data > Active data set > Remove cases with missing data (as described in Section 4.5.2).

Moreover, the default alphabetical ordering of the levels of type—“bc” , “prof” , “wc”—is not the natural ordering, and so I also reorder the levels of this factor via Data > Manage variables in active data set > Reorder factor levels to “bc”, “wc”, “prof” (see Section 3.4). This last step isn’t strictly necessary, but it makes the data analysis easier to follow.

I first fit an additive dummy regression to the Canadian prestige data, employing the model formula prestige ∼ income + education + women + type. To do so, I select Statistics > Fit models > Linear model from the R Commander menus, producing the dialog box in Figure 7.3. The automatically supplied model name is LinearModel.2, reflecting the fact that I have already fit a statistical model in the session, RegModel.1 (in Section 7.1).

Most of the structure of the Linear Model dialog is common to statistical modeling dialogs in the R Commander. If the response text box to the left of the ∼ in the model formula is empty, double-clicking on a variable name in the variable list box enters the name into the response box; thereafter, double-clicking on variable names enters the names into the right-hand side of the model formula, separated by +s (if no operator appears at the end of the partially completed formula). You can enter parentheses and operators like + and * into the formula using the toolbar in the dialog box.10 You can also type directly into the model-formula text boxes. In Figure 7.3, I simply double-clicked successively on prestige, education, income, women, and type.11 Clicking OK produces the output shown in Figure 7.4.

I already explained the general format of linear-model summary output in R. What’s new in Figure 7.4 is the way in which the factor type is handled in the linear model: Two dummy-variable regressors are automatically created for the three-level factor type. The first dummy regressor, labelled type[T.wc] in the output, is coded 1 when type is “wc”and 0 otherwise; the second dummy regressor, type[T.prof], is coded 1 when type is “prof” and 0 otherwise. The first level of type—“bc”—is therefore selected as the reference or baseline level, coded 0 for both dummy regressors.12

Consequently, the intercept in the linear-model output is the intercept for the “bc” reference level of type, and the coefficients for the other levels give differences in the intercepts between each of these levels and the reference level. Because the slope coefficients for the numeric explanatory variables education, income, and women in this additive model do not vary by levels of type, the dummy-variable coefficients are also interpretable as the average difference between each other level and “bc” for any fixed values of education, income, and women.

FIGURE 7.3: Linear Model dialog completed to fit an additive dummy-variable regression of prestige on the numeric explanatory variables education, income, and women, and the factor type.

To illustrate a structurally more complex, nonadditive model, I respecify the Canadian occupational prestige regression model to include interactions between type and education and between type and income, in the process removing women from the model—in the initial regression, the coefficient of women is small with a large p-value. 13  The Linear Model dialog (not shown) reopens in its previous state, with the model name incremented to LinearModel.3. To fit the new model, I modify the formula to read prestige ∼ type*education + type*income. Clicking OK produces the output in  Figure 7.5 .

With interactions in the model, there are different intercepts and slopes for each level of type. The intercept in the output—along with the coefficients for education and income— pertains to the baseline level “bc” of type. Other coefficients represent differences between each of the other levels and the baseline level. For example, type[T.wc] = –33.54 is the difference in intercepts between the “wc” and “bc” levels of type; 14  similarly, the interaction coefficient type[T.wc]:education = 4.291 is the difference in education slopes between the “wc” and “bc” levels. The complexity of the coefficients makes it difficult to understand what the model says about the data;  Section 7.6  shows how to visualize terms such as interactions in a complex linear model.

FIGURE 7.4: Output for the linear model prestige ∼ income + education + women + type fit to the Prestige data.

FIGURE 7.5: Output for the linear model prestige ∼ type*education + type*income fit to the Prestige data.

TABLE 7.2: Contrast-regressor codings for type generated by contr.Treatment, contr.Sum, contr.poly,, and contr.Helmert.

		Levels of type
Function	Contrast Names	“bc”	“wc”	“prof”
contr.Treatment	type[T.wc]	0	1	0
	type[T.prof]	0	0	1
contr.Sum	type[S.wc]	1	0	-1
	type[S.prof]	0	1	-1
contr.poly	type.L	−1/2	0	1/2
	type.Q	1/6	−2/6	1/6
contr.Helmert	type[H.1]	-1	1	0
	type[H.2]	-1	-1	2

7.2.4 Dummy Variables and Other Contrasts for Factors

By default in the R Commander, factors in linear-model formulas are represented by 0/1 dummy-variable regressors generated by the contr.Treatment function in the car package, picking the first level of a factor as the baseline level. 15  This contrast coding, along with some other choices, is shown in  Table 7.2 , using the factor type in the Prestige data set as an example.

The function contr.Sum from the car generates so-called “sigma-constrained” or “sum-to-zero” contrast regressors, as are used in traditional treatments of analysis of variance. 16  The standard R function contr.poly generates orthogonal-polynomial contrasts—in this case, linear and quadratic terms for the three levels of type; in the R Commander, contr.poly is the default choice for ordered factors. Finally, contr.Helmert generates Helmert contrasts, which compare each level to the average of those preceding it.

Selecting Data > Manage variables in active data set > Define contrasts for a factor produces the dialog box on the left of  Figure 7.6 . The factor type is preselected in this dialog because it’s the only factor in the data set. You can use the radio buttons to choose among treatment, sum-to-zero, Helmert, and polynomial contrasts, or define customized contrasts by selecting Other, as I’ve done here.

Clicking OK leads to the sub-dialog shown on the right of  Figure 7.6 . I change the default contrast names, .1 and .2, to [bc.v.others] and [wc.v.prof], and then fill in the contrast coefficients (i.e., the values of the contrast regressors). This choice produces contrast regressors named type[bc.v.others] and type[wc.v.prof], to be used when the factor type in the Prestige data set appears in a linear-model formula. Contrasts defined directly in this manner must be linearly independent and are simplest to interpret if they obey two additional rules: (1) The coefficients for each contrast should sum to 0, and (2) each pair of contrasts should be orthogonal (i.e., the products of corresponding coefficients for each pair of contrasts sum to 0).

FIGURE 7.6: The Set Contrasts for Factor dialog box (left) and the Specify Contrasts sub-dialog (right), creating contrasts for the factor type in the Prestige data set.

To see how these contrasts are reflected in the coefficients of the model, I refit the additive regression of prestige on education, income, women, and type, producing the output in Figure 7.7. The first contrast for type estimates the difference between “bc” and the average of the other two levels of type, holding the other explanatory variables constant, while the second contrast estimates the difference between “wc” and “prof”. This alternative contrast coding for type produces different estimates for the intercept and type coefficients from the dummy-regressor coding for typein Figure 7.4 (on page 136), but the two models have the same fit to the data (e.g., R2 = 0.8349).17

FIGURE 7.7: Output for the linear model prestige ∼ income + education + women + type fit to the Prestige data, using customized contrasts for type.

7.3 Fitting Regression Splines and Polynomials*

The second formula toolbar in the Linear Model dialog makes it easy to add nonlinear polynomial-regression and regression-spline terms to a linear model.

7.3.1 Polynomial Terms

Some simple nonlinear relationships can be represented as low-order polynomials, such as a quadratic term, using regressors x and x2 for a numeric explanatory variable x, or a cubic term, using x, x2, and x3. The resulting model is nonlinear in the explanatory variable x but linear in the parameters (the βs). R and the R Commander support both orthogonal and “raw” polynomials in linear model formulas.18

To add a polynomial term to the right-hand side of the model, single-click on a numeric variable in the Variables list box, and then press the appropriate toolbar button (either orthogonal polynomial or raw polynomial, as desired). There is a spinner in the Linear Model dialog for the degree of a polynomial term, and the default is 2 (i.e., a quadratic).

For example, inspection of the data (e.g., in a component-plus-residual plot, discussed in Section 7.8)19 suggests that there may be a quadratic partial relationship between prestige and women in the regression of prestige on education, income, and women for the Canadian occupational prestige data.20 I specify this quadratic relationship in the Linear Model dialog in Figure 7.8, using a raw second-degree polynomial, and producing the output in Figure 7.9. The quadratic coefficient in the model turns out not to be statistically significant (p = 0.15).

7.3.2 Regression Splines

Regression splines are flexible functions capable of representing a wide variety of nonlinear patterns in a model that, like a regression polynomial, is linear in the parameters. Both B-splines and natural splines are supported by the R CommanderLinear Model dialog. Adding a spline term to the right-hand side of a linear model is similar to adding a polynomial  prestige ∼ poly(women, degree=2, raw=TRUE) + ns(education, df=5) + ns(income, df=5), regressing prestige on a quadratic in women and 5-df natural splines in education and income. The output for the resulting regression model isn’t shown because the model requires graphical interpretation (see Section 7.6): The coefficient estimates for the regression splines are not simply interpretable. 21

FIGURE 7.8: Linear Model dialog with a polynomial (quadratic) term for women in the regression of prestige on education, income, and women using the Prestige data set.

FIGURE 7.9: Output from the regression of prestige on education, income, and a quadratic in women for the Prestigedata.

FIGURE 7.10: Linear Model dialog showing regression-spline and polynomial terms for the regression of prestige on education, income, and women in the Prestige data set.

7.4 Generalized Linear Models*

Briefly, generalized linear models (or GLMs), introduced in a seminal paper by Nelder and Wedderburn (1972), consist of three components:

1.  A random component specifying the distribution of the response y conditional on explanatory variables. Traditionally, the random component is a member of an exponential family—the Gaussian (normal), binomial, Poisson, gamma, or inverse Gaussian families—but both the theory of generalized linear models and their implementation in R are now more general: In addition to the traditional exponential families, R provides for quasi-binomial and quasi-Poisson families that accommodate “over-dispersed” binomial and count data.

2.  A linear predictor

ηi=β0+β1x1i+β2x2i+⋯+βkxki

on which the expectation of the response variable μi = E(yi) for the ith of n independent observations depends, where the regressors xji are prespecified functions of the explanatory variables—numeric explanatory variables, dummy regressors representing factors, interaction regressors, and so on, exactly as in the linear model.

3.  A prespecified invertible link function g(.) that transforms the expectation of the response to the linear predictor, g(μi) = ηi, and thus μi = g−1(ηi). R implements identity, inverse, log, logit, probit, complementary log-log, square root, and inverse square links, with the applicable links varying by distributional family.

The most common GLM beyond the normal linear model (i.e., the Gaussian family paired with identity link) is the binomial logit model, suitable for dichotomous (two-category) response variables. For an illustration, I’ll use data collected by Cowles and Davis (1987) on volunteering for a psychological experiment, where the subjects of the study were students in a university introductory psychology class.

The data for this example are contained in the data set Cowles in the car package, 22  which includes the following variables: neuroticism, a personality dimension with integer scores ranging from 0 to 24; extraversion, another personality dimension, also with scores from 0 to 24; sex, a factor with levels “female” and “male”; and volunteer, a factor with levels “no” and “yes” .

In analyzing the data, Cowles and Davis performed a logistic regression of volunteering on sex and the linear-by-linear interaction between neuroticism and extraversion. To fit Cowles and Davis’s model, I first read the data from the car package in the usual manner, making Cowles the active data set in the R Commander. Then I select Statistics > Fit models > Generalized linear model, producing the dialog box in  Figure 7.11 .

The Generalized Linear Model dialog is very similar to the Linear Model dialog of the preceding section: The name of the model at the top (GLM.7) is automatically generated, and you can change it if you wish. Double-clicking on a variable in the list box enters the variable into the model formula. There are toolbars for entering operators, regression splines, and polynomials into the model formula, and there are boxes for subsetting the data set and for specifying prior weights.

FIGURE 7.11: Genealized Linear Model dialog box for Cowles and Davis’s logistic regression.

What’s new in the Generalized Linear Model dialog are the Family and Link function list boxes, as are appropriate to a GLM. Families and links are coordinated: Double-clicking on a distributional family changes the available links. In each case, the canonical link for a particular family is selected by default. The initial selections are the binomial family and corresponding canonical logit link, which are coincidentally what I want for the example.

I proceed to complete the dialog by double-clicking on volunteer in the variable list, making it the response variable; then double-clicking on sex and on neuroticism; clicking the * button in the toolbar; and finally double-clicking on neuroticism—yielding the model formula volunteer ∼ sex + neuroticism*extraversion. As in the Linear Model dialog, an alternative is to type the formula directly.

Appropriate responses for a binomial logit model include two-level factors (such as volunteer in the current example), logical variables (i.e., with values FALSE and TRUE), and numeric variables with two unique values (most commonly 0 and 1). In each case, the logit model is for the probability of the second of the two values—the probability that volunteer is “yes” in the example.

Clicking the OK button produces the output in  Figure 7.12 . The Generalized Linear Model dialog uses the R glmfunction to fit the model. The summary output for a generalized linear model is very similar to that for a linear model, including a table of estimated coefficients along with their standard errors, z values (Wald statistics) for testing that the coefficients are 0, and the two-sided p-values for these tests. For a logistic regression, the R Commander also prints the exponentiated coefficients, interpretable as multiplicative effects on the odds scale—here the odds of volunteering, Pr(“yes”)/Pr(“no”).

The Wald z tests suggest a statistically significant interaction between neuroticism and extraversion, as Cowles and Davis expected, and a significant sex effect, with men less likely to volunteer than women who have equivalent scores on the personality dimensions. Because it’s hard to grasp the nature of the interaction directly from the coefficient estimates, I’ll return to this example in  Section 7.6 , where I’ll plot the fitted model.

Although I’ve developed just one example of a generalized linear model in this section—a logit model for binary data—the R Commander Generalized Linear Model dialog is more flexible:

•  The probit and complementary log-log (cloglog) link functions may also be used with binary data, as alternatives to the canonical logit link.

•  The binomial family may also be used when the value of the response variable for each case (or observation) represents the proportion of “successes” in a given number of binomial trials, which may also vary by case. In this setting, the left-hand side of the model formula should give the proportion of successes, which could be computed as successes/trials (imagining that there are variables with these names in the active data set) directly in the left-hand box of the model formula, and the variable representing the number of trials for each observation (e.g., trials)should be given in the Weights box.

•  Alternatively, for binomial data, the left-hand side of the model may be a two-column matrix specifying, respectively, the numbers of successes and failures for each observation, by typing, e.g., cbind(successes, failures) (again, imagining that these variable are in the active data set) into the left-hand-side box of the model formula.

•  Other generalized linear models are specified by choosing a different family and corresponding link. For example, a Poisson regression model, commonly employed for count data, may be fit by selecting the poisson family and canonical log link (or, to get typically more realistic coefficient standard errors, by selecting the quasipoisson family with the loglink).

FIGURE 7.12: Output from Cowles and Davis’s logistic regression (volunteer ∼ sex + neuroticism*extraversion).

7.5 Other Regression Models*

In addition to linear regression, linear models, and generalized linear models, the R Commander can fit multinomial logit models for categorical response variables with more than two categories (via Statistics > Fit models > Multinomial logit model), and ordinal regression models for ordered multi-category responses, including the proportional-odds logit modeland the ordered probit model (Statistics > Fit models > Ordinal regression model). Although I won’t illustrate these models here, many of the menu items in the Models menu apply to these classes of models. Moreover (as I will show in Chapter 9), R Commander plug-in packages can introduce additional classes of statistical models.

7.6 Visualizing Linear and Generalized Linear Models*

Introduced by Fox (1987), effect plots are graphs for visualizing complex regression models by focusing on particular explanatory variables or combinations of explanatory variables, holding other explanatory variables to typical values. One strategy is to focus successively on the explanatory variables in the high-order terms of the model—that is, terms that aren’t marginal to others (see Section 7.2.2).

In the R Commander, effect displays can be drawn for linear, generalized linear, and some other statistical models via Models > Graphs > Effect plots. Figure 7.13 shows the resulting dialog box for Cowles and Davis’s logistic regression from the previous section, GLM.7, which is the current statistical model in the R Commander. By default, the dialog offers to plot all high-order terms in the model—in this case, the sex main effect and the neuroticism-by-extraversion interaction. You may alternatively pick a subset of Predictors (explanatory variables) to plot.23 For a linear or generalized linear model, there’s also a check box for plotting partial residuals, unchecked by default, along with a slider for the span of a smoother fit to the residuals. Partial residuals and the accompanying smoother can be useful for judging departures from the functional form of the specified model, as I’ll illustrate later in this section.

Clicking OK produces the graph in Figure 7.14: The left-hand panel shows the sex main effect, with neuroticismand extraversion set to average levels. The right-hand panel shows the neuroticism-by-extraversioninteraction, for a group composed of males and females in proportion to their representation in the data set. In both graphs, the vertical volunteer axis is drawn on the logit scale but the tick-mark labels are on the estimated probability scale—that is, they represent the estimated probability of volunteering.24

In the plot of the interaction, the horizontal axis of each panel is for neuroticism, while extraversion takes on successively larger values across its range, from the lower-left panel to the upper-right panel. The value of extraversion for each panel is represented by the small vertical line in the strip labelled extraversion at the top of the panel.

FIGURE 7.13: Model Effect Plots dialog box for Cowles and Davis’s logistic regression (volunteer ∼ sex + neuroticism*extraversion).

The lines in the panels represent the combined effect of neuroticism and extraversion, and are computed using the estimated coefficients for the neuroticism:extraversion interaction along with the coefficients for the neuroticism and extraversion “main effects,” which are marginal to the interaction. It’s clear that there’s a positive relationship between volunteering and neuroticism at the lowest level of extraversion, but that this relationship becomes negative at the highest level of extraversion.

The error bars in the effect plot for sex and the gray bands in the plot of the neuroticism-by-extraversioninteraction represent point-wise 95% confidence intervals around the estimated effects. The “rug-plot” at the bottom of each panel in the display of the neuroticism-by-extraversion interaction shows the marginal distribution of neuroticism, with the lines slightly displaced to decrease over-plotting. The rug-plot isn’t terribly useful here, because neuroticism just takes on integer scores.

FIGURE 7.14: Effect plots for the high-order terms in Cowles and Davis’s logistic regression (volunteer ∼ sex + neuroticism*extraversion). The graphs are shown in monochrome; they (and the other effect plots in the chapter) were originally in color.

FIGURE 7.15: Model Effect Plots dialog box for LinearModel.2 (prestige ∼ education + income + women + type) fit to the Prestige data.

7.6.1 Partial Residuals in Effect Plots

Adding partial residuals to effect plots of numeric explanatory variables in linear and generalized linear models can be an effective tool for judging departures from the functional form (linear or otherwise) specified in the model. I’ll illustrate using the Canadian occupational prestige data. In  Sections 7.2.3  and  7.3 , I fit several models to the Prestige data, including an additive dummy-regression model (LinearModel.2 in  Figure 7.4  on  page 136 ),

prestige ∼ education + income + women + type

and a model with interactions (LinearModel.3 in  Figure 7.5  on  page 137 ),

prestige ∼ type*education + type*income

among others.

The R Commander session in this chapter is unusual in that I’ve read three data sets (Duncan, Prestige, and Cowles) and fit statistical models to each. It is much more common to work with a single data set in a session. Nevertheless, as I explained, the R Commander allows you to switch among models and data sets, and takes care of synchronizing models with the data sets to which they were fit. After making LinearModel.2 the active model, I return to the Model Effect Plots dialog, displayed in  Figure 7.15 . I check Plot partial residuals and click OK, producing  Figure 7.16 . Partial residuals are plotted for the numeric predictors but not for the factor type; this is reflected in a warning printed in the Messages pane, which I’ll simply ignore.

The solid lines in the effect plots represent the model fit to the data, while the broken lines are smooths of the partial residuals. If the lines for the fitted model and smoothed partial residuals are similar, that lends support to the specified functional form of the model. The partial residuals are computed by adding the residual for each observation to the line representing the fitted effect. It appears as if the education effect is modelled reasonably, but the income and womeneffects appear to be nonlinear.

LinearModel.3 includes interactions between type and each of education and income.  Figure 7.17  shows effect plots with partial residuals for the high-order terms in this model.

FIGURE 7.16: Effect displays with partial residuals for LinearModel.2 (prestige ∼ education + income + women + type) fit to the Prestige data.

FIGURE 7.16: Effect displays with partial residuals for LinearModel.2 (prestige ∼ education + income + women + type) fit to the Prestige data.

FIGURE 7.17: Effect displays with partial residuals for LinearModel.3 (prestige ∼ type*education + type*income) fit to the Prestige data.

Because dividing the data by type leaves relatively few points in each panel of the plots, I set the span of the smoother to a large value, 0.9. 25

The apparent nonlinearity in the relationship between prestige and income is accounted for by the interaction between income and type: 26  The right-hand display of  Figure 7.17  shows that the income slope is smaller for professional and managerial occupations (i.e., type = “prof”) than for blue-collar (“bc”) or white-collar (“wc”)occupations, and professional occupations tend to have higher incomes. The display at the left, for the education-by-type interaction, suggests that the education slope is steeper for white-collar occupations than for the other types of occupations. The smooths of the partial residuals indicate that these relationships are linear within the levels of type.

The confidence envelopes in the effect displays with partial residuals in  Figure 7.17  also make a useful pedagogical point about precision of estimation of the regression surface: Where data are sparse—or, in the extreme, absent—the regression surface is imprecisely estimated.

LinearModel.6, fit in  Section 7.3 ,

prestige ∼ poly(women, degree=2, raw=TRUE) + ns(education, df=5) + ns(income, df=5)

uses a quadratic in women along with regression splines for income and education, which should capture the unmodelled nonlinearity observed in  Figure 7.16 ; the model doesn’t include the factor type, however. I make LinearModel.6 the active model and repeat the effect plots, which are shown in  Figure 7.18 . Here, the fitted model and smoothed residuals agree well with each other.

FIGURE 7.18: Effect displays with partial residuals for LinearModel.6 (prestige ∼ ns(education, df=5) + ns(income, df=5) + poly(women, degree=2, raw=TRUE)) fit to the Prestige data.

7.7 Confidence Intervals and Hypothesis Tests for Coefficients

The Models menu includes several menu items for constructing confidence intervals and performing hypothesis tests for regression coefficients. As I explained, tests for individual coefficients in linear and generalized linear models appear in the model summaries. This section describes how to perform more elaborate tests, for example, for a related subset of coefficients.

7.7.1 Confidence Intervals

To illustrate, I’ll make Duncan’s occupational prestige regression (RegModel.1 in  Figure 7.2  on  page 131 ) the active statistical model in the R Commander, which automatically makes Duncan the active data set.

Selecting Models > Confidence intervals from the R Commander menus leads to the simple dialog box at the top of  Figure 7.19 . 27  Retaining the default 0.95 level of confidence and clicking OK produces the output at the bottom of the figure.

7.7.2 Analysis of Variance and Analysis of Deviance Tables*

You can compute an analysis of variance (ANOVA) table for a linear model or an analogous analysis of deviance table for a generalized linear model via Models > Hypothesis tests > ANOVA table. I’ll illustrate with Cowles and Davis’s logistic regression model (GLM.7 in  Figure 7.12  on  page 148 ), selecting it as the active model in the session. The ANOVA Tabledialog at the top of  Figure 7.20  offers three “types” of tests, conventionally named Types I, II, and III:

•  In addition to the intercept, there are four terms in the Cowles and Davis model: sex, neuroticism, extraversion, and the neuroticism:extraversion interaction. Type I tests are sequential, and thus test (in a short-hand terminology) sex ignoring everything else; neuroticism after sex but ignoring extraversion and the neuroticism:extraversion interaction; extraversion after sex and neuroticism ignoring the interaction; and the neuroticism:extraversion interaction after all other terms. Sequential tests are rarely sensible.

•  Type II tests are formulated in conformity with the principle of marginality ( Section 7.2.2 ): sex after all other terms, including the neuroticism:extraversion interaction; neuroticism after sex and extraversion but ignoring the neuroticism:extraversion interaction to which the neuroticism “main effect” is marginal; similarly, extraversion after sex and neuroticism but ignoring the neuroticism:extraversion interaction; and neuroticism:extraversion after all the other terms. More generally, each term is tested ignoring terms to which it is marginal (i.e., ignoring its higher-order relatives). This is generally a sensible approach and is the default in the dialog.

Type II tests for Cowles and Davis’s logistic regression are shown at the bottom of  Figure 7.20 . For a generalized linear model like Cowles and Davis’s logistic regression, the R Commander computes an analysis of deviance table with likelihood ratio tests, which entails refitting the model (in some instances twice) for each test. In this example, each likelihood ratio chi-square test has one degree of freedom because each term in the model is represented by a single coefficient. Because the neuroticism:extraversion interaction is highly statistically significant, I won’t interpret the tests for the neuroticism and extraversion main effects, which assume that the interaction is nil. The sexeffect is also statistically significant.

FIGURE 7.19: Confidence Intervals dialog and resulting output for Duncan’s occupational prestige regression (prestige ∼ education + income).

•  Type III tests are for each term after all of the others. This isn’t sensible for Cowles and Davis’s model, although it might be rendered sensible, for example, by centering each of neuroticism and extraversion at their means prior to fitting the model. Even where Type III tests can correspond to sensible hypotheses, such as in analysis of variance models, they require careful formulation in R. 28  I suggest that you avoid Type III tests unless you know what you’re doing.

FIGURE 7.20: ANOVA Table dialog and resulting Type II tests for Cowles and Davis’s logistic regression (volunteer ∼ sex + neuroticism*extraversion).

FIGURE 7.20: ANOVA Table dialog and resulting Type II tests for Cowles and Davis’s logistic regression (volunteer ∼ sex + neuroticism*extraversion).

FIGURE 7.21: Summary of a regression model fit to Duncan’s occupational prestige data forcing the education and income coefficients to be equal (prestige ∼ I(education + income)).

7.7.3 Tests Comparing Models*

The R Commander allows you to compute a likelihood ratio F test or chi-square test for two regression models, one of which is nested within the other.29 To illustrate, I return to the Duncan data set and fit a version of Duncan’s regression that sets the coefficients of education and income equal to each other, specifying the linear-model formula prestige ∼ I(income + education).30 This at least arguably makes sense in that both explanatory variables are percentages—respectively of high school graduates and of high-income earners. The resulting model, LinearModel.8, is summarized in Figure 7.21.

Picking Models > Hypothesis tests > Compare two models leads to the dialog at the top of Figure 7.22. I select the more general RegModel.1 as the first model and the more specific, constrained, LinearModel.8 as the second model, but the order of the selections is immaterial—the same F test is produced in both cases. Clicking OK in the dialog box results in the output at the bottom of Figure 7.22. The hypothesis of equal coefficients for education and income is plausible, p = 0.79—after all, the two coefficients are quite similar in the original regression (Figure 7.2 on page 131), beducation = 0.55 and bincome = 0.60.

FIGURE 7.22: Compare Models dialog and resulting output for Duncan’s occupational prestige regression, testing the equality of the education and income regression coefficients.

7.7.4 Testing Linear Hypotheses*

The menu selection Models > Hypothesis tests > Linear hypothesis allows you to formulate and test general linear hypotheses about the coefficients in a regression model. To illustrate, I’ll again use Duncan’s occupational prestige regression, RegModel.1. Figures 7.23 and 7.24 show the Test Linear Hypothesis dialog set up to test two different hypotheses, along with the corresponding output:

1.  In Figure 7.23, H0: 1 × βeducation – 1 × βincome = 0 (i.e., H0: βeducation = βincome). This is the same hypothesis that I tested by the model-comparison approach immediately above, and of course it produces the same F test.

2.  In Figure 7.24, H0: 1 × βeducation = 0, 1 × βincome = 0, which is equivalent to, and thus produces the same F test as, the omnibus null hypothesis in the linearmodel summary output, H0: βeducation = βincome = 0 (see Figure 7.2 on page 131). Because the linear hypothesis consists of two equations, the F statistic for the hypothesis has two df in the numerator.

There may be up to as many equations in a linear hypothesis as the number of coefficients in the model, with the number of equations controlled by the slider at the top of the dialog. The equations must be linearly independent of one another—that is, they may not be redundant. Initially, all of the cells in each row are 0, including the cell representing the right-hand side of the hypothesis, which is usually left at 0. The Test Linear Hypothesis dialog for a linear model provides for an optional “sandwich” estimator of the coefficient covariance matrix, which may be used to adjust statistical inference for autocorrelated or heteroscedastic errors (nonconstant error variance).

FIGURE 7.23: Testing a linear hypothesis for Duncan’s occupational prestige regression: H0: βeducation = βincome.

FIGURE 7.24: Testing a linear hypothesis for Duncan’s occupational prestige regression:H0: βeducation = βincome. = 0.

7.8 Regression Model Diagnostics*

Regression diagnostics are methods for determining whether a regression model that’s been fit to data adequately summarizes the data. For example, is a relationship that’s assumed to be linear actually linear? Do one or a small number of influential cases unduly affect the results?

Many standard methods of regression diagnostics are implemented in the R Commander Models > Numerical diagnostics and Models > Graphs menus—indeed, too many to cover in detail in this already long chapter. Luckily, most of the diagnostics dialogs are entirely straightforward, and some of the diagnostics menu items produce results directly, without invoking a dialog box. I’ll illustrate with Duncan’s occupational prestige regression (RegModel.1 in Figure 7.2on page 131). As usual, I assume that the statistical methods covered here are familiar. Regression diagnostics are taken up in many regression texts; see, in particular, Fox (2016), Weisberg (2014), Cook and Weisberg (1982), or (for a briefer treatment) Fox (1991).

The numerical diagnostics available in the R Commander include generalized variance-inflation factors (Fox and Monette, 1992) for diagnosing collinearity in linear and generalized linear models, the Breusch–Pagan test for nonconstant error variance in a linear model (Breusch and Pagan, 1979), independently proposed by Cook and Weisberg (1983), the Durbin–Watson test for autocorrelated errors in linear time-series regression (Durbin and Watson, 1950, 1951), the RESET test for nonlinearity in a linear model (Ramsey, 1969), and a Bonferonni outlier test based on the studentized residuals from a linear or generalized linear model (see, e.g., Fox, 2016, Chapter 11).

I’ll illustrate with Models > Numerical diagnostics > Breusch-Pagan test for het-eroscedasticity, leading to the dialog at the top of Figure 7.25. The default is to test for error variance that increases (or decreases) with the level of the response (through the Fitted values), but the dialog is flexible and accommodates dependence of the error variance on the explanatory variables or on a linear predictor based on any variables in the data set. I leave the dialog at its default, producing the output at the bottom of Figure 7.25. There is, therefore, no evidence that the variance of the errors in Duncan’s regression depends on the level of the response.

There are also many graphical diagnostics available through the R Commander: Basic diagnostic plots produced by the R plot function applied to a linear or generalized linear model; residual quantile-comparison plots, for example to diagnose non-normal errors in a linear model; component-plus-residual (partial-residual) plots for nonlinearity in additive linear or generalized linear models31; added-variable plots for diagnosing unusual and influential data in linear and generalized linear models; and an “influence plot”—a diagnostic graph that simultaneously displays studentized residuals, hat-values (leverages), and Cook’s distances.

I’ll selectively demonstrate these graphical diagnostics by applying an influence plot, added-variable plots, and component-plus-residual plots to Duncan’s regression (RegModel.1), making it the active model. I invite the reader to explore the other diagnostics as well.

Selecting Models > Graphs > Influence plot produces the dialog box at the top of Figure 7.26. I’ve left all of the selections in the dialog at their default values, including automatic identification of unusual points:32 Two cases are selected from each of the most extreme studentized residuals, hat-values, and Cook’s distances, potentially identifying up to six points (although this is unlikely to happen because influential points combine high leverage with a large residual). The resulting graph appears at the bottom of the figure, and, as it turns out, four relatively unusual points are identified: The occupation RR.engineer (railroad engineer) is at a high leverage point but has a small studentized residual; reporter has a relatively large (negative) studentized residual but small leverage; (railroad) conductor and, particularly, minister have comparatively large studentized residuals and moderately high leverage. The areas of the circles are proportional to Cook’s influence measure, so minister, combining a large residual with fairly high leverage, has the most influence on the regression coefficients.

FIGURE 7.25: Breusch-Pagan Test dialog and resulting output for Duncan’s occupational prestige regression (prestige ∼ education + income).

FIGURE 7.26: Influence Plot dialog and resulting graph for Duncan’s occupational prestige regression (prestige ∼ education + income).

Models > Graphs > Added-variable plots leads to the dialog box at the top of Figure 7.27; as before, the dialog allows the user to select a method for identifying noteworthy points in the resulting graphs. Once again, I retain the default automatic point identification but now increase the number of points to be identified in each graph from the default two to three.33 Clicking OK produces the graphs at the bottom of Figure 7.27. The slope of the least-squares line in each added-variable plot is the coefficient for the corresponding explanatory variable in the multiple regression, and the plot shows how the cases influence the coefficient—in effect, transforming the multiple regression into a series of simple regressions, each with other explanatory variables controlled.

FIGURE 7.28: Component+Residual Plots dialog and resulting graphs for Duncan’s occupational prestige regression (prestige ∼ education + income).

Finally, in the Models menu, Add observation statistics to data allows you to add fitted values (i.e., y), residuals, studentized residuals, hat-values, Cook’s distances, and observations indices (1, 2,…, n) to the active data set. These quantities, which (with the exception of observation indices) are named for the model to which they belong (e.g., residuals.RegModel.1), may then be used, for example, to create customized diagnostic graphs, such as an index plot of Cook’s distances versus observation indices.

7.9 Model Selection*

The R Commander Models menu includes modest facilities for comparing regression models and for automatic model selection. The menu items Models > Akaike Information Criterion (AIC) and Models > Bayesian Information Criterion (BIC) print the AIC or BIC model selection statistics for the current statistical model. Models > Stepwise model selectionperforms stepwise regression for a linear or generalized linear model, while Models > Subset model selection performs all-subsets regression for a linear model. Although I’ve never been terribly enthusiastic about automatic model selection methods, I believe that these methods do have a legitimate role, if used carefully, primarily in pure prediction applications.

I’ll illustrate model selection with the Ericksen data set in the car package. The data, described by Ericksen et al. (1989), concern the 1980 U. S. Census undercount, and pertain to 16 large cities in the United States, the remaining parts of the states to which these cities belong (e.g., New York State outside of New York City), and the other states. There are, therefore, 66 cases in all. In addition to the estimated percentage undercount in each area, the data set contains a variety of characteristics of the areas. A linear least-squares regression of undercount on the other variables in the data set reveals that some of the predictors are substantially collinear. 35  I fit this initial linear model with the formula undercount ∼ . and obtained variance-inflation factors for the regression coefficients via Models > Numerical diagnostics > Variance-inflation factors. The relevant output is in  Figure 7.29 .

Choosing Models > Subset model selection produces the dialog box at the top of  Figure 7.30 . All of the options in this dialog remain at their defaults, including using the BIC for model selection. Pressing the OK button produces the graph at the bottom of  Figure 7.30 , plotting the “best” model of each size k = 1,…,9, according to the BIC. The predictors included in each model are represented by filled-in squares, and smaller (i.e., larger-in-magnitude negative) values of the BIC represent “better” models. Notice that the regression intercept is included in all of the models. The best model overall, according to the BIC, includes the four predictors minority, crime, language, and conventional. I invite the reader to try stepwise model selection as an alternative. 36

FIGURE 7.29: Regression output and variance-inflation factors for Ericksen et al.’s Census undercount data, fitting the linear model undercount ∼ ..

FIGURE 7.29: Regression output and variance-inflation factors for Ericksen et al.’s Census undercount data, fitting the linear model undercount ∼ ..

FIGURE 7.30: Subset Model Selection dialog and resulting graph for Ericksen et al.’s Census undercount data.

1As explained in Section 7.2, Equation 7.1 also serves for more general linear models, where (some of) the xs may not be numeric explanatory variables but rather dummy regressors representing factors, interaction regressors, polynomial regressors, and so on.

2Recall that you must use the double equals sign (==) to test for equality; see Table 4.4 (page 71).

3A vector is a one-dimensional array, here of numbers.

4The sequence operator : creates an integer (whole-number) sequence, so 1:25 generates the integers 1 through 25. The c function combines its arguments into a vector, so -c(6, 16) creates a two-element vector containing the numbers –6 and –16.

5See Section 4.4.2, and in particular Table 4.4 (page 71), for information on R expressions.

6The arguments of an R function are the values given in parentheses when the function is called; if there is more than one argument, they are separated by commas.

7As formulated by Nelder (1977), the principle of marginality is deeper and more general than this characterization, but thinking of the principle in these simplified terms will do for our purposes.

8See later in this section for an explanation of how factors are handled in linear-model formulas.

9See Table 4.2 (on page 61) for an explanation of the variables in the Prestige data set.

10The second toolbar may be used to enter regression-spline and polynomial terms into the model. I’ll describe this feature in Section 7.3. If you’re unfamiliar with regression splines or polynomials, simply ignore this toolbar.

11The Linear Model dialog also includes a box for subsetting the data, and a drop-down variable list for selecting a weight variable for weighted-least-squares (as opposed to ordinary-least-squares) regression. Neither subsetting nor a weight variable is used in this example.

12The “T” in the names of the dummy-variable coefficients refers to “treatment contrasts”—a synonym for 0/1 dummy regressors—discussed further in Section 7.2.4.

13That the coefficient of women in the additive model is small doesn’t imply that women and type don’t interact. My true motive here is to simplify the example.

14Recall that a number like -3.354e+01 is expressed in scientific notation, and is written conventionally as –3.354 × 101 = –33.54.

15The function contr.Treatment is a modified version of the standard R function contr.treatment; contr.Treatment generates slightly easier to read names for the dummy variables—for example, type[T.wc] rather than typewc. Similarly, contr.Sum and contr.Helmert, discussed below, are modifications of the standard R functions contr.sum and contr.helmert.

16Multi-way analysis of variance in the R Commander, discussed in Section 6.1.3, uses contr.Sum for the factors in the ANOVA model.

17Reader: Can you see how the coefficients for type are related to each other across the two parametrizations of the model?

18The regressors of an orthogonal polynomial are uncorrelated, while those of a raw polynomial are just powers of the variable—for example, women and women2. The fit of raw and orthogonal polynomials to the data is identical: They are just alternative parametrizations of the same regression. Raw polynomials may be preferred for simplicity of interpretation of the individual regression coefficients, but orthogonal polynomials tend to produce more numerically stable computations.

19Also see the discussion of partial residuals in effect plots in Section 7.6.

20To make this example a little simpler, I’ve omitted occupational type from the regression.

21In this revised model, however, where the partial relationships of prestige to education and income are modeled more adequately, the quadratic coefficient for women is statistically significant, with p = 0.01.

22I’m grateful to Caroline Davis of York University for making the data available.

23You’re not constrained to select focal explanatory variables that correspond to high-order terms in the model, or even to terms in the model. For example, for Cowles and Davis’s logistic regression, you could select all three explanatory variables, sex, neuroticism, and extraversion, even though the three-way interaction among these variables isn’t in the model. In this case, the effect plot would be drawn for combinations of the values of the three explanatory variables.

24This strategy is used in general for effect plots of generalized linear models in the R Commander: The vertical axis is drawn on the scale of the linear predictor—the scale on which the model is linear—but labelled on the generally more interpretable scale of the response.

25Smoothing scatterplots is discussed in Section 5.4.

26Recall, however, that the explanatory variable women isn’t included in this model.

27For a generalized linear model, the Confidence Intervals dialog provides an option to base confidence intervals on either the likelihood ratio statistic or on the Wald statistic. The former requires more computation, but it is the default because confidence intervals based on the likelihood ratio statistic tend to be more accurate.

28For Type III tests to address sensible hypotheses, the contrasts used for factors in an ANOVA model must be orthogonal in the basis of the design. The functions contr.Sum, contr.Helmert, and contr.poly produce contrasts with this property, but the R Commander default dummy-coded contr.Treatment does not. For this reason, the R Commander Multi-Way Analysis of Variance dialog (Section 6.1.3) uses contr.Sum to fit ANOVA models (see, in particular, Figures 6.9 and 6.10, pages 118–119), and so the resulting models can legitimately employ Type III tests.

29F tests are computed for linear models and for generalized linear models (such as quasi-Poisson models) for which there’s an estimated dispersion parameter; chi-square tests are computed for generalized linear models (such as binomial models) for which the dispersion parameter is fixed. The same is true of analysis of variance and analysis of deviance tables.

30Recall that the identity (or inhibit) function I is needed here so that + is interpreted as addition.

31In Section 7.6, I showed how to add partial residuals to effect displays. For an additive model, that approach produces traditional component-plus-residual plots, but it is more flexible in that it can be applied as well to models with interactions.

32The default Automatic point identification has the advantage of working in the R Markdown document produced by the R Commander. As explained in Section 5.4.4, graphs that require direct interaction are not included in the R Markdown document.

33If you want to experiment with automatic identification of differing numbers of points, press the Apply button rather than OK.

34I’ll leave it as an exercise for the reader to remove the occupations minister (case 6) and conductor (case 16) from the data and refit the regression—most conveniently via the Subset expression box in the Linear Model or Linear Regression dialog; you can use the subset -c(6, 16).

35Ericksen et al. (1989) performed a more sophisticated weighted-least-squares regression.

36All-subsets regression in the R Commander is performed by the regsubsets function in the leaps package (Lumley and Miller, 2009), while stepwise regression is performed by the stepAIC function in the MASS package (Venables and Ripley, 2002). Although the distinction is not relevant to this example, where the full model is additive and all terms have one degree of freedom, stepAIC respects the structure of the model—for example, keeping dummy variables for a factor together and considering only models that obey the principle of marginality—while regsubsets does not.