Example 8.8: more Hosmer and Lemeshow

This is a special R-only entry.

In Example 8.7, we showed the Hosmer and Lemeshow goodness-of-fit test. Today we demonstrate more advanced computational approaches for the test.

If you write a function for your own use, it hardly matters what it looks like, as long as it works. But if you want to share it, you might build in some warnings or error-checking, since the user won't know its limitations the way you do. (This is likely good advice even if you are the only one to use your code!)

In R, you can add another layer of detail so that your function conforms to standards for built-in functions. This is a level of detail we don't pursue in our book, but is worth doing in many settings. Here we provided a modified version of a Hosmer-Lemeshow test sent to us by Stephen Taylor of the Auckland University of Technology. We've added a few annotations.

Note that the function accepts a glm object, rather than the two vectors our function used.

hosmerlem2 = function(obj, g=10) {
# first, check to see if we fed in the right kind of object
   stopifnot(family(obj)$family=="binomial" && family(obj)$link=="logit")
   y = obj$model[[1]]
   # the double bracket (above) gets the index of items within an object
   if (is.factor(y)) 
      y = as.numeric(y)==2
   yhat = obj$fitted.values
   cutyhat = cut(yhat, quantile(yhat, 0:g/g), include.lowest=TRUE)
   obs = xtabs(cbind(1 - y, y) ~ cutyhat)
   expect = xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
   if (any(expect < 5))
      warning("Some expected counts are less than 5. Use smaller number of groups")
   chisq = sum((obs - expect)^2/expect)
   P = 1 - pchisq(chisq, g - 2)
   # by returning an object of class "htest", the function will perform like the 
   # built-in hypothesis tests
   return(structure(list(
      method = c(paste("Hosmer and Lemeshow goodness-of-fit test with", g, "bins", sep=" ")),
      data.name = deparse(substitute(obj)),
      statistic = c(X2=chisq),
      parameter = c(df=g-2),
      p.value = P
   ), class='htest'))
}

We can run this using last entry's data from the HELP study.

ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
attach(ds)
logreg = glm(homeless ~ female + i1 + cesd + age + substance,
  family=binomial)

The results are the same as before:

> hosmerlem2(logreg)
 Hosmer and Lemeshow goodness-of-fit test with 10 bins

data:  logreg 
X2 = 8.4954, df = 8, p-value = 0.3866

Example 8.7: Hosmer and Lemeshow goodness-of-fit

The Hosmer and Lemeshow goodness of fit (GOF) test is a way to assess whether there is evidence for lack of fit in a logistic regression model. Simply put, the test compares the expected and observed number of events in bins defined by the predicted probability of the outcome. This can be calculated in R and SAS.

R

In R, we write a simple function to calculate the statistic and a p-value, based on vectors of observed and predicted probabilities. We use the cut() function (1.4.10) in concert with the quantile() function (2.1.5) to make the bins, then calculate the observed and expected counts, the chi-square statistic, and finally the associated p-value (Table 1.1). The function allows the user to define the number of bins but uses the common default of 10.

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))
}

We'll run it with some of the HELP data (available at the book web site). Note that fitted(object) returns the predicted probabilities, if the object is the result of a call to glm() with family=binomial.

ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
attach(ds)
logreg = glm(homeless ~ female + i1 + cesd + age + substance,
  family=binomial)
hosmerlem(y=homeless, yhat=fitted(logreg))

This returns the following output:

$chisq
[1] 8.495386
$p.value
[1] 0.3866328

The Design package, by Frank Harrell, includes the le Cessie and Houwelingen test (another goodness-of-fit test, Biometrics 1991 47:1267) and is also easy to run, though it requires using the package's function for logistic regression.

library(Design)
mod = lrm(homeless ~ female + i1 + cesd + age + substance,   
    x=TRUE, y=TRUE, data=ds)
resid(mod, 'gof')

Sum of squared errors     Expected value|H0                    SD
         104.1091804           103.9602955             0.1655883
                   Z                     P
           0.8991269             0.3685851

The two tests are reassuringly similar.

SAS

In SAS, the Hosmer and Lemeshow goodness of fit test is generated with the lackfit option to the model statement in proc logistic (section 4.1.1). (We select out the results using the ODS system.)

ods select lackfitpartition lackfitchisq;
proc logistic data="c:\book\help.sas7bdat";
   class substance female;
   model homeless = female i1 cesd age substance /  lackfit;
run;

This generates the following output:

            Partition for the Hosmer and Lemeshow Test

                           HOMELESS = 1            HOMELESS = 0
  Group       Total    Observed    Expected    Observed    Expected

      1          45          10       12.16          35       32.84
      2          45          12       14.60          33       30.40
      3          45          15       15.99          30       29.01
      4          45          17       17.20          28       27.80
      5          45          27       18.77          18       26.23
      6          45          20       20.28          25       24.72
      7          45          23       22.35          22       22.65
      8          45          25       25.04          20       19.96
      9          45          28       27.67          17       17.33
     10          48          32       34.95          16       13.05

Hosmer and Lemeshow Goodness-of-Fit Test

Chi-Square       DF     Pr > ChiSq

   8.4786        8         0.3882

The partition table shows the observed and expected count or events in each decile of the predicted probabilities.

The discrepancy between the SAS and R results is likely due to the odd binning SAS uses; the test is unstable in the presence of ties to the extent that some authorities suggest avoiding it. In general, with continuous predictors, the objections are not germane.