Restricted Cubic Splines with Multinomial Outcome

[ChristosK]

I had been generating spline curves for a dichotomous outcome, but now I am looking at a 3 level outcome, although then ordinal scale is not proportional. Mainly concerned with generating curves for 2 vs 0, and 1 vs 0.

This is the code I had been using:

ods select ModelANOVA ParameterEstimates SplineKnots;
proc logistic data=HEERDT.VPICU;
  effect spl = spline(Ind_DBP1 / details naturalcubic basis=tpf(noint)                 
                               knotmethod=percentiles(5) );
   model Vent24 (event ='0')= spl / selection=none;     /* fit model by using spline effects */
   output out=SplineOut predicted=Fit lower = l upper = u;    /* output predicted values for graphing */
   where CPB=1;
   Title 'Probability of Failed Extubation by 24h: Spline Adjusted Model';
run;

proc sgplot data = splineout;
loess x = ind_dbp1 y = fit/nomarkers curvelabel = "Predicted value" nolegclm nolegfit;
loess x = ind_dbp1 y = u/nomarkers   curvelabel = "Upper CI"        nolegclm nolegfit ;
loess x = ind_dbp1 y = l/nomarkers   curvelabel = "Lower CI"        nolegclm nolegfit;
run;

 This is fine for binary outcomes and I have been able to generate curves of quadratic functions of the continuous  exposure 

Ind_DBP1

I have tried entering " link=glogit" within the model statement and stated that the outcome has a reference value if 0  rather than stating a particulate value for the event.
p values come out as highly significant for the knot placements but I get the following error message when I try too plot a graph:

NOTE: PROCEDURE SGPLOT used (Total process time):

real time 0.20 seconds

cpu time 0.14 seconds

NOTE: The number of observations of the LOESS plot (5313) exceeds the limit of 5000. Specify the LOESSMAXOBS option of the ODS

GRAPHICS statement to override the limit.

NOTE: The number of observations of the LOESS plot (5313) exceeds the limit of 5000. Specify the LOESSMAXOBS option of the ODS

GRAPHICS statement to override the limit.

NOTE: The number of observations of the LOESS plot (5313) exceeds the limit of 5000. Specify the LOESSMAXOBS option of the ODS

GRAPHICS statement to override the limit.

WARNING: The LoessPlot statement named 'LOESS' will not be drawn because one or more of the required arguments were not supplied.

WARNING: The LoessPlot statement named 'LOESS1' will not be drawn because one or more of the required arguments were not supplied.

WARNING: The LoessPlot statement named 'LOESS2' will not be drawn because one or more of the required arguments were not supplied.

WARNING: A blank graph is produced. For possible causes, see the graphics template language documentation.

WARNING: The LoessPlot statement named 'LOESS' will not be drawn because one or more of the required arguments were not supplied.

WARNING: The LoessPlot statement named 'LOESS1' will not be drawn because one or more of the required arguments were not supplied.

WARNING: The LoessPlot statement named 'LOESS2' will not be drawn because one or more of the required arguments were not supplied.

WARNING: A blank graph is produced. For possible causes, see the graphics template language documentation.

WARNING: The LoessPlot statement named 'LOESS' will not be drawn because one or more of the required arguments were not supplied.

WARNING: The LoessPlot statement named 'LOESS1' will not be drawn because one or more of the required arguments were not supplied.

WARNING: The LoessPlot statement named 'LOESS2' will not be drawn because one or more of the required arguments were not supplied.

WARNING: A blank graph is produced. For possible causes, see the graphics template language documentation.

NOTE: There were 5313 observations read from the data set WORK.SPLINEOUT.

81

82 OPTIONS NONOTES NOSTIMER NOSOURCE NOSYNTAXCHECK;

95

Any ideas?

Cheers,

Christos 

[WarrenKuhfeld]

Loess is expensive, and there are default limits on the number of observations for loess fit.  Either specify the option below in the ODS GRAPHICS statement or use PBSPLINE instead of LOESS.

 

LOESSOBSMAX=n

specifies an upper limit for the number of observations that can be used with a loess plot.

[ChristosK]

I changed it to this:

proc logistic data=HEERDT.VPICU;
effect spl = spline(Ind_DBP1 / details naturalcubic basis=tpf(noint) knotmethod=percentiles(5) );
model VP24_VasoNE_012(event ='2')= spl / link=glogit ct covb clodds=wald selection=none; /* fit model by using spline effects */
output out=SplineOut predicted=Fit lower = l upper = u; /* output predicted values for graphing */
where CPB=1;
Title '';
run;
proc sgplot data = splineout;
PBSPLINE x = ind_dbp1 y = fit/nomarkers;
PBSPLINE x = ind_dbp1 y = u/nomarkers;
PBSPLINE = ind_dbp1 y = l/nomarkers;
run;

The screenshot accompanying it seems meaningless.

I'm certain that there is an effect, because a quadratic works very well at outcome level 3.

Enclosed is an additional screenshot showing the covariate matrix when reference level for the outcome is specified as '0'.

[ChristosK]

Not sure if any of this screenshots went through, trying again, but using the "add photo" option.

[ChristosK]

Not sure if any of this screenshots went through, trying again, but using the "add photo" option.

[ChristosK]

When I remove the link =glogit statement, the output looks more like a spline, but I cant tell from this what level of the outcome is specified as without link=glogit  I don't see how SAS knows what levels we are comparing to what?

[Rick_SAS]

First, you can create this plot by using the EFFECTPLOT statement. Add the statement

   EFFECTPLOT FIT;

to your PROC LOGISTIC call.

Second, in your problem, the FIT, L, and U variables are the predicted values from a regression procedure. To plot them, it is not appropriate to use LOESS or PBSPLINE statements to smooth the predicted values. Use the SERIES statement for the FIT curve and use the BAND statement for the L and U confidence limits:

ods select ModelANOVA ParameterEstimates SplineKnots;
proc logistic data=Sashelp.cars(where=(Origin in ('Asia' 'USA') & Type ^="Hybrid"));
  effect spl = spline(mpg_city / details naturalcubic basis=tpf(noint)                 
                               knotmethod=percentiles(5) );
   model Origin(event ='USA')= spl / selection=none;     /* fit model by using spline effects */
   output out=SplineOut predicted=Fit lower = l upper = u;    /* output predicted values for graphing */
   effectplot fit;  /* produce fit plot automatically */
run;

/* or produce fit plot manually by using predicted values in output data set */
proc sort data=SplineOut; by mpg_city; run;

proc sgplot data = splineout;
   band x = mpg_city lower=L upper=U;
   series x = mpg_city y = fit / curvelabel = "Predicted value";
   /* optional, uncomment if you want to overlay the upper/lower limits */
   /*
   series x = mpg_city y = U / curvelabel = "Upper CL";
   series x = mpg_city y = L / curvelabel = "Lower CL";
   */
run;

[ChristosK]

Thank you Rick:

There seems to be many ways to do this, although I'm not sure how I can tell which is"more correct" (If such a thing exists).

Why do you think this approach is better and what are the pros / cons? Was my previous approach better suited for other applications?

I was particularly interested in finding a way to run restricted cubic splines rather than just splines. Whilst I understand conceptually what this is, I have not found a way to do this in SAS yet.

Thanks,

Christos 

[Rick_SAS]

> There seems to be many ways to do this, although I'm not sure how I can tell which is"more correct" (If such a thing exists).

Yes, there are many ways to do this and other tasks. Of the two ways I showed, neither is  "more correct," but the EFFECTPLOT statement is definitely easier and requires less effort.

> Why do you think this approach is better and what are the pros / cons? Was my previous approach better suited for other applications?

This approach is better because the curves/band represent actually predicted values from a statistical model. If you use a loess or penalized b-spline to smooth the predicted values, you obtain a curve that is not the one predicted by the model. In many cases, a loess or spline curve does not pass through the predicted values.

 > I was particularly interested in finding a way to run restricted cubic splines rather than just splines. Whilst I understand conceptually what this is, I have not found a way to do this in SAS yet.

I don't know what you mean by "to do this in SAS." What is "this"? You must have used my article on regression splines to write your code. The article contains links to the SAS documentation on splines. SAS supports many kinds of splines. The documentation for PROC TRANSREG also contains details and examples. If you can specify the spline basis that you want to use, someone on this community can show you how to use that basis in a regression. 

[ChristosK]

Regarding "doing this":

Sorry about not being specific.

I was referring to "Restricted Cubic Splines". I have seen this alot in the literature in my field and am including a paper with an example  where they used Stata. (Venkatesan et al).

It seems as though they were able to:

1) Partially dismiss the effects of instability (due to fewer entries)  at the extreme ends of the independent variable of interest.

2) Assign a definitive trough value, and generate OR's in reference to this.

I don't remember reading your article. One of my friends suggested using this code a few months ago for a poster I was presenting.

I guess at the end of the day I need to be able to establish a curve that is publication worthy and that I'll be able to defend.

Will look at your article. 

Much Appreciated,

Christos

[ChristosK]

Loess output and Series output looks the same even though the shading and colours are different . Can only add one file at a time to this message so will send the series output in the next message.

[ChristosK]

This is the second  output where we used series instead of loess in the proc subplot statement .

Can anyone see whether the curve is different at all?

[Rick_SAS]

Yes, if the X values are closely spaced and the Y values are not too wild, you will often obtain a loess curve that looks very similar to a series plot. My point, however, is that it is not statistically appropriate to put a smoother through predicted values. 

For one thing, the loess curve is actually one of a family of curves that are parameterized by a smoothing parameter. The curve that is displayed is chosen to minimize a system of local regression fits, assuming that the points are data with IID errors.

The following example shows data for which the "best" loess curve does not pass anywhere near the fitted points.

data Fit;
input x y @@;
datalines;
0 0   1 5   3 5   4 0   5 5
;

proc sgplot data=Fit;
series x=x y=y / lineattrs=(color=red);
loess x=x y=y;
run;

[ChristosK]

Rick:

Hope you are enjoying the weekend.

I am struggling with getting figures for the adjusted models. I guess with out smoothing the curves, the edges will be erratic but I can't really  present these.

What if anything should I consider doing with this adjusted model in order to make it meaningful or presentable as a figure?

Christos

[Rick_SAS]

The parameter estimates indicate that you are including multiple covariates in the model (age, BMI, etc). You need to use the EFFECTPLOT SLICEFIT statement to create a sliced fit plot. See

Visualize multivariate regression models by slicing continuous variables  (see the 2nd graph, which looks like yours)

If for some reason you can't use the EFFECTPLOT statement, see

How to create a sliced fit plot in SAS