Stepwise model selection, Hosmer-Lemeshow statistics and prediction success of model in nested logistic regression in R

Question

is it possible to do stepwise (direction = both) model selection in nested binary logistic regression in R? I would also appreciate if you can teach me how to get:

• Hosmer-Lemeshow statitistic,
• Odds ratio of the predictors,
• Prediction success of the model.

I used lme4 package of R. This is the script I used to get the general model with all the independent variables:

nest.reg <- glmer(decision ~ age + education + children + (1|town), family = binomial, data = fish)

where:

• fish -- dataframe
• decision -- 1 or 0, whether the respondent exit or stay, respectively.
• age, education and children -- independent variables.
• town -- random effect (where our respondents are nested)

Now my problem is how to get the best model. I know how to do stepwise model selection but only for linear regression. (step( lm(decision ~ age + education + children, data = fish), direction +"both")). But this could not be used for binary logistic regression right? also when i add (1|town) to the formula to account for the effects of town, I get an error result.

By the way... I'm very much thankful to Manoel Galdino who provided me with the script on how to run nested logistic regression.

Thank you very much for your help.

Answer 1

I really appreciate the pointers to my book and papers and R package. Briefly, stepwise regression is invalid as it destroys all statistical properties of the result as well as faring poorly in predictive accuracy. There is no reason to use ROC curves to guide model selection (if model selection is even a good idea), because we have the optimum measure, the log-likelihood and its variants such as AIC. Thresholds for the dependent variable should be dealt with using ordinal regression instead of making a series of binary models. The Hosmer-Lemeshow test is now considered obsolete by many statisticians as well as the original authors. See the reference below (which proposes a better method, implemented in the rms package).

@ARTICLE{hos97com, author = {Hosmer, D. W. and Hosmer, T. and {le Cessie}, S. and Lemeshow, S.}, year = 1997, title = {A comparison of goodness-of-fit tests for the logistic regression model}, journal = Statistics in Medicine, volume = 16, pages = {965-980}, annote = {goodness-of-fit for binary logistic model;difficulty with Hosmer-Lemeshow statistic being dependent on how groups are defined;sum of squares test (see cop89unw);cumulative sum test;invalidity of naive test based on deviance;goodness-of-link function;simulation setup;see sta09sim} }

See also

@Article{sal09sim, author = {Stallard, Nigel}, title = {Simple tests for the external validation of mortality prediction scores}, journal = Statistics in Medicine, year = 2009, volume = 28, pages = {377-388}, annote = {low power of older Hosmer-Lemeshow test;avoiding grouping of predicted risks;logarithmic and quadratic test;scaled $\chi^2$ approximation;simulation setup; best power seems to be for the logarithmic (deviance) statistic and for the chi-square statistics that is like the sum of squared errors statistic except that each observation is weighted by $p(1-p)$} }

Answer 2

Disclaimer: This is merely a comment but it won't fit as such, so I'll leave it as a CW response.

Everything is already available in Frank Harrell's rms package (which model to choose, how to evaluate its predictive performance or how to validate it, how not to fall into the trap of overfitting or stepwise approach, etc.), with formal discussion in his textbook, Regression Modeling Strategies (Springer, 2001), and a nice set of handouts on his website.

Also, I would recommend the following papers if you're interested in predictive modeling:

1. Aliferis, C.F., Statnikov, A., Tsamardinos, I., Schildcrout, J.S., Shepherd, B.E., and Harrell, F.E. Jr (2009). Factors Influencing the Statistical Power of Complex Data Analysis Protocols for Molecular Signature Development from Microarray Data. PLoS ONE 4(3): e4922.
2. Harrell, F.E. Jr, Margolis, P.A., Gove, S., Mason, K.E., Mulholland, E.K., Lehmann, D., Muhe, L., Gatchalian, S., and Eichenwald, H.F.. (1998). Development of a clinical prediction model for an ordinal outcome: the World Health Organization Multicentre Study of Clinical Signs and Etiological agents of Pneumonia, Sepsis and Meningitis in Young Infants. WHO/ARI Young Infant Multicentre Study Group. Statistics in Medicine, 17(8): 909-44.
3. Harrell, F.E. Jr, Lee, K.L., and Mark, D.B. (1996). Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine, 15(4): 361-87.

Answer 3

I'm here again =)!

It happens that right now I'm fitting a nested logistic regression and I have to choose the better model as well. Actually I don't know how to do stepwise with lme4, nonetheless I'm not sure it is advisble to use AIC to choose the best model (better fit). Take a look at this link; is dissuss very briefly when AIC may be apropriate with hierarchical (nested) models. Just to provide some backgroud on the link, it is a FAQ about DIC, which is like AIC and is an output of Bugs.

All in all, I don't know how to answer your question, but I can provide some guidance on how I choose my own models. I take a look at two things: ROC curves and fake data simulation. ROC curves are a good graphic way to compare the fit of models regarding true positive rates and false positive rates. Depending on what is more sensitive in your case, you can choose a model with higher true positive rate or higher true negative rates.

If you use Zelig package, I think you can fit models with lme4 and plot ROC curves. The documentation has some nice examples of this.

Another way (my prefered one) is to simulate data with parameters from the model estimated. Then, you check if the simulated data resembles the real data. This is a visual comparison, but can be very effetive. If you need reference on this procedure, Gelman and Hill's book is a standard texto book on this approach. For instance, assume you fitted the model:

P(y=1) = 10 + 2*x Then you create a new variable, y.fake, which is for (i in 1:n) y.fake ~ rbinom(1,1, 10 + 2*x[i])

Then you check if y.fake resembles the true y. You can compute, for instance, the number of cases that the y.fake is equal to the true y. You can ran more than one simulation to take into acount sample variability and you can take into acount the variability of estimated coeficients as well.

Hope this helps.