The large sample property of the iterative generalized least squares estimation for hierarchical mixed effects model

Abstract

In many fields, we need to deal with hierarchically structured data. For this kind of data, hierarchical mixed effects model can show the correlation of variables in the same level by establishing a model for regression coefficients. Due to the complexity of the random part in this model, seeking an effective method to estimate the covariance matrix is an appealing issue. Iterative generalized least squares estimation method was proposed by Goldstein in 1986 and was applied in special case of hierarchical model. In this paper, we extend the method to the general hierarchical mixed effects model, derive its expressions in detail and apply it to economic examples.

Graphical abstract

Keywords

Hierarchical model / iterative generalized least squares estimation / variance-covariance components / maximum likelihood estimation

1 Introduction

For the hierarchically structured data, an important feature is the correlation between the data within the same level. For example, if we examine the math scores of students from different classes, it is clear that there should be some correlation between the scores of students in the same class, and this intra-level correlation is something that we must consider in the subsequent analysis. At this point the mutual independence of the random terms assumed by the traditional regression model no longer holds. The hierarchical model reflects the intra-level correlation of the data itself by introducing the same random term within the same level [7, 8]. As an example, take the simplest stratified model (null model):

y_{ij}=\beta +u_j+e_{ij},j=1,2,\dots ,m,i=1,2,\dots ,n_j,

where i,j are the subscripts of the first and second level, respectively, and the random terms e_{ij},u_j are mutually independent and satisfy E\left(e_{ij}\right)=E\left(u_j\right)=0, var\left(e_{ij}\right)={\sigma }_e^2, var\left(u_j\right)={\sigma }_u^2. In this model, different individuals in the same level have the same random term u_j, and the correlation coefficient between them is

\rho =\frac{\mathrm{c}\mathrm{o}\mathrm{v}\left(y_{ij},y_{i^{\prime }j}\right)}{\sqrt{\mathrm{v}\mathrm{a}\mathrm{r}\left(y_{ij}\right)\mathrm{v}\mathrm{a}\mathrm{r}\left(y_{i^{\prime }j}\right)}}=\frac{{\sigma }_u^2}{{\sigma }_u^2+{\sigma }_e^2}.

It can be seen that the weight of {\sigma }_u^2 in the total variance reflects the intra-level correlation. Thus, estimating the covariance array for the hierarchical model is an inescapable problem. The well-known and commonly used estimation methods include maximum likelihood estimation (MLE), restricted maximum likelihood estimation (RMLE), and empirical Bayesian estimation (EB). These methods have some shortcomings: MLE and RMLE involve the differentiation of vectors, which can be difficult to handle when the random part of the model becomes complex; EB has different effects on the results when estimating the prior distribution using empirical data depending on the selection of the estimator [3]. In recent years, scholars have made new explorations and attempts around the estimation method of covariance array. Lele and Taper [5] (2002) proposed the method of compound likelihood (MCL), which reduces the complexity of the algorithm; Wang and Yin [10] (2002) proposed the spectral decomposition method for the hierarchical model of a class of balanced data to obtain an explicit solution for the variance components.

The iterative generalized least squares estimation used in this paper was originally proposed by Goldstein [2] (1986) for the variance component model, i.e.,

Y_{ij}={\beta }_0+\sum _{l=1}^r{\beta }_{l,j}{x}_{l,ij}+\sum _{l=1}^p{\gamma }_l{z}_{l,j}+\left(u_j+e_{ij}\right).

In this model, the random part is independent with respect to the explanatory variable x. Here, we try to derive a specific expression for parameter estimation based on iterative generalized least squares for a more general hierarchical model containing random effects (i.e., random terms shaped as {\mu }_{0j}+u_{1j}{x}_{ij}+e_{ij}). Then we will perform the proof of the relevant properties, and show the effect of estimation through simulations.

2 Model

Consider a general two-level mixed-effects model:

(2.1)\{\begin{array}{l}{y}_{ij}={\alpha }_{0j}+{\alpha }_{1j}{x}_{ij}+e_{ij},\\ {\alpha }_{0j}={\beta }_0+u_{0j,}{\alpha }_{1j}={\beta }_1+u_{1j},\end{array}

where i is the subscript of the first level, i=1,2,\dots ,n_j, j is the subscript of the second level, j=1,2,\dots ,m. Specifically, we can consider the second level as m schools, the first level as the students in each school. The j\mathrm{t}\mathrm{h} school has n_j students. In the model, y_{ij} is the response variable, x_{ij} is the explanatory variable, and e_{ij},u_{0j},u_{1j} are the random terms that satisfies:

\begin{array}{rl}& E\left(u_{0j}\right)=E\left(u_{1j}\right)=E\left(e_{ij}\right)=0,\\ & \mathrm{v}\mathrm{a}\mathrm{r}\left(u_{0j}\right)={\sigma }_{u0}^2,\mathrm{v}\mathrm{a}\mathrm{r}\left(u_{1j}\right)={\sigma }_{u1}^2,\mathrm{c}\mathrm{o}\mathrm{v}\left(u_{0j},u_{1j}\right)={\sigma }_{u01},\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{ij}\right)={\sigma }_e^2,\end{array}

where e_{ij} and u_{0j},u_{1j} are independent. Thus, we can obtain the merged form:

(2.2)y_{ij}={\beta }_0+{\beta }_1{x}_{ij}+\left(u_{0j}+x_{ij}{u}_{1j}+e_{ij}\right).

In this way, y_{ij} is decomposed into the sum of a fixed part {\beta }_0+{\beta }_1{x}_{ij} and a random part u_{0j}+x_{ij}{u}_{1j}+e_{ij}. Thereby, given x_{ij}, the variance of y_{ij} is

\mathrm{v}\mathrm{a}\mathrm{r}\left(y_{ij}|x_{ij}\right)={\sigma }_{u0}^2+2x_{ij}{\sigma }_{u01}+x_{ij}^2{\sigma }_{u1}^2+{\sigma }_e^2.

The covariance between different individuals in the same school is

\mathrm{c}\mathrm{o}\mathrm{v}\left(y_{i_{1}j},y_{i_{2}j}|x_{ij}\right)={\sigma }_{u0}^2+\left(x_{i_{1}j}+x_{i_{2}j}\right){\sigma }_{u01}+x_{i_{1}j}{x}_{i_{2}j}{\sigma }_{u1}^2.

Thus, for the given j\mathrm{t}\mathrm{h} school, the covariance array is

(2.3)V_j=X_j\mathrm{\Omega }{X}_j^{\mathrm{T}}+{\sigma }_e^2{I}_{nj},

where

X_j=\left(\begin{array}{cc}1& x_{1j}\\ ⋮& ⋮\\ 1& x_{n_{j}j}\end{array}\right),\mathrm{\Omega }=\left(\begin{array}{cc}{\sigma }_{u0}^2& {\sigma }_{u01}\\ {\sigma }_{u01}^{}& {\sigma }_{u1}^2\end{array}\right).

In model (2.1), we are interested in the regression coefficients \beta ={\left({\beta }_0,{\beta }_1\right)}^{\mathrm{T}} and the covariance array V=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(V_1,V_2,\dots ,V_m\right). Notice that in (2.3), each piece V_j of the covariance array V can be represented by the variance-covariance component \theta ={\left({\sigma }_{u0}^2,{\sigma }_{u01},{\sigma }_{u1}^2,{\sigma }_e^2\right)}^{\mathrm{T}}. Therefore, the actual needed parameter in the model is \left({\beta }^{\mathrm{T}},{\theta }^{\mathrm{T}}\right). In the following, we give iterative generalized least squares estimates of the parameters.

3 Iterative generalized least squares estimation

Rewrite (2.2) in the form of matrix, we have

(3.1)Y=X\beta +\epsilon ,

where X is the design array, X columns are full rank, \beta ={\left({\beta }_0,{\beta }_1\right)}^{\mathrm{T}}, and \epsilon is the error vector. Given covariance array V and V is non-singular, the generalized least squares estimate of \beta is as follow:

(3.2)\hat{\beta }={\left(X^{\mathrm{T}}{V}^{-1}X\right)}^{-1}{X}^{\mathrm{T}}{V}^{-1}Y.

The estimated covariance is

\mathrm{c}\mathrm{o}\mathrm{v}(\hat{\beta })={\left(X^{\mathrm{T}}{V}^{-1}X\right)}^{-1}.

On the other hand, since E\left(\epsilon {\epsilon }^{\mathrm{T}}\right)=V,\epsilon =Y-X\beta, we can build the regression model for the covariance matrix V when \beta is known:

(3.3)\epsilon {\epsilon }^{\mathrm{T}}=V+R,

where the covariance array V=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(V_1,V_2,\dots ,V_m\right), R is the error array. Any piece V_j of V can be expressed as a linear form of the unknown parameter \theta:

(3.4)V_j=J_{n_j}{\sigma }_{u0}^2+A_j{\sigma }_{u01}^{}+B_j{\sigma }_{u1}^2+I_{n_j}{\sigma }_e^2,

where J_{n_j} is an n_j dimensional matrix with all elements equal to 1 and I_{n_j} is an n_j dimensional unit array,

\begin{array}{rl}& A_j=\left(\begin{array}{cccc}2x_{1j}& x_{1j}+x_{2j}& \cdots & x_{1j}+x_{n_{j}j}\\ x_{1j}+x_{2j}& 2x_{2j}& \cdots & x_{2j}+x_{n_{j}j}\\ ⋮& ⋮& & ⋮\\ x_{1j}+x_{n_{j}j}& x_{2j}+x_{n_{j}j}& \cdots & 2x_{n_{j}j}\end{array}\right),\\ & B_j=\left(\begin{array}{cccc}{x}_{1j}^2& x_{1j}+x_{2j}& \cdots & x_{1j}+x_{n_{j}j}\\ x_{1j}+x_{2j}& x_{2j}^2& \cdots & x_{2j}+x_{n_{j}j}\\ ⋮& ⋮& & ⋮\\ x_{1j}+x_{n_{j}j}& x_{2j}+x_{n_{j}j}& \cdots & x_{n_{j}j}^2\end{array}\right).\end{array}

Substituting (3.4) into (3.3) and doing matrix \mathrm{v}\mathrm{e}\mathrm{c}(\cdot ) operator, we get

\begin{array}{rl}{Y}^{\ast }\triangleq \left(\begin{array}{c}{\epsilon }_{11}^2\\ {\epsilon }_{11}{\epsilon }_{21}\\ ⋮\\ {\epsilon }_{n_{m}m}^2\end{array}\right)=& \left(\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right){\sigma }_{u0}^2+\left(\begin{array}{c}2x_{11}\\ x_{11}+x_{21}\\ ⋮\\ 2x_{n_{m}m}^{}\end{array}\right){\sigma }_{u01}^{}\\ & +\left(\begin{array}{c}{x}_{11}^2\\ x_{11}{x}_{21}\\ ⋮\\ x_{n_{m}m}^2\end{array}\right){\sigma }_{u1}^2+\left(\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right){\sigma }_e^2+\mathrm{v}\mathrm{e}\mathrm{c}\left(R\right).\end{array}

The above equation can be organized as the following form:

(3.5)Y^{\ast }=X^{\ast }\theta +R^{\ast },

where X^{\ast } is the design matrix, X^{\ast }columns be full rank, and R^{\ast } is the error vector. In this way, we transform the problem of estimating the covariance array into generalized least squares estimation of estimating the regression coefficient \theta ={\left({\sigma }_{\mu 0}^2,{\sigma }_{\mu 01},{\sigma }_{\mu 1}^2,{\sigma }_e^2\right)}^{\mathrm{T}} in the linear regression model (3.5). Thus, the generalized least squares estimate of \theta is

(3.6)\hat{\theta }={\left(X^{\ast \mathrm{T}}{V}^{\ast -}{X}^{\ast }\right)}^{-1}{X}^{\ast \mathrm{T}}{V}^{\ast -}{Y}^{\ast },

where V^{\ast } is the covariance array of Y^{\ast }. It can be shown that in the normal case (its proof sees the work of Searle et al. [6]). There is

(3.7)V^{\ast }=\left(V\otimes V\right)\left(I+K\right),

where \otimes denotes the Kronecker product, K is the vector permutation matrix. Further, substituting (3.7) into (3.6), using the Kronecker product and the relevant properties of vector permutation matrix [4, 9], the estimate of \theta in the normal case can be obtained:

(3.8)\hat{\theta }={\left(X^{\ast \mathrm{T}}\left(V^{-1}\otimes V^{-1}\right)X^{\ast }\right)}^{-1}{X}^{\ast \mathrm{T}}\left(V^{-1}\otimes V^{-1}\right)Y^{\ast }.

Estimated covariance array (see the appendix of Section 7 for derivation details):

(3.9)\mathrm{c}\mathrm{o}\mathrm{v}(\hat{\theta })=2{\left(X^{\ast \mathrm{T}}\left(V^{-1}\otimes V^{-1}\right)X^{\ast }\right)}^{-1}.

However, in practice, both V and \beta are unknown. Therefore, in (3.2), we replace the true V with the estimate \hat{V} of the current covariance array to find \hat{\beta }; similarly, in (3.8), we replace the true \beta, V with the current estimates \hat{\beta }, \hat{V} to find {\hat{Y}}^{\ast },{\hat{V}}^{\ast }, respectively, and then obtain \hat{\beta }. This iterative update until convergence, which is the idea of iterative generalized least squares. In summary, the steps of iterative generalized least squares estimation can be summarized as follows.

Given the initial value {\hat{\theta }}^{\left(0\right)}, {\hat{V}}^{\left(0\right)} can be obtained according to equation (3.4), and given the convergence threshold \alpha, the t+1 iteration as follows.

Step 1: {\hat{\beta }}^{(t+1)}=(X^{\mathrm{T}}{({\hat{V}}^{(t)})}^{-1}X{)}^{-1}{X}^{\mathrm{T}}({\hat{V}}^{(t)}{)}^{-1}Y;

Step 2: {\hat{\theta }}^{(t+1)}=\{X^{\ast \mathrm{T}}[{({\hat{V}}^{(t)})}^{-1}\otimes {({\hat{V}}^{(t)})}^{-1}]X^{\ast }{\}}^{-1}{X}^{\ast \mathrm{T}}[{({\hat{V}}^{(t)})}^{-1}\otimes {({\hat{V}}^{(t)})}^{-1}]{\hat{Y}}^{\ast (t+1)}, where {\hat{Y}}^{\ast (t+1)}=\mathrm{v}\mathrm{e}\mathrm{c}\{(Y-X{\hat{\beta }}^{(t+1)}){(Y-X{\hat{\beta }}^{(t+1)})}^{\mathrm{T}}\};

Step 3: The iteration ends when and hold simultaneously; otherwise, go to the step 1; where {\Vert \cdot \Vert }_{\mathrm{\infty }} is the infinite number, and for any given P-dimensional vector x, {\Vert x\Vert }_{\mathrm{\infty }}=max\left\{\left|x_1\right|,\left|x_2\right|,\dots ,\left|x_p\right|\right\}.

It can be seen that the iterative generalized least squares estimation achieves the estimation of regression coefficients as well as the covariance array, which is achieved by the iterative approach based on the generalized least squares estimation. Here, the essential idea about the estimation of the covariance array is to use the sample residual array \epsilon {\epsilon }^{\mathrm{T}} to estimate the covariance array V. For the hierarchical mixed effects model in this paper, the covariance array can be fully represented by the variance-covariance components on intra- and inter-level. Therefore, estimating the covariance array is equivalent to estimating the regression coefficients \theta in the linear regression model (3.5).

4 Properties of iterative generalized least squares estimation

Proposition 1　 For model (2.1) and its matrix form (3.1), (3.5), if Y obeys the multivariate normal distribution and the design array X, X^{\ast } are column full rank. Then for parameters \beta ,\theta, the iterative generalized least squares estimation is equivalent to the maximum likelihood estimation.

Proof　For the n-element normal distribution Y\sim N_n\left(X\beta ,V\right), the log-likelihood function is

(4.1)l=-\frac{n}{2}\ln \left(2\pi \right)-\frac{1}{2}\ln \left(V\right)-\frac{1}{2}{\left(Y-X\beta \right)}^{\mathrm{T}}{V}^{-1}\left(Y-X\beta \right).

For the regression coefficient \beta, the estimated equation is

(4.2)\frac{\mathrm{\partial }l}{\mathrm{\partial }\beta }=-\frac{1}{2}\left(2X^{\mathrm{T}}{V}^{-1}Y+2X^{\mathrm{T}}{V}^{-1}X\beta \right)=0.

Equation (3.2) is obtained by rectifying. It can be seen that the maximum likelihood estimation for \theta is equivalent to the iterative generalized least squares estimation.

The estimation of the variance component \theta is considered below. The log-likelihood function (4.1) is rewritten as the following equation:

(4.3)l=-\frac{n}{2}\ln \left(2\pi \right)-\frac{1}{2}\ln \left(V^{-1}\right)-\frac{1}{2}\mathrm{t}\mathrm{r}\left(V^{-1}S\right),

where S=\left(Y-X\beta \right){\left(Y-X\beta \right)}^{\mathrm{T}}. Since \theta is the parameter in the covariance array V, the estimated equation for \theta based on the maximum likelihood estimation is

(4.4)\frac{\mathrm{\partial }l}{\mathrm{\partial }\theta }=\frac{1}{2}\mathrm{t}\mathrm{r}\left(V\frac{\mathrm{\partial }{V}^{-1}}{\mathrm{\partial }\theta }\right)-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\theta }\mathrm{t}\mathrm{r}\left(V^{-1}S\right)=0.

In the iterative generalized least squares estimation, \hat{\theta } is obtained from equation (3.6), which can be obtained by minimizing

(4.5){\left(Y^{\ast }-X^{\ast }\theta \right)}^{\mathrm{T}}{V}^{\ast -}\left(Y^{\ast }-X^{\ast }\theta \right)={\left[\mathrm{v}\mathrm{e}\mathrm{c}\left(S-V\left(\theta \right)\right)\right]}^{\mathrm{T}}{V}^{\ast -}\mathrm{v}\mathrm{e}\mathrm{c}\left(S-V\left(\theta \right)\right).

When Y\sim N_n\left(X\beta ,V\right) holds, minimizing equation (4.5) is equivalent to minimizing

(4.6){\left[\mathrm{v}\mathrm{e}\mathrm{c}\left(S-V\left(\theta \right)\right)\right]}^{\mathrm{T}}\left(V^{-1}\otimes V^{-1}\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(S-V\left(\theta \right)\right).

It should be noted that we use different symbols V\left(\theta \right) and V to characterize whether the covariance array is related to the estimated parameter \theta. The former is the expression before the straightening operation on X\ast \theta, which is related to \theta, and the latter is the true covariance array in the model, which is treated as a constant in the derivation process below. Let equation (4.6) be G. From the relevant properties of matrix vectorization operation and trace operation, we get

(4.7)G=\mathrm{t}\mathrm{r}{\left[V^{-1}\left(S-V\left(\theta \right)\right)\right]}^2.

Taking the partial derivative of the above equation with respect to \theta, we get

(4.8)\begin{array}{rl}\frac{\mathrm{\partial }G}{\mathrm{\partial }\theta }& =\frac{\mathrm{\partial }\mathrm{t}\mathrm{r}\left[V^{-1}\left(S-V\left(\theta \right)\right)V^{-1}\left(S-V\left(\theta \right)\right)\right]}{\mathrm{\partial }\theta }\\ & =2\mathrm{t}\mathrm{r}\left[-V^{-1}S{V}^{-1}\frac{\mathrm{\partial }V\left(\theta \right)}{\mathrm{\partial }\theta }+V^{-1}V\left(\theta \right)V^{-1}\frac{\mathrm{\partial }V\left(\theta \right)}{\mathrm{\partial }\theta }\right]\\ & =2\mathrm{t}\mathrm{r}\left[S\frac{\mathrm{\partial }{V}^{-1}\left(\theta \right)}{\mathrm{\partial }\theta }+V\left(\theta \right)\frac{\mathrm{\partial }{V}^{-1}\left(\theta \right)}{\mathrm{\partial }\theta }\right].\end{array}

Setting the partial derivative of the above equation equal to 0, the estimated equation for \theta is obtained:

(4.9)\mathrm{t}\mathrm{r}\left(S\frac{\mathrm{\partial }{V}^{-1}\left(\theta \right)}{\mathrm{\partial }\theta }\right)=\mathrm{t}\mathrm{r}\left(V\left(\theta \right)\frac{\mathrm{\partial }{V}^{-1}\left(\theta \right)}{\mathrm{\partial }\theta }\right).

This is the same as (4.4). Therefore, for \theta, the maximum likelihood estimate is equivalent to the iterative generalized least squares estimate.

From the above conclusion, it can be deduced that the property of maximum likelihood estimation is appropriate for iterative generalized least squares estimation in the normal case. We know that for the variance, its maximum likelihood estimate is biased, for example, when Y\sim N\left(X\beta ,V\right), the maximum likelihood estimates for {\sigma }^2 is {\hat{\sigma }}_{ML}^2=\frac{1}{n}(Y-X\hat{\beta }{)}^{\mathrm{T}}(Y-X\hat{\beta }) and E({\hat{\sigma }}_{ML}^2)=\frac{n-\mathrm{r}\mathrm{k}\left(X\right)}{n}{\sigma }^2\ne {\sigma }^2. Therefore, the estimate of the variance obtained using iterative generalized least squares is also biased.

Regarding the large-sample property, we still consider firstly the large-sample property of maximum likelihood estimation, and then naturally transition to iterative generalized least squares estimation in the above equivalence. For the hierarchical model, since the second level samples y_1,y_2,\dots ,y_m are independent of each other, the likelihood functions are as follows:

\varphi \left(\gamma \right)=\prod _{j=1}^m{f}_j\left(y_j,\gamma \right),

where y_j=\left(y_{1j},y_{2j},\dots ,y_{n_{j}j}\right),j=1,2,\dots ,m,\gamma ={\left({\beta }^{\mathrm{T}},{\theta }^{\mathrm{T}}\right)}_{1\times k}^{\mathrm{T}}, f_j is the n_j distribution density function. It should be noted that the samples y_1,y_2,\dots ,y_m are independent but not equally distributed, depending on a common parameter vector \gamma. For the maximum likelihood estimate {\hat{\gamma }}_{ML} in this case, Bradley and Gart [1] proves that {\hat{\gamma }}_{ML} is consistent and asymptotically normal if the following assumptions hold.

Hypothesis I (i)　For any y_j\in R_j,\gamma \in \mathrm{\Omega }, there exist \frac{\mathrm{\partial }\ln f_j}{\mathrm{\partial }{\gamma }_r},\frac{{\mathrm{\partial }}^2\ln f_j}{\mathrm{\partial }{\gamma }_r\mathrm{\partial }{\gamma }_s},\frac{{\mathrm{\partial }}^3\ln f_j}{\mathrm{\partial }{\gamma }_r\mathrm{\partial }{\gamma }_s\mathrm{\partial }{\gamma }_t}, r,s,t=1,2,\dots ,k; j=1,2,\dots ,m.

Hypothesis I (ii)　For any density function f_j, and any y_j\in R_j,\gamma \in \mathrm{\Omega }, there are

where F_{jr}\left(y_j\right),F_{jrs}\left(y_j\right) are integrable in R_j, and

Hypothesis II (i)　For the density function column {\left\{f_j\right\}}_{j=1}^m, we have

\sum _{j=1}^m{\int }_{D_{1j}}{f}_j\mathrm{d}{y}_j=o\left(1\right),\sum _{j=1}^m{\int }_{D_{2j}}\left\{\frac{\mathrm{\partial }\ln f_j}{\mathrm{\partial }{\gamma }_r}\right\}{f}_j\mathrm{d}{y}_j=o\left(n^2\right),

where .

Hypothesis II (ii)　For the density function column {\left\{f_j\right\}}_{j=1}^m, we have

\sum _{j=1}^m{\int }_{D_{3j}}{f}_j\mathrm{d}{y}_j=o\left(1\right),\sum _{j=1}^m{\int }_{D_{4j}}{\left\{\frac{\mathrm{\partial }\ln f_j}{\mathrm{\partial }{\gamma }_r\mathrm{\partial }{\gamma }_s}\right\}}^2{f}_j\mathrm{d}{y}_j=o\left(n^2\right),

where ; and\underset{m\to \mathrm{\infty }}{lim}\frac{1}{m}\sum _{j=1}^m{\int }_{R_j}\frac{\mathrm{\partial }\ln f_j}{\mathrm{\partial }{\gamma }_r}\frac{\mathrm{\partial }\ln f_j}{\mathrm{\partial }{\gamma }_s}{f}_j\mathrm{d}{y}_j exists, r,s=1,2,\dots ,k.

Hypothesis II (iii)　For H_{jrst}\left(y_j\right),M_j, we have

\sum _{j=1}^m{\int }_{D_{5j}}{f}_j\mathrm{d}{y}_j=o\left(1\right),\sum _{j=1}^m{\int }_{D_{6j}}{H_{jrst}}^2\left(y_j\right)f_j\mathrm{d}{y}_j=o\left(n^2\right),

where ; and , M is a positive constant.

Hypothesis III

\underset{m\to \mathrm{\infty }}{lim}\frac{1}{m}\sum _{j=1}^m{\int }_{D_{7j}}\sum _{r=1}^k{\left\{\frac{\mathrm{\partial }\ln f_j}{\mathrm{\partial }{\gamma }_r}\right\}}^2{f}_j\mathrm{d}{y}_j=0,

where D_{7j}=\left\{y_j\in R_j:{\left[{\sum }_{r=1}^k{\left(\frac{\mathrm{\partial }\ln f_j}{\mathrm{\partial }{\gamma }_r}\right)}^2\right]}^{\frac{1}{2}}>\epsilon \sqrt{m},\mathrm{\forall }\epsilon >0\right\}.

Theorem 1　 Let Yobey a multivariate normal distribution, \hat{\gamma } be an iterative generalized least squares estimate of the parameter vector \gamma, and {\gamma }_0 be the true value. Then

(1) Under the conditions that hypothesis I, II hold, \hat{\gamma } has consistence: \hat{\gamma }\stackrel{P}{\to }{\gamma }_0,m\to \mathrm{\infty };

(2) Under the conditions that hypothesis I, II, and III, \hat{\gamma } is asymptotically normal: \sqrt{m}\left(\hat{\gamma }-{\gamma }_0\right)\stackrel{L}{\to }N\left(0,J_0^{-1}\right), where J_0=\left(\underset{m\to \mathrm{\infty }}{lim}\frac{1}{m}\sum _{j=1}^{m}E{\left(\frac{\mathrm{\partial }\ln f_j}{\mathrm{\partial }{\gamma }_r}\frac{\mathrm{\partial }\ln f_j}{\mathrm{\partial }{\gamma }_s}\right)|}_{\gamma ={\gamma }_0}\right).

Proof　The proof of the theorem can be obtained directly from the large-sample property of {\hat{\gamma }}_{ML} given by Bradley and Gart [1] (1962) under the guarantee of the above equivalence proposition.

5 Simulation

To test the estimation effect of the iterative generalized least squares method, we first simulate model (2.1). Given the regression coefficients \beta ={\left({\beta }_0,{\beta }_1\right)}^{\mathrm{T}}={\left(10,2\right)}^{\mathrm{T}}, variance-covariance components \theta =({\sigma }_{u0}^2,{\sigma }_{u01},{\sigma }_{u1}^2,{\sigma }_e^2{)}^{\mathrm{T}}=(1.2,0.5,1,1{)}^{\mathrm{T}}, define \mathrm{\Omega } the same as in equation (2.3). x_{ij} is generated from the standardized normal distribution. The random term {\left(u_{0j},u_{1j}\right)}^{\mathrm{T}} of the second level is generated from two-dimensional normal distribution N\left(0,\mathrm{\Omega }\right), i.e., {\left(u_{0j},u_{1j}\right)}^{\mathrm{T}}\stackrel{iid}{\sim }N\left(0,\mathrm{\Omega }\right), j=1,2,\dots ,m. The random term e_{ij} of the first level is generated from one-dimensional normal distribution N\left(0,{\sigma }_e^2\right), i.e., e_{ij}\stackrel{iid}{\sim }N\left(0,{\sigma }_e^2\right), i=1,2,\dots ,n_j. Taking m=100, n_j is generated from the Poisson distribution with average of 25. Thus, y_{ij} is obtained according to y_{ij}={\beta }_0+{\beta }_1{x}_{ij}+\left(u_{0j}+u_{1j}{x}_{ij}+e_{ij}\right). For illustration purposes, as above, we consider the first and second levels of the model here as the student level and the school level, respectively, so that y_{ij} is the response variable corresponding to the ith student in jth school. More intuitively, we plot the scatter plots for two of the schools as follows.

As can be seen in Fig.1, the intercept and regression coefficients are not the same for each school because different schools correspond to different school-level random terms, and are perturbed based on {\beta }_0,{\beta }_1. We estimated the regression coefficients \beta and variance components \theta using iterative generalized least squares, and the estimation results are shown in Tab.1.

It can be seen that the estimated values of the parameters deviate less from the true values, and the standard deviations of the estimates are also small, which indicates that the parameter estimates obtained by using the iterative generalized least squares method fluctuate in a smaller range around the true values and have better estimation results for model (2.1).

Further, we generalize the model by introducing explanatory variables in the second level, which gives the following hierarchical model:

(5.1)\{\begin{array}{l}{y}_{ij}={\alpha }_{0j}+{\alpha }_{1j}{x}_{ij}+e_{ij},\\ {\alpha }_{0j}={\beta }_{00}+{\beta }_{01}{\omega }_{0j}+u_{0j},{\alpha }_{1j}={\beta }_{10}+{\beta }_{11}{\omega }_{1j}+u_{1j}.\end{array}

The assumptions about the distribution and independence of the random terms in the model are the same as in model (2.1). The two levels are combined to obtain the following form:

(5.2)y_{ij}={\beta }_{00}+{\beta }_{01}{\omega }_{0j}+{\beta }_{10}{x}_{ij}+{\beta }_{11}{\omega }_{1j}{x}_{ij}+\left(u_{0j}+x_{ij}{u}_{1j}+e_{ij}\right).

Given the true values of parameters \left({\beta }_{00},{\beta }_{01},{\beta }_{10},{\beta }_{11},{\sigma }_{u0}^2,{\sigma }_{u01},{\sigma }_{u1}^2,{\sigma }_e^2\right), the data x_{ij},{\omega }_{0j},{\omega }_{1j} are generated independently from the standard normal distribution. Multiple simulations are performed and the estimated results are averaged to obtain Tab.2.

Tab.2 shows that the estimates of the parameters are close to the true values and the standard deviations are also small, which shows that the generalized least squares method is also applicable for the extended model iterations.

6 Example

In this section, we try to develop the hierarchical mixed-effects model to analyze the relationship between GDP and total foreign trade (TMX) for several major cities with export-oriented economies. Seven coastal cities with export-oriented economies, namely Tianjin, Shanghai, Dalian, Qingdao, Guangzhou, Shenzhen and Ningbo, are selected as the study subjects, and the GDP indicators and TMX indicators of the seven cities in recent years are recorded based on the data released by the National Bureau of Statistics. For Tianjin and Shanghai, we record the data for a total 25 years from 1990‒2014, and for the remaining five cities, we record the data for a total of 9 years from 2006‒2014. It should be noted that, the above data imbalance is caused by the fact that data of municipalities directly under the central government are more accessible and have longer coverage time than those of other cities. However, the model and estimation method involved in this paper do not have any requirement on whether the data are balanced or not, and can be used to analyze this data set. Using the economic indicators of different cities in different years as the first level, the first level model is developed as follows:

(6.1){\mathrm{G}\mathrm{D}\mathrm{P}}_{ij}={\alpha }_{0j}+{\alpha }_{1j}{\mathrm{T}\mathrm{M}\mathrm{X}}_{ij}+e_{ij},j=1,2,\dots ,7,i=1,2,\dots ,n_j,

where n_j=9 or 25. Since the data are centralized, the intercept term {\alpha }_{0j} represents the GDP value corresponding to the j\mathrm{t}\mathrm{h} city when its TMX is at the total average level {\overline{\mathrm{T}\mathrm{M}\mathrm{X}}}_{ij}, which should be related to the city's own average {\overline{\mathrm{T}\mathrm{M}\mathrm{X}}}_j. Therefore, when building a linear model for the intercept term, we introduce{\overline{\mathrm{T}\mathrm{M}\mathrm{X}}}_j as an explanatory variable to build the second level model:

(6.2)\{\begin{array}{l}{\alpha }_{0j}={\beta }_{00}+{\beta }_{01}{\overline{\mathrm{T}\mathrm{M}\mathrm{X}}}_j+u_{0j},\\ {\alpha }_{1j}={\beta }_{10}+u_{1j}.\end{array}

Combining (6.1) and (6.2), we get

(6.3){\mathrm{G}\mathrm{D}\mathrm{P}}_{ij}={\beta }_{00}+{\beta }_{01}{\overline{\mathrm{T}\mathrm{M}\mathrm{X}}}_{ij}+{\beta }_{10}{\mathrm{T}\mathrm{M}\mathrm{X}}_{ij}+\left(u_{0j}+{\mathrm{T}\mathrm{M}\mathrm{X}}_{ij}{u}_{1j}+e_{ij}\right).

The regression coefficients \left({\beta }_{00},{\beta }_{01},{\beta }_{10}\right) and variance-covariance components \left({\sigma }_{u0}^2,{\sigma }_{u01},{\sigma }_{u1}^2,{\sigma }_e^2\right) in model (6.3) were estimated using the iterative generalized least squares method, and the results were shown in Tab.3.

From Tab.3, we can see that , which means that when {\mathrm{T}\mathrm{M}\mathrm{X}}_{ij}={\overline{\mathrm{T}\mathrm{M}\mathrm{X}}}_{ij}, if the average level of foreign trade of the city is high in recent years, it will lower the GDP of the city in that year. This means the economy of a city depends excessively on foreign trade, which will have a negative impact on economic development. In addition, the average effect of TMX on GDP, i.e., the slope’s expectation {\beta }_{10}=0.78, and the slope’s standard deviation \sqrt{{\sigma }_{u1}}=0.46. So, the effect {\alpha }_{1j} of TMX on GDP is roughly in the interval \left({\beta }_{10}-2\sqrt{{\sigma }_{u1}^2},{\beta }_{10}+2\sqrt{{\sigma }_{u1}^2}\right)=\left(0.32,1.24\right). This indicates that the effect of TMX on GDP is always positive for different cities. However, the effect of TMX on GDP is higher in some cities and can reach about four times the low effect. Furthermore, the correlation coefficient \rho =\frac{{\sigma }_{u01}}{\sqrt{{\sigma }_{u0}^2}\sqrt{{\sigma }_{u1}^2}}>0 between the city-level random terms u_{0j},u_{1j}. Therefore, if a city adds the positive perturbation to the intercept term, then the effect of TMX on GDP tends to be larger.

7 Appendix

Conclusions related to the micro-quotient of matrices

\frac{\mathrm{\partial }}{\mathrm{\partial }t}\ln \left|A\left(t\right)\right|=\mathrm{t}\mathrm{r}\left(A^{-1}\left(t\right)\frac{\mathrm{\partial }A\left(t\right)}{\mathrm{\partial }t}\right),\frac{\mathrm{\partial }A\left(t\right)}{\mathrm{\partial }t}=-A^{-1}\left(t\right)\frac{\mathrm{\partial }A\left(t\right)}{\mathrm{\partial }t}{A}^{-1}\left(t\right).

Conclusion related to matrix vectorization operations, trace operations

\mathrm{v}\mathrm{e}\mathrm{c}(ABC)=\left(C^{\mathrm{T}}\otimes A\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(B\right),\mathrm{t}\mathrm{r}\left(AB\right)={\left(\mathrm{v}\mathrm{e}\mathrm{c}\left(A^{\mathrm{T}}\right)\right)}^{\mathrm{T}}\mathrm{v}\mathrm{e}\mathrm{c}\left(B\right).

Vector permutation matrix

Definition: If for any given m\times n dimensional matrix, there is always\mathrm{v}\mathrm{e}\mathrm{c}\left(A\right)=K\mathrm{v}\mathrm{e}\mathrm{c}\left(A^{\mathrm{T}}\right), then the mn\times mn dimensional matrix K is said to be vector permutation matrix.

Properties: (1) K=K^{\prime }, (2) K^2=I, (3) {\left(I+K\right)}^2=2\left(I+K\right), (4) K(A\otimes B)=\left(B\otimes A\right)K.

From the above properties, the following conclusions can be drawn:

\left(V\otimes V\right){\left(I+K\right)}^{-}=\frac{1}{4}\left(I+K\right)\left(V^{-1}\otimes V^{-1}\right).

The derivation of Eq. (3.9) Since

\hat{\theta }={\left(X^{\ast \mathrm{T}}\left(V^{-1}\otimes V^{-1}\right)X^{\ast }\right)}^{-1}{X}^{\ast \mathrm{T}}\left(V^{-1}\otimes V^{-1}\right)Y^{\ast },

note thatM={\left(X^{\ast \mathrm{T}}\left(V^{-1}\otimes V^{-1}\right)X^{\ast }\right)}^{-1}, then

(7.1)\begin{array}{rl}\mathrm{c}\mathrm{o}\mathrm{v}(\hat{\theta })& =M{X}^{\ast \mathrm{T}}\left(V^{-1}\otimes V^{-1}\right)\mathrm{c}\mathrm{o}\mathrm{v}\left(Y^{\ast }\right){\left(V^{-1}\otimes V^{-1}\right)}^{\mathrm{T}}{X}^{\ast }{M}^{\mathrm{T}}\\ & =M{X}^{\ast T}\left(V^{-1}\otimes V^{-1}\right)\left(V\otimes V\right)\left(I+K\right){\left(V^{-1}\otimes V^{-1}\right)}^{\mathrm{T}}{X}^{\ast }{M}^{\mathrm{T}},\end{array}

where \left(V\otimes V\right)\left(I+K\right){\left(V^{-1}\otimes V^{-1}\right)}^{\mathrm{T}}{X}^{\ast }=X^{\ast }+\left(V\otimes V\right)K{\left(V^{-1}\otimes V^{-1}\right)}^{\mathrm{T}}{X}^{\ast }.

Here X^{\ast }=\left(\mathrm{v}\mathrm{e}\mathrm{c}\left(A_1\right)|\mathrm{v}\mathrm{e}\mathrm{c}\left(A_2\right)|\mathrm{v}\mathrm{e}\mathrm{c}\left(A_3\right)|\mathrm{v}\mathrm{e}\mathrm{c}\left(A_4\right)\right), A_1,A_2,A_3,A_4 correspond to design arrays of {\sigma }_{u0}^2, {\sigma }_{u01}, {\sigma }_{u1}^2, {\sigma }_e^2, respectively, and all are symmetric arrays. Thus:

(7.2)\begin{array}{rl}& (V\otimes V)(I+K)(V^{-1}\otimes V^{-1}{)}^{\mathrm{T}}{X}^{\ast }\\ & =X^{\ast }+(V\otimes V)K(V^{-1}\otimes V^{-1}{)}^{\mathrm{T}}(\mathrm{v}\mathrm{e}\mathrm{c}(A_1)|\mathrm{v}\mathrm{e}\mathrm{c}(A_2)|\mathrm{v}\mathrm{e}\mathrm{c}(A_3)|\mathrm{v}\mathrm{e}\mathrm{c}(A_4))\\ & =X^{\ast }+(V\otimes V)K(\mathrm{v}\mathrm{e}\mathrm{c}(V^{-1}{A}_1{V}^{-1})|\mathrm{v}\mathrm{e}\mathrm{c}(V^{-1}{A}_2{V}^{-1})\\ & |\mathrm{v}\mathrm{e}\mathrm{c}(V^{-1}{A}_3{V}^{-1})|\mathrm{v}\mathrm{e}\mathrm{c}(V^{-1}{A}_4{V}^{-1}))\\ & =X^{\ast }+(V\otimes V)(\mathrm{v}\mathrm{e}\mathrm{c}(V^{-1}{A}_1{V}^{-1})|\mathrm{v}\mathrm{e}\mathrm{c}(V^{-1}{A}_2{V}^{-1})\\ & |\mathrm{v}\mathrm{e}\mathrm{c}(V^{-1}{A}_3{V}^{-1})|\mathrm{v}\mathrm{e}\mathrm{c}(V^{-1}{A}_4{V}^{-1}))\\ & =X^{\ast }+(V\otimes V)(V^{-1}\otimes V^{-1}{)}^{\mathrm{T}}{X}^{\ast }\\ & =2X^{\ast }.\end{array}

Substituting (7.2) into (7.1), we get

\mathrm{c}\mathrm{o}\mathrm{v}(\hat{\theta })=2{\left(X^{\ast \mathrm{T}}\left(V^{-1}\otimes V^{-1}\right)X^{\ast }\right)}^{-1}.◻

References

[1] Bradley R A, Gart J J. The asymptotic properties of ML estimators when sampling from associate populations. Biometrika 1962; 49(1/2): 205–214

[2] Goldstein H. Hierarchical mixed linear model analysis using iterative generalized least squares. Biometrika 1986; 73(1): 43–56

[3] GoldsteinH. Hierarchical Statistical Models, 4th ed. Chichester: John Wiley & Sons, 2011

[4] Henderson H V, Searle S R. The vec-permutation matrix the vec operator and Kronecker products: A review. Linear Multilinear Algebra 1981; 9(4): 271–288

[5] Lele S, Taper M L. A composite likelihood approach to (co)variance components estimation. J Statist Plann Inference 2002; 103(1/2): 117–135

[6] SearleS RCasellaGMcCullochC E. Variance Components. Hoboken, NJ: John Wiley & Sons, 2006

[7] SnijdersT A BBoskerR J. Hierarchical Analysis: An Introduction to Basic and Advanced Hierarchical Modeling. London Sage, 1999

[8] TianM Z. Higher hierarchical Quantile Regression Modeling Theory. Beijing: Science Press, 2015

[9] WangS GShiJ HYinS JWuM X. Introduction to linear models. Beijing: Science Press, 2004

[10] Wang S G, Yin S J. A new estimate of the parameters in linear mixed models. Sci Sin Math 2002; 32(5): 434–443

RIGHTS & PERMISSIONS

Higher Education Press 2023