Changing group number in Hosmer Lemeshow Test

[matthewdye14]

I have created a logistic regression model from 13K observations that models a dichotomous outcome (0 or 1). I wanted to run the Hosmer Lemeshow Goodness of Fit test to determine whether the model was a good fit to the data. The HL test, however, is very sensitive to small differences in large datasets and the literature suggests changing the number of groups constructed by the HL test from its default of 10 to 35. Does anyone know how to change the number of groups used by the HL test in SAS Enterprise Guide?

Thank you!

[TomKari]

It looks to me like you can't change the number of groups either in Enterprise Guide, or in the underlying SAS procedure.

If you don't get a good reply from the communities, I suggest you ask SAS technical support. They are really good with questions like this, and might be able to suggest a workaround.

Tom

[Ksharp]

HL test has low power I think, maybe that is reason why SAS don't offer you change number of groups.

I have done it before for SAS-L. Here are SAS Data step code and IML code . Pick up one you like .

data have;
 set sashelp.cars(keep=origin type invoice);
call streaminit(1234);
y=rand('bern',.75);
run;

ods select none;
proc logistic data=have ;
class origin type;
model y(event='1')=origin type invoice;
score data=have out=pred_have(keep= y  p_1);
run;
ods select all;

proc sort data=pred_have;by p_1;run;



%let groups=20;
proc rank data=pred_have out=group_have groups=&groups ;
 var p_1;
 ranks rank_p_1;
run;
proc summary data=group_have ;
by rank_p_1;
var y p_1 ;
output out=summary(drop=_type_) sum=observed1 expected1 ;
run;
data summary;
retain rank_p_1 _freq_ observed1 expected1 observed0 expected0;
 set summary;
 observed0=_freq_-observed1 ;
 expected0=_freq_-expected1;
 label rank_p_1='Group' _freq_='Total' 
      observed1 ='Observed(P=1)' expected1='Expected(P=1)' 
      expected0='Expected(P=0)' observed0='Observed(P=0)';
run;

proc print label noobs ;run;
proc sql noprint;
select sum((observed1-expected1 )**2/expected1 +(observed0-expected0)**2/expected0)
 into : chi
 from summary;
quit;
data result;
 chi=χ
 df=&groups-2;
 p=1-cdf('CHISQ',chi,df); 
label chi='Chi-Square' df='DF' p='Pr > ChiSq';
run;
proc print label noobs ;run;

data have;
 set sashelp.cars(keep=origin type invoice);
call streaminit(1234);
y=rand('bern',.75);
run;

ods select none;
proc logistic data=have ;
class origin type;
model y(event='1')=origin type invoice;
score data=have out=pred_have;
run;
ods select all;






proc iml;
use pred_have;
read all var {y p_0 p_1} into x;
close;

groups=20;

call sort(x,3);
nrow=nrow(x);
id=1:nrow;
cutpoints=do(0,nrow+1,(nrow+1)/groups);
g=t(bin(id,cutpoints));

level=unique(g);
n=ncol(level);
all=j(n,6,0); /*Group N O1 E1 O0 E0*/
do i=1 to n;
 idx=loc(g=level[i]);
 all[i,1]=level[i];
 all[i,2]=nrow(g[idx]);
 all[i,3]= sum(x[idx,1]=1);
 all[i,4]= sum(x[idx,3]); 
 all[i,5]= sum(x[idx,1]=0);
 all[i,6]= sum(x[idx,2]);
end;
print all[l='' c={Group Total Observed1 Expected1 Observed0 Expected0}];
chi=sum((all[,3]-all[,4])##2/all[,4]) + sum((all[,5]-all[,6])##2/all[,6]);
df=groups-2;
p=1-cdf('CHISQ',chi,df); 
print chi[l='Chi-Square'] df[l='DF'] p[l='Pr > ChiSq'];
quit;