[R] Hosmer-Lemeshow 'goodness of fit'

viostorm rob.schutt at gmail.com
Mon May 9 01:50:39 CEST 2011


I'm trying to do a Hosmer-Lemeshow 'goodness of fit' test on my logistic
regression model.

I found some code here:

http://sas-and-r.blogspot.com/2010/09/example-87-hosmer-and-lemeshow-goodness.html

The R code is above is a little complicated for me but I'm having trouble
with my answer:

Hosmer-Lemeshow: p=0.6163585
le Cessie and Houwelingen test (Design library): p=0.2843620

The above link indicated they should be approximately equal which in my case
they are not, any suggestions or is there a package function people would
recommend in R for use with a logistic regression model?

Thanks in advance,

-Rob Schutt

--------------------------------
Robert Schutt, MD, MCS 
Resident - Department of Internal Medicine
University of Virginia, Charlottesville, Virginia 


------------------


########################################################
# Compute the Hosmer-Lemeshow 'goodness-of-fit' test

cd.full_model = glm(formula = Collaterals ~ CHF + Age + CABG 
	+ relevel (as.factor (num.obst.vessels),"one") 
	+ Current.smoker + DM + HTN + ace.inhibitor 
	+ MI, family = binomial(link = "logit"))

hosmerlem = function(y, yhat, g=10) {
  cutyhat = cut(yhat,
     breaks = quantile(yhat, probs=seq(0,
       1, 1/g)), include.lowest=TRUE)
  obs = xtabs(cbind(1 - y, y) ~ cutyhat)
  expect = xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
  chisq = sum((obs - expect)^2/expect)
  P = 1 - pchisq(chisq, g - 2)
  return(list(chisq=chisq,p.value=P))
}

hosmerlem(y=Collaterals, yhat=fitted(cd.full_model))

########################################################33
# Compute the le Cessie and Houwelingen test 

f <- lrm(Collaterals ~ CHF + Age + CABG 
	+ relevel (as.factor (num.obst.vessels),"one") 
	+ Current.smoker + DM + HTN + ace.inhibitor 
	+ MI, x=TRUE, y=TRUE)

library(Design)  
resid(f, 'gof')

Output:

> hosmerlem(y=Collaterals, yhat=fitted(cd.full_model))
$chisq
[1] 6.275889

$p.value
[1] 0.6163585

> resid(f, 'gof')
Sum of squared errors     Expected value|H0                    SD 
          118.5308396           118.3771115             0.1435944 
                    Z                     P 
            1.0705717             0.2843620 


--
View this message in context: http://r.789695.n4.nabble.com/Hosmer-Lemeshow-goodness-of-fit-tp3508127p3508127.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list