Assessing the goodness of fit of logistic regression models in large samples: A modification of the Hosmer-Lemeshow test

Evaluating the goodness of fit of logistic regression models is crucial to ensure the accuracy of the estimated probabilities. Unfortunately, such evaluation is problematic in large samples. Because the power of traditional goodness of fit tests increases with the sample size, practically irrelevant discrepancies between estimated and true probabilities are increasingly likely to cause the rejection of the hypothesis of perfect fit in larger and larger samples. This phenomenon has been widely documented for popular goodness of fit tests, such as the Hosmer-Lemeshow test. To address this limitation, we propose a modification of the Hosmer-Lemeshow approach. By standardizing the noncentrality parameter that characterizes the alternative distribution of the Hosmer-Lemeshow statistic, we introduce a parameter that measures the goodness of fit of a model but does not depend on the sample size. We provide the methodology to estimate this parameter and construct confidence intervals for it. Finally, we propose a formal statistical test to rigorously assess whether the fit of a model, albeit not perfect, is acceptable for practical purposes. The proposed method is compared in a simulation study with a competing modification of the Hosmer-Lemeshow test, based on repeated subsampling. We provide a step-by-step illustration of our method using a model for postneonatal mortality developed in a large cohort of more than 300 000 observations.     

Tests for Hazard Transformation

The semiparametric Cox proportional hazards model is routinely adopted to model time-to-event data. Proportionality is a strong assumption, especially when follow-up time, or study duration, is long. Zeng and Lin (J. R. Stat. Soc., Ser. B, 69:1–30, 2007) proposed a useful generalisation through a family of transformation models which allow hazard ratios to vary over time. In this paper we explore a variety of tests for the need for transformation, arguing that the Cox model is so ubiquitous that it should be considered as the default model, to be discarded only if there is good evidence against the model assumptions. Since fitting an alternative transformation model is more complicated than fitting the Cox model, especially as procedures are not yet incorporated in standard software, we focus mainly on tests which require a Cox fit only. A score test is derived, and we also consider performance of omnibus goodness-of-fit tests based on Schoenfeld residuals. These tests can be extended to compare different transformation models. In addition we explore the consequences of fitting a misspecified Cox model to data generated under a true transformation model. Data on survival of 1043 leukaemia patients are used for illustration.     

Testing goodness-of-fit in logistic case-control studies

We present a goodness-of-fit test for the logistic regressionmodel under case-control sampling. The test statistic is constructedvia a discrepancy between two competing kernel density estimatorsof the underlying conditional distributions given case-controlstatus. The proposed goodness-of-fit test is shown to comparevery favourably with previously proposed tests for case-controlsampling in terms of power. The test statistic can be easilycomputed as a quadratic form in the residuals from a prospectivelogistic regression maximum likelihood fit. In addition, theproposed test is affine invariant and has an alternative representationin terms of empirical characteristic functions.     

Score Test for Conditional Independence Between Longitudinal Outcome and Time to Event Given the Classes in the Joint Latent Class Model

Summary .  Latent class models have been recently developed for the joint analysis of a longitudinal quantitative outcome and a time to event. These models assume that the population is divided in  G  latent classes characterized by different risk functions for the event, and different profiles of evolution for the markers that are described by a mixed model for each class. However, the key assumption of conditional independence between the marker and the event given the latent classes is difficult to evaluate because the latent classes are not observed. Using a joint model with latent classes and shared random effects, we propose a score test for the null hypothesis of independence between the marker and the outcome given the latent classes versus the alternative hypothesis that the risk of event depends on one or several random effects from the mixed model in addition to the latent classes. A simulation study was performed to compare the behavior of the score test to other previously proposed tests, including situations where the alternative hypothesis or the baseline risk function are misspecified. In all the investigated situations, the score test was the most powerful. The methodology was applied to develop a prognostic model for recurrence of prostate cancer given the evolution of prostate-specific antigen in a cohort of patients treated by radiation therapy.     

Testing inflated zeros in binomial regression models

Binomial regression models are commonly applied to proportion data such as those relating to the mortality and infection rates of diseases. However, it is often the case that the responses may exhibit excessive zeros; in such cases a zero‐inflated binomial (ZIB) regression model can be applied instead. In practice, it is essential to test if there are excessive zeros in the outcome to help choose an appropriate model. The binomial models can yield biased inference if there are excessive zeros, while ZIB models may be unnecessarily complex and hard to interpret, and even face convergence issues, if there are no excessive zeros. In this paper, we develop a new test for testing zero inflation in binomial regression models by directly comparing the amount of observed zeros with what would be expected under the binomial regression model. A closed form of the test statistic, as well as the asymptotic properties of the test, is derived based on estimating equations. Our systematic simulation studies show that the new test performs very well in most cases, and outperforms the classical Wald, likelihood ratio, and score tests, especially in controlling type I errors. Two real data examples are also included for illustrative purpose.     

Goodness of Fit Tests for Proportional Hazards Regression Models

We present a test of goodness of fit for the proportional hazard regression model. The test is based on a score statistic for testing against local mixture alternatives. Contrary to the findings of several other authors, we detect a significant lack of fit in Freireich's leukemia data.     

Use of the score test as a goodness-of-fit measure of the covariance structure in genetic analysis of longitudinal data

Model selection is an essential issue in longitudinal data analysis since many different models have been proposed to fit the covariance structure. The likelihood criterion is commonly used and allows to compare the fit of alternative models. Its value does not reflect, however, the potential improvement that can still be reached in fitting the data unless a reference model with the actual covariance structure is available. The score test approach does not require the knowledge of a reference model, and the score statistic has a meaningful interpretation in itself as a goodness-of-fit measure. The aim of this paper was to show how the score statistic may be separated into the genetic and environmental parts, which is difficult with the likelihood criterion, and how it can be used to check parametric assumptions made on variance and correlation parameters. Selection of models for genetic analysis was applied to a dairy cattle example for milk production.     

Neighbourhood effects influence drought‐induced mortality of savanna trees in Australia

Questions: The following hypotheses of neighbourhood effects on drought‐induced mortality are evaluated: (A) drought‐induced stem death is randomly distributed in space, (B) stems are predisposed to drought‐induced death through negative density‐dependent effects and (C) stems are predisposed to drought‐induced death due to local deficits in plant available resources. Location: Central Queensland, Australia. Methods: Recent mass mortality of woody stems was surveyed and mapped in three 1.21‐ha quadrats within Eucalyptus melanophloia‐dominated savanna. A multi‐faceted analytical approach was adopted including spatial pattern analyses, two logistic regressions of neighbourhood density effects on survival and spatial autocorrelation analyses of model residuals. Results: Mortality was concentrated in stems ≤15‐cm diameter at breast height (DBH). Survival was aggregated or random in quadrats 1 and 3 and random o regular in quadrat 2. Small neighbour density had a negative effect on survival in all quadrats. In addition, the second model identified a positive relationship between survival and living neighbour density in quadrat 3 (indicating a resource patch effect), but a negative relationship in quadrat 2 (density effect). Analysis of model residuals showed that neighbour density explained mortality equally well across quadrat 2, but not across quadrats 1 and 3. Conclusions: There was evidence in support of hypotheses B (neighbour density) and C (resource heterogeneity). We found strong support for an interaction between microsite quality and neighbourhood stem densities, and suggest that this interaction is driven by plant available water.     

Functional form diagnostics for Cox's proportional hazards model

We propose a new type of residual and an easily computed functional form test for the Cox proportional hazards model. The proposed test is a modification of the omnibus test for testing the overall fit of a parametric regression model, developed by Stute, González Manteiga, and Presedo Quindimil (1998, Journal of the American Statistical Association93, 141-149), and is based on what we call censoring consistent residuals. In addition, we develop residual plots that can be used to identify the correct functional forms of covariates. We compare our test with the functional form test of Lin, Wei, and Ying (1993, Biometrika80, 557-572) in a simulation study. The practical application of the proposed residuals and functional form test is illustrated using both a simulated data set and a real data set.     

Permutation tests for random effects in linear mixed models

Inference regarding the inclusion or exclusion of random effects in linear mixed models is challenging because the variance components are located on the boundary of their parameter space under the usual null hypothesis. As a result, the asymptotic null distribution of the Wald, score, and likelihood ratio tests will not have the typical χ(2) distribution. Although it has been proved that the correct asymptotic distribution is a mixture of χ(2) distributions, the appropriate mixture distribution is rather cumbersome and nonintuitive when the null and alternative hypotheses differ by more than one random effect. As alternatives, we present two permutation tests, one that is based on the best linear unbiased predictors and one that is based on the restricted likelihood ratio test statistic. Both methods involve weighted residuals, with the weights determined by the among- and within-subject variance components. The null permutation distributions of our statistics are computed by permuting the residuals both within and among subjects and are valid both asymptotically and in small samples. We examine the size and power of our tests via simulation under a variety of settings and apply our test to a published data set of chronic myelogenous leukemia patients.     

The Mizon–Richard Encompassing Test for the Cox and Aalen Additive Hazards Models

Summary .   The Cox hazards model ( Cox, 1972 , Journal of the Royal Statistical Society, Series B 34, 187–220) for survival data is routinely used in many applied fields, sometimes, however, with too little emphasis on the fit of the model. A useful alternative to the Cox model is the Aalen additive hazards model ( Aalen, 1980 , in Lecture Notes in Statistics-2 , 1–25) that can easily accommodate time changing covariate effects. It is of interest to decide which of the two models that are most appropriate to apply in a given application. This is a nontrivial problem as these two classes of models are nonnested except only for special cases. In this article we explore the Mizon–Richard encompassing test for this particular problem. It turns out that it corresponds to fitting of the Aalen model to the martingale residuals obtained from the Cox regression analysis. We also consider a variant of this method, which relates to the proportional excess model ( Martinussen and Scheike, 2002 , Biometrika 89, 283–298). Large sample properties of the suggested methods under the two rival models are derived. The finite-sample properties of the proposed procedures are assessed through a simulation study. The methods are further applied to the well-known primary biliary cirrhosis data set.     

One‐Sided Test to Assess Correlation in Linear Logistic Models using Estimating Equations

A score‐type test is proposed for testing the hypothesis of independent binary random variables against positive correlation in linear logistic models with sparse data and cluster specific covariates. The test is developed for univariate and multivariate one‐sided alternatives. The main advantage of using score test is that it requires estimation of the model only under the null hypothesis, that in this case corresponds to the binomial maximum likelihood fit. The score‐type test is developed from a class of estimating equations with block‐diagonal structure in which the coefficients of the linear logistic model are estimated simultaneously with the correlation. The simplicity of the score test is illustrated in two particular examples.     

Spending degrees of freedom in a poor economy: A case study of building a sightability model for moose in northeastern Minnesota

Sightability models are binary logistic-regression models used to estimate and adjust for visibility bias in wildlife-population surveys. Like many models in wildlife and ecology, sightability models are typically developed from small observational datasets with many candidate predictors. Aggressive model-selection methods are often employed to choose a best model for prediction and effect estimation, despite evidence that such methods can lead to overfitting (i.e., selected models may describe random error or noise rather than true predictor–response curves) and poor predictive ability. We used moose (Alces alces) sightability data from northeastern Minnesota (2005–2007) as a case study to illustrate an alternative approach, which we refer to as degrees-of-freedom (df) spending: sample-size guidelines are used to determine an acceptable level of model complexity and then a pre-specified model is fit to the data and used for inference. For comparison, we also constructed sightability models using Akaike's Information Criterion (AIC) step-down procedures and model averaging (based on a small set of models developed using df-spending guidelines). We used bootstrap procedures to mimic the process of model fitting and prediction, and to compute an index of overfitting, expected predictive accuracy, and model-selection uncertainty. The index of overfitting increased 13% when the number of candidate predictors was increased from three to eight and a best model was selected using step-down procedures. Likewise, model-selection uncertainty increased when the number of candidate predictors increased. Model averaging (based on R = 30 models with 1–3 predictors) effectively shrunk regression coefficients toward zero and produced similar estimates of precision to our 3-df pre-specified model. As such, model averaging may help to guard against overfitting when too many predictors are considered (relative to available sample size). The set of candidate models will influence the extent to which coefficients are shrunk toward zero, which has implications for how one might apply model averaging to problems traditionally approached using variable-selection methods. We often recommend the df-spending approach in our consulting work because it is easy to implement and it naturally forces investigators to think carefully about their models and predictors. Nonetheless, similar concepts should apply whether one is fitting 1 model or using multi-model inference. For example, model-building decisions should consider the effective sample size, and potential predictors should be screened (without looking at their relationship to the response) for missing data, narrow distributions, collinearity, potentially overly influential observations, and measurement errors (e.g., via logical error checks). © 2011 The Wildlife Society.     

Estimation of rates-across-sites distributions in phylogenetic substitution models

Previous work has shown that it is often essential to account for the variation in rates at different sites in phylogenetic models in order to avoid phylogenetic artifacts such as long branch attraction. In most current models, the gamma distribution is used for the rates-across-sites distributions and is implemented as an equal-probability discrete gamma. In this article, we introduce discrete distribution estimates with large numbers of equally spaced rate categories allowing us to investigate the appropriateness of the gamma model. With large numbers of rate categories, these discrete estimates are flexible enough to approximate the shape of almost any distribution. Likelihood ratio statistical tests and a nonparametric bootstrap confidence-bound estimation procedure based on the discrete estimates are presented that can be used to test the fit of a parametric family. We applied the methodology to several different protein data sets, and found that although the gamma model often provides a good parametric model for this type of data, rate estimates from an equal-probability discrete gamma model with a small number of categories will tend to underestimate the largest rates. In cases when the gamma model assumption is in doubt, rate estimates coming from the discrete rate distribution estimate with a large number of rate categories provide a robust alternative to gamma estimates. An alternative implementation of the gamma distribution is proposed that, for equal numbers of rate categories, is computationally more efficient during optimization than the standard gamma implementation and can provide more accurate estimates of site rates.     

Fit and Frailties in Proportional Hazards Regression

We exploit a conjectured equivalence between proportional hazards models with frailties and a particular subclass of non proportional hazards models, specifically those with declining effects, to address the question of fit. A goodness of fit test of the proportional hazards assumption against an alternative of declining regression effect is equivalent to a test for the presence of frailties. Such tests are now widely available in standard software. Although a number of tests of the proportional hazards assumption have been developed there is no test that directly formulates the alternative in terms of a non‐specified monotonic decline in regression effect and that enables a quantification of this in terms of a simple index. The index we obtain lies between zero and one such that, for any given set of covariates, values of the index close to one indicate that the fit cannot essentially be improved by allowing the possibility of regression effects to decline. Values closer to zero and away from one indicate that the fit can be improved by relaxing the proportional hazards constraint in this particular direction. (© 2004 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

Tests for Trend in Binary Response

Tests for trend are important in analyzing data where the binary response in ordered categories is of interest. An example is in toxicology where the response in various dose groups is observed. For testing an association between the dose and the response the approach from Cochran and Armitage is widely used. However the result of this test is highly dependent on the scores assigned to the dose groups. Various dose assignments can lead to different outcomes. As an alternative the isotonic regression, a nonparametric method, is proposed. The outcome of this approach is independent of the quantification of the dose. Both methods (Cochran‐Armitage test and isotonic regression) are compared within a simulation study to an isotonic version of the Pearson's Chi‐squared test and the Wilcoxon rank sum test.     

A simple hierarchical approach to modeling distributions of substitution rates

Genetic sequence data typically exhibit variability in substitution rates across sites. In practice, there is often too little variation to fit a different rate for each site in the alignment, but the distribution of rates across sites may not be well modeled using simple parametric families. Mixtures of different distributions can capture more complex patterns of rate variation, but are often parameter-rich and difficult to fit. We present a simple hierarchical model in which a baseline rate distribution, such as a gamma distribution, is discretized into several categories, the quantiles of which are estimated using a discretized beta distribution. Although this approach involves adding only two extra parameters to a standard distribution, a wide range of rate distributions can be captured. Using simulated data, we demonstrate that a "beta-" model can reproduce the moments of the rate distribution more accurately than the distribution used to simulate the data, even when the baseline rate distribution is misspecified. Using hepatitis C virus and mammalian mitochondrial sequences, we show that a beta- model can fit as well or better than a model with multiple discrete rate categories, and compares favorably with a model which fits a separate rate category to each site. We also demonstrate this discretization scheme in the context of codon models specifically aimed at identifying individual sites undergoing adaptive or purifying evolution.     

Monitoring habitat types by the mixed multinomial logit model using panel data

Habitats in the Wadden Sea, a world heritage area, are affected by land subsidence resulting from natural gas extraction and by sea level rise. Here we describe a method to monitor changes in habitat types by producing sequential maps based on point information followed by mapping using a multinomial logit regression model with abiotic variables of which maps are available as predictors.In a 70 ha study area a total of 904 vegetation samples has been collected in seven sampling rounds with an interval of 2–3 years. Half of the vegetation plots was permanent, violating the assumption of independent data in multinomial logistic regression. This paper shows how this dependency can be accounted for by adding a random effect to the multinomial logit (MLN) model, thus becoming a mixed multinomial logit (MMNL) model. In principle all regression coefficients can be taken as random, but in this study only the intercepts are treated as location-specific random variables (random intercepts model). With six habitat types we have five intercepts, so that the number of extra model parameters becomes 15, 5 variances and 10 covariances.The likelihood ratio test showed that the MMNL model fitted significantly better than the MNL model with the same fixed effects. McFadden-R² for the MMNL model was 0.467, versus 0.395 for the MNL model. The estimated coefficients of the MMNL and MNL model were comparable; those of altitude, the most important predictor, differed most. The MMNL model accounts for pseudo-replication at the permanent plots, which explains the larger standard errors of the MMNL coefficients. The habitat type at a given location-year combination was predicted by the habitat type with the largest predicted probability. The series of maps shows local trends in habitat types most likely driven by sea-level rise, soil subsidence, and a restoration project.We conclude that in environmental modeling of categorical variables using panel data, dependency of repeated observations at permanent plots should be accounted for. This will affect the estimated probabilities of the categories, and even stronger the standard errors of the regression coefficients.     

A Generalized Logistic Model for Quantal Response Bioassay

In bioassay, where different levels of the stimulus may represent different doses of a drug, the binary response is the death or survival of an individual receiving a specified dose. In such applications, it is common to model the probability of a positive response P at the stimulus level x by P = F(x′β), where F is a cumulative distribution function and β is a vector of unknown parameters which characterize the response function. The two most popular models used for modelling binary response bioassay involve the probit model [BLISS (1935), FINNEY (1978)], and the logistic model [BERKSON (1944), BROWN (1982)]. However, these models have some limitations. The use of the probit model involves the inverse of the standard normal distribution function, making it rather intractable. The logistic model has a simple form and a closed expression for the inverse distribution function, however, neither the logistic nor the probit can provide a good fit to response functions which are not symmetric or are symmetric but have a steeper or gentler incline in the central probability region. In this paper we introduce a more realistic model for the analysis of quantal response bioassay. The proposed model, which we refer to it as the generalized logistic model, is a family of response curves indexed by shape parameters m₁ and m₂. This family is rich enough to include the probit and logistic models as well as many others as special cases or limiting distributions. In particular, we consider the generalized logistic three parameter model where we assume that m₁ = m, m is a positive real number, and m₂ = 1. We apply this model to various sets of data, comparing the fit results to those obtained previously by other dose-response curves such as the logistic and probit, and showing that the fit can be improved by using the generalized logistic.