Panel data sets have recently been developed in various areas, and many recent studies have analyzed panel, or longitudinal data sets. Often a dichotomous dependent variable occur in survival analysis, biomedical and epidemiological studies that is analyzed by a generalized linear mixed effects model (GLMM). The most common estimation method for the binary panel data may be the maximum likelihood (ML). Many statistical packages provide ML estimates; however, the estimates are computed from numerically approximated likelihood function. For instance, R packages,
The panel regression model with individual specific effects has the following specification:
where
The probit-normal panel regression model assumes that
The likelihood function of the GLMM is often quite complex in that it includes multidimensional integrals; in addition, it is also difficult to maximize the likelihood function directly. Many authors proposed to use an integral approximation method for the computation of likelihood function, and then maximize the approximated likelihood function. Various integral approximation methods, such as Gauss–Hermite quadratures, Markov chain Monte Carlo, have been employed for this purpose. These approximation methods enable the ML estimation. Theoretically, these approximations can be improved to any precision; however, such integral approximation methods are inadequate for high dimensional integrals in practice. It may be difficult to get a good approximation when
Currently, there are many statistical packages to get the ML estimate of a GLMM model. Most of these packages employ the first approach. For instance,
This paper demonstrate the effect of approximation on the precision of ML estimate and its standard error. For this purpose, we first calculate the ML estimate base on theoretical log-likelihood function, and then compare this ML estimate with the results of
In this section, we wish to discuss the computation of ML estimate in the random effects panel probit model. Because the theory of ML in the GLMM is quite standard, most of this section may be well known, but we include this section for completeness.
Let
Each element of
where
To maximize the likelihood, it needs to evaluate multidimensional integrals.
The theoretical ML estimate is the solution of likelihood equations obtained by differentiating the log of (
where tr(
where
Thus, the likelihood equations are the function of conditional expectations which are related to up to the second order moments of a multivariate truncated normal random vector. One may use the EM algorithm given in McCulloch (1994) or that of Chan and Kuk (1997) for another version of the EM algorithm. Even though we may use one of these algorithms, the computation of the second order moments is essential.
To avoid the calculation of moments, many researchers approximate the conditional moments at this stage, and then apply the EM or variations of EM algorithms. This kind of approach includes the SEM (Celeux and Diebolt, 1985), the SAEM (Celeux
There is a history of the moment calculation for a multivariate truncated normal random vector. Using the moment generating function or recurrence relationships, many moment calculation methods have been proposed under various conditions, see Arismendi (2013). Among them, Kan and Robotti (2017) provide a suitable method for our needs. That is, using their algorithm, we can compute
The Hessian matrix
where
and
To compute Var(
Kennedy and Gentle (1980) and Searle
is satisfied, where
To demonstrate the main feature of ML estimates from two differently approximated likelihood functions, we use Union Wage data in
Table 1 shows the ML estimates and their standard errors computed by two different approaches. We borrowed
The score values of theoretical ML are close to zero. Thus, we may consider the reported value as the computationally attainable ML estimate. Note that the two approximation approaches give similar values of estimate and standard error to ML for slope parameters, but give different estimates of variance component. It seems that the main difference between two approaches occurs in the estimation of the variance component. If we look at the value of the score calculated by (
To confirm the conjecture based on an example, we perform a simulation study under the design used in Harris
The first regressor
We computed the differences between the approximated estimate and the ML estimate during the simulation to see how well the Gauss–Hermite and simulated maximum likelihood methods approximate the ML estimate. The first row of Figure 1 is the bean plot of these quantities when
It is observable the estimates of slope parameter given by the Gauss–Hermite method are almost same as the ML method. The Gauss–Hermite method provides better approximation than the simulated maximum likelihood method for slope parameter; however, the approximation of the latter method is also good in that the maximum bias is less than 0.02, which may be negligible in practice. Note also the simulated maximum likelihood method is better for the variance component. The Gauss–Hermite method has a tendency to give inflated estimate of variance component. It seems that this bias would be more critical. Thus, it may be concluded that the simulated maximum likelihood method is preferable to the Gauss–Hermit method in this study. The second row of Figure 1 is presented to see the approximation of the standard error of estimate. Both methods provide suitable standard error for the slope parameter, but the Gauss–Hermit method is unable to compute the proper standard error of variance component. It usually gives large value of estimate and standard error for the variance component. We also examined the case,
We wish to examine the behavior of two approaches, when
As before, two approaches behave very similarly in the estimation of slope parameters for all cases. In particular, they are essentially same as the ML, when
The recently developed algorithm for the high order moments of a truncated multivariate normal distribution enables the exact computation of the ML estimate in various generalized linear mixed effects models such as the random-effects panel probit model; however, it is difficult to use the algorithm because of computational burden. Most statistical packages employ numerical integral methods, rather than the exact method, to approximate the ML estimate. Theoretically the approximation based on a numerical method can be improved to any precision, but it is also known that numerical integral approximation is inadequate for high dimensional integration in practice. As a result, different packages may give quite different outputs. Thus, we examine the adequacy of two most popular approximation methods, the Gauss–Hermite quadrature and the simulated maximum likelihood method, in the random-effects panel probit model. Based on a simulation study, we found that both methods have the ability of approximating the exact ML estimate and standard error quite well for regression coefficients when the time period is moderate. However, unlike the simulated maximum likelihood, the Gauss–Hermite quadrature does not adequately approximate both quantities for variance component. It has a tendency to give inflated estimates and standard errors of variance component. The Gauss–Hermite quadrature method shows slightly better performance than the other for the estimation of regression coefficient; however, we recommend the simulated maximum likelihood method with a sufficiently large number of Halton draws for practice.
This work was supported by Hanshin University research grant.
Estimates of UnionWage data
Gauss–Hermite ( |
Simulated ML (40 and 1600 Halton draws) | Theoretical ML | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
Estimate | Std.err | Score | Estimate | Std.err | Score | Estimate | Std.err | Score | Estimate | Std.err | Score | |
(Intercept) | −0.9210 | 0.2374 | −0.2646 | −0.9211 | 0.2381 | −0.1730 | −0.9210 | 0.2374 | −0.0014 | −0.9210 | 0.2374 | 2.29e−05 |
wage | 0.5272 | 0.1363 | −0.2995 | 0.5288 | 0.1367 | −0.2892 | 0.5272 | 0.1363 | −0.0030 | 0.5272 | 0.1363 | 1.15e−04 |
expr | −0.0409 | 0.0297 | −2.2501 | −0.0414 | 0.0297 | −0.7111 | −0.0409 | 0.0297 | 0.0026 | −0.0409 | 0.0297 | 8.31e−06 |
rural | 0.2333 | 0.1468 | 0.0378 | 0.2339 | 0.1472 | −0.0082 | 0.2333 | 0.1468 | −0.0015 | 0.2333 | 0.1468 | 3.27e−05 |
sigma | 0.1395 | 0.3811 | −3.4369 | 0.1151 | 0.2085 | −1.2273 | 0.0983 | 0.2712 | 0.0284 | 0.0986 | 0.2695 | −1.07e−04 |
Summary of simulation, (
Gauss–Hermite | Simulated ML | Maximum Likelihood | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Ave. Est | MSE | Ave. Std | Cover. | Ave. Est | MSE | Ave. Std | Cover. | Ave. Est | MSE | Ave. Std | Cover. | |||
30 | 5 | 1.5563 | 0.2121 | 0.4503 | 0.952 | 1.5564 | 0.2121 | 0.4504 | 0.952 | 1.5563 | 0.2121 | 0.4503 | 0.952 | |
−1.0579 | 0.0730 | 0.2457 | 0.954 | −1.0578 | 0.0730 | 0.2457 | 0.954 | −1.0578 | 0.0729 | 0.2457 | 0.954 | |||
1.0396 | 0.3136 | 0.5040 | 0.940 | 1.0397 | 0.3139 | 0.5041 | 0.940 | 1.0395 | 0.3136 | 0.5040 | 0.940 | |||
1.3356 | 0.3715 | 0.4190 | 0.813 | 0.9621 | 0.0972 | 0.2961 | 0.986 | 0.9621 | 0.0972 | 0.2963 | 0.986 | |||
10 | 1.5187 | 0.1461 | 0.3847 | 0.948 | 1.5188 | 0.1461 | 0.3847 | 0.948 | ||||||
−1.0148 | 0.0321 | 0.1769 | 0.940 | −1.0148 | 0.0321 | 0.1769 | 0.940 | |||||||
1.0186 | 0.1926 | 0.4171 | 0.932 | 1.0185 | 0.1925 | 0.4172 | 0.932 | |||||||
1.3554 | 0.2138 | 0.2891 | 0.714 | 0.9584 | 0.0454 | 0.2044 | 0.948 | |||||||
50 | 5 | 1.5256 | 0.1212 | 0.3443 | 0.940 | 1.5251 | 0.1212 | 0.3443 | 0.938 | 1.5256 | 0.1212 | 0.3443 | 0.940 | |
−1.0250 | 0.0304 | 0.1823 | 0.972 | −1.0248 | 0.0304 | 0.1822 | 0.974 | −1.0249 | 0.0304 | 0.1823 | 0.972 | |||
1.0565 | 0.1724 | 0.3863 | 0.936 | 1.0574 | 0.1725 | 0.3863 | 0.936 | 1.0565 | 0.1724 | 0.3863 | 0.936 | |||
1.3670 | 0.3463 | 0.3248 | 0.691 | 0.9896 | 0.0609 | 0.2297 | 0.980 | 0.9896 | 0.0609 | 0.2297 | 0.980 | |||
10 | 1.5036 | 0.1027 | 0.3014 | 0.936 | 1.5037 | 0.1027 | 0.3014 | 0.936 | ||||||
−1.0075 | 0.0203 | 0.1433 | 0.960 | −1.0074 | 0.0203 | 0.1433 | 0.960 | |||||||
1.0063 | 0.1118 | 0.3255 | 0.940 | 1.0061 | 0.1117 | 0.3255 | 0.940 | |||||||
1.3805 | 0.1946 | 0.2244 | 0.516 | 0.9762 | 0.0255 | 0.1587 | 0.964 |