基于SPSRs准则下的分类混合模型预测方法研究

The Study of Classified Mixed Model Prediction Method Based on SPSRs Rules

摘要: 分类混合模型预测(CMMP)方法是近年来小区域估计领域中提出的一种新方法，该方法是在待预测效应识别后的基础上形成的方法，较传统的混合效应预测方法有更高的预测精度，得到许多统计学者的关注。最早的分类混合模型预测方法是基于均方预测误差(MSPE)准则进行分类识别构造最佳预测。MSPE准则虽然是一个具有较好数学性质(对称性和平滑性)的不确定性度量准则，但是其不是一个严格适当的评分准则(SPSRs)。因此，提出了基于SPSRs准则(即对数评分)进行分类识别，构造最佳预测的方法。首先，在最佳预测的基础上构造了SPSRs分类器，并进行识别预测；其次分析了该预测的渐近性质，并通过数值模拟证明了该方法较经典的回归预测方法具有更高的准确度；最后，给出实例进一步论证了我们的理论结果。

Abstract: The Classified Mixed Model Prediction (CMMP) method is a newly proposed method in the field of small area estimation in recent years. The prediction accuracy of CMMP has attracted the attention of many statisticians. The earliest Classified Mixed Model Prediction is based on the Mean Square Prediction Error (MSPE) criterion for classification, identification and construction of the best pre-diction. Although the MSPE criterion is an uncertainty measurement criterion with good mathe-matical properties (symmetry and smoothness), it is not a Strictly Proper Scoring Rules (SPSRs). Therefore, we propose a method for classification and identification based on SPSRs criterion (i.e. logarithmic score) to construct the best prediction. Firstly, on the basis of the best prediction, the SPSRs classifier is constructed, identified and predicted. Secondly, the asymptotic properties of the prediction are analyzed, and the numerical simulation proves that this method has higher accuracy than the classical regression prediction method. Finally, an example is given to further demonstrate our theoretical results.

1. 引言

在现实生活中许多领域都涉及对相关特征预测的问题，而且其问题通常在个体层次或亚群体层次发生。在对特征预测问题分析时，可将新观测数据(含有特征信息的数据)在训练集(已有)数据中，找到一个与其相匹配的群组，即借用这个群组中训练集数据的相关信息来提高预测的准确性。基于这思想，Jiang等 [1] 提出了分类混合模型预测(Classified Mixed Model Prediction, CMMP)，通过数值模拟显示，在预测准确性方面上，CMMP方法显著优于回归预测(Regression Prediction, RP)方法。CMMP方法得到统计学者的认同，此外该方法在小区域估计领域存在广泛的应用，如刘育孜等 [2] 通过用CMMP方法对农作物面积的估算；王婕雯等 [3] 通过用CMMP方法对小麦面积估算。同时，基于CMMP思想，有许多文献对该方法进行了改进和推广，如Sun等 [4] 将适用于连续响应变量的线性混合模型的CMMP方法拓展到集群二分类数据(离散响应变量)，从而提出了分类混合逻辑模型预测方法。此外，Sun等 [5] 又对该方法进行改进，结合协变量数据中信息进行分类匹配，提出新的分类混合模型预测方法，其预测准确性也优于(Mixed Model Prediction, MMP)方法 [6]。

CMMP方法的匹配判别策略度量是均方预测误差(Mean Squared Prediction Error, MSPE)，研究发现，在这种准则下，即使训练集数据中的某组数据和新观测数据的确来自同一群组，但是在判别匹配关系即 值时，仍有可能存在不正确的情况，由于这种匹配的错误率，从而影响CMMP方法预测的准确性。因此通过改善CMMP方法的匹配准则，提高匹配的正确率，就有可能提高CMMP方法预测的正确性。Gneiting等 [7] 提出严格适当评分准则(Strictly Proper Scoring Rules, SPSRs)。SPSRs优点在于分布预测，比点预测更加实用，同时分布预测能提供更多有用的信息。SPSRs准则也得到各学者的相应研究，比如Merkle和Steyvers [8]，Landes [9]，Du, H.L. [10] 等。根据SPSRs的定义，MSPE准则满足“适当的(Proper)”的条件，但不满足“严格适当的(Strictly Proper)”条件，然而“严格适当的(Strictly Proper)”的评分规则优于“适当的(Proper)”的评分规则。从这个角度看，将CMMP方法中的匹配准则(MSPE)换成SPSRs准则，则匹配的正确率可能会提高，从而提高预测的准确性，为此我们提出SPSRs准则下分类混合效应预测记为CMMP-SPSRs。在SPSRs中，关于分类变量的评分规则提出多种方法，如Brier评分(Brier score)，对数评分(Logarithmic score)和0~1评分(Zero-one score)等，本章选择其中的对数评分(Logs)进行相关分析。

本文内容安排如下：第二节将在嵌套误差回归模型中，提出CMMP-SPSRs方法的混合效应预测，并验证该方法预测的优良性。在第三节通过数值模拟分析，将CMMP-SPSRs方法与RP方法比较。第四节用一个真实数据对CMMP-SPSRs方法进行验证。

2. 基于SPSRs准则下分类混合模型预测

假设已知一组训练集数据为， $y_{ij},i=1,\cdots ,m;j=1,\cdots ,n_i$ 并且知其分组情况，即 $y_{ij}$ 属于第i组的第j个数据。假定训练集数据可采用嵌套误差回归(Nested-Error Regression, NER)模型来建模，如下所示

$y_{ij}=x_{ij}^{\text{T}}\beta +{\alpha }_i+{ϵ}_{ij},$ (1)

其中 $x_{ij}$ 是已知协变量向量， $\beta$ 为未知的回归系数向量，即固定效应， ${\alpha }_i$ 是随机效应， ${ϵ}_{ij}$ 为误差项。同时也假设随机效应 ${\alpha }_i$ 和误差项 ${ϵ}_{ij}$ 相互独立，且有 ${\alpha }_i~N\left(0,{\sigma }_{\alpha }^2\right)$ 和 ${ϵ}_{ij}~N\left(0,{\sigma }_{ϵ}^2\right)$。

假设现有一组新观测数据为 $y_{n,1},\cdots ,y_{n,j},1\le j\le n_{new}$ 其下标n与上文含义相同，只做符号区分并没有其他含义。假定新观测数据也可以嵌套误差回归模型来建模

$y_{n,j}=x_n^{\text{T}}\beta +{\alpha }_I+{ϵ}_{n,j},1\le j\le n_{new},$ (2)

其中 $x_n$ 为不依赖j的协变量；同时参数 $\beta$ 和(1)式中的 $\beta$ 相同； ${\alpha }_I$ 可能为：第一种情况(匹配 $I\in \left\{1,\cdots ,m\right\}$ )，与 ${\alpha }_i,1\le i\le m$ 中的某一个相同，但并不知道具体哪一个，需要在训练集数据中寻找与新观测数据相匹配的群组；第二种情况(不匹配 $I\notin \left\{1,\cdots ,m\right\}$ )，一个全新的随机效应。假设 $E\left({\alpha }_I\right)=0$, $var\left({\alpha }_I\right)={\sigma }_{n,\alpha }^2$ 有界。 ${ϵ}_{n,j},1\le j\le n_{new}$ 是新观测误差项且相互独立，并假设 $E\left({ϵ}_{n,j}\right)=0$, $var\left({ϵ}_{n,j}\right)={\sigma }_{n,ϵ}^2$ 有界，且和训练集数据中的 ${\alpha }_i$， ${ϵ}_{ij}$ 也相互独立。此外， ${\alpha }_I$ 和 ${ϵ}_{n,j}$ 并不要求服从正态分布，同时也不要求 ${\sigma }_{n,\alpha }^2={\sigma }_{\alpha }^2$， ${\sigma }_{n,ϵ}^2={\sigma }_{ϵ}^2$。则需要预测的混合效应为

$\theta =E\left(y_{n,j}|{\alpha }_I\right)=x_n^{\text{T}}\beta +{\alpha }_I,$ (3)

其中此时的参数 $\beta ,{\sigma }_{\alpha }^2,{\sigma }_{ϵ}^2$ 表示已知的；同时记 $\hat{\beta },{\hat{\sigma }}_{\alpha }^2,{\hat{\sigma }}_{ϵ}^2$ 分别表示其对应的一致估计，可以利用训练集数据得到 $\hat{\beta },{\hat{\sigma }}_{\alpha }^2,{\hat{\sigma }}_{ϵ}^2$。理由如下：由于新观测数据中的样本量通常不是很大，如果仅仅通过新观测数据提供的信息来对参数进行估计是不够的。虽然新观测数据不足，但是训练集数据的样本量通常是很多的。因此，通过训练集数据来对相关参数进行估计远优于新观测数据。很明显，如果参数估计值更准确，此时的经验最佳预测和最佳预测两者的预测也会更接近，同时得到的预测结果也更好。本章参数 $\beta ,{\sigma }_{\alpha }^2,{\sigma }_{ϵ}^2$ 统一采用极大似然估计。

2.1. 新观测数据与训练集存在匹配关系

假设新观测数据与训练集数据中的某一群组来自同一组群，即 $I=i\in \left\{1,\cdots ,m\right\}$，但并不知道属于哪一组，故需要寻找具体的参数I，可根据已知的训练集数据提供的相关信息来估计出这个具体的I值，即 $\hat{I}\in \left\{1,\cdots ,m\right\}$，则新观测数据的混合效应(3)式可写成 $\theta =x_n^{\text{T}}\beta +{\alpha }_i$，通过此表达式，可以看出此时新观测数据与训练集数据中的第i群组相匹配，故 $y_{ij}$ 与(3)式混合效应的关系为： $y_1,\cdots ,y_{i-1},{\left(y_i^{\text{T}},\theta \right)}^{\text{T}},y_{i+1},\cdots ,y_m$。混合效应 $\theta =x_n^{\text{T}}\beta +{\alpha }_i$ 的最佳预测，即 $\theta |y$，给定训练集数据后的条件期望，

${\tilde{\theta }}_{\left(i\right)}=E\left(\theta |y\right)=E\left(\theta |y_1,\cdots ,y_m\right)=E\left(\theta |y_i\right).$

根据模型的正态性假定，则有,

$\left(\begin{array}{c}{y}_i\\ {\alpha }_i\end{array}\right)~N\left(\begin{array}{cc}\left(\begin{array}{c}{x}_i^{\text{T}}\beta \\ 0\end{array}\right)& \left(\begin{array}{cc}{\sigma }_{\alpha }^2{J}_{n_i}+{\sigma }_{ϵ}^2{I}_{n_i}& {\sigma }_{\alpha }^2{1}_{n_i}\\ {\sigma }_{\alpha }^2{1}_{n_i}^{\text{T}}& {\sigma }_{\alpha }^2\end{array}\right)\end{array}\right),$

其中训练集数据 并且 $1_n$ 记为n维全1列向量， $I_n$ 记为n阶单位矩阵， $J_n$ 表示元素全为1的n阶矩阵。那么待预测新观测数据的混合效应最佳预测为

$\begin{array}{l}{\mu }_i={\tilde{\theta }}_{\left(i\right)}=E\left(x_n^{\text{T}}\beta +{\alpha }_i|y_i\right)\\ =x_n^{\text{T}}\beta +E\left({\alpha }_i|y_i\right)\\ =x_n^{\text{T}}\beta +\frac{n_i{\sigma }_{\alpha }^2}{n_i{\sigma }_{\alpha }^2+{\sigma }_{ϵ}^2}\left({\bar{y}}_{i.}-{\bar{x}}_{i.}^{\text{T}}\beta \right).\end{array}$ (4)

同时，其 $\theta |y$，给定训练集数据后的条件的方差为，

$\begin{array}{l}{\sigma }_i^2=Var\left(\theta |y\right)=Var\left(\theta |y_1,\cdots ,y_m\right)=Var\left(\theta |y_i\right)\\ =Var\left(x_n^{\text{T}}\beta +{\alpha }_i|y_i\right)\\ \text{}\text{}\text{}=Var\left(x_n^{\text{T}}\beta \right)+Var\left({\alpha }_i|y_i\right)+2Cov\left(x_n^{\text{T}}\beta ,{\alpha }_i|y_i\right)\\ =\frac{{\sigma }_{\alpha }^2{\sigma }_{ϵ}^2}{n_i{\sigma }_{\alpha }^2+{\sigma }_{ϵ}^2}.\end{array}$ (5)

将(4)式中的参数 $\beta ,{\sigma }_{\alpha }^2,{\sigma }_{ϵ}^2$ 分别用它们的一致估计 $\hat{\beta },{\hat{\sigma }}_{\alpha }^2,{\hat{\sigma }}_{ϵ}^2$ 代替，则可得到经验最佳预测，

${\hat{\mu }}_i=x_n^{\text{T}}\hat{\beta }+\frac{n_i{\hat{\sigma }}_{\alpha }^2}{{\hat{\sigma }}_{ϵ}^2+n_i{\hat{\sigma }}_{\alpha }^2}\left({\bar{y}}_{i.}-{\bar{x}}_{i.}^{\text{T}}\hat{\beta }\right),$ (6)

同样，将(5)式中的参数换成其对应的估计，则 $\theta |y$，给定训练集数据后的条件方差估计表达式为，

${\hat{\sigma }}_i^2=\frac{{\hat{\sigma }}_{\alpha }^2{\hat{\sigma }}_{ϵ}^2}{{\hat{\sigma }}_{ϵ}^2+n_i{\hat{\sigma }}_{\alpha }^2},$ (7)

综上所述，既有 $\theta |y$ 的期望也有方差，故其表达形式为： $\theta |y~N\left({\mu }_i,{\sigma }_i^2\right)$。从而 $\theta |y$ 的密度函数为，

$f\left(\theta |y\right)=f_i\left(\theta |y\right)=\frac{1}{\sqrt{2\pi }{\sigma }_i^2}\exp \left\{-\frac{{\left(\theta -{\mu }_i\right)}^2}{2{\sigma }_i^2}\right\}.$ (8)

根据 $\theta |y$ 的密度函数(8)式，又根据SPSRs的定义，可得对数评分(Logs)表达式为，

$\log \left\{f_i\left(\theta |y\right)\right\}\propto -\left\{\log {\sigma }_i^2+\frac{{\left({\mu }_i-\theta \right)}^2}{{\sigma }_i^2}\right\}.$ (9)

接下来对参数I进行估计。采用SPSRs匹配准则来估计参数I进行估计，根据其定义有

$\begin{array}{c}{\text{SPSR}}_i=E\left\{\log {\hat{\sigma }}_i^2+\frac{{\left({\hat{\mu }}_i-\theta \right)}^2}{{\hat{\sigma }}_i^2}\right\}\\ =E\left\{\log {\hat{\sigma }}_i^2+\frac{{\left({\hat{\mu }}_i-{\bar{y}}_n+{\overline{ϵ}}_n\right)}^2}{{\hat{\sigma }}_i^2}\right\}\\ =E\left\{\log {\hat{\sigma }}_i^2+\frac{{\left({\hat{\mu }}_i-{\bar{y}}_n\right)}^2}{{\hat{\sigma }}_i^2}+\frac{2\left({\hat{\mu }}_i-{\bar{y}}_n\right){\overline{ϵ}}_n}{{\hat{\sigma }}_i^2}+\frac{{\left({\overline{ϵ}}_n\right)}^2}{{\hat{\sigma }}_i^2}\right\}\\ =E\left\{\log {\hat{\sigma }}_i^2+\frac{{\left({\hat{\mu }}_i-{\bar{y}}_n\right)}^2}{{\hat{\sigma }}_i^2}+\frac{{\sigma }_{ϵ}^2}{\left(n_{new}{\hat{\sigma }}_i^2\right)}\right\}.\end{array}$ (10)

那么与(10)式相对应的SPSRs即为期望符号内的表达式，从而提出基于SPSRs的匹配准则

${\hat{I}}_{LogS}=\arg {\min }_{1\le i\le m}\left\{\log {\hat{\sigma }}_i^2+\frac{{\left({\hat{\mu }}_i-{\bar{y}}_n\right)}^2}{{\hat{\sigma }}_i^2}+\frac{{\hat{\sigma }}_{ϵ}^2}{\left(n_{new}{\hat{\sigma }}_i^2\right)}\right\},$ (11)

其中 $\theta ={\bar{y}}_n-{\overline{ϵ}}_n$， ${\bar{y}}_n=n_{new}^{-1}{\displaystyle {\sum }_{j=1}^{n_{new}}{y}_{n,j}}$， ${\overline{ϵ}}_n=n_{new}^{-1}{\displaystyle {\sum }_{j=1}^{n_{new}}{ϵ}_{n,j}}$。

最终，通过将 ${\hat{I}}_{LogS}$ 替换(6)式中的i，则 $\theta$ 的分类混合效应预测(Classified Mixed-Effect Predictor, CMEP)为 $\hat{\theta }={\hat{\mu }}_{\hat{I}}$。

下面的定理阐述了当给定一些合理条件下，则可保持分类混合效应预测的渐进性质。

记 $N=\left(m,n_i,1\le i\le m,n_{new}\right)$ 是本文中所涉及的全部样本容量。假设方差参数 $\psi$ 的参数空间为 $\Psi$，而参数 $\beta$ 的参数空间是p维的欧几里得空间 $R^p$。 $\Vert A\Vert =\sqrt{{\lambda }_{\max }\left(A^{\prime }A\right)}$ 表示矩阵A的范数。让 ${\mu }_{\left(0i^{\prime }\right)}$ 表示(4)式等号右边的部分，但其中i被 $i^{\prime }$ 替换，并且参数 $\beta ,\psi$ 为真参数向量。当参数 $\beta ,\psi$ 换成其一致估计 $\hat{\beta },\hat{\psi }$ 时，则对应其经验最佳预测，记为 ${\hat{\mu }}_{\left(i^{\prime }\right)}$。做如下假设：

A1. 当 $I=i\in \left\{1,\cdots ,m\right\}$ 情形下，将i换成 $i^{\prime }$ 时，模型(1)依然成立；下标由I换成i时，模型(2)也成立。

A2. 参数空间 $\Psi$ 的内点为参数 $\psi$ 的真值 ${\psi }_0$。

A3. $x_{ijk},1\le i\le m,1\le j\le n_i$ 是有界的，其中 $x_{ijk}$ 为 $x_{ij}$ 的第k个元素。

A4. 记 $R_{\max }=\underset{1\le i\le m}{\max }{n}_i$, $R_{\min }=\underset{1\le i\le m}{\min }{n}_i$，若 $n_{new}\to \infty ,R_{\min }\to \infty$ 时， $R_{\max }-R_{\min }\to 0$，且 $\hat{\beta }-{\beta }_0=O_p\left(a_N\right)$, $\hat{\psi }-{\psi }_0=O_p\left(b_N\right)$， ${\max }_{1\le i\le m}\left|{\alpha }_i\right|=O_p\left(c_N\right)$，而且当 $\psi ={\psi }_0$ 时， ${\max }_{1\le i\le m}\left|I_{n_i}{ϵ}_i\right|=O_p\left(d_N\right)$ 其中 $a_N,b_N$ 等为正常数，且使得 $d_N\left(a_N+b_N+c_N{b}_N\right)\to 0$， $d_N^2\to 0$。为了不失一般性，假定 $c_N\ge 1$。

定理1. 如果假设A1~A4都成立，那么有 ${\hat{\mu }}_{\hat{I}}\stackrel{p}{\to }\theta$，即SPSRs准则下的分类混合效应预测满足渐进性。

证明：首先考虑当 $I=i\in \left\{1,\cdots ,m\right\}$ 存在匹配关系时的情形。通过本文上述理论的推导与证明，从而可知新观测数据的分类混合效应预测(CMEP)可表示为 ${\hat{\mu }}_{\hat{I}}=x_n^{\text{T}}\hat{\beta }+n_{\hat{I}}{\hat{\sigma }}_{\alpha }^2/{\hat{\sigma }}_{ϵ}^2+n_{\hat{I}}{\hat{\sigma }}_{\alpha }^2\left({\bar{y}}_{\hat{I}.}-{\bar{x}}_{\hat{I}.}^{\text{T}}\hat{\beta }\right)$。接下来证明 ${\hat{\mu }}_{\hat{I}}\stackrel{p}{\to }\theta$。

记 $b={\sigma }_{ϵ}^2/{\sigma }_{\alpha }^2$，则有 ${\sigma }_i^2={\sigma }_{ϵ}^2\cdot {\sigma }_{\alpha }^2/{\sigma }_{ϵ}^2+n_i{\sigma }_{\alpha }^2={\sigma }_{ϵ}^2/b+n_i$，同时也有 $B_i=n_i{\hat{\sigma }}_{\alpha }^2/{\hat{\sigma }}_{ϵ}^2+n_i{\hat{\sigma }}_{\alpha }^2=n_i/b+n_i\le 1$。设 $B_{0i^{\prime }}$ 表示 $B_{i^{\prime }}$ 中 ${\sigma }_{\alpha }^2={\sigma }_{\alpha 0}^2$， ${\sigma }_{ϵ}^2={\sigma }_{ϵ0}^2$，而有

$\begin{array}{c}{\mu }_{\left(0i^{\prime }\right)}=x_n^{\text{T}}{\beta }_0+B_{0i^{\prime }}\left({\bar{y}}_{i^{\prime }.}-{\bar{x}}_{i^{\prime }.}^{\text{T}}{\beta }_0\right)\\ =x_n^{\text{T}}{\beta }_0+\frac{1}{n_{i^{\prime }}+b_0}{1}_{n_{i^{\prime }}}^{\text{T}}\left(y_{i^{\prime }}-x_{i^{\prime }}{\beta }_0\right)\\ =x_n^{\text{T}}{\beta }_0+\frac{1}{n_{i^{\prime }}+b_0}{1}_{n_{i^{\prime }}}^{\text{T}}\left(1_{n_{i^{\prime }}}{\alpha }_{i^{\prime }}+{ϵ}_{i^{\prime }}\right),\end{array}$

$\theta =x_n^{\text{T}}\beta +{\alpha }_I=x_n^{\text{T}}{\beta }_0+{\alpha }_i,$

因此，有

$\begin{array}{c}{\mu }_{\left(0i^{\prime }\right)}-\theta =x_n^{\text{T}}{\beta }_0+\frac{1}{n_{i^{\prime }}+b_0}{1}_{n_{i^{\prime }}}^{\text{T}}\left(1_{n_{i^{\prime }}}{\alpha }_{i^{\prime }}+{ϵ}_{i^{\prime }}\right)-x_n^{\text{T}}{\beta }_0-{\alpha }_i\\ ={\alpha }_{i^{\prime }}-{\alpha }_i+\left(\frac{n_{i^{\prime }}}{n_{i^{\prime }}+b_0}-I_q\right){\alpha }_{i^{\prime }}+\frac{1}{n_{i^{\prime }}+b_0}{1}_{n_{i^{\prime }}}^{\text{T}}{ϵ}_{i^{\prime }}\\ ={\alpha }_{i^{\prime }}-{\alpha }_i-\frac{b_0}{n_{i^{\prime }}+b_0}{\alpha }_{i^{\prime }}+\frac{1}{n_{i^{\prime }}+b_0}{1}_{n_{i^{\prime }}}^{\text{T}}{ϵ}_{i^{\prime }},\end{array}$ (12)

另外，通过泰勒公式展开有

$\begin{array}{c}{\hat{\mu }}_{\left(i^{\prime }\right)}-{\mu }_{\left(0i^{\prime }\right)}=x_n^{\text{T}}\hat{\beta }+\frac{1}{n_{i^{\prime }}+\hat{b}}{1}_{n_{i^{\prime }}}^{\text{T}}\left(y_{i^{\prime }}-x_{i^{\prime }}\hat{\beta }\right)-x_n^{\text{T}}{\beta }_0-\frac{1}{n_{i^{\prime }}+b_0}{1}_{n_{i^{\prime }}}^{\text{T}}\left(y_{i^{\prime }}-x_{i^{\prime }}{\beta }_0\right)\\ =\frac{\partial {\hat{\mu }}_{\left(i^{\prime }\right)}}{\partial {\beta }^{\prime }}\left(\hat{\beta }-{\beta }_0\right)+\frac{\partial {\hat{\mu }}_{\left(i^{\prime }\right)}}{\partial {\psi }^{\prime }}\left(\hat{\psi }-{\psi }_0\right),\end{array}$ (13)

其中 $\partial {\hat{\mu }}_{\left(i\right)}/\partial {\beta }^{\prime }$ 表示 $\partial {\hat{\mu }}_{\left(i\right)}/\partial \beta$ 在点 ${\left({\beta }^{\prime },{\psi }^{\prime }\right)}^{\prime }$ 上的微分，其中 ${\left({\beta }^{\prime },{\psi }^{\prime }\right)}^{\prime }$ 介于 ${\left({{\beta }^{\prime }}_0,{{\psi }^{\prime }}_0\right)}^{\prime }$ 与 ${\left(\widehat{{\beta }^{\prime }},\widehat{{\psi }^{\prime }}\right)}^{\prime }$ 之间。

通过结合上述的(12)式和(13)式，同时又根据上面理论可知 ${\bar{y}}_n=\theta +{\bar{\epsilon }}_n$。从而有

$\begin{array}{c}{\hat{\mu }}_{\left(i^{\prime }\right)}-{\bar{y}}_n={\hat{\mu }}_{\left(i^{\prime }\right)}-\theta -{\overline{ϵ}}_n={\hat{\mu }}_{\left(i^{\prime }\right)}-{\mu }_{\left(0i^{\prime }\right)}+{\mu }_{\left(0i^{\prime }\right)}-\theta -{\overline{ϵ}}_n\\ ={\alpha }_{i^{\prime }}-{\alpha }_i-{\overline{ϵ}}_n-\frac{b_0}{n_{i^{\prime }}+b_0}{\alpha }_{i^{\prime }}+\frac{1}{n_{i^{\prime }}+b_0}{1}_{n_{i^{\prime }}}^{\text{T}}{ϵ}_{i^{\prime }}+\frac{\partial {\hat{\mu }}_{\left(i^{\prime }\right)}}{\partial {\beta }^{\prime }}\left(\hat{\beta }-{\beta }_0\right)+\frac{\partial {\hat{\mu }}_{\left(i^{\prime }\right)}}{\partial {\psi }^{\prime }}\left(\hat{\psi }-{\psi }_0\right)\\ ={\alpha }_{i^{\prime }}-{\alpha }_i-{\overline{ϵ}}_n+{\xi }_{i^{\prime }},\end{array}$ (14)

其中根据上述的假设有 ${\xi }_{i^{\prime }}\le O_p\left(a_N\right)+O_p\left[\left(c_N+1\right)\cdot b_N\right]+O_p\left(d_N\right)$。

另一方面 $\hat{I}$ 使 $\log {\hat{\sigma }}_{i^{\prime }}^2+{\left({\hat{\mu }}_{i^{\prime }}-{\bar{y}}_n\right)}^2/{\hat{\sigma }}_{i^{\prime }}^2+{\hat{\sigma }}_{ϵ}^2/\left(n_{new}{\hat{\sigma }}_{i^{\prime }}^2\right)$ 在 $1\le i^{\prime }\le m$ 中达到最小，即可表示为

$\begin{array}{c}\hat{I}=\arg {\min }_{1\le i^{\prime }\le m}\left\{\log {\hat{\sigma }}_{i^{\prime }}^2+\frac{{\left({\hat{\mu }}_{i^{\prime }}-{\bar{y}}_n\right)}^2}{{\hat{\sigma }}_{i^{\prime }}^2}+\frac{{\hat{\sigma }}_{ϵ}^2}{n_{new}{\hat{\sigma }}_{i^{\prime }}^2}\right\}\\ =\arg {\min }_{1\le i^{\prime }\le m}\left\{\log {\hat{\sigma }}_{i^{\prime }}^2+\frac{{\left({\hat{\mu }}_{i^{\prime }}-{\bar{y}}_n\right)}^2}{{\hat{\sigma }}_{i^{\prime }}^2}+\frac{\hat{b}+n_{i^{\prime }}}{n_{new}}\right\}\\ =\arg {\min }_{1\le i^{\prime }\le m}\left\{-\log \left(\hat{b}+n_{i^{\prime }}\right)+\frac{{\left({\hat{\mu }}_{i^{\prime }}-{\bar{y}}_n\right)}^2\left(\hat{b}+n_{i^{\prime }}\right)}{{\hat{\sigma }}_{ϵ}^2}+\frac{\hat{b}+n_{i^{\prime }}}{n_{new}}\right\},\end{array}$

因此，有

$\begin{array}{l}-\log \left(\hat{b}+n_{\hat{I}}\right)+\frac{{\left({\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right)}^2\left(\hat{b}+n_{\hat{I}}\right)}{{\hat{\sigma }}_{\epsilon }^2}+\frac{\hat{b}+n_{\hat{I}}}{n_{new}}\\ \le -\log \left(\hat{b}+n_I\right)+\frac{{\left({\hat{\mu }}_I-{\bar{y}}_n\right)}^2\left(\hat{b}+n_I\right)}{{\hat{\sigma }}_{ϵ}^2}+\frac{\hat{b}+n_I}{n_{new}},\end{array}$ (15)

整理得

${\left({\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right)}^2\left(\hat{b}+n_{\hat{I}}\right)-{\left({\hat{\mu }}_I-{\bar{y}}_n\right)}^2\left(\hat{b}+n_I\right)\le {\hat{\sigma }}_{ϵ}^2\left[\log \frac{\hat{b}+n_{\hat{I}}}{\hat{b}+n_I}+\frac{n_I-n_{\hat{I}}}{n_{new}}\right],$

当 $n_i=n_j\left(i=j\right)$ (表示 $n_i$ 是一个固定值，只与参数i相关)时，显然有 ${\left({\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right)}^2\le {\left({\hat{\mu }}_I-{\bar{y}}_n\right)}^2$ ；如果当 $n_i\ne n_j\left(i\ne j\right)$ 时，记 $R_{\max }=\underset{1\le i\le m}{\max }{n}_i$， $R_{\min }=\underset{1\le i\le m}{\min }{n}_i$，如果有 $n_{new}\to \infty$， $\underset{1\le i\le m}{\min }{n}_i\to \infty$， $R_{\max }-R_{\min }\to 0$，则有 ${\left({\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right)}^2\le {\left({\hat{\mu }}_I-{\bar{y}}_n\right)}^2$。故，综上所述有，

${\left({\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right)}^2\le {\left({\hat{\mu }}_I-{\bar{y}}_n\right)}^2.$ (16)

令 ${\xi }_{\hat{I}}$ 和 ${\xi }_I$ 分别表示 ${\xi }_i$ 项中参数i换成其对应参数 $\hat{I}$ 和I；同时为了方便区分，故将 ${\xi }_{i^{\prime }}$ 用 $s_N$ 来表示。则一方面通过(14)式有

$\begin{array}{c}{\left({\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right)}^2={\left({\alpha }_{\hat{I}}-{\alpha }_I-{\overline{ϵ}}_n+{\xi }_{\hat{I}}\right)}^2\\ ={\left({\alpha }_{\hat{I}}-{\alpha }_I-{\overline{ϵ}}_n\right)}^2+{\xi }_{\hat{I}}^2+2\left({\alpha }_{\hat{I}}-{\alpha }_I-{\overline{ϵ}}_n\right)\cdot {\xi }_{\hat{I}}\\ ={\left({\alpha }_{\hat{I}}-{\alpha }_I-{\overline{ϵ}}_n\right)}^2-2{\xi }_{\hat{I}}\cdot {\overline{ϵ}}_n+{\xi }_{\hat{I}}^2+2\left({\alpha }_{\hat{I}}-{\alpha }_I\right){\xi }_{\hat{I}}\\ \ge {\left({\alpha }_{\hat{I}}-{\alpha }_I-{\overline{ϵ}}_n\right)}^2-2{\xi }_{\hat{I}}\cdot {\overline{ϵ}}_n-O_p\left(c_N\right)s_N,\end{array}$ (17)

另一方面，同理，通过(14)式有

$\begin{array}{c}{\left({\hat{\mu }}_I-{\bar{y}}_n\right)}^2={\left({\alpha }_I-{\alpha }_I-{\overline{ϵ}}_n+{\xi }_{\hat{I}}\right)}^2\\ ={\left(-{\overline{ϵ}}_n+{\xi }_I\right)}^2\\ ={\overline{ϵ}}_n^2-2{\xi }_I\cdot {\overline{ϵ}}_n+{\xi }_I^2\\ \le {\overline{ϵ}}_n^2-2{\xi }_I\cdot {\overline{ϵ}}_n+s_N^2,\end{array}$ (18)

结合(16)，(17)和(18)式，有

$\begin{array}{l}{\left({\alpha }_{\hat{I}}-{\alpha }_I-{\bar{\epsilon }}_n\right)}^2-2{\xi }_{\hat{I}}\cdot {\overline{ϵ}}_n-O_p\left(c_N\right)s_N\le {\overline{ϵ}}_n^2-2{\xi }_I\cdot {\overline{ϵ}}_n+s_N^2\\ {\left({\alpha }_{\hat{I}}-{\alpha }_I\right)}^2-2\left({\alpha }_{\hat{I}}-{\alpha }_I\right)\cdot {\overline{ϵ}}_n+{\overline{ϵ}}_n^2\le 2{\xi }_{\hat{I}}\cdot {\overline{ϵ}}_n-2{\xi }_I\cdot {\bar{\epsilon }}_n+{\overline{ϵ}}_n^2+O_p\left(c_N\right)s_N+s_N^2\\ {\left({\alpha }_{\hat{I}}-{\alpha }_I\right)}^2\le 2\left({\alpha }_{\hat{I}}-{\alpha }_I+{\xi }_{\hat{I}}-{\xi }_I\right)\cdot {\overline{ϵ}}_n+O_p\left(c_N\right)s_N+s_N^2,\end{array}$ (19)

根据 ${\hat{\mu }}_{\hat{I}}-\theta ={\hat{\mu }}_{\hat{I}}-{\bar{y}}_n+{\bar{y}}_n-\theta =\left({\alpha }_{\hat{I}}-{\alpha }_I\right)-{\overline{ϵ}}_n+{\xi }_{\hat{I}}+{\overline{ϵ}}_n=\left({\alpha }_{\hat{I}}-{\alpha }_I\right)+{\xi }_{\hat{I}}$，再结合(19)式可得

$\begin{array}{c}{\left({\hat{\mu }}_{\hat{I}}-\theta \right)}^2={\left(\left({\alpha }_{\hat{I}}-{\alpha }_I\right)+{\xi }_{\hat{I}}\right)}^2\\ ={\left({\alpha }_{\hat{I}}-{\alpha }_I\right)}^2+2\left({\alpha }_{\hat{I}}-{\alpha }_I\right)\cdot {\xi }_{\hat{I}}+{\xi }_{\hat{I}}^2\\ \le 2\left({\alpha }_{\hat{I}}-{\alpha }_I+{\xi }_{\hat{I}}-{\xi }_I\right)\cdot {\overline{ϵ}}_n+O_p\left(c_N\right)s_N+s_N^2+2\left({\alpha }_{\hat{I}}-{\alpha }_I\right)\cdot {\xi }_{\hat{I}}+{\xi }_{\hat{I}}^2\\ \le 2\left({\alpha }_{\hat{I}}-{\alpha }_I+{\xi }_{\hat{I}}-{\xi }_I\right)\cdot {\overline{ϵ}}_n+O_p\left(c_N\right)s_N+2s_N^2．\end{array}$ (20)

设 $\hat{\eta }={\xi }_{\hat{I}}-{\overline{ϵ}}_n$ 和 ${\eta }_{\text{1}}={\xi }_{i^{\prime }}-{\overline{ϵ}}_n$，则有 $\left|\hat{\eta }\right|\vee \left|{\eta }_{\text{1}}\right|\le s_N+\left|{\overline{ϵ}}_n\right|=o_p\left(1\right)$。又假设存在 $\mathcal{A}=\left\{\left|\left({\alpha }_{\hat{I}}-{\alpha }_I\right)\right|\ge 1,\left|\hat{\eta }\right|\le 1/4\right\}$，则一方面有

${\left\{{\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right\}}^2={\left\{\left({\alpha }_{\hat{I}}-{\alpha }_I\right)+\hat{\eta }\right\}}^2\ge \left(1/2\right){\left\{\left({\alpha }_{\hat{I}}-{\alpha }_I\right)\right\}}^2.$

另一方面有

${\left\{{\hat{\mu }}_{i^{\prime }}-{\bar{y}}_n\right\}}^2={\left\{\left({\alpha }_{i^{\prime }}-{\alpha }_I\right)+{\eta }_1\right\}}^2\le 2\left[{\left({\alpha }_{i^{\prime }}-{\alpha }_I\right)}^2+{\eta }_1^2\right].$

由此结合(16)式得到

$\begin{array}{l}{\left\{{\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right\}}^2\le {\left\{{\hat{\mu }}_{i^{\prime }}-{\bar{y}}_n\right\}}^2\\ \frac{1}{2}{\left\{\left({\alpha }_{\hat{I}}-{\alpha }_I\right)\right\}}^2\le 2\left[{\left({\alpha }_{i^{\prime }}-{\alpha }_I\right)}^2+{\eta }_1^2\right]\\ {\left\{\left({\alpha }_{\hat{I}}-{\alpha }_I\right)\right\}}^2\le 4\left[{\left({\alpha }_{i^{\prime }}-{\alpha }_I\right)}^2+{\eta }_1^2\right]\\ \left|{\alpha }_{\hat{I}}-{\alpha }_I\right|\le 2\sqrt{{\left({\alpha }_{i^{\prime }}-{\alpha }_I\right)}^2+{\eta }_1^2},\end{array}$

由此可见，当 $\left|\hat{\eta }\right|\le 1/4$ 情况成立时，则有 $\left|\left({\alpha }_{\hat{I}}-{\alpha }_I\right)\right|=O_p\left(1\right)$， $\left|{\alpha }_{\hat{I}}-{\alpha }_I\right|\le 1\vee \left[2\sqrt{{\left({\alpha }_{i^{\prime }}-{\alpha }_I\right)}^2+{\eta }_1^2}\right]=O_p\left(1\right)$。因此，回到最初的问题中，即由(20)式可知 ${\left({\hat{\mu }}_{\hat{I}}-\theta \right)}^2\le 2\left[O_p\left(1\right)+o_p\left(1\right)\right]o_p\left(1\right)+O_p\left(c_N\right)s_N+2s_N^2=o_p\left(1\right)$。从而，定理得证。

2.2. 新观测数据与训练集数据不匹配

在不假设匹配关系存在的情况下，则匹配关系可能存在也有可能不存在，即不存在新观测数据与训练集数据中的某一组群相匹配的前提，这比较符合实际。如果匹配关系存在则可按2.1节推导；当匹配关系不存在时，即 $I\ne i\in \left\{1,\cdots ,m\right\}$，则新观测数据与训练集数据是相互独立的。接下来讨论匹配关系不存在的情况。

如果不匹配，则此时的混合效应为(3)式，从而有 $\left(\theta ,{ϵ}_n\right)$ 与训练集数据相互独立，故可得的最佳预测为，

$\tilde{\theta }=E\left(\theta |y\right)=E\left(\theta |y_1,\cdots ,y_m\right)=E\left(\theta \right),$

同理，与匹配情况一样，BP为，

$\tilde{\mu }=\tilde{\theta }=x_n^{\text{T}}\beta ,$ (21)

其 $\theta |y$，方差为

$\begin{array}{l}{\sigma }^2=Var\left(\theta |y\right)=Var\left(\theta |y_1,\cdots ,y_m\right)=Var\left(\theta \right)\\ =Var\left(x_n^{\text{T}}\beta +{\alpha }_I\right)\\ =Var\left(x_n^{\text{T}}\beta \right)+Var\left({\alpha }_I\right)+2Cov\left(x_n^{\text{T}}\beta ,{\alpha }_I\right)\\ \text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}={\sigma }_{n,\alpha }^2.\end{array}$ (22)

将(21)的参数 $\beta$ 用它的一致估计 $\hat{\beta }$ 代替，则可 $\theta$ 得到相应的经验最佳预测 $\hat{\mu }=x_n^{\text{T}}\hat{\beta }$。

同样，将(22)式中的参数 ${\sigma }_{n,\alpha }^2$ 换成其对应的一致估计 ${\hat{\sigma }}_{n,\alpha }^2$，则有 $\theta |y$，的方差估计为： ${\hat{\sigma }}^2={\hat{\sigma }}_{n,\alpha }^2$。

综上所述，可得 $\theta |y$ 表达形式为： $\theta |y~N\left(\mu ,{\sigma }^2\right)$。从而， $\theta |y$ 的密度函数可以表示为

$f\left(\theta |y\right)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{{\left(\theta -\mu \right)}^2}{2{\sigma }^2}\right\}.$ (23)

同理，这种情况下，则对数评分(Logs)表达式为

$\log \left\{f\left(\theta |y\right)\right\}\propto -\left\{\log {\sigma }^2+\frac{{\left(\mu -\theta \right)}^2}{{\sigma }^2}\right\}.$ (24)

与匹配情况中的(10)式的推导过程相同，可得出此时的SPSRs为：

$\begin{array}{c}\text{SPSRs}=E\left\{\log {\hat{\sigma }}^2+\frac{{\left(\hat{\mu }-\theta \right)}^2}{{\hat{\sigma }}^2}\right\}\\ =E\left\{\log {\hat{\sigma }}^2+\frac{{\left(\hat{\mu }-{\bar{y}}_n+{\overline{ϵ}}_n\right)}^2}{{\hat{\sigma }}^2}\right\}\\ =E\left\{\log {\hat{\sigma }}^2+\frac{{\left(\hat{\mu }-{\bar{y}}_n\right)}^2}{{\hat{\sigma }}^2}+\frac{2\left(\hat{\mu }-{\bar{y}}_n\right){\overline{ϵ}}_n}{{\hat{\sigma }}^2}+\frac{{\left({\overline{ϵ}}_n\right)}^2}{{\hat{\sigma }}^2}\right\}\\ =E\left\{\log {\hat{\sigma }}^2+\frac{{\left(\hat{\mu }-{\bar{y}}_n\right)}^2}{{\hat{\sigma }}^2}+\frac{{\sigma }_{ϵ}^2}{\left(n_{new}{\hat{\sigma }}^2\right)}\right\}.\end{array}$ (25)

将(10)式和(16)式比较，区别是：将 ${\hat{\mu }}_i$ 换成 $\hat{\mu }$ ； ${\hat{\sigma }}_i^2$ 换成 ${\hat{\sigma }}^2$。故将之前的方法做如下延伸：令 $\hat{I}$ 由(3.11)式给出，比较 $\log {\hat{\sigma }}_{\hat{I}}^2+{\left({\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right)}^2/{\hat{\sigma }}_{\hat{I}}^2+{\sigma }_{ϵ}^2/\left(n_{new}{\hat{\sigma }}_{\hat{I}}^2\right)$ 和 $\log {\hat{\sigma }}^2+{\left(\hat{\mu }-{\bar{y}}_n\right)}^2/{\hat{\sigma }}^2+{\sigma }_{ϵ}^2/\left(n_{new}{\hat{\sigma }}^2\right)$ 的大小。如果前者更小，则 $\theta$ 的CMEP为 ${\hat{\mu }}_{\hat{I}}$ ；否则， $\theta$ 的CMEP是 $\hat{\mu }$。

3. 数值模拟分析

本节主要对SPSRs准则下CMEP的预测效果进行分析，主要将CMEP-SPSRs方法与PR方法的MSPE结果进行比较，讨论其预测效果。针对新观测数据与训练数据是否匹配两种情形分别讨论，具体步骤如下：

1) 新观测数据与训练集数据存在匹配关系

步骤1：训练集数据 $y_{ij}$ 按照(1)式生成，其中 $n_i=5,1\le i\le m$。给定 $\beta ={\left(5,1\right)}^{\text{T}}$， $x_{ij}={\left(1,x_{ij2}\right)}^{\text{T}}$，协变量 $x_{ij2}$ 服从 $N\left(0,1\right)$ 的随机数；群组特定的随机效应 ${\alpha }_i$ 服从 $N\left(0,{\sigma }_{\alpha }^2\right)$ 的随机数，其中给定 ${\sigma }_{\alpha }^2$ ；误差项 ${ϵ}_{ij}$ 服从 $N\left(0,{\sigma }_{ϵ}^2\right)$ 的随机数，其中给定 ${\sigma }_{ϵ}^2$。

步骤2：新观测数据 $y_{n,j}$ 按照(2)式生成，其中给定 $\beta ={\left(5,1\right)}^{\text{T}}$， $x_n={\left(1,x_{n2}\right)}^{\text{T}}$，协变量 $x_{n2}$ 服从 $N\left(0,1\right)$ 的随机数；将(1)式中 ${\alpha }_i$， ${ϵ}_{ij}$ 分别替换为 ${\alpha }_I$， ${ϵ}_{n,j}$。故需要预测的混合效应为(3)式。

步骤3：先基于(4)式可获得 $\theta$ 的最佳预测 ${\tilde{\theta }}_{\left(i\right)}$，同时又根据(5)式可得到 $\theta |y$ 的方差为 ${\sigma }_i^2$ ；然后得到参数的最大似然估计， $\hat{\beta },{\hat{\sigma }}_{\alpha }^2,{\hat{\sigma }}_{ϵ}^2$。最后将参数估计代回(4)与(5)式，则分别可获得 $\theta$ 的经验最佳预测 ${\hat{\mu }}_i$ 和 $\theta |y$ 方差估计 ${\hat{\sigma }}_i^2$。

步骤4：通过(11)式获得 ${\hat{I}}_{LogS}$。接着将 ${\hat{I}}_{LogS}$ 替换(6)式中的i，则得到 $\theta$ 的分类混合效应预测(CMEP)： $\hat{\theta }={\hat{\mu }}_{\hat{I}}$。

步骤5：记 ${\hat{\theta }}_{n,r}=x_n^{\text{T}}{\hat{\beta }}_{LS}$ 表示 ${\theta }_n$ 的标准回归预测(RP)，其中 ${\hat{\beta }}_{LS}$ 表示为 $\beta$ 的最小二乘(Least-Square, LS)估计。因此要计算出 ${\hat{\beta }}_{LS}$，从而得到RP。

步骤6：分别计算CMMP-SPSRs方法和RP方法的MSPE： ${\text{MSPE}}_1={\left(\hat{\theta }-\theta \right)}^2$， ${\text{MSPE}}_2={\left({\hat{\theta }}_{n,r}-\theta \right)}^2$。

通过Matlab软件按以上步骤进行编程，除步骤2以外，所有步骤重复100次，然后获得CMMP-SPSRs方法和RP方法预测混合效应的MSPE平均值，这样有利于降低误差，最后将MSPE的平均值作为评价预测准确性的指标。结果如表1~3中所示。%Improve为CMMP-SPSRs方法和RP方法MSPE的相对大小。

Table 1. Comparison of two MSPE methods with matching relationship, fixed at m = 50 , n n e w = 5 , σ ϵ 2 = 1

表1. 存在匹配关系两种方法MSPE的比较，固定 $m=50,n_{new}=5,{\sigma }_{ϵ}^2=1$

Table 2. Comparison of two MSPE methods with matching relationship, fixed at m = 50 , n n e w = 5 , σ α 2 = 1

表2. 存在匹配关系两种方法MSPE的比较，固定 $m=50,n_{new}=5,{\sigma }_{\alpha }^2=1$

Table 3. Comparison of two MSPE methods with matching relationship, fixed at m = 50 , σ ϵ 2 = 1 , σ α 2 = 1

表3. 存在匹配关系两种方法MSPE的比较，固定 $m=50,{\sigma }_{ϵ}^2=1,{\sigma }_{\alpha }^2=1$

$\text{\%Improve}=\left(\frac{\text{MSPEofRP}-\text{MSPEofCMMP}}{\text{MSPEofCMMP}}\right)\times 100\text{\%}\text{.}$

从表1~3中数据可以发现RP方法的MSPE都大于CMMP-SPSRs方法的MSPE。具体来说，表1：当只有随机效应方差 ${\sigma }_{\alpha }^2$ 改变时，随着 ${\sigma }_{\alpha }^2$ 值增加，两者方法的MSPE都有增大，但通过%Improve可以发现RP方法的MSPE增大的速度高于CMMP-SPSRs方法；同时在MSPE方面上，发现CMMP-SPSRs方法小于RP方法，故CMMP-SPSRs方法优于RP方法。表2：当只有误差项方差 ${\sigma }_{ϵ}^2$ 改变时，随着 ${\sigma }_{ϵ}^2$ 值增加，CMMP-SPSRs方法和RP方法的MSPE也增加，但CMMP-SPSRs方法始终小于RP方法，即CMMP-SPSRs方法更优。表3：当只有 $n_{new}$ 改变时，从数据中可以看出，CMMP-SPSRs方法混合效应预测的准确性随着 $n_{new}$ 值增加而提高，这符合实际，当新观测数据增多，从而预测的准确性也会增加。RP方法并没有因为 $n_{new}$ 的改变而发生明显的变化，且MSPE值一直比CMMP-SPSRs方法大，因此CMMP-SPSRs方法优于RP方法。

2) 新观测数据与训练集数据不匹配

步骤1：训练集数据 $y_{ij}$ 可以按照模型(1)式生成，其中 $n_i=5,1\le i\le m$。同时给定 $\beta ={\left(1,2,3\right)}^{\text{T}}$， $x_{ij}={\left(1,x_{1,ij},x_{2,ij}\right)}^{\text{T}}$，协变量 $x_{1,ij},x_{2,ij}$ 由 $N\left(0,1\right)$ 分布产生的随机数；群组特定的随机效应 ${\alpha }_i$ 服从 $N\left(0,{\sigma }_{\alpha }^2\right)$ 的随机数，其中给定 ${\sigma }_{\alpha }^2$ ；误差项 ${ϵ}_{ij}$ 服从 $N\left(0,{\sigma }_{ϵ}^2\right)$ 的随机数，其中给定 ${\sigma }_{ϵ}^2$。

步骤2：新观测数据 $y_{n,j}$ 可以按照模型(2)式生成，其中同时给定 $\beta ={\left(1,2,3\right)}^{\text{T}}$， $x_n={\left(1,x_{1,n},x_{2,n}\right)}^{\text{T}}$，协变量 $x_{1,n},x_{2,n}$ 由 $N\left(0,1\right)$ 分布产生的随机数；将(1)式中 ${\alpha }_i$， ${ϵ}_{ij}$ 分别替换为 ${\alpha }_I$， ${ϵ}_{n,j}$。故需要预测的混合效应为(3)式。

步骤3：首先基于(21)式获得 $\theta$ 的最佳预测 $\tilde{\theta }$，同时通过(22)式得到 ${\sigma }^2$。然后得到参数的MLE。最后可获得 $\theta$ 的经验最佳预测 $\hat{\mu }$ 和 ${\hat{\sigma }}^2$。

步骤4：与匹配关系中的步骤5相同，得到RP。

步骤5：基于(1)匹配情况下得到 $\hat{I}$ 和 $\theta$ 的分类混合效应预测： $\hat{\theta }={\hat{\mu }}_{\hat{I}}$。同时计算 $\log {\hat{\sigma }}_{\hat{I}}^2+{\left({\hat{\mu }}_{\hat{I}}-{\bar{y}}_n\right)}^2/{\hat{\sigma }}_{\hat{I}}^2+{\sigma }_{ϵ}^2/n_{new}{\hat{\sigma }}_{\hat{I}}^2$ 和 $\log {\hat{\sigma }}^2+{\left(\hat{\mu }-{\bar{y}}_n\right)}^2/{\hat{\sigma }}^2+{\sigma }_{ϵ}^2/n_{new}{\hat{\sigma }}^2$ 的大小。如果前者较小，则 $\theta$ 的分类混合效应预测为 ${\hat{\mu }}_{\hat{I}}$ ；否则， $\theta$ 的分类混合效应预测是 $\hat{\mu }$。

步骤6：与(1)匹配中的步骤6一样。

通过表4和表5中的数据比较可以发现，在不匹配情况下，两者的MSPE的变化趋势和匹配情况相类似，即CMMP-SPSRs方法相对RP方法依旧存在优越，同时%Improve的变化趋势也是相似的。由此，可以得出不管新观测数据和训练集数据中是否存在匹配关系，CMMP-SPSRs方法相对于RP方法预测都更具有优越性。

4. 实例应用分析

本节以电视学校和家庭预防与戒烟项目(Television School and Family Smoking Prevention and Cessation Project, TVSFP) [11] 的数据来对CMMP-SPSRs方法验证，其中研究对象是来自美国加利福尼亚州洛杉矶和圣地亚哥学校的七年级学生。最初该项目主要以学校为基础的社会抵抗课程和以电视为基础的预防和戒烟方面的独立与联合效果研究。为了对嵌套误差回归模型的演示，选择TVSFP数据中的一个子集，即取洛杉矶中的28所学校，这些学校被随机分配为四种研究条件中的一种：1) 抵制社会的课程(social-resistance classroom curriculum, CC)；2) 电视干预(television intervention, TV)；3) 结合CC和TV；4) 不做任何处理。烟草与健康知识量表(Tobacco and Health knowledge scale, THKS)得分是主要研究结果变量之一，作为本次实验的响应变量。THKS是由七个问卷项目组成，用于评估学生的烟草和健康知识。目前数据只涉及干预前和干预后时间点完成的THKS。该数据是来自28所学校的1600名学生，每所学校有1至13间教室，每个教室有2至28名学生。该数据可在http://www.hsph.harvard.edu/fitzmaur/ala/tvsfp.txt上查找。

Table 4. Comparison of two MSPE methods without matching relationship, fixed at m = 50 , n n e w = 5 , σ ϵ 2 = 1

表4. 不存在匹配关系两种方法MSPE的比较，固定 $m=50,n_{new}=5,{\sigma }_{ϵ}^2=1$

Table 5. Comparison of two MSPE methods without matching relationship, fixed at m = 50 , n n e w = 5 , σ α 2 = 1

表5. 不存在匹配关系两种方法MSPE的比较，固定 $m=50,n_{new}=5,{\sigma }_{\alpha }^2=1$

将这28所学校的数据作为训练集数据的一个子集。因为并不知道训练集(总体)数据，所以假设总体数据为相应学校数据的10倍(重复)，故而，相应学校的总体均值与样本均值相同的。在28所学校中选择其中一所学校作为新观测数据，但并不知道这新观测数据来自具体哪一个学校，故而需要对其进行估计。接下来分别采用CMMP-SPSRs方法和RP方法对新观测的混合效应(均值)进行预测。最后为28所学校中的每所都重复以上操作。

假设训练集数据符合嵌套误差回归模型为

$y_{ij}={\beta }_0+{\beta }_1{x}_{ij1}+{\beta }_2{x}_{ij2}+{\beta }_3{x}_{ij3}+{\beta }_4{x}_{ij4}+{\alpha }_i+{\epsilon }_{ij},$ (26)

$i=1,\cdots ,m$， $j=1,\cdots ,n_i$， $m=28$， $n_i$ 为2至28不等。其中 $y_{ij}$ 是干预后的THKS得分， $x_{ij1}$ 为干预前的THKS得分，即不做任何处理的数据； $x_{ij2}$ 为CC的数据； $x_{ij3}$ 为TV的数据； $x_{ij4}$ 为结合CC和TV的数据。 $\beta ={\left({\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4\right)}^{\text{T}}$ 是固定效应且为未知参数，其中 ${\alpha }_i,{\epsilon }_{ij}$ 和(1)式中的假设相同。参数 $\psi ={\left(\beta ,{\sigma }_{\alpha }^2,{\sigma }_{\epsilon }^2\right)}^{\text{T}}$ 的最大似然估计可由(26)式得到，记为 $\hat{\psi }={\left(\hat{\beta },{\hat{\sigma }}_{\alpha }^2,{\hat{\sigma }}_{\epsilon }^2\right)}^{\text{T}}$，同时通过计算也可获得参数 $\beta$ 的最小二乘估计记为 ${\beta }_{LS}$。依次选择一所学校的模拟样本作为新观测数据，28所学校则作为训练集数据。假设新观测数据也符合嵌套误差回归模型(NER)，则需要预测的混合效应为 $\theta ={\beta }_0+{\beta }_1{x}_{ij1}+{\beta }_2{x}_{ij2}+{\beta }_3{x}_{ij3}+{\beta }_4{x}_{ij4}+{\alpha }_{Ｉ}$。

接下来分别用CMMP-SPSRs方法和RP方法对新观测数据的混合效应进行预测，在预测准确性方面比较了这两种方法。结果如表6所示，其中CMMP-SPSRs和RP列表示为绝对预测误差，即，

$\left|\hat{\theta }-n_i^{-1}{\displaystyle {\sum }_{j=i}^{n_i}{y}_{ij}}\right|$

$\left|{\hat{\theta }}_{n,r}-n_i^{-1}{\displaystyle {\sum }_{j=i}^{n_i}{y}_{ij}}\right|$

其中 $\hat{\theta }$ 表示CMMP-SPSRs方法的混合效应预测总体均值， ${\hat{\theta }}_{n,r}$ 表示RP方法下的混合效应预测总体均值， ${\bar{y}}_i=n_i^{-1}{\displaystyle {\sum }_{j=i}^{n_i}{y}_{ij}}$ 表示为总体均值。

Table 6. CMMP vs RP for TVSFP data

表6. TVSFP数据的CMMP vs RP

通过表6中的数据可以看出，在这28所学校中，CMMP-SPSRs方法的预测误差都比RP方法小，换句话讲，CMMP-SPSRs方法的预测准确性比RP方法更优。改进的百分比由4.6%到329.8%不等。新观测数据和训练集数据间的匹配关系可能存在也有可能并不存在，本文通过匹配关系判别策略，利用训练集数据中尽量相似的群组中信息来提高预测的准确性。

5. 总结

本文主要讨论了嵌套误差模型中分类混合效应预测问题，主要对CMMP方法在识别准则的预测理论和方法进行改进，总结如下：对参数I的匹配准则方法的改进，在已有的分类混合模型CMMP方法，其匹配准则采用的方法为均方预测误差。为了优化CMMP方法，故提出了新的匹配准则方法：SPSRs。本文主要对CMMP-SPSRs方法和RP方法进行了比较，通过大量数值模拟可看出，当匹配关系存在或者不存在，CMMP-SPSRs方法对于回归预测方法依然保持着很好的预测效果。

下一步的工作是在对待预测效应进行识别时，只讨论了两种可能，即新的样本属于某个小域或者不属于任何小域。结合模型结构误定和SPSRs准则，我们可以进一步考虑待预测样本按一定概率(或者模糊隶属度)落在小域中的相应预测方法。