Goodness of fit statistics (Hosmer-Lemeshow) after -glm-

[Joshua Nealon]

Dear all,
Despite looking, I cannot find an answer to what seems a simple problem and wonder if someone here can help.

I'm trying to assess the goodness of fit of a complimentary log-log binomial regression, parameterized using the -glm- command. I would like to use the Hosmer-Lemeshow test to explore a hypothesis that my model is statistically similar to a saturated model.

Following logistic, etc., the respective command in STATA is:

-estat gof -

I would like to perform the equivalent, following -glm-. The full command is:

glm var1, family (binomial) link(cloglog) offset (var2) cluster (var3)

A side note considering the objective - I have read Paul Allison's and others' caution to the use of the HL test. If others have experience and could direct me to code/literature relating to alternative tests within STATA to explore the "acceptability" of this model fit (again, I have sought but not found them), any guidance would be highly appreciated.

I am using STATAv 15.0, with many thanks in advance.

Josh.

[Clyde Schechter]

You can emulate an H-L like calibration statistic following -glm- as follows:

predict phat, mu
xtile decile = phat, nq(10)
collapse (sum) observed = var1 predicted = phat, by(decile)
gen chi2_component = (observed-predicted)^2/predicted
summ chi2_component, meanonly
display "H-L chi2 = " as result %3.2f `r(mean)'

As this code overwrites the data in memory, you may prefer to -preserve- first and then -restore-.

For my part, I tend not to rely on up or down "tests" of calibration in my work. My preference is to actually look at the observed and predicted values in each decile. This enables me to see, for example, if my model fits well in one range but poorly in another, and such things often suggest ways to improve the model. Also, in larger data sets, I will usually split the range of predicted values into a larger number of groups; typically targeting a group size of 50-100. Just my practice.

[Joshua Nealon]

Hi Clyde,
Thanks so much, really appreciate the feedback. This works perfectly.

Just FYI, in this instance, we have an underlying assumption about the distribution of the data because the outcome (historical infectious disease exposure) is a direct function of age. If the data do not fit this distribution, it's interesting that this assumption of constant exposure risk has been broken, a finding which warrants investigation/discussion.

Appreciate your advice to look at the observed/predicted probabilities per decile; doing this will contribute to the discussion above.

Thanks again, will let you know how this goes.

Josh.

[Joshua Nealon]

Hi Clyde,
Sorry to come back to you here - I guess the code you kindly shared calculates the "mean" chi2 per decile but the chi2 statistic should be the sum of the "(o-e)^2/e" for each decile.

Of course this can be altered by multiplying by the number of groups, and I have done this to calculate a p-value (sampling from the chi2 distribution with 8 degrees of freedom).

This isn't exactly the result I expected (in Excel, the p-value is different...), I wonder if you have any advice on this?

Thanks again for your thoughts, if any; highly appreciated.

Josh.

glm var1, family (binomial) link(cloglog) offset (var2) cluster (var3) eform

predict phat, mu
xtile decile = phat, nq(10)
collapse (sum) observed = var1 predicted = phat, by(decile)
gen chi2_component = (observed-predicted)^2/predicted
summ chi2_component, meanonly
display "H-L chi2 = " as result %3.2f `r(mean)'
display "p-value = "  chi2(8,`r(mean)'*10)

[Clyde Schechter]

Sorry, yes, I should have said `r(sum)' instead of `r(mean)'. Or your solution of multiplying by 10 accomplishes the same thing, and either should be applied to the H-L chi2 statistic itself, not just calculating the p-value.

Also, your calculation of the p-value is incorrect. The chi2() function gives you the left tail area, but what you need is the right tail area. So, putting it altogether, it should be:

glm var1, family (binomial) link(cloglog) offset (var2) cluster (var3) eformpredict phat, mu
xtile decile = phat, nq(10)
collapse (sum) observed = var1 predicted = phat, by(decile)
gen chi2_component = (observed-predicted)^2/predicted
summ chi2_component, meanonly
display "H-L chi2 = " as result %3.2f `r(sum)'
display "p-value = " chi2tail(8,`r(sum)')

As for the discrepancy from Excel, I have no idea. Excel's statistical functions are not reliable. But perhaps the discrepancy you are seeing reflects having used chi2() when chi2tail() is needed.

[Joshua Nealon]

Hi Clyde,
I cannot thank you enough for your help - I now have graphs of the predicted vs observed probabilities, associated p-values and data to describe the variance for discussion.
Really appreciate it.
Josh.

[Etwal Bou Raad]

Hello. I was trying to test gof form my glm model. Is there any way to retain my data after collapse?
Thank you so much.

[Matthew Williams]

You can use preserve (before the collapse command) and restore after you have done with L-H code.