When a continuous predictor in a model has an unknown, complex association with the response, one common approach is to include higher-order effects of the predictor (quadratic, cubic, and so on) in the model. These polynomial effects can capture aspects of the complexity of the predictor's association, but the shapes that they can accommodate are limited. A more flexible approach is to use a spline effect to represent the predictor. A spline is a series of piecewise polynomials joined at points called knots. You can use the EFFECT statement to define splines for one or more predictors in your model. The EFFECT statement offers many ways to define a spline. See "EFFECT Statement: Spline Effects" in the "Shared Concepts and Topics" chapter of the SAS/STAT User's Guide. See also SAS Note 57975 for more on using spline effects.
Tools for assessing a continuous variable involved in interactions, with and without splines, are discussed further in SAS Note 67024. As described there, the Margins macro (SAS Note 63038) is generally the most useful tool. However, that macro does not support some models including models with constructed effects such as splines.
After fitting a model that includes a spline effect for a continuous predictor, it is often of interest to plot the association of the spline and the estimated response mean. This can be done using either the PLOTS= option or the EFFECTPLOT statement in the modeling procedure, if available. Alternatively, if the procedure offers the STORE statement, you can use the STORE statement followed by the EFFECTPLOT statement in PROC PLM.
In addition, you might want to estimate the effect of changing the original predictor by one or more units and plot how that effect changes across the predictor's range. For a categorical predictor, the effect of changing it can be assessed using the LSMEANS and related statements, but those statements do not operate on continuous predictors. When the effect measure is the response mean difference over a one unit change, it approximates the slope of a line drawn tangent to the spline curve at each setting of the original predictor and any other predictors in the model. This is the marginal effect of the predictor. For a binary response model such as a logistic model, the odds ratio or relative risk (prevalence ratio) might be the preferred effect measure to be estimated and plotted as it changes over the original predictor.
The following illustrates how these plots can be obtained.
Plot the Change in the Odds Ratio over a Splined Predictor
For a binary response model, such as a logistic model, the odds ratio is typically a statistic of interest that assesses how a unit change in a predictor changes the odds of the response event. If the predictor appears as a main effect in the model, its estimated odds ratio is a single value. But if the predictor is represented in the model by a spline effect, the odds ratio changes at each level of the original predictor.
The ODDSRATIO statement in PROC LOGISTIC can present a plot of the changing odds ratio estimates along with confidence limits at each requested value. Without the PLOTS= option, it presents the plot in a vertical orientation. To get the more usual horizontal orientation, add the PLOTS= option as shown below.
The following example uses the Kyphosis data in the example titled "Generalized Additive Model with Binary Data" in the PROC GAM documentation. The EFFECT statement defines a natural (or restricted) cubic spline on the StartVert predictor and includes it in the model along with another predictor, NumVert. The first EFFECTPLOT statement plots the spline against the predicted event probability over the range of StartVert. The second EFFECTPLOT statement with the LINK option does the same but with the log odds (or logit, which is the linear predictor) on the vertical axis. The ODDSRATIO statement provides a table and plot of the odds ratio estimates for StartVert over the range from 2 to 15. The PLOTS= option uses the TYPE= suboption to change the orientation of the plot. The STORE statement saves the fitted model for use in the next example.
proc logistic data=kyphosis plots(only)=oddsratio(type=vertical); effect SplSV=spline(StartVert/naturalcubic basis=tpf(noint)); model Kyphosis(event="1") = SplSV NumVert; effectplot fit(x=StartVert); effectplot fit(x=StartVert) / link; oddsratio StartVert / at(StartVert=2 to 15); store kmod; ods output orplot=or; run;
The first two plots from the EFFECTPLOT statements show how the spline on StartVert is associated with the event probability and the log odds of the event (the linear predictor). As noted at the bottom of the plots from the EFFECTPLOT statements, the plots control for NumVert by holding NumVert at its mean. The plot from the ODDSRATIO statement shows how the event odds change for a unit increase in StartVert for unit increases starting at each of the specified values. Notice in the second plot how a unit increase in StartVert has almost zero change on the log odds from 2 to 6. Since the difference in log odds is a log odds ratio, a zero log odds ratio translates to an odds ratio estimate of 1 by exponentiation and that is what the ODDSRATIO plot shows for StartVert values from 2 to 6 where the log odds plot begins to drop, which translates to the odds ratio becoming less than 1. From about StartVert=12, that drop is constant implying that the odds ratio stays approximately the same - at about 0.6:
Instead of using bars for the individual odds ratios as above, the odds ratios can be displayed as a curved line fitted to the odds ratios and a confidence band. This can be done by adding an ODS OUTPUT statement to save the data that produces the plot and then using the capabilities of PROC SGPLOT to redraw the plot as desired. This is one of the methods for modifying ODS graphs discussed in SAS Note 24529.
Using the odds ratio plot data saved by the ODS OUTPUT statement above, the following statements create a numeric variable, XVAL, containing the value of StartVert that appears in the text labels on the horizontal axis in the above odds ratio plot. Then, in PROC SGPLOT, the BAND and SERIES statements are used with the appropriate variables in the saved plot data to create the desired plot.
data or;
set or(where=(displaylabel ? 'StartVert')); length xval 8; xval=scan(displaylabel,-1,'='); run; proc sgplot noautolegend; band upper=uppercldisplay lower=lowercldisplay x=xval / transparency=.5; series y=oddsratioestdisplay x=xval; xaxis label='Startvert start point' grid; yaxis label='Odds Ratio' grid; refline 1 / axis=y; title "Odds Ratios with 95% Wald Confidence Limits"; run;
Estimate and Plot the Risk Difference over a Splined Predictor
In the previous model, suppose that the effect of unit increases in StartVert on the probability of kyphosis is of interest rather than the odds ratio. Note that the probability of event (kyphosis in this example) is the response mean in a logistic model. You can use the following method that employs the NLMeans macro (SAS Note 62362). This method can similarly be used for estimating and plotting differences in means for other generalized linear models.
Using the NLMeans macro and PROC LOGISTIC
These statements refit the model in PROC LOGISTIC and use the ESTIMATE statement to estimate the means (kyphosis probabilities) over a desired range of StartVert values (2 to 15) and with NumVert fixed at its mean (4.20482). Nonpositional syntax is used to specify the value of the splined variable, StartVert, in each case. The E option produces the coefficients tables, and the ODS OUTPUT statement saves them in data set C.
proc logistic data=kyphosis;
effect SplSV = spline(StartVert/naturalcubic basis=tpf(noint));
model Kyphosis(event="1") = SplSV NumVert;
estimate
'2' intercept 1 SplSV [1,2] NumVert 4.20482,
'3' intercept 1 SplSV [1,3] NumVert 4.20482,
'4' intercept 1 SplSV [1,4] NumVert 4.20482,
'5' intercept 1 SplSV [1,5] NumVert 4.20482,
'6' intercept 1 SplSV [1,6] NumVert 4.20482,
'7' intercept 1 SplSV [1,7] NumVert 4.20482,
'8' intercept 1 SplSV [1,8] NumVert 4.20482,
'9' intercept 1 SplSV [1,9] NumVert 4.20482,
'10' intercept 1 SplSV [1,10] NumVert 4.20482,
'11' intercept 1 SplSV [1,11] NumVert 4.20482,
'12' intercept 1 SplSV [1,12] NumVert 4.20482,
'13' intercept 1 SplSV [1,13] NumVert 4.20482,
'14' intercept 1 SplSV [1,14] NumVert 4.20482,
'15' intercept 1 SplSV [1,15] NumVert 4.20482 /
ilink e;
ods output coef=c;
store kmod;
run;
Results from the ESTIMATE statement appear below.
The following NLMeans macro call reads the saved model (kmod) and the coefficients used by the ESTIMATE statement to compute the differences between successive pairs of means as specified by diff=seq. Since the original model is a logistic model, link=logit is specified. options=reverse is also specified so that differences in increasing StartVert are computed.
%nlmeans(instore=kmod, coef=c, link=logit, diff=seq,
options=reverse, title=One unit risk differences)
In the Label column in the NLMeans results, "1" in the kth position and "-1" in the (k-1)st position indicate that the (k-1)st mean in the previous table of means is subtracted from the kth mean forming the risk difference estimate for a one-unit increase appearing in the Estimate column. So, 0.001730 is the change in kyphosis probability from increasing StartVert from 2 to 3 and -0.02851 is the probability change from an increase from 14 to 15.
The NLMeans macro automatically saves its results in a data set named EST. The following DATA step adds the StartVert values that begin each one unit increase from the previous data set of predicted means. The differences in response mean (probability of kyphosis) are then plotted using PROC SGPLOT.
data SVrange; do StartVert=2 to 14; output; end; run; data kplot; merge SVrange est; run; proc sgplot data=kplot noautolegend; band upper=upper lower=lower x=StartVert; series y=estimate x=StartVert / markers; refline 0; title "Kyphosis risk difference for unit increases in StartVert"; run;
The plot shows that the effect of increasing StartVert by one unit decreases the probability of kyphosis most near StartVert=10.
Using the RISK statement in PROC LOGSELECT
You can easily produce the same risk difference estimates using the RISK statement in PROC LOGSELECT beginning in SAS® Viya® 2026.05. See SAS KB0046630, which discusses possible differences in results from the RISK statement and the NLMeans macro. After starting a CAS session and creating a CAS libref called SASCAS1, these statements fit the same model as above. The RISK statement produces the one-unit change risk differences at values of StartVert from 2 to 14.
data sascas1.kyphosis; set kyphosis; run; proc logselect data=sascas1.kyphosis; effect SplSV=spline(StartVert/naturalcubic basis=tpf(noint)); model Kyphosis(event="1") = SplSV NumVert; risk startvert / at(startvert=2 to 14); displayout risks=risks; run;
Using the data set of risk differences created by the DISPLAYOUT statement, the following statement produces the plot of the risk differences.
data risks; set sascas1.risks; length StartVert 8; StartVert=scan(at,-1,"="); run; proc sgplot data=risks noautolegend; band upper=upper lower=lower x=StartVert; series y=estimate x=StartVert / markers; refline 0; title "Kyphosis risk difference for unit increase in StartVert"; run;
Estimate and Plot the Relative Risk over a Splined Predictor
Now suppose that the effect of changing StartVert on the relative risk – that is, the ratio of kyphosis probabilities – is of interest rather than the odds ratio or risk difference.
Using the NLMeans macro and PROC LOGISTIC
The NLMeans macro can again be used just by adding ratio in options=. This method can be used for estimating and plotting ratios of means for any generalized linear model that does not use the identity link:
%nlmeans(instore=kmod, coef=c, link=logit, diff=seq,
options=ratio reverse)
data kplot;
merge SVrange est;
run;
proc sgplot data=kplot noautolegend;
band upper=upper lower=lower x=StartVert;
series y=estimate x=StartVert / markers;
refline 1;
title "Relative risk for unit increase in StartVert";
run;
In the Label column, similar to the table of differences above, "1" in the kth position and "/1" in the (k-1)st position indicate that the kth mean in the previous table of means is the numerator and the (k-1)st mean is the denominator forming the relative risk estimate in the Estimate column:
Using the RISK statement in PROC LOGSELECT
As done above for risk differences, PROC LOGSELECT can produce the relative risk estimates and plot. The TYPE=RELATIVE option request risk ratios rather than differences.
proc logselect data=sascas1.kyphosis; effect SplSV=spline(StartVert/naturalcubic basis=tpf(noint)); model Kyphosis(event="1") = SplSV NumVert; risk startvert / type=relative at(startvert=2 to 14); displayout risks=risks; run; data risks; set sascas1.risks; length StartVert 8; StartVert=scan(at,-1,"="); run; proc sgplot data=risks noautolegend; band upper=upper lower=lower x=StartVert; series y=estimate x=StartVert / markers; refline 0; title "Kyphosis risk ratios for unit increase in StartVert"; run;
Estimate and Plot the Mean Difference in an Ordinary Regression Model
In the following example, the ORTHOREG procedure is used to model the logSalary variable in the SASHELP.BASEBALL data set. While the NLMeans macro can be used as above, specifying link=identity, the mean differences can be estimated directly in ORTHOREG using the ESTIMATE statement. When the model includes a spline effect, it is easiest to use the nonpositional syntax available in that statement, particularly when many estimates are needed. To obtain a single or very few estimates, you can use the HAZARDRATIO approach shown in the spline example in SAS Note 67024 to show the coefficients needed in the ESTIMATE statement using the traditional, positional syntax.
The PROC ORTHOREG statements below fit the model. The EFFECTPLOT statement plots the curve showing the association of the YrMajor predictor on the logSalary mean while holding the covariates at their means. The ESTIMATE statement estimates the effects of changing YrMajor by one unit starting at values from YrMajor=6 to 18. The nonpositional syntax allows you to specify coefficients that apply to levels of YrMajor by referring to values of YrMajor directly without having to deal with the individual spline components. The statement then translates that specification into the coefficients that apply to the individual spline components. So, a specification such as splYr [-1,6] [1,7] computes the YrMajor=7 - YrMajor=6 difference in the same way that a specification such as GROUP -1 1 differences the parameters of a binary GROUP predictor in a model that has GROUP as a CLASS variable.
The ODS OUTPUT statement saves the results from the ESTIMATE statement in data set DIFFS. The RENAME option changes the name of the variable that contains the labels in the ESTIMATE statements from LABEL to YrMajor. This data set is then used in PROC SGPLOT to plot the mean difference estimates along with their confidence limits. A reference line is drawn at a mean difference of zero.
proc orthoreg data=sashelp.baseball;
effect splYr=spline(YrMajor / naturalcubic basis=tpf(noint));
model logSalary = splYr nAtBat nHits nBB crBB nOuts;
effectplot fit(x=YrMajor);
estimate
'6' splYr [-1,6] [1,7] ,
'7' splYr [-1,7] [1,8] ,
'8' splYr [-1,8] [1,9] ,
'9' splYr [-1,9] [1,10] ,
'10' splYr [-1,10] [1,11],
'11' splYr [-1,11] [1,12],
'12' splYr [-1,12] [1,13],
'13' splYr [-1,13] [1,14],
'14' splYr [-1,14] [1,15],
'15' splYr [-1,15] [1,16],
'16' splYr [-1,16] [1,17],
'17' splYr [-1,17] [1,18],
'18' splYr [-1,18] [1,19] / e cl;
ods output estimates=diffs(rename=(label=YrMajor));
run;
proc sgplot data=diffs noautolegend;
band upper=upper lower=lower x=YrMajor;
series y=estimate x=YrMajor / markers;
refline 0;
title "Unit change effect on mean(logSalary)";
run;
As with the logistic model above, if you imagine lines drawn tangent to the spline curve in the first plot at YrMajor values from 6 to 18, you can see that those lines have positive slope beginning at 6, then the slope decreases to zero at the peak of the curve near 11, and then becomes negative beyond 11. This trend can be seen in the table of response mean differences and is depicted in the plot from PROC SGPLOT:
Estimating a Splined Predictor's Effect when the Spline Is Involved in Interaction
If the spline effect is involved in interactions with other predictors, then the effect of a unit change in the splined predictor is a function of more than just the spline parameters. As a result, the ESTIMATE statement specification is more complex. But the nonpositional syntax in the ESTIMATE statement makes it relatively easy to write. As noted in the previous example, you can use the HAZARDRATIO approach shown in the spline example in SAS Note 67024 to show the coefficients needed in the ESTIMATE statement using the traditional, positional syntax when only a few estimates are needed.
Consider the spline example discussed in SAS Note 67024 in which the splined BaseDeficit predictor interacts with age and the square of age. As shown there, the contrast of model parameters that estimates the change in response mean for a unit change in BaseDeficit involves only model effects that contain the spline. The parameters for the intercept, age, and the square of age drop out when estimating the difference in two BaseDeficit settings. Then, for each model term, the two sets of brackets that follow first contain a 1 for the higher BaseDeficit level in the difference or a -1 for the lower level. This is followed by a comma and then the level(s) needed for that term. For the spline term, only the BaseDeficit value is needed. For example, s [-1,-10] [1,-9] is the spline term contribution to the change from -10 to -9 in BaseDeficit. Similarly, s*age [-1, 10 -10] [1,10 -9] is the contribution for the s*age interaction term when age=10. Note that the value of the spline always follows ordinary variable values. Finally, s*age*age [-1, 10 10 -10] [1,10 10 -9] is the contribution for s*age*age. Note that age appears twice in the term, so its value, 10, must also appear twice in the brackets before the BaseDeficit value.
Rather than writing multiple specifications for the one unit differences at multiple levels of BaseDeficit, a simple macro can be written to generate them. The BD macro below generates a separate ESTIMATE statement as described above. The resulting ESTIMATE statements generate estimates of the mean response change for unit increases in BaseDeficit starting at BaseDeficit values from -20 to -10 and with age fixed at 10:
%macro bd();
%do i=-20 %to -10;
estimate "&i"
s [-1,&i] [1,%eval(&i+1)]
s*age [-1,10 &i] [1,10 %eval(&i+1)]
s*age*age [-1,10 10 &i] [1,10 10 %eval(&i+1)] / cl;
%end;
%mend;
proc orthoreg data=diabetes; effect s = spline(BaseDeficit / naturalcubic); model logCP = s|age|age; effectplot fit(x=BaseDeficit) / at (age=10); %bd() ods select fitplot estimates; ods output estimates=diffs(rename=(label=BaseDeficit)); run;
proc sgplot data=diffs noautolegend; band upper=upper lower=lower x=BaseDeficit; scatter y=estimate x=BaseDeficit; series y=estimate x=BaseDeficit; refline 0; title "Mean difference for one unit change"; run;
The first plot shows how the BaseDeficit spline is associated with the response mean. Note the downward trend, leveling at about BaseDeficit=-14, and then increasing. This is reflected in the difference estimates that follow and in the plot of the differences: