Plot the Change in the Odds Ratio over a Splined Predictor using PROC SURVEYLOGISTIC?

[jkromine]

I am using complex survey data and want to use restricted cubic splines to plot the change in the Odds Ratio over a Splined Predictor using PROC SURVEYLOGISTIC.

I have tried using PROC PLM and the EFFECTPLOT statement, however, the Y axis displays "predicted probabilities", rather than displaying odds ratios. 

How do I use PROC SURVEYLOGISTIC to display splines that can easily be interpreted as odds ratios?

Code that I have tried:

TITLE "Splines of effect of age on opioid prescribing";
PROC SURVEYLOGISTIC DATA=work.ed2014_2021_v4;
EFFECT spl_age = spline(age / DETAILS NATURALCUBIC BASIS=TPF KNOTMETHOD=PERCENTILELIST(5 25 50 75 95));
MODEL opioids1ormore (EVENT='1') = spl_age/NOINT CLPARM;
CLUSTER cpsum;
STRATA cstratm;
WEIGHT patwt1k;
STORE splineModel;
OUTPUT OUT=SplineOut PREDICTED=fit;
RUN;

PROC PLM RESTORE=splineModel PLOTS=ALL NOCLPRINT;
EFFECTPLOT FIT(X=age);
ODS OUTPUT fitplot=fitplot;
SCORE DATA=work.ed2014_2021_v4 OUT=plm_out PREDICTED=_predicted LCLM=lower UCLM=upper;
RUN;

[StatDave]

This note shows how the odds ratio plot can be done for non-survey data using PROC LOGISTIC. The note also shows how you could similarly plot the changing risk differences or relative risks. However, since PROC SURVEYLOGISTIC does not have a PLOTS= option as in PROC LOGISTIC, the approach shown in the note must be used which introduces a dummy CLASS variable so that LSMEANS statements can be used. To obtain the odds ratios, the NLEST macro must be used rather than the NLMEANS macro as shown in the note.

The idea is to use pairwise differences of the coefficient vectors from the LSMEANS statements to yield log odds ratio estimates for each one-unit change in the original predictor and then exponentiate each difference to yield the odds ratio estimates. In the code below, the coefficient vectors for the original model are displayed (the introduced dummy variable is removed) after transposing and adding names as used by the NLEST macro. This makes it easy to see what each one-unit difference reduces to when differencing the two vectors involved. The FD data set contains the resulting expressions that exponentiate each difference. Note that the results are similar to the plot produced by PROC LOGISTIC in the beginning section of the note.

proc surveylogistic data=kyphosis;
  effect SplSV=spline(StartVert/naturalcubic);
  model Kyphosis(event="1") = SplSV NumVert;
  store kmod;
  run;
data k2;
   set kyphosis(in=in1) kyphosis;  
   dummy=in1; 
   run;
ods exclude all;
proc surveylogistic data=k2;
  class dummy / param=glm;
  effect SplSV=spline(StartVert/naturalcubic);
  model Kyphosis(event="1") = dummy SplSV NumVert;
  lsmeans dummy / at StartVert=2 e;
  lsmeans dummy / at StartVert=3 e;
  lsmeans dummy / at StartVert=4 e;
  lsmeans dummy / at StartVert=5 e;
  lsmeans dummy / at StartVert=6 e;
  lsmeans dummy / at StartVert=7 e;
  lsmeans dummy / at StartVert=8 e;
  lsmeans dummy / at StartVert=9 e;
  lsmeans dummy / at StartVert=10 e;
  lsmeans dummy / at StartVert=11 e;
  lsmeans dummy / at StartVert=12 e;
  lsmeans dummy / at StartVert=13 e;
  lsmeans dummy / at StartVert=14 e;
  lsmeans dummy / at StartVert=15 e;
  ods output coef=c(drop=row2 where=(parameter ^? "dummy"));
  run;
ods select all;
proc transpose data=c out=c2 prefix=B_P; 
   by lmatrix; var row1; 
   run;
proc print data=c2; run;
data fd; 
   length label f $32767; 
   infile datalines delimiter=',';
   input label f; 
   datalines;
StartVert=2 to 3 , exp(b_p3)
StartVert=3 to 4 , exp(b_p3)
StartVert=4 to 5 , exp(b_p3)
StartVert=5 to 6 , exp(b_p3 + 0.0496*b_p4)
StartVert=6 to 7 , exp(b_p3 + (0.6305-0.0496)*b_p4)
StartVert=7 to 8 , exp(b_p3 + (2.4467-0.6305)*b_p4)
StartVert=8 to 9 , exp(b_p3 + (6.204-2.4467)*b_p4)
StartVert=9 to 10 , exp(b_p3 + (12.579-6.204)*b_p4)
StartVert=10 to 11 , exp(b_p3 + (21.572-12.579)*b_p4)
StartVert=11 to 12 , exp(b_p3 + (32.506-21.572)*b_p4)
StartVert=12 to 13 , exp(b_p3 + (44.675-32.506)*b_p4)
StartVert=13 to 14 , exp(b_p3 + (57.375-44.675)*b_p4)
StartVert=14 to 15 , exp(b_p3 + (70.125-57.375)*b_p4)
;
%nlest(instore=kmod, fdata=fd, title=1 Unit Odds Ratios)

Then this code produces the plot from the EST data set that is automatically created containing the results from the macro.

proc sgplot data=est noautolegend;
  band upper=upper lower=lower x=label;
  scatter y=estimate x=label;
  series y=estimate x=label;
  refline 1;
  run;