Pearson-type goodness-of-ﬁt test with bootstrap maximum likelihood estimation

: The Pearson test statistic is constructed by partitioning the data into bins and computing the diﬀerence between the observed and expected counts in these bins. If the maximum likelihood estimator (MLE) of the original data is used, the statistic generally does not follow a chi-squared distribution or any explicit distribution. We propose a bootstrap- based modiﬁcation of the Pearson test statistic to recover the chi-squared distribution. We compute the observed and expected counts in the par- titioned bins by using the MLE obtained from a bootstrap sample. This bootstrap-sample MLE adjusts exactly the right amount of randomness to the test statistic, and recovers the chi-squared distribution. The bootstrap chi-squared test is easy to implement, as it only requires ﬁtting exactly the same model to the bootstrap data to obtain the corresponding MLE, and then constructs the bin counts based on the original data. We examine the test size and power of the new model diagnostic procedure using simulation studies and illustrate it with a real data set.


Introduction
The model goodness-of-fit test is an important component of model fitting, because model misspecification may cause severe bias and even lead to incorrect inference. The classical Pearson chi-squared test can be traced back to the pioneering work of Pearson (1900). Since then, various model selection and diagnostic tests have been proposed in the literature (Claeskens and Hjort, 2008). In contrast to model selection which concerns multiple models under consideration and eventually selects the best fitting model among them, model diagnostic tests are constructed for a single model, and the goal is to examine whether the model fits the data adequately. The commonly used criterion-based model selection procedures include the Akaike information criterion (AIC) by Akaike (1973) and Bayesian information criterion (BIC), which, however, cannot be used for testing the fit of a single model. For model diagnostics, a common practice is to plot the model residuals versus the predictive outcomes. If the model fits the data adequately, we expect the residuals would be fluctuating around the zero axis, which can thus be used as a graphical checking tool for model misspecification. More sophisticated statistical tests may be constructed based on the partial or cumulative sum of residuals (for example, see Su and Wei, 1991;Stute, Manteiga and Quindimilm, 1998;and Stute and Zhu, 2002).
The classical Pearson chi-squared test statistic is constructed by computing the expected and observed counts in the partitioned bins (Pearson, 1900). More specifically, let (y 1 , . . . , y n ) denote a random sample from the distribution F β0 (y), where β 0 is the true parameter characterizing the distribution function. We are interested in examining whether the sample is from F β0 (y); that is, the null hypothesis is H 0 : F β0 (y) is the true distribution for the observed data, and the alternative is H 1 : F β0 (y) is not the true distribution for the observed data. In the Pearson test, we first partition the sample space into K nonoverlapping bins, and let p k denote the probability assigned to bin k, for k = 1, . . . , K. When the true parameter value β 0 is known, we can easily count the number of observations falling into each prespecified bin. We denote the observed count for bin k as m k . The Pearson goodness-of-fit test statistic takes the form of which asymptotically follows the χ 2 (K−1) distribution under the null hypothesis. We may replace the expected counts np k in the denominator of (1.1) by the observed counts m k , which is asymptotically equivalent to Q 1 (β 0 ), and also follows the χ 2 (K−1) distribution.
However, the true parameter β 0 is often unknown in practice. As a consequence, we need to estimate β 0 in order to construct the bin probabilities or bin counts. For non-regression settings with independent and identically distributed (i.i.d.) data, Chernoff and Lehmann (1954) showed that using the maximum likelihood estimator (MLE) of β 0 based on the original data,β, the test statistic does not follow a χ 2 distribution or any explicit known distribution. In particular, we denote the corresponding estimates for the bin probabilities by p k (β), and define Generally speaking, Q(β) does not follow a χ 2 distribution asymptotically, but it stochastically lies between two χ 2 distributions with different degrees of freedom. This feature of the Pearson-type χ 2 test weakens its generality and limits its applicability to a variety of regression models for which the maximum likelihood estimation procedure dominates. Although some numerical procedures can be used to approximate the null distribution, but they are typically quite computationally intensive (e.g., see Imhof, 1961;and Ali, 1984). If we apply the maximum likelihood estimation to the grouped data and denote the corresponding MLE asβ g , then asymptotically follows a χ 2 K−r−1 distribution with r indicating the dimensionality of β. More recently, Johnson (2004) took a Bayesian approach to constructing a χ 2 test statistic in the form of whereβ is a sample from the posterior distribution of β. In the Bayesian χ 2 test, the partition is constructed as follows. We prespecify 0 ≡ s 0 < s 1 < · · · < s K ≡ 1, and let p k = s k − s k−1 , and let m k (β) be the count of y i 's satisfying Fβ(y i ) ∈ [s k−1 , s k ), for i = 1, . . . , n. Johnson (2004) showed that Q Bayes (β) is asymptotically distributed as χ 2 (K−1) regardless of the dimensionality of β. Intuitively, by generating a posterior sampleβ, Q Bayes (β) recovers the χ 2 distribution and the degrees of freedom that are lost due to computing the MLE of β. However, the Bayesian χ 2 test requires implementation of the usual Monte Carlo Markov chain (MCMC) procedure, which is computationally intensive and also depends on the prior distribution of β. In particular, the prior distribution on β must be noninformative. A major class of noninformative prior distributions are improper priors, which, however, may lead to improper posteriors. If some informative prior distribution is used for β, the asymptotic χ 2 distribution of Q Bayes (β) may be distorted, i.e., Q Bayes (β) is sensitive to the prior distribution of β. In addition, the Pearson-type statistic is largely based on the frequentist maximum likelihood approach, and thus combining a Bayesian posterior sample with the Pearson test is not natural. As a result, Q Bayes (β) cannot be generally used in the classical maximum likelihood framework. Johnson (2007) further developed Bayesian model assessment using pivotal quantities along the similar direction in the Bayesian paradigm.
Our goal is to overcome the dependence of Q Bayes (β) on the prior distribution and further expand the Pearson-type goodness-of-fit test to regression models in the classical maximum likelihood paradigm. We propose a bootstrap χ 2 test to evaluate model fitting, which is easy to implement, and does not require tedious computations other than calculating the MLE of the model parameter by fitting exactly the same model to a bootstrap sample of the original data. The new test statistic maintains the elegance of the Pearson-type formulation, as the right amount of randomness is produced as a whole set through a bootstrap sample to recover the classical χ 2 test. The proposed bootstrap χ 2 test does not require intensive MCMC sampling, and also it is more objective because it does not depend on any prior distribution. Moreover, it is more natural to combine the bootstrap procedure with the classical maximum likelihood estimation in the Pearson test, in contrast to using a posterior sample in the Bayesian paradigm.
The rest of this article is organized as follows. In Section 2, we propose the bootstrap χ 2 goodness-of-fit test using the MLE of the model parameter obtained from a bootstrap sample of the data, and derive the asymptotic distribution for the test statistic. In Section 3, we conduct simulation studies to examine the bootstrap χ 2 test in terms of the test size and statistical power, and also illustrate the proposed method using a real data example. Section 4 gives concluding remarks, and technical details are outlined in the appendix.

Pearson χ 2 test with bootstrap
Let (y i , Z i ) denote the i.i.d. data for i = 1, . . . , n, where y i is the outcome of interest and Z i is the r-dimensional covariate vector for subject i. For ease of exposition, we take the generalized linear model (GLM) to characterize the association between y i and Z i (McCullagh and Nelder, 1989). It is well known that GLMs are suitable for modeling a broad range of data structures, including both continuous and categorical data (e.g., binary or Poisson count data). We assume that the density function of y i is from an exponential family in the form of where θ i is a location parameter, φ is a scalar dispersion parameter, and a i (·), b(·) and c(·) are known functions. The linear predictor and b ′′ (·) represent the first and second derivatives, respectively.
We are interested in testing whether the model in (2.1) fit the observed data adequately. We illustrate the bootstrap χ 2 test under the GLM framework as follows. We first take a simple random sample with replacement from the observed data {(y i , Z i ), i = 1, . . . , n}, and denote the bootstrap sample as We then fit the original regression model to the bootstrap sample and obtain the MLE of β, denoted as β * . We partition the range . . , n} and β * , we then compute the Pearson-type bin counts for each partition. Let m k (β * ) denote the number of subjects satisfying and then we define The proposed bootstrap chi-squared statistic Q Boot (β * ) has the following asymptotic property.
Theorem 2.1. Under the regularity conditions in the appendix, Q Boot (β * ) asymptotically converges to a chi-squared distribution with K − 1 degrees of freedom, χ 2 (K−1) , under the null hypothesis. We outline the key steps of the proof in the appendix. For continuous distributions, m k (β * ) in (2.2) can be obtained in a straightforward way. However, if the data are from a discrete distribution, the corresponding distribution function F (·) is a step function. In this case, we replace the step function with a piecewise linear function that connects the jump points, and redefine F β * (y i |Z i ) to be a uniform distribution between the two adjacent endpoints of the line segment. In particular, for binary data we define where β * is the MLE for a bootstrap sample under the logistic regression. If y i = 0, then we take F β * (y i |Z i ) to be a uniform draw from (0, π * i ); and if y i = 1, we take F β * (y i |Z i ) to be a uniform draw from (π * i , 1). In the Poisson regression, for each given subject with (y i , Z i ), we can calculate the Poisson mean µ * i = exp(β * T Z i ) based on the bootstrap sample MLE β * . We then take In the proposed goodness-of-fit test, the MLE of β needs to be calculated only once based on one bootstrap sample, and thus computation is not heavier than the classical Pearson chi-squared test. On the other hand, the test result depends on one particular bootstrap sample, which can be different for different bootstrap samples. Ideally, we may eliminate the randomness by calculating E{Q Boot (β * b )|data}, where the expectation is taken over all the bootstrap samples conditional on the original data. In practice, we may take a large number of bootstrap samples, and for each of them we construct a chi-squared test statistic. Although these chi-squared values are correlated, the averaged chi-squared test statistic may provide an approximation to E{Q Boot (β * b )|data}. In terms of empirical distribution functions, Durbin (1973) and Stephens (1978) studied the half-sample method and random substitution for goodness-offit tests for distributional assumptions. In particular, using the randomly chosen half of the samples without replacement, the same distribution can be obtained as if the true parameters are known. Nevertheless, our bootstrap procedure not only examines the distributional assumptions, but it also checks the mean structure of the model.

Simulations
We carried out simulation studies to examine the finite sample properties of the proposed bootstrap χ 2 goodness-of-fit test. We focused on the GLMs by simulating data from the linear model, the Poisson regression model, and the logistic model, respectively. We took the number of partitions K = 5 and the sample sizes n = 50, 100, and 200. For each model, we independently generated two covariates: the first covariate Z 1 was a continuous variable from the standard normal distribution and the second Z 2 was a Bernoulli variable taking a value of 0 or 1 with an equal probability of 0.5. We set the intercept β 0 = 0.2, and the two slopes corresponding to Z 1 and Z 2 , β 1 = 0.5 and β 2 = −0.5. Under the linear regression model, we simulated the error term from a normal distribution with mean zero and variance 0.01 under the null hypothesis. The Poisson log-linear regression model took the form of where Z 1 and Z 2 were generated in the same way as those in the linear model. The logistic model assumed the success probability p in the form of and all the rest of setups are the same as before. We conducted 1,000 simulations under each configuration. The simulation results evaluating the test levels are summarized in Table 1. We can see that for each of the five prespecified significance levels of α = 0.01 up to 0.5, the bootstrap χ 2 test clearly maintains the type I error rate under each model. As the sample size increases, the test sizes become closer to the corresponding nominal levels. Figure 1 exhibits the quantile-quantile (Q-Q) plots under each modeling structure with n = 100. Clearly, the proposed bootstrap χ 2 test recovers the χ 2 distribution, as all of the Q-Q plots using the MLE from a bootstrap sample closely match the straight diagonal lines. This demonstrates that the proposed bootstrap χ 2 test performed well with finite sample sizes. We also computed the classical Pearson test statistic when using the MLE calculated from the original data. The corresponding Q-Q plots are presented in Figure 1 as well. The Pearson test statistics using the original data MLE are lower than the expected χ 2 (4) quantiles. This confirms the findings by Chernoff and Lehmann (1954) and also extends their conclusions for the i.i.d. case to the general regression models.
We further examined the power of the proposed bootstrap χ 2 test by simulating data from the alternative hypothesis. Under the linear model, we simulated the error terms from a student t (2) distribution with two degrees of freedom, i.e., ǫ ∼ t (2) in model (3.1). The covariates were generated similarly to those in the null case. We took the number of partitions K = 5, the sample size n = 150, and conducted 1,000 simulations. Under the linear model with the t (2) error, the power of our χ 2 test was 0.893. In another simulation with the linear model, we generated data from an alternative model with an extra quadratic term of covariate Z 1 , that is, y = β 0 + β 1 Z 1 + β 2 Z 2 + γZ 2 1 + ǫ, while the null model is still given by (3.1). The power of the proposed χ 2 test was 0.817 for γ = 0.15, and 0.940 for γ = 0.2.
Similarly, for the Poisson regression model, we added an extra quadratic term in the Poisson mean function under the alternative model, that is, µ = exp(β 0 + β 1 Z 1 + β 2 Z 2 + γZ 2 1 ). If γ = 0.5, the power of our χ 2 test was 0.829, and if γ = 0.6, the power increased to 0.962. We also examined the case where the alternative model was from a negative binomial distribution but with the same mean as that of the Poisson mean. In particular, we took the mean of the negative binomial distribution µ = exp(β 0 + β 1 Z 1 + β 2 Z 2 ) and the negative binomial parameter p = r/(r + µ). The probability mass function of the negative binomial distribution is given by which converges to a Poisson distribution (the null model), as r → ∞. When r = 0.7, the power of our chi-squared test was 0.838, and when r = 0.8, the corresponding power was 0.783. Finally, we examined the test power for the logistic regression model using the proposed χ 2 test. Under the alternative hypothesis, we added a quadratic term in the logistic model, where Z 1 was simulated from a uniform distribution on (1, 2) and Z 2 was still a binary covariate. As γ = 0 corresponded to the null model, we took γ = 0.4 to yield a power of 0.897 for our test, and γ = 0.5 to have a power of 0.976. We also examined a different modeling structure by taking p = Ψ(β 0 + β 1 Z 1 + β 2 Z 2 ), where Ψ(·) is the cumulative distribution function of an exponential distribution. Under this alternative, the power of the proposed test was 0.929. Our test uses the bootstrap data MLE, which recovers the chi-squared distribution. In contrast, the Chernoff and Lehmann test statistic does not follow any explicit distribution.

Application
As an illustration, we applied the proposed goodness-of-fit test to a well-known steam data set described in Draper and Smith (1998). The steam study contained n = 25 observations measured at intervals from a steam plant. The outcome variable was the monthly use of steam, and the covariates of interest included the operating days per month and the average atmospheric temperature. The steam data set was analyzed using a linear regression model, which involved three unknown regression parameters and the variance of the errors. The linear model was claimed to be of adequate fit based on the plot of residuals versus the predicted outcomes. This was also confirmed by the Durbin-Watson test (Draper and Smith, 1998). To quantify the model fit in a more objective way, we applied the proposed bootstrap χ 2 test to examine how well the linear model fit the data from the steam study. Because the sample size was quite small, we partitioned the range of [0, 1] into 3 or 4 intervals, i.e., K = 3 or 4. We took 10,000 bootstrap samples from the original data, and for each of them, we computed the MLEs of the model parameters. Based on these MLEs, we constructed our χ 2 (K−1) test statistics by plugging the bootstrap sample MLEs in the Pearson-type statistic. In Figure 2, we show the histograms of the proposed χ 2 (2) and χ 2 (3) statistics for K = 3 and 4, respectively. We can see that among 10,000 bootstrap χ 2 test statistics only 2.92% of the test statistics exceed the critical value at the significance level of α = 0.05 for K = 3, while 2.75% for K = 4. Our findings provided strong evidence for the model fit, and thus confirmed that the linear regression model adequately fit the steam data.

Discussion
We have proposed a bootstrap-based modification to the classical Pearson χ 2 goodness-of-fit test for regression models, which is a major extension of the work of Chernoff and Lehmann (1954) and Johnson (2004). The new procedure replaces the classical MLE from the original data by the MLE from a bootstrap sample. Using the MLE of a bootstrap sample adjusts the right amount of randomness to the test statistic. Not only does the proposed method restore the degrees of freedom, but also the χ 2 distribution itself, which would have been a nonstandard distribution lying between two χ 2 distributions with different degrees of freedom. Our simulation studies have shown that the proposed test statistic performs well with small sample sizes, and increasingly so as the sample size increases.
Compared with the well-known Akaike information criterion (AIC) and Bayesian information criterion (BIC), we may use the averaged value of the chi-squared statistics computed from a large number of bootstrap samples for model selection or comparison. A smaller value of the averaged chi-squared statistic indicates a better fitting model. It is worth noting that there is no scale associated with the AIC and BIC statistics, thus they are not meaningful alone. In other words, the AIC and the BIC by themselves do not provide any information on the goodness-of-fit of a single model, and they are only interpretable when comparing two or more competing models. In contrast, not only can our averaged bootstrap χ 2 statistic be used for model comparison or model selection, but also it is closely related to the χ 2 distribution, and as an approximation, one would know how well a model fits the data based on the corresponding χ 2 distribution. That is, the proposed test can be used for both model diagnostic and model selection at the same time. For example, a very large value of the averaged χ 2 K−1 value for a small K may shed doubt on the model fit. For the i.i.d. data, the minimum χ 2 statistic estimates the unknown parameter β by minimizing the χ 2 statistic or maximizing the grouped-data likelihood (Cramér, 1946). The minimum χ 2 statistic may not be directly applicable in regression settings due to difficulties involved in grouping the data with regression models. Also, it is challenging to generalize the proposed bootstrap Pearson-type statistic to censored data with commonly used semiparametric Cox proportional hazards model in survival analysis (Cox, 1972;Akritas, 1988;and Akritas and Torbeyns, 1997). Future research is warranted along these directions.

Appendix A: Proof of Theorem 2.1
We assume the conditions (a)-(d) in Cramér (1946, pp. 426-427), and the regularity conditions in Chernoff and Lehmann (1954, p. 581). The conditions in Cramér (1946) are sufficient to prove the χ 2 distribution when using the grouped data MLE. We essentially require the likelihood to be a smooth function of the parameter, the information in the sample increases with the sample size, and the third-order (partial) derivatives of the density function exist. Let β 0 be the true value of the parameter β, letβ be the MLE of β based on the original observations, and let β * be the MLE of β based on the bootstrap sample. Denotê We have that If we follow the notation of (5) in Chernoff and Lehmann (1954), The remaining term is of the same order as the standard deviation of I{s k−1 ≤ Fβ(y i |Z i ) < s k } − I{s k−1 ≤ F β0 (y i |Z i ) < s k }, which takes the value of 0 with probability 1 − O(β − β 0 ), and the value of 1 or −1 with probability O(β − β 0 ) = O p (n −1/2 ). Thus, we can further writê We now show that b k − p k can be approximated by p k −p k , in the classical MLE construction. Note that s k = G(β 0 , β 0 , s k ) = G(β,β, s k ). Denoting G 1 (α, γ, s) = ∂G(α, γ, s) ∂α , G 2 (α, γ, s) = ∂ 2 G(α, γ, s) ∂α∂γ T , we have that = G(β, β 0 , s k ) − G(β, β 0 , s k−1 ) − G(β 0 , β 0 , s k ) + G(β 0 , β 0 , s k−1 ) − G(β,β, s k ) + G(β,β, s k−1 ) + G(β 0 ,β, s k ) − G(β 0 ,β, s k−1 ) wherev k is defined in (7) of Chernoff and Lehmann (1954). We now consider the first term in (A.1). Following the bootstrap principle, the conditional distribution of this term should be the same as that of the second one. We show that in fact the two terms are identically distributed as n → ∞, and they are independent. As an intermediate result, we have already established that m k − m k √ n = {G 1 (β 0 , β 0 , s k ) − G 1 (β 0 , β 0 , s k−1 )} T √ n(β − β 0 ) + o p (1).
Following a similar derivation, we have that