In a binary-response model in which the predictor of primary interest interacts with a continuous covariate, you might want to plot the change in the estimated effect of the primary predictor across the range of the covariate. You could assess the effect of the predictor's effect using one of several statistics such as its odds ratio or the ratio or difference in event probabilities comparing two of its levels (the relative risk and risk difference, respectively).
Note that when the primary predictor is not involved in an interaction, then its odds ratio, relative risk, and risk difference are constant. This odds ratio is automatically provided by PROC LOGISTIC. You can estimate the relative risk as described in SAS Note 23003, and you can estimate the risk difference as shown in SAS Note 37228.
In the following example, n individuals are given various doses of one of two treatments, and the number surviving, r, is recorded. Primary interest is in the effect of the treatment variable, Treat, on the event probability, but that effect is expected to change depending on the level of the Dose variable. Note that in these data, different doses are used in the treatment levels. A logistic model that allows for TREAT and DOSE to interact is fit to the data. This SAS Note shows two ways of fitting the model and saving the effect estimates.
As described above, these data are counts of individual responses aggregated within combinations of the treatment and dose levels. The analyses and plots that follow could similarly be done on data in which observations represent individual responses as shown in Individual level data below.
These statements create the data set and compute the observed event probability, Obsp.
data epidemic; input Treat $ Dose n r Sex @@; Obsp=r/n; datalines; A 2.17 142 142 0 A .57 132 47 1 A 1.68 128 105 1
A 1.08 126 100 0 A 1.79 125 118 0
B 1.66 117 115 1 B 1.49 127 114 0 B 1.17 51 44 1 B 2.00 127 126 0 B .80 129 100 1 ;
These statements plot the change in the observed event probability across the dose for each of the treatment groups.
proc sort data=epidemic; by treat dose; run; proc sgplot data=epidemic; series y=Obsp x=Dose / group=Treat; yaxis label="Observed probability" min=0 max=1 grid; xaxis grid; run;
Using the LOGISTIC Procedure and the NLEST Macro
The logistic model allowing for interaction of the Treat and Dose predictors is fit using the LOGISTIC procedure. The CLASS statement declares Treat as categorical and represents it in the model using a 0,1-coded variable with Treat='B' as the reference level. The change in the predicted event probability, separately for each treatment group, is plotted by the EFFECTPLOT statement. You can use the ODDSRATIO statement to estimate and plot the odds ratio for Treat at a specified set of Dose levels. However, the relative risk (risk ratio) cannot be estimated by PROC LOGISTIC. Estimates and a plot of the risk difference across Dose are also not easily available. The STORE statement saves the fitted model in an item store named LogMod.
proc logistic data=epidemic; class Treat(ref='B') / param=ref; model r/n=Treat|Dose; effectplot slicefit(x=Dose sliceby=Treat) / noobs clm; oddsratio treat / at (dose=0 to 2 by 0.5); store LogMod; run;
While the interaction is not significant (p=0.653), it is left in the model and the analyses and plots are produced for illustrative purposes. The lack of substantial interaction can be seen by adding the LINK option in the EFFECTPLOT statement.
Using the saved model, the NLEST macro can estimate any function of the model parameters. See NLEST macro documentation for details and additional examples of its use. In f=, specify the function to be estimated using B_P1, B_P2, and so on to represent the parameters in the order shown in the PROC LOGISTIC results. The specified function is estimated for each observation in the score= data set and saved in a new data set named in outscore=. By default, 95% confidence intervals are produced. If a different confidence level is desired, specify alpha= in the NLEST macro call.
The following calls of the NLEST macro estimate the relative risk, risk difference, and odds ratio at each combination of Treat and Dose. Those statistics are saved in separate data sets. No displayed results are produced. In each macro call, the f= specification is a difference or ratio of two expressions, the first for Treat='A' and the second for Treat='B'. The LOGISTIC function is used when the event probability is needed (as for the relative risk and risk difference). The EXP function is used to obtain the odds in each treatment for the odds ratio. Note that the specification within each EXP or LOGISTIC function is a log odds for one of the two treatments. Note also that the definition of the odds ratio in the last NLEST call could be further simplified as exp(b_p2+b_p4*dose).
%nlest(instore=logmod,
f=logistic(b_p1+b_p2+b_p3*Dose+b_p4*Dose)/logistic(b_p1+b_p3*Dose),
score=epidemic, outscore=RR)
%nlest(instore=logmod,
f=logistic(b_p1+b_p2+b_p3*Dose+b_p4*Dose)-logistic(b_p1+b_p3*Dose),
score=epidemic, outscore=Diff)
%nlest(instore=logmod,
f=exp(b_p1+b_p2+b_p3*Dose+b_p4*Dose)/exp(b_p1+b_p3*Dose),
score=epidemic, outscore=OR)
In addition to the estimated statistic (in the Predicted Value column) and confidence limits comparing Treat='A' to Treat='B' at each Dose, each data set contains a test that the statistic equals zero. For the odds ratios and relative risks, you might prefer to test the equality to one rather than zero. To do that, add null=1 in the NLEST call. The following statements display the relevant variables in the OR data set. Note that any differences in confidence limits on the odds ratios between the results from PROC LOGISTIC and from the NLEST macro, or from PROC NLMIXED below, are caused by differing estimation methods.
proc sort data=OR; by Dose; run; proc print data=OR label; id Dose Obsp; var Pred--Pr; title "Odds Ratio Estimates"; run;
Using the NLMIXED Procedure
The same logistic model can be fit using the NLMIXED procedure. The PREDICT statement in NLMIXED allows you to estimate and save the odds ratio, relative risk, or risk difference for plotting. So, this is an alternative to the approach above.
The following statements refit the logistic model. The LOGISTIC function is used to define the model on the event probability, p, as 1/(1+exp(-xβ)), where β is the set of model parameters, b1, b2, b3, and b4, and xβ is some function of those parameters. The MODEL statement indicates that the number surviving, r, follows the binomial distribution with parameters n and p.
The PREDICT statements estimate the event probability, relative risk, risk difference, and odds ratio at each combination of Treat and Dose and save those statistics in separate data sets. The first PREDICT statement saves the estimated event probabilities, p, in data set Pred. As done in the NLEST macro above, the specified expression in each of the following statements is a difference or ratio of two expressions, the first for Treat='A' and the second for Treat='B'. The LOGISTIC function is used when the event probability is needed (as for the relative risk and risk difference). The EXP function is used to obtain the odds in each treatment for the odds ratio. Note that xβ alone is the log odds. Note also that the definition of the odds ratio in the last PREDICT statement could be further simplified as exp(b2+b4*dose).
proc nlmixed data=epidemic; p=logistic(b1 + b2*(Treat='A') + b3*Dose + b4*Dose*(Treat='A')); model r ~ binomial(n,p); predict p out=Pred; predict logistic(b1+b2+b3*Dose+b4*Dose) / logistic(b1+b3*Dose) out=RR; predict logistic(b1+b2+b3*Dose+b4*Dose) - logistic(b1+b3*Dose) out=Diff; predict exp(b1+b2+b3*Dose+b4*Dose) / exp(b1+b3*Dose) out=OR; run;
The resulting model fit is very similar to the fit from PROC LOGISTIC above.
Plotting the Treatment Effect
With the statistics saved in separate data sets, you can plot them, with confidence bands, showing the change in each statistic across the range of Dose. The first plot uses the Pred data set saved by PROC NLMIXED and shows the predicted event probabilities for each treatment, effectively replicating the plot from the EFFECTPLOT statement in PROC LOGISTIC above.
proc sgplot data=Pred; band upper=Upper lower=Lower x=Dose / group=Treat transparency=.8; series y=Pred x=Dose/group=treat; yaxis label="Predicted probability" min=0 max=1 grid; xaxis grid; run;
The next plot shows the change in relative risk, which is not readily obtained from PROC LOGISTIC.
proc sort data=RR; by Dose; run; proc sgplot data=RR noautolegend; band upper=Upper lower=Lower x=Dose / transparency=.8; series y=Pred x=Dose; yaxis label="Relative risk" min=0 grid; xaxis grid; run;
The plot of the change in the odds ratio follows. Note that this plot is shown in the more typical horizontal orientation as opposed to the vertical orientation provided by the ODDSRATIO statement in PROC LOGISTIC. Also, the odds ratio at every distinct Dose level is used resulting in a plot with connected odds ratio estimates and continuous confidence bands.
proc sort data=OR; by Dose; run; proc sgplot data=OR noautolegend; band upper=Upper lower=Lower x=Dose / transparency=.8; series y=Pred x=Dose; yaxis label="Odds ratio" grid; xaxis grid; run;
Finally, a plot of the risk difference is provided, which is also not easily available from PROC LOGISTIC.
proc sort data=Diff; by Dose; run; proc sgplot data=Diff noautolegend; band upper=Upper lower=Lower x=Dose / transparency=.8; series y=Pred x=Dose; yaxis label="Risk difference" grid; xaxis grid; run;
Individual Level Data
The same analyses and plots as above can be done on individual response data with minimal change to the syntax. The following statements create individual response data from the original data set of aggregated data and produce the same results. Note that the NLEST macro calls and the code in the Plotting the Treatment Effect section can be used unchanged.
data ind;
set epidemic; do i=1 to r; y=1; output; end; do i=1 to n-r; y=0; output; end; run; proc summary data=ind nway; class Treat Dose; var y; output out=Agg mean=Obsp; run; proc sort data=agg; by treat dose; run; proc sgplot data=agg; series y=Obsp x=Dose / group=Treat; yaxis label="Observed probability" min=0 max=1 grid; xaxis grid; run; proc logistic data=ind; class Treat(ref='B') / param=ref; model y(event="1")=Treat|Dose; effectplot slicefit(x=Dose sliceby=Treat) / noobs clm; oddsratio treat / at (dose=0 to 2 by 0.5); store LogMod; run; proc nlmixed data=ind; p=logistic(b1 + b2*(Treat='A') + b3*Dose + b4*Dose*(Treat='A')); model y ~ binary(p); predict p out=Pred2; predict logistic(b1+b2+b3*Dose+b4*Dose) / logistic(b1+b3*Dose) out=RR; predict logistic(b1+b2+b3*Dose+b4*Dose) - logistic(b1+b3*Dose) out=Diff; predict exp(b1+b2+b3*Dose+b4*Dose) / exp(b1+b3*Dose) out=OR; run;