Plotting Calibration curve/line

[rajdesai223]

Hi, 

Do you mind posting the entire code from the start please? I don't understand what is phat_mean? Is it the mean of the predicted probability? And what is ob_risk ?

Thanks

[Rick_SAS]

The two blog posts that I linked to have complete code for analyzing a simulated data set.

[rajdesai223]

Hi Rick, 

I have a binary outcome variable called adherence(codes as 0 and 1). I ran the codes you suggested.

Here is my code.

proc logistic data=cases plots=effect;
class sex1 race group age1 hmo ER FmsScore;
model adherence(event="0")=sex1 race group age1 hmo ER FmsScore;
output out=LogiOut predicted=PredProb; /* save predicted probabilities in data set */
run;

proc rank data=LogiOut out=LogiDecile groups=10;
var PredProb;
ranks Decile;
run;

/* Then compute the mean predicted prob and the empirical proportions (and CI) for each decile */
proc means data=LogiDecile noprint;
class Decile;
types Decile;
var adherence PredProb;
output out=file mean=obs_risk PredProbMean
lclm=Lower_CI uclm=Upper_CI;
run;

proc sgplot data=file noautolegend aspect=1;
lineparm x=0 y=0 slope=1 / lineattrs=(color=grey pattern=dash);
loess x=PredProbMean y=obs_risk; /* if you want a smoother based on deciles */
scatter x=PredProbMean y=ob_risk / yerrorlower=Lower_CI yerrorupper=Upper_CI;
yaxis label="Observed Probability of Outcome";
xaxis label="Predicted Probability of Outcome";
run;

Here are the issues I am facing(I have attached a csv file for you)

1) There are unequal number of observations in each decile 

2) my observed risk and predicted risk are very different in each group because of which I am not getting the plot correct.

Please help. 

Thank You again.

[Rick_SAS]

The main issue is that you are modeling 

adherence(event="0")

whereas in my blog I model adherence(event="1"). Thus your graph does not fall along the line y=x.

That's not an important concern. You can either model adherence(event="1") or you can change the LINEPARM statement to conform to the situation that you have, as follows. (Note that I changed some variable names to conform to the names in your CVS file.):

proc sgplot data=A noautolegend aspect=1;
   lineparm x=0 y=1 slope=-1 / lineattrs=(color=grey pattern=dash) clip;
   loess x=PredProbMean y=obs_risk;
   scatter x=PredProbMean y=obs_risk / yerrorlower=yLower yerrorupper=yUpper;
   yaxis label="Observed Probability of Outcome";
   xaxis label="Predicted Probability of Outcome";
run;

[rajdesai223]

Thank You Rick. This was very helpful.