A General Framework for Comparing Predictions and Marginal Effects across Models

Abstract

Many research questions involve comparing predictions or effects across multiple models. For example, it may be of interest whether an independent variable’s effect changes after adding variables to a model. Or, it could be important to compare a variable’s effect on different outcomes or across different types of models. When doing this, marginal effects are a useful method for quantifying effects because they are in the natural metric of the dependent variable and they avoid identification problems when comparing regression coefficients across logit and probit models. Despite advances that make it possible to compute marginal effects for almost any model, there is no general method for comparing these effects across models. In this article, the authors provide a general framework for comparing predictions and marginal effects across models using seemingly unrelated estimation to combine estimates from multiple models, which allows tests of the equality of predictions and effects across models. The authors illustrate their method to compare nested models, to compare effects on different dependent or independent variables, to compare results from different samples or groups within one sample, and to assess results from different types of models.

Keywords

marginal effects predictions cross-model comparisons nonlinear models seemingly unrelated estimation

1. Introduction

Many research questions require comparing predictions or effects across models. For example, one might compare whether the predicted number of poor mental health days for a 70-year-old differ from the number of poor physical health days for a similar person. Or, one could test whether the probability of tenure for a female scientist with a given set of characteristics is the same as that for a male scientist with the same characteristics. In terms of effects, a researcher might expect that the relationship between two variables can be explained by the addition of another variable. In this case, a model is fit to determine the initial relationship, followed by a second model that accounts for additional covariates. As another example, consider the situation in which a reviewer asks whether a different operationalization of the focal independent variable would change the effect of interest. In this case, two models are fit in which one uses the original operationalization and the other uses the alternative. As a last example, a researcher may want to compare an independent variable’s effect across different samples. For each of these examples, multiple models are fit, and tests of the equality of predictions or effects across models are used to answer the substantive question.

Fitting multiple models is straightforward, but testing the equality of effects across models is more challenging. Often, researchers informally assess cross-model differences by examining only the statistical significance of an effect within each model; this approach has been criticized because it is necessary to directly test cross-model differences to determine their statistical significance (Gelman and Stern 2006; Mustillo, Lizardo, and McVeigh 2018). Much methodological work has considered this topic and provides solutions for particular applications. For example, multiple solutions have been proposed to test for differences in effects across groups for linear (Chow 1960) and nonlinear (Allison 1999; Breen, Holm, and Karlson 2014; Long and Mustillo forthcoming; Williams 2009) models, for mediation in linear and nonlinear models (Karlson, Holm, and Breen 2012; MacKinnon et al. 2002), and for testing the ordinality assumption in ordinal logit (Brant 1990; Wolfe and Gould 1998). However, these solutions are sometimes limited in their applicability. A given solution might only work with mediation, be limited to tests of regression coefficients for the log-odds rather than marginal effects on the probability, or be unable to handle applications in which multiple regression parameters are needed to compute the effect, such as when a model includes both age and age squared.

We propose a general and flexible framework for comparing predictions and marginal effects across models.¹ Our method uses seemingly unrelated estimation (SUEST) to combine estimates from multiple models, which allows cross-model tests of predictions and marginal effects (Weesie 1999). This approach can be used for almost any linear or nonlinear regression model (e.g., binary, ordinal, nominal, count) and allows us to quantify predictions and effects in the natural metric of the dependent variable (e.g., probabilities instead of log-odds). By combining a wide range of applications with the flexibility of SUEST, we subsume a wide range of cross-model comparisons under a single framework. Our goal is not to critique existing methods for specific applications but rather to provide a general and flexible approach that can be applied to many statistical problems.

We begin by establishing our notation for models and predictions, followed by a review of marginal effects. We then show how SUEST can be used to combine estimates from multiple models, making it possible to test the equality of predictions and effects across models. We then discuss the implications of various methods of testing effects across models. The second half of the article applies our framework to applications that are commonly encountered in applied data analysis using models for binary, continuous, count, ordinal, and nominal dependent variables. We end with a discussion of software for implementing our methods.

2. Models and Predictions

Predictions are often of fundamental interest for understanding and interpreting the results of a regression model. Consider the regression model

η (x) = G (x β),

where $x$ is a vector of covariates, $β$ is a vector of regression coefficients, and $η$ is the predicted outcome. $G (x β)$ shows that the prediction is a function of $x β$ . For example, for the linear regression model $η = \hat{y} = x β$ , for the binary logit model we can express the predictions in the predicted probability metric with the formula $η = \Pr (y = 1) = \frac{\exp (x β)}{1 + \exp (x β)}$ . We use models for binary and continuous outcomes to explain our methods, but our results can be extended easily to models for nominal, ordinal, and count outcomes (see Agresti 2013; Long and Freese 2014), as illustrated by our examples at the end of this article.

Although most of our article focuses on comparing marginal effects, cross-model comparisons of predictions may also be of interest. For example, after fitting separate models for men and women, we might compare a prediction for a 60-year-old man with a prediction for a 60-year-old woman. Predictions can also be made over the range of a variable and then compared across models, as we illustrate in example 6.1. In either case, when making predictions, important decisions must be made regarding the specific values of the independent variables. To make these values explicit, we use the more precise notation $η (x_{k} = x_{k}^{*}, x_{- k} = x^{*})$ to highlight the specific values at which a prediction is made. The focal independent variable $x_{k}$ is set to the value $x_{k}^{*}$ . For example, $x_{k} = {\bar{x}}_{k}$ or $x_{k} = 1$ . The notation $x_{- k}$ is the vector created from $x$ by removing the focal variable $x_{k}$ . This vector of control variables is set equal to the values in $x^{*}$ . We discuss decisions about what values to use for $x_{k}^{*}$ and $x^{*}$ in more detail below.

3. Marginal Effects

Marginal effects summarize how changes in a focal independent variable affect the predicted value of the outcome, holding other variables at specific values. The marginal effect of $x_{k}$ is the change in $η$ for a change in $x_{k}$ holding the variables in $x_{- k}$ at the specific values $x^{*}$ . There are two varieties of marginal effects (for further details, see Long and Freese 2014:270–80). A discrete change is the change in $η$ for a discrete or finite change in $x_{k}$ . An instantaneous change is the change in $η$ for an infinitely small change in $x_{k}$ (i.e., the first derivative at a specified point). In linear regression, marginal effects can be identical to regression coefficients. However, when nonlinearities are introduced into the model, say by adding powers of independent variables, marginal effects no longer equal the regression coefficients. In models such as logit and probit, the regression coefficients are not equal to marginal effects on the probability of the outcome. In this article, we primarily focus on the discrete change, but our methods can be used with an instantaneous rate of change as well. The critical idea is that a focal independent variable is changing while other covariates are not.

We begin by reviewing two types of marginal effects: (1) marginal effects at representative values (MER), in which covariates are held at theoretically interesting or representative values, and (2) average marginal effects (AMEs) that average the marginal effects computed at the observed values of the covariates for each observation. After defining these effects, we explain how to test if two or more effects are equal.

3.1. Discrete Changes at Representative Values

The discrete change with respect to $x_{k}$ at representative values (DCR) is the change in $η$ as $x_{k}$ changes from start to end while holding other variables at specified values:

DC R_{x_{k}} = η (x_{k} = e n d, x_{- k} = x^{*}) - η (x_{k} = s t a r t, x_{- k} = x^{*}) .

We use start and end to indicate the specific values of interest for $x_{k}$ with other covariates held at the values in $x^{*}$ . Commonly, the variables in $x_{- k}$ are held at their means (referred to as the discrete change at the mean), so the discrete change can be interpreted as the effect of a change in $x_{k}$ for an “average person” in the sample. However, other values could be used, such as means for specific groups of interest or values that define substantively interesting observations. For a continuous focal variable $x_{k}$ , the start and end values can be any values of interest. For example, $x_{k}$ could increase from a start value of the mean to an end value of the the mean plus one standard deviation $s_{k}$ :

DC R_{x_{k} + s_{k}} = η (x_{k} = {\bar{x}}_{k} + s_{k}, x_{- k} = x^{*}) - η (x_{k} = {\bar{x}}_{k}, x_{- k} = x^{*}) .

Increases of one unit or any other value could also be used. If the model includes polynomials or interactions, all associated variables must change in tandem. For example, if $x_{agesq} = x_{age} \times x_{age}$ , then $x_{agesq}$ must change from $100$ to $121$ when $x_{age}$ changes from $10$ to $11$ . If $x_{k}$ is binary, start is 0 and end is 1:

DC R_{binary x_{k}} = η (x_{k} = 1, x_{- k} = x^{*}) - η (x_{k} = 0, x_{- k} = x^{*}) .

(1)

3.2. Average Discrete Changes

The discrete change for $x_{k}$ could also be computed at the observed values for observation $i$ :

D C_{x_{ik}} = η (x_{k} = e n d_{i}, x_{- k} = x_{- k, i}) - η (x_{k} = star t_{i}, x_{- k} = x_{- k, i}) .

The subscript i for $e n d_{i}$ and $star t_{i}$ indicates that these values can vary by observation. $x_{- k} = x_{- k, i}$ indicates the control variables are equal to their observed values. As with the DCR, the discrete change for the $i^{th}$ observation can be computed for different amounts of change in $x_{k}$ . For example, for a continuous focal variable we can compute the effect when $x_{k}$ increases by $δ$ from its observed value $x_{ik}$ :

D C_{x_{ik} + δ} = η (x_{k} = x_{ik} + δ, x_{- k} = x_{- k, i}) - η (x_{k} = x_{ik}, x_{- k} = x_{- k, i}) .

One or a standard deviation are common values of $δ$ , but other values (e.g., 10) can be used. When $x_{k}$ is binary, the only reasonable change is from 0 to 1:

D C_{binary x_{ik}} = η (x_{k} = 1, x_{- k} = x_{- k, i}) - η (x_{k} = 0, x_{- k} = x_{- k, i}) .

It is also possible to change $x_{k}$ between the same two values for every observation, such as increasing income from 10,000 to 20,000 while holding the variables in $x_{- k}$ at their observed values.

The average discrete change (ADC) for $x_{k}$ is the average of the discrete changes for each observation:

AD C_{x_{k}} = \frac{1}{N} \sum_{i}^{N} η (x_{k} = {e n d}_{i}, x_{- k} = x_{- k, i}) - η (x_{k} = {s t a r t}_{i}, x_{- k} = x_{- k, i}) .

(2)

Often, it is of interest to average the effect across the entire sample of $N$ observations, as shown in equation (2). In this case, the ADC is the average effect in the sample. In some situations, it is more useful to compute the average over a subsample. For example, one may wish to estimate the average effect separately for different groups, such as for men and for women. In this case, only the relevant observations (e.g., those for women) are used to calculate the ADC for that group.

3.3. Risk Ratios

Discrete changes quantify effects in terms of the difference between predictions at ending and starting values. For example, equation (1) shows the calculation of a DCR for a binary independent variable. For models interpreted in terms of predicted probabilities (e.g., logit and probit), the ratio of these two predictions can also be of interest, which is referred to as the risk ratio (RR; for details, see Norton, Miller, and Kleinman 2013). The RR when $x_{k}$ changes from $0$ to $1$ (e.g., the prediction for men compared with the prediction for women) is defined as follows:

R R_{binary x_{k}} = \frac{η (x_{k} = 1, x_{- k} = x^{*})}{η (x_{k} = 0, x_{- k} = x^{*})} .

As with discrete changes, the values of the covariates in $x^{*}$ can held at whatever values are of interest. One could also make predictions at the observed values for each observation and then average the ratios across the sample, as proposed by Norton et al. (2013). Although we focus here on discrete changes, the methods we develop can be used with RRs or other measures defined as functions of predictions.

3.4. Marginal Effects for Linear and Nonlinear Relationships

Marginal effects can be estimated for almost any regression model, although how they are interpreted depends on whether the relationship between a predictor and the outcome is approximately linear over the region of interest. For example, if the model is $l n (y) = α + β x + ϵ$ , then the relationship between $x$ and $l n (y)$ is linear, but the relationship between $x$ and $y$ is nonlinear. If a model is $y = α + β x + ϵ$ , the relationship between $x$ and $y$ is linear throughout the range of $x$ . If one adds $x^{2}$ to the model, $y = α + β x + γ x^{2} + ϵ$ , the model is still linear in the sense that the right-hand side is a linear combination of the regressors, but there is a nonlinear relationship between $x$ and $y$ . Over a limited range of $x$ , however, the relationship could be approximately linear.

Whether the relationship between a predictor and the outcome is linear also depends on how the outcome is modeled. For example, the logit model is typically thought of as a nonlinear model in the sense that the relationship between predictors and the probability $π$ of the outcome is nonlinear. But, if the model is written in terms of the logit as $l n [π / (1 - π)] = α + β x$ , the model is linear in the logit, although the relationship between $x$ and $π$ is still nonlinear. However, if the observed data occur in a region where $π$ is between $. 4$ and $. 6$ , then the relationship between $x$ and $π$ will be essentially linear (Long 1997). When computing effects and comparing effects across models, it is important to consider whether the relationship being modeled is linear or nonlinear over the data space of substantive interest. This point is illustrated in Section 6.

4. Comparing Models Using SUEST

To test whether predictions or marginal effects from multiple models are equal, the models often need to be fit simultaneously to compute cross-model covariances. This section explains why this is necessary and introduces SUEST as a method for estimating cross-model covariances.

4.1. Testing the Equality of Predictions and Marginal Effects

We begin by reviewing tests for predictions or effects from a single model. We then consider the case of comparing predictions or effects across models, which is more complicated because covariances among the estimates across models are not available if the models are fit separately. In Section 4.2 we explain how SUEST makes it possible to compute these covariances and why it is important to have estimates of these covariances.

Suppose we are predicting individuals’ overall life satisfaction on the basis of their age, accounting for whether they have a college degree (0 or 1) and whether they earn more than $50,000 per year (0 or 1). The dependent variable is binary, with 1 = high satisfaction and 0 = low satisfaction. We fit a regression model for life satisfaction that includes age, education, and income as independent variables. Let $π_{1} = π (age = 60, educ = 0, inc = 0)$ be the probability of high life satisfaction for someone who is 60 years old, did not go to college, and has lower income. Similarly, $π_{2} = π (age = 70, educ = 0, inc = 0)$ is a prediction for a similar person who is 70 years old. The Wald statistic to test $H_{0} : π_{1} = π_{2}$ is

z = \frac{{\hat{π}}_{1} - {\hat{π}}_{2}}{\sqrt{{\hat{σ}}_{1}^{2} + {\hat{σ}}_{2}^{2} - 2 {\hat{σ}}_{1, 2}}},

(3)

where ${\hat{σ}}_{k}$ is the standard error of ${\hat{π}}_{k}$ and ${\hat{σ}}_{1, 2}$ is the covariance between the estimates ${\hat{π}}_{1}$ and ${\hat{π}}_{2}$ . The numerator is the estimated difference in probabilities, and the denominator is the standard error of this difference. The null hypothesis is rejected if $z$ exceeds the critical value. This test can be generalized to test more complex hypotheses, such as the equality of more than two effects (Cameron and Trivedi 2005:135–39). To test if the effects of education and income are equal, let $D C_{educ}$ be the discrete change of education on satisfaction and $D C_{inc}$ be the discrete change of income. The Wald statistic to test $H_{0} : D C_{educ} = D C_{inc}$ is

z = \frac{{\hat{DC}}_{educ} - {\hat{DC}}_{inc}}{\sqrt{{\hat{σ}}_{educ}^{2} + {\hat{σ}}_{inc}^{2} - 2 {\hat{σ}}_{educ, inc}}} .

(4)

Calculation of the $z$ statistics in equations (3) and (4) requires two estimates, the standard error of each estimate, and the covariance between the estimates.

When testing the equality of regression coefficients from a single model, the standard errors and covariances are obtained from the covariance matrix of the estimates. When comparing estimates that are not regression coefficients, such as discrete changes for a logit model, the estimates of effects are computed from the regression coefficients, and the covariance matrix for these effects is obtained using the delta method (Agresti 2013:72–77; Bishop, Fienberg, and Holland 1975:486–97), bootstrapping (Efron and Tibshirani 1993), or simulation (King, Tomz, and Wittenberg 2000).²

The Wald tests in equations (3) and (4) can also be used to test predictions and marginal effects across models. The practical problem is that if the models are fit separately, there is no estimate of the covariance between the estimates being compared. As discussed in Section 5, in some applications we might know that these covariances are zero, but in most applications they need to be estimated. We now explain how this can be done using SUEST.

4.2. SUEST

SUEST combines covariance matrices across multiple models and computes the cross-model covariances often needed for tests of predictions and effects across models.³ SUEST can be implemented either by stacking data and then fitting models simultaneously, or by fitting the models separately and then combining the results. Weesie (1999) developed the method, which is based on results from Hausman (1978) and White (1982). Because SUEST combines covariance matrices of multiple models, both within- and cross-model variances and covariances of the regression coefficients are available. Predictions and marginal effects can then be estimated along with the necessary variances and covariances for testing the equality of estimates across models.

To demonstrate the intuition behind the SUEST approach, we adapt an example from Weesie (1999) that illustrates SUEST using the method of stacking. Suppose we want to determine whether the number of social roles a person holds has the same effect on having depression (depress) as on feeling one has bad health (sick). The data set contains the following information:

(\begin{matrix} depress & sick & roles & id \\ 0 & 0 & 5 & 1 \\ 1 & 1 & 1 & 2 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 1 & 8 & N \end{matrix}) (unstacked dataset)

where id is a unique identification number for each respondent. First, we regress depress on roles (model 1), we then regress sick on roles (model 2). After fitting these models separately, we cannot use the formula in equation (4) to test if the effect of roles is the same across the two models, because the covariance between the two estimates has not been computed. To estimate this covariance, we use the unstacked data set to create a stacked data set with $2 N$ observations:

(\begin{matrix} y & r o l e s m 1 & r o l e s m 2 & i s m o d e l 1 & i d \\ d e p r e s s = 0 & 5 & 0 & 1 & 1 \\ d e p r e s s = 1 & 1 & 0 & 1 & 2 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ d e p r e s s = 1 & 8 & 0 & 1 & N \\ s i c k = 0 & 0 & 5 & 0 & 1 \\ s i c k = 1 & 0 & 1 & 0 & 2 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ s i c k = 1 & 0 & 8 & 0 & N \end{matrix}) (stacked dataset) .

The first $N$ rows of y contain the variable depress from the unstacked data set; the last $N$ rows contain the variable sick. The first $N$ rows for rolesm1 are the same as roles in the unstacked dataset with the last $N$ rows fixed to 0, for reasons we will explain shortly. rolesm2 is 0 in the first $N$ rows and then reproduces roles. The first half of the data are used to fit model 1 for depress; the second half are used to fit model 2 for sick. For example, if we regress y on rolesm1 using the first N rows of the stacked data set, we obtain exactly the same results as using the unstacked data set.

Two other variables are included in the stacked data set that are necessary to jointly fit the two models, which allows us to estimate the needed covariances. Variable ismodel1 is 1 if the observation is needed to fit model 1 and 0 if the observation is used to fit model 2. Variable id keeps track of the clustering caused by having two observations for each participant. It shows, for example, that row 1 of the stacked data set is the same observation as row $N$ + 1.

To jointly fit the two models, we regress y on rolesm1, rolesm2, and ismodel1. Panel A of Table 1 presents estimates of the binary logit coefficients and their standard errors from three models. Model 1 is the logit for depress fit using the unstacked data set; model 2 is the logit for sick fit using the unstacked data set; and model 3 fits a logit for y regressed on rolesm1, molesm2, and ismodel1 using the stacked data set.⁴ Model 3 was fit using the cluster-robust variance estimator (Liang and Zeger 1986; White 1984) to adjust for the duplication of each observation in the stacked data set (details below). The coefficient for roles in model 1 is identical to the coefficient for rolesm1 in model 3. The $N$ observations for which rolesm1 was fixed to 0 had no effect on the coefficient estimates. The standard errors differ slightly because of the adjustment for clustering that is applied when the models are fit jointly; the standard errors in model 3 are identical to those that would be obtained if models 1 and 2 were fit using the robust sandwich estimator of the variance (described in the next section). Because model 3 fits models 1 and 2 jointly, we obtain ${\hat{σ}}_{{\hat{β}}_{rolesm 1}, {\hat{β}}_{rolesm 2}}$ , which is the covariance between ${\hat{β}}_{roles}$ from model 1 and ${\hat{β}}_{roles}$ from model 2.

Table 1.

Regression Estimates and Average Discrete Changes from Logit Models 1 and 2 Fit Separately and Model 3 Fit Jointly with Stacked Data

	Model 1	Model 2	Model 3
	depress	sick	Stacked
Panel A: Regression coefficients
roles	−.201*** (.025)	−.261*** (.037)
rolesm1			−.201*** (.026)
rolesm2			−.261*** (.037)
ismodel1			.811*** (.231)
Constant	−.040(.148)	−.852*** (.206)	−.040(.150)
Panel B: Average discrete changes
roles	−.035*** (.004)	−.021*** (.003)
rolesm1			−.035*** (.004)
rolesm2			−.021*** (.002)
$N$	4,307	4,307	8,614

Note: Standard errors are in parentheses.

***

p < .001 (two-tailed test).

It would be tempting to test the equality of the effect of roles across the two models by comparing the estimated logit coefficients. However, the regression coefficients are not appropriate for comparing effects across logit and probit models because a change in the size of the coefficient across models can reflect both confounding and rescaling of the model (Karlson et al. 2012). Marginal effects do not have this limitation and are thus appropriate for comparing effects across logit and probit models (Breen, Karlson, and Holm 2018). Panel B of Table 1 presents AMEs for an instantaneous change in roles. As with the coefficient estimates, the AME sizes are the same in the individual models 1 and 2 as in the stacked data set model 3. Using equation (4), the test statistic for the hypothesis that the AMEs for social roles are equal for the two outcomes is $z = 2.99$ , which has a two-tailed p value of .003. Therefore, we reject the null hypothesis that the AMEs for social roles are equal across the two outcomes.

Although we use stacking to illustrate the logic of SUEST, in practice stacking data sets has several limitations. First, stacking only works when the models being jointly fit are of the same type (e.g., two binary logits; two linear regressions). Stacking does not work if we want to combine, for example, estimates from a logit model with those from a probit model. Second, stacking constrains ancillary parameters, such as the cutpoints in an ordinal logit, to be the the same across models. Next, we turn to a more formal method of implementing SUEST that does not have these limitations.

4.3. Obtaining Cross-Model Covariances with Equation-Level Scores

In most applications of SUEST, models are fit individually and the estimates are combined post hoc. This method of implementing SUEST, developed by Weesie (1999), builds on the work of Hausman (1978) and Clogg, Petkova, and Haritou (1995). Cross-model covariances are estimated using an application of the robust sandwich estimator of the variance (Arellano 1987; White 1982). For models $m = 1, . . ., n$ , the robust sandwich estimator is

Var ({\hat{β}}_{m}) = D_{m}^{- 1} \sum_{i} w_{mi} u_{mi} u_{mi}^{'} D_{m}^{- 1},

where $D_{m}$ is the Jacobian of the $m^{th}$ model, $w_{mi}$ are weights that can sometimes be set to zero if an observation is not in both models, and $u_{mi}$ are the equation-level scores from the individually estimated models. For linear regression, the equation-level score is the residual for each observation; with maximum likelihood estimation, the equation-level score is the first derivative of the log-likelihood. Using the sandwich estimator, the cross-model covariance for models $p$ and $q$ is the following:

Cov ({\hat{β}}_{p}, {\hat{β}}_{q}) = D_{p}^{- 1} \sum_{i} w_{i} u_{pi} u_{qi}^{'} D_{q}^{- 1} .

The estimates of the variance of individual parameters are equivalent to those obtained from separate models that apply the robust sandwich estimator. However, by combining the covariance matrices we are also able to obtain the covariances among the estimates across and within each model. This is accomplished by concatenating the scores at the level of each observation, which are treated as clusters. Recall that in the earlier stacked data example, we introduced clustering by repeating each observation twice. Here, too, if the models are estimated on the same observations, we handle the clustering and obtain estimates of cross-model covariances using the cluster-robust sandwich estimator.

In theory, estimates from almost any models can be combined using SUEST. In practice, current implementations of SUEST are limited to models that can produce equation-level scores (see Weesie 1999), although alternative implementations of the approach may allow for a broader range of models (Canette 2014). In our examples in Section 6, we demonstrate the SUEST approach for linear regression, binary logit, negative binomial regression, ordinal logit, and multinomial logit.

5. Comparing Tests of Cross-Model Differences

Section 4.1 presented Wald tests of the hypothesis that predictions or marginal effects are equal across models. More generally, we can test the equality of any function $τ$ of predictions (e.g., marginal effects, probabilities, rates, RRs). The Wald test of $H_{0} : τ_{1} = τ_{2}$ is

z = \frac{{\hat{τ}}_{1} - {\hat{τ}}_{2}}{\sqrt{{\hat{σ}}_{1}^{2} + {\hat{σ}}_{2}^{2} - 2 {\hat{σ}}_{1, 2}}},

(5)

where ${\hat{σ}}_{k}$ is the standard error of ${\hat{τ}}_{k}$ and ${\hat{σ}}_{1, 2}$ is the covariance between the estimates ${\hat{τ}}_{1}$ and ${\hat{τ}}_{2}$ . The denominator is the standard error of the difference ${\hat{τ}}_{1} - {\hat{τ}}_{2}$ , where ${\hat{σ}}_{1}^{2} + {\hat{σ}}_{2}^{2}$ quantifies the uncertainty of the individual estimates. The covariance ${\hat{σ}}_{1, 2}$ adds information about sampling covariation between ${\hat{τ}}_{1}$ and ${\hat{τ}}_{2}$ . For example, a positive covariance between ${\hat{τ}}_{1}$ and ${\hat{τ}}_{2}$ indicates that in repeated sampling when the estimate ${\hat{τ}}_{1}$ from the first model is large and positive, the estimate ${\hat{τ}}_{2}$ from the second model is also likely to be large and positive. A positive covariance results in a larger $z$ statistic because subtracting $2 {\hat{σ}}_{1, 2}$ decreases the variability introduced by ${\hat{σ}}_{1}^{2} + {\hat{σ}}_{2}^{2}$ . A negative covariance produces a smaller $z$ statistic, reflecting greater uncertainty in our estimates of the difference. As we will show, both positive and negative covariances are possible.

Some tests for cross-model differences assume that $σ_{1, 2} = 0$ in equation (5) (e.g., Paternoster et al. 1998). These tests use the statistic $z'$ defined as

z' = \frac{{\hat{τ}}_{1} - {\hat{τ}}_{2}}{\sqrt{{\hat{σ}}_{1}^{2} + {\hat{σ}}_{2}^{2}}},

(6)

which implicitly assumes $σ_{1, 2} = 0$ . In the next section we consider when it is appropriate to assume $σ_{1, 2} = 0$ . We then consider the effect of incorrectly assuming that the covariance is zero, which is sometimes done to simplify computation of the test.

5.1. When Do Cross-Model Covariances Need to Be Estimated?

As a general rule, cross-model covariances ${\hat{σ}}_{1, 2}$ need to be estimated when the same observations are used across models. For example, tests of mediation or tests of effects across nested models require an estimate of the cross-model covariances. Accordingly, $z'$ in equation (6) is an incorrect test that can lead to incorrect conclusions. In other cases, the covariance across models is known to be zero and $z'$ is appropriate. However, there is no downside in using $z$ in these cases as $z$ equals $z'$ . In contrast, when the cross-model covariances are nonzero, $z$ does not equal $z'$ and $z'$ will be incorrect.

Cross-model covariances can sometimes be zero when the models are fit using mutually exclusive observations. For example, with a simple random sample, the cross-model covariances are zero if model 1 includes only men and model 2 includes only women. Cross-model covariances might also be zero in meta-analyses that compare effects from separate samples. Keep in mind, however, that even when different observations are used across models, cross-model covariances can be nonzero. For example, consider data on middle school students nested within different schools. Even if separate models are fit for boys and girls, corrections for clustering, such as the cluster-robust variance estimator, produce nonzero cross-model covariances. Similarly, corrections for stratification in sampling produce nonzero covariances even with mutually exclusive observations across models (Long and Mustillo forthcoming). The critical point is that using mutually exclusive observations across models does not imply that the cross-model covariances are zero. Below, we demonstrate the various consequences that non-zero cross-model covariances can have on inferences for cross-model differences.

5.2. The Impact of Accounting for Cross-Model Covariances

To provide an intuitive understanding of the impact of the covariance between effects on tests of their difference, we simulated a data set of 2,000 observations for two examples. The first example compares the effects of $x_{1}$ across models for the binary dependent variables $y_{1}$ and $y_{2}$ . The second example compares the effect of $x_{1}$ on the binary dependent variable $y_{1}$ in a model with just $x_{1}$ and a second model that also includes a potential mediator $x_{2}$ . We use the bootstrap method to obtain 1,000 estimates of the AMEs from each model and the cross-model difference in effects. The bootstrap test statistic for the hypothesis that the difference is zero is $z_{bs}$ (Efron and Tibshirani 1993). We compare these results to those obtained with SUEST using both $z$ , which accounts for cross-model covariance, and $z'$ , which assumes the covariance is zero.

For the first example, Figure 1A shows the distribution of 1,000 bootstrapped estimates of the two effects. These estimates have a correlation of −0.056. Figure 1B shows the distribution of the estimates of the difference in effects, leading to $z_{bs} = 3.12$ for the test that the difference is zero. Using SUEST, the test statistic is $z' = 3.16$ when assuming the covariance is zero, compared with $z = 3.09$ when the covariance is included. The test statistics are similar because in the formula for the standard error of the difference in equation (5), the effect of the covariance ${\hat{σ}}_{1, 2} * 2 = - 0.000008$ is small compared with the variance of the estimates ${\hat{σ}}_{1}^{2} + {\hat{σ}}_{2}^{2} = 0.000193$ . $z$ properly incorporates information that we have slightly more uncertainty in the estimate of the difference than we would expect if we only examined the uncertainty of each individual estimate, but in this example the results are almost identical to those using $z'$ where we assume $σ_{1, 2} = 0$ .

Figure 1.

Bootstrapped estimates of average marginal effects (A) and cross-model difference (B): effect of $x_{1}$ across different dependent variables ( $N = 2, 000$ ).

For the second example, Figure 2A shows the distribution of bootstrapped estimates of the two AMEs, where the bootstrapped estimates now have a large correlation of 0.857. Intuitively, this correlation suggests that even though the values of the individual estimates vary across the bootstrapped samples, the difference between these estimates is relatively stable. Figure 2B shows the distribution of the estimates of the difference in effects. On the basis of these estimates, the test statistic is $z_{bs} = 4.68$ . Using SUEST, the z statistics computed with and without the constraint that the covariance is zero lead to different conclusions regarding the hypothesis that the effects are equal: $z' = 1.78$ versus $z = 4.69$ . These results differ because when computing the standard error of the difference, the effect of subtracting twice the covariance ${\hat{σ}}_{1, 2} * 2 = 0.000182$ from ${\hat{σ}}_{1}^{2} + {\hat{σ}}_{2}^{2} = 0.000213$ is substantial. As Figure 2B makes clear, we have considerably less uncertainty, shown by the tighter confidence intervals, about the difference in the estimates than would be suggested by looking only at the variability of the individual estimates.

Figure 2.

Bootstrapped estimates of average marginal effects (A) and cross-model difference (B): effect of $x_{1}$ across base and mediation models ( $N = 2, 000$ ).

The key point is that assuming that cross-model covariances are zero can lead to incorrect estimates of the standard error of the difference. Our framework provides a general method for obtaining cross-model covariances of the estimates, which are often critical for determining the correct test statistic when comparing estimates across models.

5.3. Statistical Significance and Substantive Importance

Thus far, we have discussed the necessary steps to determine the correct test for the equality of estimates across models, focusing on issues related to obtaining the correct estimate of the uncertainty in estimates of the difference (i.e., the denominator in equation 5). Note that no test of statistical significance can determine whether a cross-model difference is substantively important. Thus, we should also consider whether the difference in the estimates in the numerator of equation (5) is substantively interesting. Across different topics, applications, and literatures, various effect sizes may or may not be substantively meaningful.

We believe our focus on predictions and marginal effects has important benefits for determining the substantive importance of effect sizes, as they can be expressed in terms of the natural metric of a dependent variable (e.g., probabilities instead of log-odds). Still, a substantive decision must be made. In our second example above, the marginal effect in the base model without the mediator is 0.083, which reduces to 0.057 after adding the mediating variable $x_{2}$ . The cross-model difference of 0.026 is statistically significant ( $p < . 001$ ), but determining the substantive importance of a 0.026 reduction in the effect is the purview of researchers engaged in the relevant substantive field.⁵ Although the result is statistically significant, no statistical test, including those we suggest, can determine the substantive importance of an effect. However, expressing effects in the natural metric of the dependent variable aids researchers’ ability to make decisions about substantive importance.

6. Examples and Applications

In this section, we present a series of examples that illustrate how our framework can be used to compare predictions and marginal effects across models in a variety of common situations. Our first two examples compare effects across nested models. This situation often arises in the context of mediation, where researchers informally assess the impact of a mediating variable by examining whether the effect of interest changes after accounting for other variables. This approach has been criticized for not providing a formal statistical test (Gelman and Stern 2006; Mustillo et al. 2018), something our method addresses.⁶

Example 6.1 uses linear regression to compare predictions and effects when nonlinearities are introduced by adding polynomial terms. Example 6.2 is a more traditional mediation analysis that fits a series of nested logit models and examines how the marginal effect of a variable changes as controls are added to the model. The next two examples show how to compare predictions and marginal effects for different variables. Example 6.3 compares different operationalizations of a predictor, and example 6.4 compares effects of the same predictor on different outcome variables. These examples compare effects computed from the same type of model. Example 6.3 uses two binary logits, and example 6.4 uses negative binomial regression for two count outcomes. Example 6.5 shows how SUEST can combine estimates from different types of models, comparing predictions and marginal effects from a multinomial logit to those from ordinal logit. Finally, example 6.6 illustrates comparing effects across different samples. Our framework is not limited to these examples; it can be applied generally to applications involving comparison of results from multiple models.

Analyses were computed in Stata 15.1 using SPost13 (Long and Freese 2014). Section 7 discusses software issues. The data sets and script files to reproduce our examples are available at https://www.trentonmize.com/software/mecompare.

6.1. Predictions and Marginal Effects to Summarize Curvilinear Relationships

Our first example examines a nonlinear relationship using linear regression. In this context, nonlinearities most often arise when powers of continuous independent variables are included in the model, such as income and income squared. Comparing the effects of such variables is complicated by two issues. First, the magnitude of the effect depends on the value of the independent variable where it is evaluated. Second, the effect depends on multiple regression coefficients.

Consider the relationship between income and depressive symptoms. Earning more should protect against depression, but we expect the relationship will be nonlinear because at some point increases in income will likely produce diminishing benefits. Data for this example come from the National Longitudinal Study of Adolescent to Adult Health Wave IV (Harris 2009). The dependent variable is a scale of depressive symptoms that we treat as continuous. The focal independent variable is income measured in thousands of dollars. Model 1 includes income and income squared as predictors to allow a nonlinear relationship between income and depression. Age, gender, and race are included as control variables.

The lefthand panel of Figure 3 plots the predictions from model 1, showing a nonlinear relationship between income and depressive symptoms. As shown in Table 2, the coefficient for income squared is statistically significant at the $0.001$ level. Even though the relationship is nonlinear, it is often useful to summarize the effect. To do this, we calculate the ADC for a standard deviation increase in income. We find that, on average, a standard deviation increase in income, about $27,000, is associated with 0.816 fewer depressive symptoms ( $p < . 001$ ).

Figure 3.

Predicted depressive symptoms on the basis of personal income in linear regression models including income and income squared ( $N = 4, 307$ ).

Table 2.

Effect of Income on Depressive Symptoms before and after Controlling for Job Satisfaction (n = 4,307)

	Model 1	Model 2	Cross-Model
	Base Model	Mediation Model	Difference
Regression coefficients^a
Income	–.515***(.065)	–.432***(.063)
Income squared	.002***(.001)	.002***(.001)
Average discrete change^b
Income (+ SD)	–.816***(.076)	–.648***(.074)	.168***(.022)

Note: All models control for age, gender, and race. Standard errors are in parentheses.

The regression coefficients are rescaled to represent a $10,000 increase in income.

Average discrete change is for a standard deviation increase in income.

***

p < .001 (two-tailed test).

It is unlikely that income has only a direct effect on depression. Indirect effects are likely underlying the relationship shown in the lefthand panel of Figure 3. One possibility is that people with more income work in more satisfying jobs. To test whether the effect of income changes after accounting for this factor, model 2 adds a measure of job satisfaction (see the righthand panel of Figure 3). Although not a dramatic change, the relationship between income and depression is attenuated from the relationship in model 1.

Figure 4 plots the differences in predictions across the models. At about $50,000 and above, model 2 predicts higher depressive symptoms than model 1. For example, our conclusions about the predicted number of depressive symptoms for someone who earns $100,000 would differ across the two models. Figure 4 shows differences in predictions at specific income values, but the ADC is helpful to summarize differences across the two models. In model 1, the ADC equals $- 0.816$ , which reduces to $- 0.648$ in model 2 because of accounting for job satisfaction (column 2 of Table 2). A test of the difference in ADCs across models (column 3 of Table 2) shows that accounting for job satisfaction decreases the effect of income by $0.168$ , a significant reduction in effect size ( $z = 7.590, p < . 001$ ). Note the importance of accounting for the cross-model covariance in the calculation of $z$ ; for this example, $z' = 1.606$ , which highlights the misleading conclusions that are possible when the covariance is incorrectly assumed to be zero.

Figure 4.

Difference in predicted depressive symptoms across models ( $N = 4, 307$ ).

6.2. Marginal Effects across Nested Logit Models

This example compares marginal effects across a series of nested logit models. Consider that people with college degrees tend to report happier lives than those without college degrees. A college degree confers myriad possible benefits that might partially explain this relationship. Theoretically, we suspect that receiving a college degree first and foremost leads to better and higher paying jobs, which has a direct effect on life happiness. Our analyses examine whether the effect of a college degree changes after accounting for earnings, occupational prestige, and other variables. In linear regression, this question could be answered by fitting a series of nested models and determining whether the regression coefficient for college is significantly attenuated. In logit models, however, we cannot compare regression coefficients because they are only identified up to a scale factor (Amemiya 1981:1489; Karlson et al. 2012; McKelvey and Zavoina 1975). In contrast, marginal effects do not have this limitation and can be used to compare effects across models.

Data for this example come from the General Social Survey (GSS; Smith et al. 2016), pooling data from 2000 to 2016 and including all employed individuals. Happiness is a binary variable equal to 1 if a respondent reports being “very happy” and 0 if not. Education is measured as a binary variable indicating receipt of a college degree. Control variables include age, age squared, marital status, parental status, gender, political views, religious affiliation, and year in which the survey was answered. The two variables we expect may explain some of the relationship between college degree status and happiness are hourly wages and occupational prestige, both of which are continuous.

Panel A of Table 3 presents the ADC of a college degree on happiness for a series of four nested models. Results in panel B show tests that the discrete changes are equal across models. The only independent variable in model 1 is having a college degree. The ADC indicates that, on average, the probability of being happy is 0.072 higher for individuals with a college degree ( $p < . 001$ ). Model 2 adds demographic control variables, reducing the ADC for college to 0.060, which remains significant. A direct test of the difference in the ADC from model 1 to model 2 shows that adding the controls significantly decreases the effect of college by 0.012 ( $p < . 01$ ; see Panel B of Table 3). Model 3 considers whether the effect of college can be further explained by the association between education and hourly wages. After accounting for wages, the effect of college continues to be significant (ADC $= . 036; p < . 01$ ) but is again reduced, here by 0.024 ( $p < . 001$ ). Finally, model 4 adds occupational prestige, decreasing the ADC for college to a nonsignificant 0.019, a reduction from model 3 of 0.017 (cross-model difference $p < . 001$ ). Thus, the results show significant reductions in the effect of college with each addition of independent variables across models.

Table 3.

Effect of Having a College Education on Happiness Using Average Discrete Changes from Binary Logit Model (n = 9,216)

	Model 1	Model 2	Model 3	Model 4
	College Only	+ Controls	+ Wages	+ Prestige
Panel A: Average discrete change
College degree	.072***(.010)	.060***(.011)	.036**(.011)	.019(.012)
Panel B: Cross-model differences
$AD C_{Model 1} - AD C_{Model 2}$	.072–.060 =	.012**(.004)
$AD C_{Model 2} - AD C_{Model 3}$	.060–.036 =	.024***(.004)
$AD C_{Model 4} - AD C_{Model 3}$	.036–.019 =	.017***(.004)

Note: Controls include age, age squared, marital status, parental status, gender, political views, religious affiliation, and survey year. Standard errors are in parentheses.

p < .01 and ***p < .001 (two-tailed test).

6.3. Variable Comparisons: Predictions and Marginal Effects Using Alternative Predictors

Researchers often perform specification checks on important variables to ensure their findings are robust. Our framework provides a way to systematically compare the magnitude of effects across alternative operationalizations of a predictor. For example, when comparing sexual orientation groups, a researcher may want to see if different ways of measuring sexual orientation influence substantive conclusions. How someone identifies and how someone behaves are distinct dimensions of sexual orientation (Mize 2015), and these factors can be the basis of different measurement strategies. We consider whether classifying someone’s sexual orientation on the basis of sexual orientation identity versus history of sexual behavior leads to different conclusions.

Data for this example come from the GSS pooled from 2008 to 2016 (Smith et al. 2016). The outcome of interest is whether a person views same-sex relationships as “not wrong at all” (coded 0) or as “sometimes,”“almost always,” or “always” wrong (coded 1). For simplicity, category 1 is referred to as viewing same-sex relationships as wrong. Both sexual orientation identity and sexual orientation behavior are coded into three categories: heterosexual, bisexual, and gay/lesbian. We fit two binary logit models predicting views of same-sex relationships on the basis of a person’s sexual orientation identity (model 1) or sexual orientation behavior (model 2) along with some control variables.

Figure 5 shows predicted probabilities from the models. Views of same-sex relationships are nearly identical across models for individuals classified as heterosexual and are similar for those classified as gay/lesbian. For individuals classified as bisexual, however, the probability of viewing same-sex relationships as wrong is substantially lower when they are classified on the basis of their sexual identity rather than their past sexual behavior. Panel A of Table 4 presents the same information shown in Figure 5 (the predicted probabilities of thinking same-sex relationships are wrong). Column 1 shows probabilities from model 1 using sexual orientation identity, column 2 from model 2 using sexual behavior, and column 3 shows tests of the cross-model difference. The formal tests support the general inferences described earlier: there is a significant difference in opinions of same-sex relationships for bisexual individuals depending on whether classification is based on identity or behavior.

Figure 5.

Predicted probabilities of thinking same-sex relationships are sometimes wrong on the basis of sexual orientation identity and behavior ( $N = 4, 921$ ).

Table 4.

Predicted Probabilities and Average Discrete Changes for Effects of Sexual Orientation on Thinking Same-Sex Relationships Are Sometimes Wrong (n = 4,921)

	Model 1	Model 2
	Sexual Identity	Sexual Behavior	Cross-Model Difference
Panel A: Predicted probabilities
Heterosexual	.531	.533	.002(.002)
Bisexual	.257	.436	.179***(.041)
Gay/lesbian	.104	.171	.067(.040)
Panel B: Average discrete changes
Heterosexual – bisexual	.274***(.041)	.097***(.023)	.177***(.043)
Heterosexual – gay/lesbian	.428***(.034)	.362***(.038)	.066(.041)
Bisexual – gay/lesbian	.153**(.052)	.265***(.044)	–.111(.059)

Note: All models include controls for gender, education, age, race, and year of the survey. Standard errors are in parentheses.

p < .01 and ***p < .001 (two-tailed test).

Columns 1 and 2 of panel B in Table 4 present pairwise comparisons among the three sexual orientation groups in the probability of viewing same-sex relationships as wrong. Column 3 shows the cross-model difference in these effects. The contrasts between heterosexual and gay/lesbian and between bisexual and gay/lesbian do not differ significantly across the two models. However, there is a significantly larger difference between heterosexual and bisexual individuals when using the identity measure, compared with using the sexual behavior measure. In this example, different ways of measuring sexual orientation change the effect size.

6.4. Variable Comparisons: Marginal Effects on Different Outcomes

This example uses negative binomial regression to model two count outcomes: (1) the number of days in the past month a person’s mental health was not good and (2) the number of days in the past month a person’s physical health was not good. Independent variables include standard demographic characteristics. The data come from the 2002, 2006, 2010, and 2014 GSS (Smith et al. 2016). Typically, to compare effects across outcomes, the outcomes need to have the same level of measurement. For example, it would be difficult to compare a marginal effect of a predictor on a continuous outcome to the marginal effect for the same predictor on a binary outcome. For our example, days of poor mental and days of physical health are counts measured with identical response options (select a number between 0 and 30 days). Therefore, we are comfortable comparing whether the effect of a predictor differs across the outcomes. For example, we can compare whether the effect of age on the number of days reporting poor mental health is different from the effect of age on days of poor physical health.

Coefficients from a negative binomial regression model are in the metric of the log of the expected count (when exponentiated these coefficients are multiplicative effects on the expected count or rate). By exponentiating the predictions using the formula $E (y | x) = \exp (x b)$ , the predictions and marginal effects are expressed in terms of the rate or expected count. We use the equations in Section 3 to compute the ADC in the number of days of poor health for each outcome.

The first column of Table 5 presents ADCs for the effects of each predictor on days of poor mental health, column 2 contains the effects on physical health, and column 3 shows the difference in ADCs across the two outcomes. We find, for example, that women report $0.99$ more days of poor mental health, and 0.77 more days of poor physical health, per month, than do men. Although the effect of gender is $0.22$ larger for mental health, the difference is not statistically significant. In contrast, being married significantly reduces the days of poor mental health by $1.01$ , and the effect of marriage on physical health is nonsignificant. A test of the cross-model difference shows that the effect of marriage on mental health is significantly larger than the effect of marriage on physical health (cross-model difference $= - . 851$ , $p < . 01$ ). Similarly, the effect of age differs significantly across the two outcomes, with aging associated with fewer poor mental health days but more poor physical health days.

Table 5.

Average Discrete Changes for the Number of Days Reporting Poor Mental and Physical Health (n = 5,062)

	Poor Mental Health Days	Poor Physical Health Days	Cross-Model Difference
Woman	.993***(.208)	.773***(.174)	.220(.229)
Married	–1.010***(.230)	–.159(.192)	–.851**(.250)
Parent	.274(.248)	–.269(.219)	.543*(.274)
College	–.879***(.229)	–.542**(.189)	–.337(.254)
Age + SD	–.462***(.108)	.491***(.118)	–.953***(.132)
Family income + SD	–.445***(.118)	–.381***(.102)	–.064(.133)
Black vs. White	–1.016***(.258)	–.531*(.217)	–.485(.279)
Other vs. White	–.438(.356)	.145(.343)	–.583(.385)
Other vs. Black	.578(.400)	.676(.374)	–.098(.431)

Note: Average discrete changes for continuous variables are for a standard deviation increase. All models also include a control for survey year. Standard errors are in parentheses.

p < .05, **p < .01, and ***p < .001 (two-tailed test).

6.5. Model Comparison: Predictions and Marginal Effects across Nominal and Ordinal Models

A researcher might prefer an ordinal model instead of a nominal model because the ordinal model requires estimating far fewer parameters. However, if the relationship between an independent variable and the outcome does not conform to the constraints imposed by an ordinal model, an ordinal model can produce results that do not accurately describe the relationships in the data (Long 2014). Several statistical tests can be used to examine the assumptions of ordinal logistic regression (e.g., Brant 1990; Wolfe and Gould 1998). Here we focus on whether the substantive conclusions one would draw from predicted probabilities and marginal effects differ across the two models. The dependent variable is party affiliation, coded 0 = strong Democrat, 1 = Democrat, 2 = independent, 3 = Republican, and 4 = strong Republican. Age is the focal independent variable, with age and age-squared included in the models. Controls for education, parental status, marital status, family income, employment status, and region of the United States are also included. This example uses pooled GSS data from 2010 to 2016 (Smith et al. 2016).

To compare how age is associated with party affiliation across the ordinal logit model and the multinomial logit model, we begin by computing the predicted probabilities of each political party affiliation as age changes from 20 to 80, holding other variables at their means. For both the ordinal and multinomial models, we make predictions for all $j$ outcome categories. For the ordinal logit model, the probability for outcome $m$ is

\Pr (y = m | x) = F (τ_{m} - x β) - F (τ_{m - 1} - x β),

where $F$ is the cumulative distribution function for the logistic distribution with $Var (ε) = π^{2} / 3$ , and $τ_{m}$ are the estimated cutpoints. For the multinomial logit model, the probability for outcome $m$ is

\Pr (y = m | x) = \frac{\exp (x β_{m | 1})}{\sum_{j = 1}^{J} \exp (x β_{j | 1})} .

Plotting the predictions, Figure 6 shows clear differences in the effect of age across the two models. For example, age has almost no effect on the probability of being independent in the ordinal model (left panel), but it has a large negative effect in the multinomial model (right panel). As another example, in the ordinal model the effect of age on being a strong Republican is slightly negative, but it has a stronger and positive effect in the multinomial model. Similarly, age has a much stronger effect on being a strong Democrat in the nominal model than in the ordinal model. Overall, the effects of age are weaker in the ordinal logit model than in the multinomial logit model.

Figure 6.

Comparing predicted probabilities of political party affiliation by age from ordinal logit and multinomial logit models ( $N = 8, 179$ ).

Figure 6 also shows that the effects of age vary across the age range. To illustrate this point, we estimate discrete changes for a 10-year increase in age for someone who is 20 years old and for someone who is 60 years old. Other variables are held at their means. Table 6 shows these effects along with tests of the equality of the effects across models. For someone who is 20 years old, the effect of a 10-year increase in age is significantly different across the ordinal and nominal models for three of the five outcome categories (top panel of Table 6). For example, the effect of aging on the probability of identifying as a strong Democrat is positive in both models, but it is significantly larger in the nominal model. Even more striking, the effects of age on being a strong Republican are in opposite directions across the two models.

Table 6.

Comparing Marginal Effects of Age on Political Party Affiliation from Ordinal Logit and Multinomial Logit Models (n = 8,179)

Discrete Change of Years in Age	Ordinal Logit Model	Multinomial Logit Model	Cross-Model Difference
From age 20 to 30
Strong Democrat	.020***(.004)	.032***(.003)	–.012***(.003)
Democrat	.023***(.005)	–.010(.011)	.034***(.009)
Independent	–.003***(<.001)	<–.001(.010)	–.003(.010)
Republican	–.024***(.005)	–.031**(.011)	.007(.009)
Strong Republican	–.016***(.004)	.010**(.003)	–.026***(.004)
From age 60 to 70
Strong Democrat	–.003(.004)	.029***(.007)	–.032***(.005)
Democrat	–.002(.003)	–.018**(.007)	.016*(.006)
Independent	.001(.001)	–.037***(.004)	.037***(.004)
Republican	.003(.004)	.009(.005)	–.006(.005)
Strong Republican	.002(.002)	.017***(.004)	–.016***(.003)

Note: All models control for race, income, education, marital status, parental status, and employment status. Standard errors are in parentheses.

p < .05, **p < .01, and ***p < .001 (two-tailed test).

For someone who is 60 years old, the effects of age are nonsignificant for all party affiliations in the ordinal model (bottom panel of Table 6). In stark contrast, in the nominal model, the effects of age are significant for four of the five outcomes. Cross-model differences of the effects are significant for these same four outcome categories. It is clear that the ordinal and multinomial models lead to different substantive conclusions: Figure 6 shows quite different patterns across the two models, and Table 6 shows these differences are often statistically significant. Long (2014) explained how the constraints implicit in ordinal models can lead to these differences.

6.6. Group Differences: Marginal Effects on Nonoverlapping Observations

A special case of our general framework is comparing effects across nonoverlapping or mutually exclusive sets of observations. For example, suppose a question is asked in different years of a repeated, cross-sectional survey. A researcher may wish to fit separate models for each year, with different respondents being interviewed at each time point. Another common situation in which observations do not overlap is when fitting models across groups, such as separate models for men and women or separate models for different religious groups. In these cases, we can fit the models separately for each group/sample or fit a single model with interactions between the group/sample indicator and each independent variable. Although these approaches may seem different, they are statistically equivalent (see Long and Mustillo forthcoming). It is tempting to assume that the cross-model covariances are zero when models use nonoverlapping observations, which is often done in the literature. However, as explained in Section 5.1, there are common situations in which mutually exclusive observations across models do not imply that the cross-model covariances are zero. To illustrate this point, we treat the region that respondents were sampled from as a source of clustering and use the cluster-robust variance estimator. Because of this, the cross-model covariances are not zero, and thus $z$ and $z'$ are not equivalent.

For our example, we are interested in whether political polarization surrounding the role government should play in providing health care increased over the past 30 years. Our outcome is whether a respondent thinks the government is responsible for providing health care (coded 1) or believes it is an individual’s responsibility (coded 0). Our focal independent variable is a binary measure of political views indicating if a respondent identifies as conservative or not. All models include controls for gender, race, age, income, education, marital status, parental status, and employment status. Although this question has been asked regularly in the GSS, our example compares responses from only the 1986 and 2016 surveys.

Table 7 presents results from binary logit models fit using the 1986 sample and the 2016 sample. All predictions and marginal effects are computed using only the sample-specific observations. The first column presents the predicted probabilities for respondents who do not identify as conservative. From 1986 to 2016, the probability of thinking health care is the government’s responsibility increased by 0.073 for respondents who do not identify as conservative ( $p < . 01$ ). The second column, in contrast, shows a significant decrease of 0.094 for conservatives ( $p < . 01$ ), indicating that conservatives in 2016 were less likely to view health care as the government’s responsibility. The third column presents the discrete change for being conservative computed using only sample-specific observations. Conservatives had significantly lower predicted probabilities of thinking the government is responsible for providing health care in 1986 (ADC $= - 0.092$ ) and in 2016 (ADC $= - 0.258$ ). However, the polarization between conservatives and nonconservatives increased over time, as the effect of being conservative is significantly larger in 2016 than in 1986 (cross-model difference of ADCs $= 0.166; p < . 001$ ).

Table 7.

Views of Government’s Responsibility in Providing Health Care Using the 1986 (n = 1,254) and 2016 (n = 1,670) General Social Survey Samples

	Pr(Government Responsible)
	Not Conservative	Conservative	ADC
1986 sample	.516	.425	–.092*(.044)
2016 sample	.589	.331	–.258***(.022)
2016 – 1986 difference	.073**(.027)	–.094**(.028)	.166***(.042)

Note: Standard errors are in parentheses. Models control for gender, race, age, income, education, marital status, parental status, and employment status. Standard errors are corrected for clustering within geographic regions. Predictions and average discrete changes (ADCs) are calculated using only the sample-specific observations.

p < .05, **p < .01, and ***p < .001 (two-tailed test).

7. Software Implementation

To test the equality of effects across models, it is often necessary to compute covariances across models. Even in the special case where cross-model covariances are zero, the software needed to test the equality of estimates generally requires combined covariance matrices for both models, including explicit zeros for covariances constrained to be 0. In Section 4 we explained how this can be done either by stacking data sets or by using equation-level scores as Stata’s suest command does. Canette (2014) explained that the SUEST method can also be implemented using generalized structural equation modeling software, such as Stata’s gsem command (see also Lindsey 2016).

As discussed earlier, stacking has some limitations. First, it can be used only when the models being jointly fit are of the same type (e.g., both are binary logit models). Second, stacking forces ancillary parameters to be equal across models, which is too limiting for many applications. Although identical results for most of our examples can be obtained from either suest or gsem, Stata’s generalized structural equation modeling command gsem makes postestimation calculation of predictions and marginal effects much simpler. In addition, gsem allows predictions and AMEs to be computed separately across different groups or samples. For instance, in example 6.6, using gsem we computed AMEs that were averaged only over the 1986 sample or only the 2016 sample, which would not have been possible if the pair of models were fit using suest.⁷

Overall, the specific software used to implement our approach is a matter of convenience, because identical results can, with few exceptions, be obtained with either approach. For those who wish to use gsem, as we do in the example script files available at https://www.trentonmize.com/software/mecompare, it is important to run update all, because Stata has recently made changes that make some postestimation commands run after gsem much faster. We are currently working on a Stata package to simplify the computation of the type of effect comparisons shown throughout our examples; interested readers are encouraged to visit the Web site mentioned above for updates.

Finally, as demonstrated in Section 5, bootstrapping might be another way to compute tests of cross-model differences that account for nonzero covariances across models. However, we caution that more work needs to be done to confirm that bootstrapping returns results similar to those from the SUEST method we advocate. One advantage of the SUEST approach is that it is much less computationally intensive than bootstrapping. For example, for the estimates shown in Figure 2, we used 1,000 bootstrap replications, which took roughly 150 times longer than obtaining the comparable estimates and cross-model comparisons using the SUEST approach.

8. Conclusions

Researchers almost always estimate more than one model in the course of examining a research question. Our article presents a general framework for comparing predictions and effects across models. We emphasize the utility of this approach for comparing marginal effects, as it provides a useful way to summarize an independent variable’s effect across many types of linear and nonlinear models and across outcomes measured in different metrics. In addition, marginal effects avoid the rescaling problems of the regression coefficients in logit- and probit-based models, making them appropriate summary measures to compare effect sizes across a range of models.

When cross-model comparisons are of interest, it is important to provide a statistical test of whether the effect size is the same across models. This is easily accomplished using SUEST. SUEST is an extremely general method for obtaining cross-model covariances for any model that produces equation-level scores. General methods for estimating marginal effects are also available for most common regression models. Our examples illustrate a variety of situations in which comparing effects across models provides valuable substantive insights. Our approach is not, however, limited to the applications we considered; we suspect many other cases of comparing predictions and effects across models could fruitfully use the general framework we advance in this article.

Footnotes

Acknowledgements

We thank David Dalenberg and Josh Doyle for their suggestions and the editor and reviewers at Sociological Methodology for their instrumental feedback. Jeff Pitblado and David Drukker provided valuable help on using Stata to implement our framework.

Authors’ Note

All data and script files are available on the first author’s Web site ().

ORCID iD

Trenton D. Mize

Notes

Author Biographies

Trenton D. Mize is an assistant professor of sociology and a core faculty member for the Advanced Methodologies cluster at Purdue University. His research primarily focuses on (1) gender and sexuality, with a particular focus on labor market outcomes; (2) categorical data analysis and experimental methodology; and (3) health behavior and health inequalities. Recent publications examine sexual orientation inequalities in the labor market (American Sociological Review), how status influences sexual orientation categorization (American Sociological Review), gender inequalities in leadership positions (Social Problems), and nonlinear interaction effects (Sociological Science).

Long Doan is an assistant professor of sociology at the University of Maryland. His research focuses on social psychological processes underpinning patterns of inequality. Current projects examine the emotional consequences of time use, responses to identity threats, and attitudes toward gender and sexuality.

J. Scott Long is a distinguished professor emeritus of sociology and statistics at Indiana University. His methodological research deals with categorical data analysis and methods for reproducible research. He is a fellow of the American Statistical Association and the Sociological Research Association. He is the 2002 recipient of the Paul F. Lazarsfeld Memorial Award for Distinguished Contributions to Sociological Methodology.

References

Agresti

Alan

. 2013. Categorical Data Analysis. 3rd ed. Hoboken, NJ: John Wiley.

Allison

Paul. D.

1999. “Comparing Logit and Probit Coefficients across Groups.” Sociological Methods & Research 28(2):186–208.

Amemiya

Takeshi

. 1981. “Qualitative Response Models: A Survey.” Journal of Economic Literature 19(4):1483–1536.

Arellano

Manuel

. 1987. “Practitioners’ Corner: Computing Robust Standard Errors for Within-Groups Estimators.” Oxford Bulletin of Economics and Statistics 49(4):431–34.

Bishop

Yvonne

Fienberg

Stephen E.

Holland

Paul W.

1975. Discrete Multivariate Analysis: Theory and Practice. Cambridge, MA: MIT Press.

Brant

Rollin

. 1990. “Assessing Proportionality in the Proportional Odds Model for Ordinal Logistic Regression.” Biometrics 46(4):1171–78.

Breen

Richard

Holm

Anders

Karlson

Kristian Bernt

. 2014. “Correlations and Nonlinear Probability Models.” Sociological Methods & Research 43(4):571–605.

Breen

Richard

Karlson

Kristian Bernt

Holm

Anders

. 2018. “Interpreting and Understanding Logits, Probits, and Other Nonlinear Probability Models.” Annual Review of Sociology 44:39–54.

Cameron

A. Colin

Trivedi

Pravin K.

2005. Microeconometrics: Methods and Applications. New York: Cambridge University Press.

10.

Canette

Isabel

. 2014. “Using gsem to Combine Estimation Results.” The Stata Blog. Retrieved May 16, 2019. http://blog.stata.com/2014/08/18/using-gsem-to-combine-estimation-results/.

11.

Chow

Gregory C.

1960. “Tests of Equality between Sets of Coefficients in Two Linear Regressions.” Econometrica 28:591–605.

12.

Clogg

Clifford C.

Petkova

Eva

Haritou

Adamantios

. 1995. “Statistical Methods for Comparing Regression Coefficients between Models.” American Journal of Sociology 100(5):1261–93.

13.

Dowd

Bryan E.

Greene

William H.

Norton

Edward C.

2014. “Computation of Standard Errors.” Health Services Research 49(2):731–50.

14.

Efron

Bradley

Tibshirani

Robert J.

1993. An Introduction to the Bootstrap. New York: Chapman & Hall.

15.

Gelman

Andrew

Stern

Hal

. 2006. “The Difference between Significant and Not Significant Is Not Itself Statistically Significant.” American Statistician 60(4):328–31.

16.

Gould

William

. 1996. “crc43: Wald Test of Nonlinear Hypotheses after Model Estimation.” Stata Technical Bulletin Reports 5:15–18.

17.

Harris

Kathleen M.

2009. “The National Longitudinal Study of Adolescent to Adult Health.”Chapel Hill: Carolina Population Center, University of North Carolina at Chapel Hill.

18.

Hausman

Jerry A.

1978. “Specification Tests in Econometrics.” Econometrica 46(6):1251–72.

19.

Karlson

Kristian Bernt

Holm

Anders

Breen

Richard

. 2012. “Comparing Regression Coefficients between Models Using Logit and Probit: A New Method.” Sociological Methodology 42:286–313.

20.

King

Gary

Tomz

Michael

Wittenberg

Jason

. 2000. “Making the Most of Statistical Analyses: Improving Interpretation and Presentation.” American Journal of Political Science 44(2):347–61.

21.

Liang

Kung-Yee

Zeger

Scott L.

1986. “Longitudinal Data Analysis Using Generalized Linear Models.” Biometrika 73(1):13–22.

22.

Lindsey

Charles

. 2016. “Multiple Equation Models: Estimation and Marginal Effects Using gsem.” The Stata Blog. Retrieved May 17, 2019. https://blog.stata.com/2016/06/07/multiple-equation-models-estimation-and-marginal-effects-using-gsem/.

23.

Long

J. Scott

. 1997. Regression Models for Categorical and Limited Dependent Variables. Thousand Oaks, CA: Sage.

24.

Long

J. Scott

. 2014. “Regression Models for Nominal and Ordinal Outcomes.” Pp. 173–203 in Regression Analysis and Causal Inference, edited by Best

Wolf

London: Sage Ltd.

25.

Long

J. Scott

Freese

Jeremy

. 2014. Regression Models for Categorical Dependent Variables Using Stata. 3rd ed. College Station, TX: Stata Press.

26.

Long

J. Scott

Mustillo

Sarah A.

Forthcoming. “Using Predictions and Marginal Effects to Compare Groups in Regression Models for Binary Outcomes.” Sociological Methods and Research.

27.

MacKinnon

David P.

Lockwood

Chondra M.

Hoffman

Jeanne M.

West

Stephen G.

Sheets

Virgil

. 2002. “A Comparison of Methods to Test Mediation and Other Intervening Variable Effects.” Psychological Methods 7(1):83–104.

28.

McKelvey

Richard D.

Zavoina

William

. 1975. “A Statistical Model for the Analysis of Ordinal Level Dependent Variables.” Journal of Mathematical Sociology 4:103–20.

29.

Mize

Trenton D.

2015. “What Social Psychology Can Contribute to the Study of Sex, Gender, and Sexual Orientation.” Sociology Compass 9(12):1066–81.

30.

Mustillo

Sarah A.

Lizardo

Omar A.

McVeigh

Rory M.

2018. “Editors’ Comment: A Few Guidelines for Quantitative Submissions.” American Sociological Review 83(6):1281–83.

31.

Norton

Edward C.

Miller

Morgen M.

Kleinman

Lawrence C.

2013. “Computing Adjusted Risk Ratios and Risk Differences in Stata.” Stata Journal 13(3):492–509.

32.

Paternoster

Raymond

Brame

Robert

Mazerolle

Paul

Piquero

Alex

. 1998. “Using the Correct Statistical Test for the Equality of Regression Coefficients.” Criminology 36(4):859–66.

33.

Phillips

P.C.B.

Park

Joon Y.

1988. “On the Formulation of Wald Tests of Nonlinear Restrictions.” Econometrica 56(5):1065–83.

34.

Preacher

Kristopher J.

Hayes

Andrew F.

2008. “Asymptotic and Resampling Strategies for Assessing and Comparing Indirect Effects in Multiple Mediator Models.” Behavior Research Methods 40(3):879–91.

35.

Preacher

Kristopher J.

Kelley

Ken

. 2011. “Effect Size Measures for Mediation Models: Quantitative Strategies for Communicating Indirect Effects.” Psychological Methods 16(2):93–115.

36.

Smith

Tom W.

Davern

Michael

Freese

Jeremy

Morgan

Stephen L.

2016. “General Social Surveys: 1972–2016.” Chicago: National Opinion Research Center, University of Chicago.

37.

Weesie

Jeroen

. 1999. “Seemingly Unrelated Estimation and Cluster-Adjusted Sandwich Estimator.” Stata Technical Bulletin 9:231–48.

38.

White

Halbert

. 1982. “Maximum Likelihood Estimation of Misspecified Models.” Econometrica 50(1):1–24.

39.

White

Halbert

. 1984. Asymptotic Theory for Econometricians. San Diego, CA: Academic Press.

40.

Williams

Richard

. 2009. “Using Heterogeneous Choice Models to Compare Logit and Probit Coefficients across Groups.” Sociological Methods & Research 37(4):531–59.

41.

Wolfe

Rory

Gould

William

. 1998. “sg76: An Approximate Likelihood-Ratio Test for Ordinal Response Models.” Stata Technical Bulletin 42:24–27.