There are two ways of estimating how much the event probability changes when a given predictor is changed by one unit. The marginal effect of a predictor is defined as the partial derivative of the event probability with respect to the predictor of interest. A more direct measure is the change in predicted probability for a unit change in the predictor.
Being a derivative, the marginal effect is the slope of the line that is drawn tangent to the fitted probability curve at the selected point. Note that the marginal effects depend on the variable settings that correspond to the selected point at which this tangent line is drawn, so the marginal effect of a variable is not constant. The slope of the tangent line is the change in event probability, p, measured at two points one unit apart along this straight line. If the fitted probability curve is approximately linear (as it is near p=0.5) at the selected point, then the tangent line will closely approximate the fitted curve and the marginal effect will closely approximate the probability change when changing the predictor by one unit. But in areas where the curve is nonlinear (near the smallest and largest values of p), the marginal effect might deviate substantially from the true change. Also, note that the derivative is not strictly defined for discrete explanatory variables. So, in these cases, directly computing the change in predicted probability, as discussed below, might be preferred.
For a binary logistic main-effects model, logit(p)=Σixiβi , the marginal effect of xi is equal to p(1p)bi , where p is the event probability at the chosen setting of the predictors and bi is the parameter estimate for xi . The binary probit main-effects model is Φ-1(p)=Σixiβi , where Φ-1 is the inverse of the cumulative normal distribution function, or probit. The marginal effect of xi in the probit model is equal to φ(x'b)bi , where φ(x'b) is the density function of the standard normal distribution evaluated at x'b, x'b is the product of the row vector of chosen covariate values, x, and the column vector of parameter estimates, b, and bi is the parameter estimate for xi . For a general definition of marginal effects that includes ordinal responses, see the section below on the ordinal logistic model and the Details section of the QLIM documentation.
This example illustrates three methods of estimating marginal effects in a binary logistic model. It also illustrates the change in the predicted probabilities method. Beginning in SAS 9.1, marginal effects for binary or ordinal responses, computed as partial derivatives, are available from SAS/ETS PROC QLIM by specifying the MARGINAL option in the OUTPUT statement. While no procedure in SAS/STAT software directly computes marginal effects, the partial derivative can be computed using results from the procedure used to fit the model. Note that many SAS procedures can fit the binary logistic model as discussed in this note on the kinds of logistic models available in SAS.
This example uses the cancer remission data presented in the example titled "Stepwise Logistic Regression and Predicted Values" in the PROC LOGISTIC documentation. The following statements create the data set.
data Remiss; input remiss cell smear infil li blast temp; datalines; 1 .8 .83 .66 1.9 1.1 .996 1 .9 .36 .32 1.4 .74 .992 0 .8 .88 .7 .8 .176 .982 0 1 .87 .87 .7 1.053 .986 1 .9 .75 .68 1.3 .519 .98 0 1 .65 .65 .6 .519 .982 1 .95 .97 .92 1 1.23 .992 0 .95 .87 .83 1.9 1.354 1.02 0 1 .45 .45 .8 .322 .999 0 .95 .36 .34 .5 0 1.038 0 .85 .39 .33 .7 .279 .988 0 .7 .76 .53 1.2 .146 .982 0 .8 .46 .37 .4 .38 1.006 0 .2 .39 .08 .8 .114 .99 0 1 .9 .9 1.1 1.037 .99 1 1 .84 .84 1.9 2.064 1.02 0 .65 .42 .27 .5 .114 1.014 0 1 .75 .75 1 1.322 1.004 0 .5 .44 .22 .6 .114 .99 1 1 .63 .63 1.1 1.072 .986 0 1 .33 .33 .4 .176 1.01 0 .9 .93 .84 .6 1.591 1.02 1 1 .58 .58 1 .531 1.002 0 .95 .32 .3 1.6 .886 .988 1 1 .6 .6 1.7 .964 .99 1 1 .69 .69 .9 .398 .986 0 1 .73 .73 .7 .398 .986 ;
These statements fit a logistic model to the binary REMISS response using predictors BLAST and SMEAR. The MARGINAL option in the OUTPUT statement adds marginal effect estimates to the OUT= data set. The resulting data set contains a marginal effect estimate for each observation using the predictor values in that observation. PROC PRINT displays the OUT= data set and all marginal effects.
proc qlim data=remiss; model remiss=blast smear / discrete(d=logistic); output out=outqlim marginal; run; proc print data=outqlim noobs; var smear blast meff:; run;
The MARGINAL option provides the marginal effect on the probabilities of both response levels and for both predictors. So, variable Meff_P1_blast is the marginal effect of BLAST on the first level (0) of the response.
Similarly, for the probit model:
proc qlim data=remiss; model remiss=blast smear / discrete; output out=outqlim marginal; run; proc print data=outqlim noobs; var smear blast meff:; run;
To compute the marginal effects using results from a model fit with PROC LOGISTIC (or other modeling procedure), include an ODS OUTPUT statement to write the ParameterEstimates table to a data set. Also specify the OUTPUT statement to save the predicted probabilities for a logistic model (use the P= option), or the x'b values for a probit model (use the XBETA= option) to a data set.
proc logistic data=remiss; model remiss(event="1")=blast smear; ods output parameterestimates=logparms; output out=outlog p=p; run;
Use PROC TRANSPOSE to arrange the parameter estimates as a row and rename them so as not to conflict with the original variable names.
proc transpose data=logparms out=tlog (rename=(blast=tblast smear=tsmear)); var estimate; id variable; run;
Then use a DATA step to combine the parameter estimates and OUT= data sets and compute the marginal effects for each observation in the original data. Only the marginal effects for the response level representing the event of interest (REMISS=1) are computed below. The marginal effect for REMISS=0 could be similarly computed.
data outlog; if _n_=1 then set tlog; set outlog; MEffBlast = p*(1-p)*tblast; MEffSmear = p*(1-p)*tsmear; run; proc print noobs; var smear blast MEff:; run;
Notice that the computed marginal effects match the P2 marginal effects from PROC QLIM. The highlighted values are discussed in the next section.
For the probit model, use the LINK=PROBIT option in PROC LOGISTIC (or use PROC PROBIT) to fit the model. After transposing the parameter estimates data set, use the PDF function to compute the marginal effect for the probit model.
proc logistic data=remiss; model remiss(event="1")=blast smear / link=probit tech=newton; ods output parameterestimates=prbparms; output out=outprb xbeta=xb; run; proc transpose data=prbparms out=tprb (rename=(blast=tblast smear=tsmear)); var estimate; id variable; run; data outprb; if _n_=1 then set tprb; set outprb; MEffBlast = pdf('NORMAL',xb)*tblast; MEffSmear = pdf('NORMAL',xb)*tsmear; run; proc print noobs; var smear blast MEff:; run;
The effect of changing a predictor from one level to another can be directly computed by estimating pxipxj , the difference in event probabilities at levels i and j of the predictor. For a categorical predictor, xj is often an adjacent level (for ordinal predictors) or a reference level (for nominal predictors). For continuous predictors, it is common to look at the effect of a unit change in the predictor: px+1px . But changes of more or less than one unit may be of interest.
To estimate the difference in probabilities, create a data set containing one or more pairs of observations in which the predictor of interest changes by the desired amount. Set all other predictors to the same desired values, such as their means or reference levels, in both observations of the pair. Specify this data set in the DATA= option of the SCORE statement in PROC LOGISTIC. The estimated change in event probability for the given change in the predictor is the difference in predicted values between the two observations of each pair.
For example, data set A below contains three pairs of observations. The PAIR=1 observations will investigate the change in probability associated with changing SMEAR by one unit, from 0.3 to 1.3. The other predictor, BLAST, is held at 1.072. Pair 2 looks at the effect of a unit increase in BLAST from 0.5 to 1.5. In the plot of the fitted model below, note that this is in the relatively straight part of the logistic curve. The third pair computes the probability change for the same unit amount, but this time in the 0 to 1 range where the logistic model is more curved. For the last two pairs, SMEAR is held at 0.63 which is near its mean.
data a; input pair blast smear; datalines; 1 1.072 0.3 1 1.072 1.3 2 0.5 0.63 2 1.5 0.63 3 0 0.63 3 1 0.63 ;
These statements fit the logistic model and use the SCORE statement to compute predicted probabilities for the observations in data set A.
proc logistic data=remiss; model remiss(event="1")=blast smear; score data=a out=preda; run;
PROC TRANSPOSE is used to rearrange the predicted values so that differences can be more easily computed.
proc transpose data=preda out=tpreda; by pair; var p_1; copy blast smear; run;
Finally, the differences in predicted probabilities are computed using a DATA step and displayed by PROC PRINT.
data tpreda; set tpreda; pchange=col2-col1; run; proc print noobs; by pair; var smear blast pchange; run;
The same steps can be used for the probit model.
To compare the marginal effects and the changes in the predicted probabilities, some plots of the logistic curve are helpful. These statements refit the logistic model and create two plots showing the fitted probability curve as a function of BLAST at SMEAR=0.63 and as a function of SMEAR at BLAST=1.072.
ods graphics on; proc logistic data=remiss; model remiss(event="1")=blast smear; score data=a out=preda; effectplot fit(x=smear) / at(blast=1.072) extend=1 noobs nolimits; effectplot fit(x=blast) / at(smear=0.63) noobs nolimits; run;
In the plot below, the logistic curve is very linear. Consequently, the marginal effect (-0.106, highlighted above) closely approximates the change in predicted probabilities (-0.105 = PCHANGE for PAIR=1).
Note that the marginal effect most closely approximates the unit change in predicted probability in the linear part of the logistic curve shown below. The logistic curve is approximately linear for BLAST between 0.5 and 1.5. Over this range, the change in predicted probability is estimated to be 0.38 (PCHANGE for PAIR=2). The marginal effect, highlighted in the table above, is approximately 0.41. For BLAST in the range of 0 to 1, the curve is more nonlinear and the change in predicted probability is about 0.30 (PCHANGE for PAIR=3). A tangent line drawn at (BLAST=1.072, SMEAR=0.63) is steeper than a line connecting the points (BLAST=0, SMEAR=0.63) and (BLAST=1, SMEAR=0.63), so the marginal effect noticeably exceeds the actual change in predicted probability.
The standard error and confidence limits for the marginal effect can be obtained by fitting the model with PROC NLMIXED and using the PREDICT statement to estimate the marginal effect. The following statements fit the same logistic model as above. Note that when you use the BINARY or BINOMIAL distribution in the MODEL statement, the response variable must have only two values: 0 or 1. If you encounter model fitting problems such as convergence failure, you may need to include the PARMS statement to specify the parameter estimates obtained from PROC LOGISTIC or PROC PROBIT as the initial values for the model parameters in PROC NLMIXED.
proc nlmixed data=remiss; p=1/(1+exp(-(Intercept+bb*blast+bs*smear))); model remiss ~ binomial(1,p); predict p*(1-p)*bb out=b; predict p*(1-p)*bs out=s; run;
The PREDICT statements create two data sets containing marginal effect estimates for the two predictors — data set B containing the BLAST estimates and data set S containing the SMEAR estimates. Both are displayed by PROC PRINT.
proc print data=b noobs; var smear blast pred stderrpred lower upper; title "BLAST marginal effects"; run; proc print data=s noobs; var smear blast pred stderrpred lower upper; title "SMEAR marginal effects"; run;
Each data set contains the marginal effect estimate for each observation, its standard error, and a 95% confidence interval. The marginal effect estimates match those found using PROC QLIM and PROC LOGISTIC.
Estimates of marginal effects at predictor settings other than those in the original data can be obtained by simply adding observations with the desired settings to the input data set with the response set to missing. Such observations are ignored when fitting the model but marginal effects are computed as for the other observations.
For example, you may want to obtain marginal effect estimates when all predictors are at their means. The MEANS and DATA steps below compute the means of the predictors and add that observation to the data. Since the data set from PROC MEANS does not contain the response variable, REMISS, it is set to missing when added to the REMISS data set. PROC QLIM fits the model as before.
proc means data=remiss; var blast smear; output out=means mean=; run; data WithMean; set means remiss; run; proc qlim data=WithMean; model remiss=blast smear / discrete(d=logistic); output out=outqlim marginal; run; proc print data=outqlim noobs; where remiss=.; var smear blast meff:; run;
The PRINT step displays the marginal effects for the added observation.
The same approach can be used with PROC LOGISTIC, PROC PROBIT, or PROC NLMIXED. Simply add the above MEANS and DATA steps before the modeling step and specify the new data set in the DATA= option of the modeling procedure.
As noted above, the marginal effect is the partial derivative of the event probability with respect to the variable of interest, xi:
For the case of simple main-effects models as discussed above, logit(p)=Σixiβi , the final partial derivative is just βi yielding p(1-p)βi as the marginal effect of xi as before. For a higher-order model, such as a model involving xi in an interaction or quadratic effect, the marginal effect is slightly more complex. Consider the cancer remission data and a model that includes the main effects of BLAST and SMEAR as well as their interaction:
logit(p) = β0 + βsSMEAR + βbBLAST + βsbSMEAR·BLAST
For this model, the partial derivative of x'β with respect to SMEAR is βs+βsbBLAST, so the marginal effect for SMEAR is p(1-p)(βs+βsbBLAST). Similarly for BLAST.
The marginal effects for this model can be obtained using the following statements.
proc nlmixed data=remiss; p=1/(1+exp(-(Intercept+bb*blast+bs*smear+bsb*smear*blast))); model remiss ~ binomial(1,p); predict p*(1-p)*(bb+bsb*smear) out=b; predict p*(1-p)*(bs+bsb*blast) out=s; run;
Suppose the possible response values are ordered with levels i=1, 2, ... , k. Under the ordinal logistic model (proportional odds model), the probability of response level i is the difference in the cumulative probabilities at level i and level i-1.
pi = F(αi+x'β) - F(αi-1+x'β) ,
where αi is the ith intercept, β contains all non-intercept parameters, and F(x) is the logistic cumulative distribution function F(x)=exp(x)/(1+exp(x)). Then the marginal effect of the jth predictor, xj, on pi is
For a model containing only main-effects, = βj as in the binary logistic model discussed above. For more complex models, replace with the resulting function.
The ordinal model can be fit in many procedures including LOGISTIC, PROBIT, QLIM, and NLMIXED. However, only PROC QLIM can directly provide marginal effect estimates. In the following example, the response is the severity of symptoms with ordered levels: none, mild, severe. These statements create the data set with a numerically coded response Y with levels 1, 2, and 3 corresponding to increasing severity of the symptoms.
data multi; input Prep $ Dose Symptoms $ N; if symptoms='None' then y=1; else if symptoms='Mild' then y=2; else y=3; LDose=log10(Dose); datalines; stand 10 None 33 stand 10 Mild 7 stand 10 Severe 10 stand 20 None 17 stand 20 Mild 13 stand 20 Severe 17 stand 30 None 14 stand 30 Mild 3 stand 30 Severe 28 stand 40 None 9 stand 40 Mild 8 stand 40 Severe 32 test 10 None 44 test 10 Mild 6 test 10 Severe 0 test 20 None 32 test 20 Mild 10 test 20 Severe 12 test 30 None 23 test 30 Mild 7 test 30 Severe 21 test 40 None 16 test 40 Mild 6 test 40 Severe 19 ;
These statements fit the ordinal logistic model and display the marginal effect estimates. The ordinal probit model can be fit using the DISCRETE(DIST=NORMAL) option. PROC QLIM models the probabilities at the high end of the response scale and cumulates the probabilities over increasingly lower response levels.
proc qlim data=multi; freq N; model y=LDose / discrete(dist=logit); output out=outqlim marginal; run; proc print noobs; var y symptoms ldose meff:; run;
Notice that marginal effect estimates are provided for each response level for each predictor in the model.
Alternatively, you can fit the model in PROC LOGISTIC or PROC PROBIT and compute the marginal effect estimates from the cumulative predicted probabilities. These statements fit the same ordinal model as in PROC QLIM above, save the parameter estimates of the model and the cumulative predicted probabilities to data sets, and compute the marginal effects. Note that after transposing the parameter estimates data set, the variable containing the parameter estimate for the LDose predictor is named COL3. The DESCENDING response variable option causes the probabilities at the high end of the response scale to be modeled as above. The resulting probabilities are cumulated over the response levels 3, 2, and 1 and are named CP_3, CP_2, and CP_1 so that CP_1 = 1. The EFFECTPLOT statement produces a plot of the predicted probabilities for the individual (not cumulative) response levels.
proc logistic data=multi; freq N; model y(descending)=LDose; effectplot fit(x=ldose) / individual; output out=outlog predprobs=cumulative; ods output ParameterEstimates=logparms; run; proc transpose data=logparms out=tlog; var estimate; run; data margeff; if _n_=1 then set tlog; set outlog; Meff3=(cp_3*(1-cp_3))*col3; Meff2=(cp_2*(1-cp_2)-cp_3*(1-cp_3))*col3; Meff1=(cp_1*(1-cp_1)-cp_2*(1-cp_2))*col3; run; proc sort; by ldose; run; proc print noobs; var y symptoms ldose meff:; run;
Examining the plot of the individual probabilities, the slope of a line tangent to the Pr(Y=3) curve is positive and increasingly so as LDose increases. This is reflected in the marginal effects for level 3 (Meff3). Apparently an inflection point is crossed before the highest LDose level resulting in a slight decrease in the marginal effect at that level. The Pr(Y=2) curve has a maximum at about LDose = 1.4 so the slope of a tangent line should be decreasingly positive as LDose increases to that maximum and then be increasingly negative beyond. The marginal effects show this pattern. Finally, the decreasing Pr(Y=1) curve seems to have an inflection point near LDose = 1.3 and the marginal effects show this with an increasingly negative value up to the inflection point and decreasingly negative beyond.
The ordinal model can also be fit using PROC NLMIXED as shown in this note, and PREDICT statements can be added to provide an estimate of the marginal effect on the probability of each response level as well as its standard error and a confidence interval.
Suppose the possible response values are unordered with levels i=1, 2, ... , k. Under the generalized logit model commonly used for nominal responses, the probability of response level i is
pi = exp(x'βi)/Σj(exp(x'βj))
Then the marginal effect of the jth predictor, xj, on pi is
For a model containing only main-effects, = βij and = βkj. For more complex models, replace these partial derivatives with the resulting functions.
The generalized logit model can be fit by the LOGISTIC, GLIMMIX, CATMOD, and NLMIXED procedures. Marginal effects are not directly available, but can be computed using the parameter estimates and individual predicted probabilities from any of these procedures.
The following example uses the remote-sensing data presented in the example titled "Linear Discriminant Analysis of Remote-Sensing Data on Crops" in the DISCRIM documentation. The response is type of crop with five possible levels. X1 is one of four variables used to predict the type of crop. The following statements fit a generalized logit model with X1 as predictor and saves the parameter estimates and individual predicted probabilities to data sets. Marginal effects are computed using the above formula for each of the crops using the values of X1 in each of the observations. Note that four generalized logits can be defined on the five crop types. Consequently, the parameter for the last crop type (Sugarbeets) is constrained to zero. After transposing the data set of parameter estimates, the four parameter estimates for X1 are named COL5, COL6, COL7, and COL8. The EFFECTPLOT statement produces a plot of the predicted probabilities for the individual response levels.
proc logistic data=crops; model crop = x1 / link=glogit; effectplot fit(x=x1) / noobs nolimits; output out=preds predprobs=individual; ods output ParameterEstimates=betas; run; proc transpose data=betas out=rowbetas; var estimate; run; data margeff; if _n_=1 then set rowbetas; set preds; SumBetaPred=col5*IP_Clover + col6*IP_Corn + col7*IP_Cotton + col8*IP_Soybeans; MEClover=IP_Clover*(col5-SumBetaPred); MECorn=IP_Corn*(col6-SumBetaPred); MECotton=IP_Cotton*(col7-SumBetaPred); MESoybeans=IP_Soybeans*(col8-SumBetaPred); MESugarbeets=IP_Sugarbeets*(-SumBetaPred); run; proc sort nodupkey; by x1; run; proc print; id x1; var me:; run;
The values of the marginal effects reflect the slopes of lines tangent to each of the crop curves at each X1 setting. For instance, lines tangent to the Soybeans curve have positive slopes up to about 19, then become negative after 20, and essentially zero beyond 50.
Again, the model could be fit in PROC NLMIXED including PREDICT statements to estimate the marginal effects and provide confidence limits.
|Product Family||Product||System||SAS Release|
|Topic:||SAS Reference ==> Procedures ==> PROBIT|
Analytics ==> Regression
SAS Reference ==> Procedures ==> LOGISTIC
SAS Reference ==> Procedures ==> QLIM
Analytics ==> Econometrics
Analytics ==> Categorical Data Analysis
SAS Reference ==> Procedures ==> NLMIXED
|Date Modified:||2005-05-19 10:40:54|
|Date Created:||2002-12-16 10:56:39|