topic Re: Anyone could help me to solve this program code for natural cubic spline? in New SAS User

Anyone could help me to solve this program code for natural cubic spline?

cjchung — Mon, 30 Dec 2019 06:52:58 GMT

Hello, everyone, I want to explore the association between PM2.5 (X) and COPD (Y).

I defined the variable, PM2.5 as non-linear and treated with a natural cubic spline.

COPD.match_newcopd is my data file. Here is the code....

proc sort data=COPD.match_newcopd 
out=cubic (keep=copd pm25a mari_g diab hear Arthritis asthma cancer smk_g drk_g nut_g);
by pm25a;
run;
ods select Model ANOVA ParameterEstimates SplineKnots;
proc logistic data=cubic outdesign=SplineBasis;
effect spl = spline(pm25a / naturalcubic basis=tpf(noint) details knotmethod=EQUAL(5)) ; 
model COPD (event ='1')= spl /SELECTION = none;
output out=cubicOut5 predicted=Fit upper=u lower=l; /* output predicted values for graphing */
effectplot fit;
quit;

data a; set cubicout;
Fitexp=exp(Fit);
uexp=exp(u);
lexp=exp(l);
data copd.multipm; set work.a;run;

proc sort data=work.a; by pm25a; 
title "Restricted Cubic Spline Regression";
ods graphics/height=6.5in width=10in imagename="Figure2" border=off;
proc sgplot data=a noautolegend;
where Fitexp ^=.;
band x = pm25a lower=Lexp upper=Uexp;
xaxis label='PM2.5 (ug/m3)';
yaxis label='Risk of COPD' ;
series x=pm25a y=Fitexp / curvelabel = "Predicted value";
series x=pm25a y=uexp / curvelabel = "Upper CL";
series x=pm25a y=lexp / curvelabel = "Lower CL";
run;

Question 1: When I add the confounders (such as marriage) in the logistic model, I will get the terrible figure......

like: model COPD (event ='1')= spl marriage /SELECTION = none;

Is anything wrong with the above code?

Question 2: How I choose the best model with knots=3 or 5 by comparing the different AIC values? How I write the SAS code?

Thank you

Re: Anyone could help me to solve this program code for natural cubic spline?

PGStats — Mon, 30 Dec 2019 20:31:54 GMT

Just a guess. You might have more than two levels for COPD. What happens if you add the statement

where COPD in (0, 1);

in your call to proc logistic?

Re: Anyone could help me to solve this program code for natural cubic spline?

cjchung — Tue, 31 Dec 2019 01:25:50 GMT

Hi, thanks for your reply. I checked it and there were just two levels of COPD, Yes and No. Therefore I still did not know the problem.

Thank you.

Re: Anyone could help me to solve this program code for natural cubic spline?

PGStats — Tue, 31 Dec 2019 04:54:19 GMT

OK, I think the problem is the superposition of predicted values for diffrent marriage groups. Try sgpanel:

proc sort data=work.a; 
by marriage pm25a; 
run;

title "Restricted Cubic Spline Regression";
ods graphics/height=6.5in width=10in imagename="Figure2" border=off;
proc sgpanel data=a noautolegend;
where Fitexp ^=.;
panelby marriage;
band x = pm25a lower=Lexp upper=Uexp;
colaxis label='PM2.5 (ug/m3)';
rowaxis label='Risk of COPD' ;
series x=pm25a y=Fitexp / curvelabel = "Predicted value";
series x=pm25a y=uexp / curvelabel = "Upper CL";
series x=pm25a y=lexp / curvelabel = "Lower CL";
run;

or call proc sgplot with a BY marriage statement.