To better ask my question, I have provided some of the outputs from both a 16 variable model (fit) and a 17 variable model (fit2) below (all predictor variables in these models are continuous, where the only difference between these models is that fit does not contain variable 17 (var17)):

fit                    Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       
 Obs        102849    LR chi2   13602.84    R2       0.173    C       0.703    
  0          69833    d.f.            17    g        1.150    Dxy     0.407    
  1          33016    Pr(> chi2) <0.0001    gr       3.160    gamma   0.416    
 max |deriv| 3e-05                          gp       0.180    tau-a   0.177    
                                            Brier    0.190       

fit2                 Model Likelihood       Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       
 Obs        102849    LR chi2   13639.70    R2       0.174    C       0.703    
  0          69833    d.f.            18    g        1.154    Dxy     0.407    
  1          33016    Pr(> chi2) <0.0001    gr       3.170    gamma   0.412    
 max |deriv| 3e-05                          gp       0.180    tau-a   0.177    
                                            Brier    0.190          

I used Frank Harrell's rms package to build these lrm models. As you can see, these models do not appear to vary much, if at all, across Discrimination Indexes and Rank Discrim. Indexes; however, using lrtest(fit,fit2), I was provided with the following results:

 L.R. Chisq         d.f.            P 
3.685374e+01     1.000000e+00    1.273315e-09 

As such, we would reject the null hypothesis of this likelihood ratio test; however, I would assume this is likely due to the large sample size (n = 102849) as these models appear to perform in a similar fashion. Furthermore, I am interested in finding a better way of formally comparing nested binary logistic regression models when n is large.

I greatly appreciate any feedback, R scripts, or documentation that can steer me in the right direction in terms of comparing these types of nested models! Thanks!

  • $\begingroup$ What is the purpose of removing variable 17? $\endgroup$ – Michael M Nov 27 '13 at 5:45
  • $\begingroup$ This is a toy example; however, I am typically asked to build models with 8-12 variables, and removing variables that do not contribute to a model is a primary interest for me. Variable 17 just appears to mean very little to the model as a whole (in terms of predictability), yet the likelihood ratio test tells us there is a significant difference between the two models (likely due to large n rather than an actual difference in these two models). As such, I am hoping to find a way to compare these two models (finding a method that does not indicate a difference between these two models) $\endgroup$ – Matt Reichenbach Nov 27 '13 at 13:33
  • $\begingroup$ (1) not sure I quite understand what you're looking for. But in medicine the issue with using discrimination like c-statistic is well established, c-statics can be unchanged even with addition of significant variable, and has lead to development of reclassification indices (circ.ahajournals.org/content/121/15/1768.full) (2) are the AIC/BIC similar? a number of variable information gain criteria may be more useful than discrimination criteria. $\endgroup$ – charles Nov 28 '13 at 19:15
  • 1
    $\begingroup$ I think there is a typo in your 1st paragraph. It is stated that fit2 is a 17 variable model, but it is also the model that omits V17. You may want to edith this. $\endgroup$ – tomka Dec 2 '13 at 10:07
  • 1
    $\begingroup$ @tomka, I changed fit2 to fit in the above example as per your correction. Thanks! $\endgroup$ – Matt Reichenbach Dec 2 '13 at 19:49

(1) There is an extensive literature on why one should prefer full models to restricted/parsimonious models. My understanding are few reasons to prefer the parsimonious model. However, larger models may not be feasible for many clinical applications.

(2) As far as I know, Discrimination/Discrimination indexes aren’t (?should not be) used as a model/variable selection parameter. They aren’t intended for this use and as a result there may not be much of a literature on why they shouldn’t be used for model building.

(3) Parsimonious models may have limitations that aren’t readily apparent. They may be less well calibrated than larger models, external/internal validity may be reduced.

(4) The c statistic may not be optimal in assessing models that predict future risk or stratify individuals into risk categories. In this setting, calibration is as important to the accurate assessment of risk. For example, a biomarker with an odds ratio of 3 may have little effect on the cstatistic, yet an increased level could shift estimated 10-year cardiovascular risk for an individual patient from 8% to 24%

Cook N.R.; Use and misuse of the ROC curve in the medical literature. Circulation. 115 2007:928-935.

(5) AUC/c-statistic/discrimination is known to be insensitive to significant predictor variables. This is discussed in the Cook reference above, and the motivating force behind the development of net reclassification index. Also discussed in Cook above.

(6) Large datasets can still lead to larger models than desired if standard variable selection methods are used. In stepwise selection procedures often a p-value cut-off of 0.05 is used. But there is nothing intrinsic about this value that means you should choose this value. With smaller datasets a larger p-value (0.2) may be more appropriate, in larger datasets a smaller p-value may be appropriate (0.01 was used for the GUSTO I dataset for this reason).

(7) While AIC is often use for model selection, and is better supported by the literature, BIC may be a valid alternative in larger datasets. For BIC model selection the chi-squared must exceed log(n), thus it will result in smaller models in larger datasets. (Mallow’s may have similar characteristics)

(8) But if you just want a max of 10 or 12 variables, the easier solution is something like bestglm or leaps packages were you just set the maximum number of variables you want to consider.

(9) if you just want a test that will make the two models look the same, and aren't too worried about the details, you could likely compare the AUC of the two models. Some packages will even give you a p-value for the comparison. Doesn't seem advisable.

Ambler G (2002) Simplifying a prognostic model: a simulation study based on clinical data
Cook N.R.; Use and misuse of the ROC curve in the medical literature. Circulation. 115 2007:928-935.
Gail M.H., Pfeiffer R.M.; On criteria for evaluating models of absolute risk. Biostat. 6 2005:227-239.

(10) Once the model has been build, c-statistics/decimation indexes may not be the best approach to comparing models and have well documented limitations. Comparisons should likely also at the minimum include calibration, reclassification index.

Steyerber (2010) Assessing the performance of prediction models: a framework for some traditional and novel measures

(11) It may be a good idea to go beyond above and use decision analytic measures.

Vickers AJ, Elkin EB. Decision curve analysis: a novel method for evaluating prediction models. Med Decis Making. 2006;26:565-74.
Baker SG, Cook NR, Vickers A, Kramer BS. Using relative utility curves to evaluate risk prediction. J R Stat Soc A. 2009;172:729-48.
Van Calster B, Vickers AJ, Pencina MJ, Baker SG, Timmerman D, Steyerberg EW. Evaluation of Markers and Risk Prediction Models: Overview of Relationships between NRI and Decision-Analytic Measures. Med Decis Making. 2013;33:490-501

---Update--- I find the Vickers article the most interesting. But this still hasn't been widely accepted despite many editorials. So may not be of much practical use. The Cook and Steyerberg articles are much more practical.

No one likes stepwise selection. I'm certainly not going to advocate for it. I might emphasize that most of the criticisms of stepwise assumes EPV<50 and a choice between a full or pre-specified model and a reduced model. If EPV>50 and there is a commitment to a reduce model the cost-benefit analysis may be different.

The weak thought behind comparing c-statistics is that they may not be different and I seem to remember this test being significantly underpowered. But now I can't find the reference, so might be way off base on that.

  • $\begingroup$ (1) I am aware that full models are preferred, but I have over 1k vars to choose from and I am required to build these smaller models due to industry-specific requirements. (2) That makes sense! (3) Agreed! (4) True (5) Interesting $\endgroup$ – Matt Reichenbach Dec 3 '13 at 13:40
  • $\begingroup$ (6) Agreed; however, stepwise procedures are very questionable as is, and the lower the p value cutoff, the more biased these types of models become, regardless of sample size. (7) “ For BIC model selection the chi-squared must exceed log(n)”, this seems very useful. Thanks! (8) bestglm and leaps packages are very computationally expensive and take days to run with datasets like the ones that I work with, but thank you for the potential ideas. $\endgroup$ – Matt Reichenbach Dec 3 '13 at 13:40
  • $\begingroup$ (9) These p-values would be significant even if the models were almost exactly the same due to the large sample size alone. (10) I need to research calibration and reclassification index more, thanks! (11) I am very interested to read these articles, would you recommend that I start with Vickers? Thanks! $\endgroup$ – Matt Reichenbach Dec 3 '13 at 13:41

One option is to use pseudo R-square measures for both models. A strong difference in pseudo R-square would suggest that the model fit strongly decreases by omitting V17.

There are different kinds of Pseudo R-squares available. An overview can be found here, for example:


A popular measure is Nagelkerke R-square. It varies between 0 and 1 and, with care, can be interpreted like R-squared from a simple linear regression model. It is based on a transformed ratio of estimated likelihoods of the full model to the intercept-only model.

You could estimate it for fit and fit2, respectively, and compare the relative size to get an indication on your problem. A substantially higher Nagelkerke R-square for fit would suggest that fit2 looses a lot of predictive power by omission of V17.

In lrm the stats value provides Nagelkerke's R-squared. So giving fit$stats should provide you with an estimate. See also ?lrm.

  • $\begingroup$ I am familiar with Nagelkerke's R-square; however, my question lies in what is "a substantially higher Nagelkerke R-Square for fit"? In the above example, there is a difference of 0.001 as we have Nagelkerke's R-square of 0.173 and 0.174 for fit and fit2, respectively. Do you have any references as to what "a substantially higher Nagelkerke R-Square" is? Thanks! $\endgroup$ – Matt Reichenbach Dec 2 '13 at 15:46
  • $\begingroup$ @Matt: I think there are no generalized guidelines on interpreting Nagelkerke's R² or other pseudo R² measures. However, note that it is a transformed measure of 'likelihood reduction' by inclusion of covariates to the intercept only model, which makes it similar to 'explained variance' indicated by standard R² in linear regression. In that sense, I would interpret the .173/.174 difference as very small. A stronger difference would be sth. bridging deciles. However, I suggest to estimate other pseudo R² measures, e.g. McFadden or Cox/Snell to check robustness of this conclusion. $\endgroup$ – tomka Dec 2 '13 at 20:57
  • $\begingroup$ I agree that is difference is very small, but I wish I could find a reference that states what a "small" difference is... I do appreciate your thoughts. Thanks again! $\endgroup$ – Matt Reichenbach Dec 3 '13 at 12:47
  • 1
    $\begingroup$ no problem! sorry for not upvoting you sooner! I will post again, if I find an answer as to what a "small" difference is in terms of pseduo R-squared! thanks! $\endgroup$ – Matt Reichenbach Dec 3 '13 at 18:42

I just read about this. The proper way to do this is use R's glm's finalmodel output and look for "Residual deviance: " and derive the delta between the two models and use this value in a chi-squared test using df equal to the # of predictor terms dropped. And that is your p value.

Applied Regression Modeling Iaian Pardoe 2nd edition 2012 pg 270


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.