Modified generalized method of moments for a robust estimation of polytomous logistic model

The maximum likelihood estimation (MLE) method, typically used for polytomous logistic regression, is prone to bias due to both misclassification in outcome and contamination in the design matrix. Hence, robust estimators are needed. In this study, we propose such a method for nominal response data with continuous covariates. A generalized method of weighted moments (GMWM) approach is developed for dealing with contaminated polytomous response data. In this approach, distances are calculated based on individual sample moments. And Huber weights are applied to those observations with large distances. Mellow-type weights are also used to downplay leverage points. We describe theoretical properties of the proposed approach. Simulations suggest that the GMWM performs very well in correcting contamination-caused biases. An empirical application of the GMWM estimator on data from a survey demonstrates its usefulness.


INTRODUCTION
Polytomous logistic regression models for multinomial data are a powerful technique for relating dependent categorical responses to both categorical and continuous explanatory covariates (McCullagh & Nelder, 1989;Liu & Agresti, 2005). In practice, however, the model building process can be highly influenced by peculiarities in the data. The maximum likelihood estimation (MLE) method, typically used for the polytomous logistic regression model (PLRM), is prone to bias due to both misclassification in outcome and contamination in the design matrix (Pregibon, 1982;Copas, 1988). Hence, robust estimators are needed.
A generalized method of moments (GMM) estimation can be formed as a substitute of MLE. The GMM is particularly useful when the moment conditions are relatively easy to obtain. GMM has been extensively studied in econometrics (Hansen, 1982;Newey & West, 1987;Pakes & Pollard, 1989;Hansen, Heaton & Yaron, 1996;Newey & McFadden, 1994). Under some regularity conditions, the GMM estimator is consistent (Hansen, 1982). With an appropriately chosen weight matrix, GMM achieves the same efficiency as the MLE (Hayashi, 2000). Furthermore, under certain circumstances, GMM provides more flexibility, such as dealing with endogeneity through instrumental variables (Baum, Schaffer & Stillman, 2002).
Like MLE, GMM estimation can be easily corrupted by aberrant observations (Ronchetti & Trojani, 2001). Such observations can bring up disastrous bias on standard parameter estimates if they are not properly accounted for, see Huber (1981), Hampel et al. (2005), and Rousseeuw & Leroy (2003). So we propose a modified estimation method based on an outlier robust variant of GMM. The method is different from the kernel-weighted GMM developed for linear time-series data by Kuersteiner (2012) in that this is a data-driven method for defining weights. The new approach is evaluated using asymptotic theory, simulations, and an empirical example.
The robust GMM estimator is motivated by the data from a 2006 study on hypertension in a sample of the Chinese population. 520 people completed the survey. Observed variables included demographics, social-economic status, weight, height, blood pressure, and food consumption. Sodium intakes were calculated based on overall food consumption. Among those covariates, age, body mass index (BMI), and sodium intakes are all continuous. Based on blood pressure measurements, subjects were classified into 4 categories: Normal, Pre-hypertension, Stage 1 and Stage 2 hypertension. Table 1 lists the summary statistics of the sample. One of the research objectives is to examine the association between hypertension and risk factors in the population. Since the proportional odds assumption is violated (Score test for the proportional odds assumption gives χ 2 = 182.27 with a degree of freedom of 8, p < 0.0001), we apply the polytomous logistic model, using the normal category as the reference level. In the case of J category, the polytomous logit model have J − 1 comparisons. Each comparison have a set of parameters for all covariates in the model. Therefore, the generalized logit model is not parsimonious when comparing with the proportional odds model. But the simultaneous estimation of all parameters is more efficient than separate models for each comparison. It is another option for ordinal response data, especially when a proportional odds model does not fit the data well. Table 2 lists the output from the model estimated by MLE. It is obvious that, if MLE is used, the estimates is inconsistent for sodium intakes, particularly the negative coefficient of sodium intake for the odds between the Stage 2 hypertension and the Normal categories. The inconsistency is more obvious when we plot the odds with respect to the sodium intake, the downward trend of the odds in Fig. 2A. This result contradicts the previous finding that there is a strong relationship between sodium intake and hypertension, see for example National Research Council (2005), He & MacGregor (2004) and references therein. Besides, Fig. 2A also shows another strange situation: the higher starting points  for the odds between the Pre-hypertension and the Normal categories. The scatter plot ( Fig. 1) between distances and leverages suggests some observations are possible outliers: Observations 21,33,85,92,194,274,336,414,459,483,and 489 have large distances,and Observations 37,83,263,459,483,485,and 490 have large leverages, which are red-colored. The paper is set up as follows. In the next section we presents the basic notations, model, and standard GMM. "A robust GMM" introduces the outlier robust GMM  estimator, and gives a detailed exposition of its implementation. In "Results", we compares the performance of the standard MLE with the new estimator using a Monte-Carlo experiment. And we apply both estimators to real epidemiological data, and illustrate the usefulness of the robust estimator for application oriented researchers. We conclude with a discussion of advantages and limitations of the approach. The supporting document gathers the proofs of the asymptotic property.

The baseline-category logit model
Assume a random sample of size n from a large population. Each element in the population may be classified into one of J categories, denoted by y i = (y i1 ,y i2 ,...,y iJ ) the multinomial trial for subject i, where y ij = 1 when the response is in category j and y ij = 0 otherwise, i = 1,...,n, j = 1,...,J. Thus,  j y ij = 1. Suppose p explanatory covariates, with at least one of them being continuous, are observed. Define , denote the probability that the observation of Y belongs to category j, given covariates x i , we assume the relationship between the probability π j and x can be modeled as: where β T j = (β j0 ,β j1 ,...,β jp ). Here we set the first category as reference class. This model is called a baseline-category logit model (Agresti, 2012) or generalized logit model (Stokes, Davis & Koch, 2009). MLE is usually used for obtaining parameter estimation of this model. Here we present an alternative estimation method formed with the GMM.

Estimation using GMM
The baseline-category logit model can be viewed as a multivariate model. Define In the GMM framework, we define where with the corresponding sample moment condition The GMM estimation ofβ M can be obtained by minimizing the following quadratic objective function where Σ n (β) can be the empirical variance-covariance matrix given by Or, for the best efficiency of the GMM estimation, we can take the information matrix of the polytomous logit model (PLRM), that is, In general,β M can be computed via an iterative procedure (Hansen, Heaton & Yaron, 1996). Under standard regularity conditions, the GMM estimatorβ M exists and converges in probability to the true parameter β 0 (Hansen, 1982). A proof of asymptotic normality of GMM can be found on p. 2148 of Newey & McFadden (1994).

A robust GMM
In this section we introduce the outlier robust GMM estimator. In the following subsection, we specify moment conditions used for robust estimation. And the details on the implementation of the estimator follows.

The generalized method of weighted moments
The main principle used in the robust GMM estimator is that we replace moment conditions by a set of observation weighted moment conditions. Instead of Eq. (3), we define where Then the estimation can be based on the moment conditions Consequently, the generalized method of weighted moments (GMWM) estimates can be defined bŷ where Here we take the summation as the sample moment condition. The advantage of using the summation is that it can lead us to a direct estimation of covariance matrix. It is clear to see that this definition is analogous to the standard GMM. If we choose w i = 1 and c i = 0 for all observations, the moment conditions in (6) are reduced to the standard moment conditions. Therefore, the standard GMM is a special case of the GMWM.
In order to specify the weights for the robust GMM estimator, we need the following definition of a distance, which is based on individual moment conditions: The weight is assigned based on There are several alternative specifications of weight functions available in the literature (Huber, 1981;Hampel et al., 2005). In this study, the Huber's weights are applied: The above specification of weight function requires a value of the tuning constant c d . Both the outlier sensitivity and the efficiency of the estimator are determined by the constant. On the one hand, the estimator should be reasonably efficient if the sample contains no outlier. On the other hand, the estimator should be insensitive to outliers. To determine c d , understanding the distribution of d i (β) is critical. Clearly, u w i (β) is a column vector, and d i (β) is a scalar quadratic distance, so we set c d = χ −2 1 (0.975)/n, where χ −2 p (·) is the quantile of the χ 2 distribution with p degrees of freedom.
If we take the information matrix (5) of the PLRM as Σ w n (β), we can compute leverage for each observation: where σ w i is the ith component of Σ w n (β). Then, a Mallows-type weight can be defined based on trace(H i ); that is, w x = w(trace(H i )), to downplay the observations with high leverages. Lesaffre & Albert (1989) suggest that the practical rule for isolating leverage points might set c x = 2(p + 1)(J − 1)/n. In this study, we give observations with large leverages 0 weights, An approach often used to combine the two weights is w i = w d · w x (Heritier et al., 2009 The consistency correction vector c i is defined as

Implementation of the estimator
The continuous updating estimation method is applied in this study for estimating the regression coefficients and corresponding variance. The procedure is detailed as follows: 1. Apply an initial value β (0) for computing Σ n  β  .
6. Continue this procedure until convergence criteria are met.
For the starting value β (0) , a reasonable choice is the MLE estimation based on the original data.
In the appendix, we proved that, under some regularity assumptions, we can have that β w is consistent for β 0 . And by studying the behavior of the weighted moment equations in a neighborhood of β 0 , we showed that the asymptotic linearity ensures the applicability of the central limit theorem for the asymptotic normality of GMWM.

Monte Carlo simulations
In this section we investigate the properties of the GMWM estimator using a Monte-Carlo study. We generate data with three response categories and two covariates which are from multivariate normal distribution with 0 mean and identity covariance. The true coefficient matrix β 0 is Based on the specified coefficient values and using the probability based on the model (1), we compute the category-specific probabilities for each subject. Then, using the computed probabilities, we determine the most likely category to which each subject belongs. This decision is made through random generation from the multinomial distribution with the probability vector as a parameter. For instance, multinomial categories in R-Language is the probability vector, n i is the number of random vectors to draw, and N i is the total number of objects that are put into J-categories. In our case, n i = N i = 1 for all subjects and J = 3. Two sample sizes, 100 and 1000, are examined. For each sample size, we run the simulation 1000 times. Average biases and MSEs are calculated and tabulated. Table 3 shows the results from randomly generated data with no outliers added. When the sample size is small, GMWM will give greater biases on β 30 and β 31 compared to the MLE method. For the sample size 1000, biases on these two parameters increase too, but not so obviously. Variances will also be inflated due to the weights we applied.
Outliers are generated from a multivariate normal distribution with the mean vector = (2,3) and identity covariance I 2 . For these outliers, their responses are intentionally misclassified, that is, they are placed within a different category from those predicted categories based on the true parameters. Table 4 lists simulation results with outliers added. For estimations from datasets with 5% outliers, bias correction from the GMWM is excellent. However, when the datasets have 10% outliers, biases on estimations of some parameters (β 21 and β 22 in this simulation) are decreased, but not completely corrected.

Application
For the hypertension data, the criterion for identifying observations with large distances is c d = 0.22, and the criterion for identifying leverage points is c x = 0.12. Applying the GMWM estimator, those blue-colored points in Fig. 1 are automatically downweighted, and red-colored points have 0 weight. The GMWM method indeed eliminates those inconsistencies: the coefficient of sodium intake for the odds model between the Stage 2 hypertension and the Normal categories is no longer negative, see the right side of Table 2. As the results indicate, age, gender, and BMI all had significant impact on hypertension status. For example, one unit increase in BMI resulted in an increase of 1.26 (95% confidence interval [1.16-1.35]) times in likelihood to have Stage 2 hypertension when compared with the normal status. And with one year age increase, a subject was 1.07 (95% CI [1.06-1.10]) times more likely to have Stage 2 hypertension than to stay at the normal healthy status. Contrary to the MLE results for sodium intakes, which were difficult to make a conclusion due to inconsistent estimate, we now find that sodium intakes were statistically significant. When a daily intake of sodium increased one gram, a subject were 1.26 (95% CI [1.15-1.37]) times more likely to have Stage 1 hypertension, and 1.25 (95% CI [1.17-1.35]) times more likely to have Stage 2 hypertension. These results are consistent with the findings from previous studies (National Research Council, 2005;He & MacGregor, 2004).

DISCUSSION
A reasonable choice to fit ordinal response data is the proportional odds model if the proportional odds assumption is not violated. Proportional odds models can take the ordinal information into modeling. And it reduces the number of parameters which is needed by the generalized logit model. Unfortunately, our data does not met the fundamental assumption of proportional odds models, which makes us choose to treat the outcome as a nominal response.
A datum with a nominal response and some continuous covariates is commonly seen in many scientific areas, such as sociology, economy, and biomedical studies. In order to be able to deal with outliers, we modified the GMM estimator to replace the standard moment conditions with weighted moment conditions, so that aberrant observations automatically receive less weight. We proved that the proposed method has good asymptotic behavior. When outliers are present, the GMWM estimator give much smaller biases than the estimations derived from the traditional MLE method. This method can be adapted to check whether results obtained with the traditional MLE approach are driven only by a few outlying observations. The weights produced from the robust procedure can be used to diagnose the cause of the differences and to indicate routes for model re-specification.

APPENDIX: CONSISTENCY AND ASYMPTOTIC NORMALITY
In this appendix, we introduce the assumptions for the asymptotic analysis of GMWM, and outline the derivations on the main asymptotic properties of GMWM.
We make the following sets of regularity assumptions regarding properties of the moment functions and identification assumptions. Assumption I I1. B is a compact parametric space.
I2. Σ is a positive definite matrix.
When the above assumptions are met, we can prove thatβ w is consistent for β 0 . We begin with studying the behavior of the weighted moment equations in a neighborhood of β 0 . And proving their asymptotic linearity is followed. The linearity ensures the applicability of the central limit theorem for the asymptotic normality of GMWM. Theorem 1. Let the assumptions F and I hold, then the GMWM estimatorβ w is asymptotically normal, that is,