Bootstrap模型平均置信区间

Model-Averaged Bootstrap Confidence Interval

摘要: 模型平均作为能够解决模型不确定性的统计推断方法越来越受欢迎，模型平均应用也从点估计扩展到区间估计。目前模型平均置信区间常用方法包括Wald区间和MATA-Wald区间，这些区间都是在参数估计的分布渐近正态的情况下构建的，所以这两种方法在拟合正态数据时效果是非常好的。在本文中，我们将参数估计的Bootstrap分布来近似其真实分布，进而提出了Bootstrap模型平均置信区间，并对这3种方法分别在正态数据和偏态数据的模型下的表现进行模拟比较。模拟结果表明在正态数据的模型下，Bootstrap模型平均置信区间在上错误率和覆盖率方面仅略差于MATA区间。在偏态数据的模型下，百分位数Bootstrap模型平均置信区间的覆盖率要优于Wald区间和MATA区间，并且提供了较好的错误率，特别是上错率，都更接近名义水平。所以Boorstrap区间表现出很好的适应性。

Abstract: Model averaging is becoming more and more popular as a statistical inference method to solve model uncertainty, and the application of model averaging is also extended from point estimation to interval estimation. At present, the common methods of model averaged confidence interval include Wald interval and MATA-Wald interval. These intervals are constructed under the condition that the distribution of parameter estimation is asymptotically normal. Therefore, these two methods are very good in fitting normal data. In this paper, the bootstrap distribution of parameter estimation was approximated to the real distribution, and then the model-averaged bootstrap interval was proposed. The performance of the three methods under normal and skewered data models was simulated and compared respectively. The simulation results show that the model-averaged bootstrap interval is only slightly worse than that of the MATA model in the upper error rate and coverage under the normal data model, Under the model with skewed data, the coverage of the model-averaged bootstrap interval is better than that of the Wald interval and the MATA interval, and it also provides a better error rate, especially the upper error rate, which is closer to the nominal level. So the bootstrap interval shows good adaptability.

1. 引言

传统上，参数估计是基于单一模型的假设，该模型的选择是一个模型选择的过程。这种方法忽视了模型的不确定性 [1]，导致参数估计的过拟合。

在许多应用领域进行模型参数估计时，考虑模型不确定性的情况越来越普遍。这种不确定性之所以会出现，是因为我们不能确保模型选择过程总是会选择相同的最佳模型。传统上，经过模型选择过程挑选出最佳模型之后，再使用挑选出的最佳模型来进行参数估计。该参数的估计依赖于这个最佳模型的选择。由于模型选择的方法一般依赖于挑选出的最佳模型，这样就会丢弃其他模型所特有的信息。这种方法忽视了模型的不确定性。为了解决模型选择带来的这些问题，需要考虑在某种程度上接受所有的候选模型。考虑到这个问题的重要性，人们开始深入学习并研究模型平均方法，其作为一种简单、通用、易理解的方法，应用于多个模型的研究过程。这解决了出现的模型不确定性问题。模型平均方法的出现起初更多应用在点估计上，后来随着理论研究的深入，区间估计也开始使用这种方法。

模型平均有Beyes模型平均 [2] 和频率模型平均 [3] 两方面，本文主要介绍频率模型平均。在频率方法

的设置中，模型平均估计 $\hat{\theta }$ 定义为单个模型估计的加权和： $\hat{\theta }={\displaystyle {\sum }_{i=1}^R{\omega }_i\widehat{{\theta }_i}}$ 其中 $\widehat{{\theta }_i}$ 是模型 $M_i$ 下的估计，模型权重 ${\omega }_i$ 是由信息准则决定的，例如AIC [4]，和式是基于R个候选模型的集合。

目前已经提出了一些构造频率模型平均置信区间的方法。形式为 $\hat{\theta }\pm z_{\alpha }\cdot \widehat{se}\left(\hat{\theta }\right)$ 的Wald区间，其中 $z_{\alpha }$ 是标准正态分布的 $1-\alpha$ 分位数。其依赖于 $\hat{\theta }$ 的标准差 $se\left(\hat{\theta }\right)$ 的精确估计，Hjort和Claeskens [5] 表明 $se\left(\hat{\theta }\right)$ 估计的困难性，以及Wald区间在覆盖率方面表现很差。Buckland [6] [7] 等人提出了 $\widehat{se}\left(\hat{\theta }\right)$ 的一些选择。

另一个可选的构造频率模型平均置信区间的方法是MATA-Wald [8] 方法，由Turek和Fletcher [7] 提出。在MATA方法中，每个置信限度定义为满足单个模型Wald区间错误率的和等于期望错误率的结果，MATA在小样本时的覆盖率方面要强于现有的模型平均的方法。Fletcher和Turek [9] 将MATA的构造方式应用到轮廓似然区间中，提出了模型平均轮廓似然区间，轮廓似然区间在偏态数据的拟合上表现出明显的优势。Kabaila [10] [11] [12] 等人研究了MATA置信区间的一些性质。Yu [13] 等人提出了MATA区间的转换版本。

本文利用单模型下参数的Bootastrap置信区间以及类比MATA-Wald模型平均置信区间的构造原理，构造了Bootstrap模型平均置信区间，通过模拟计算来说明其优势。第二部分，介绍了Wald区间的构造方法以及相应的标准差的估计、MATA-Wald模型平均置信区间的方法及Bootstrap模型平均置信区间的构造方法。第三部分，分别在响应变量是正态和偏态的情形下进行了模拟研究，从上下错误率、覆盖率和平均区间长度分别比较了Wald区间、MATA区间和Bootstrap区间的表现性能。第四部分，对本文的研究进行了总结。

2. 方法及思想

假设考虑R个候选模型 $\left\{M_1,M_2,\cdots ,M_R\right\}$，感兴趣的参数 $\theta$ 对所有模型都有相同的解释意义。假设模型平均估计 $\hat{\theta }={\displaystyle {\sum }_i{\omega }_i\widehat{{\theta }_i}}$ 近似服从正态分布，其中 ${\omega }_i$ 和 $\widehat{{\theta }_i}$ 分别对应模型 $M_i$ 的权重和参数 $\theta$ 的估计。计算参数 $\theta$ 的 $100\left(1-\alpha \right)\%$ Wald置信区间形式为

$\hat{\theta }\pm z_{\alpha }\cdot \widehat{se}(\; \theta \; ^\; )$

其中 $z_{\alpha }$ 是标准正态分布的 $100\left(1-\alpha \right)\%$ 分位数， $\widehat{se}\left(\hat{\theta }\right)$ 是 $\hat{\theta }$ 方差的估计。尽管这种方法很简单，但是计算一个合适的 $se\left(\hat{\theta }\right)$ 的估计是很困难的。由Burnham [5] 和Anderson等人提出的两种 $\widehat{se}\left(\hat{\theta }\right)$ 估计为

$\widehat{se}\left(\hat{\theta }\right)={\displaystyle {\sum }_i{\omega }_i}\sqrt{{\left(\widehat{se}\left(\hat{\theta }\right)\cdot t_{{\nu }_i}/z_{\alpha }\right)}^2+{\left(\widehat{{\theta }_i}-\hat{\theta }\right)}^2}$ (2.1)

$\widehat{se}\left(\hat{\theta }\right)={\displaystyle {\sum }_i\sqrt{{\omega }_i\left[{\left(\widehat{se}\left(\hat{\theta }\right)\cdot t_{{\nu }_i}/z_{\alpha }\right)}^2+{\left(\widehat{{\theta }_i}-\hat{\theta }\right)}^2\right]}}$ (2.2)

$t_{{\nu }_i}$ 是自由度为 ${\nu }_i$ 的t分布的分位数， ${\nu }_i$ 是与模型 $M_i$ 相关的残差的自由度。由Cauchy-Schwarz不等式可知， ${\widehat{se}}_2\ge {\widehat{se}}_1$，则根据(2.2)式得到的Wald置信区间比由(2.1)式得到的Wald置信区间长度宽，故其覆盖率大，本文我们只考虑(2.2)式得到的Wald置信区间。Wald区间在覆盖率方面表现得很差，并且忽视了 $\hat{\theta }$ 的非正态性 [8]，下面基于单个模型来获取参数 $\theta$ 的模型平均置信区间。

假设恰有一个候选模型是真实的，记 $c_i$ 为未知的指示变量，如果模型 $M_i$ 为真，则 $c_i=1$ 否则 $c_i=0$，假设我们知道 $\widehat{{\theta }_i}$ 的真实分布记为 $G_i(\cdot )$，参数 $\theta$ 的置信区间 $\left[{\theta }_L,{\theta }_U\right]$ 满足下列方程

$\{\begin{array}{l}{\displaystyle {\sum }_{i=1}^R{c}_i{G}_i\left({\theta }_L\right)}=\alpha \\ {\displaystyle {\sum }_{i=1}^R{c}_i{G}_i\left({\theta }_U\right)}=1-\alpha \end{array}$ (2.3)

因为 $c_i$ 是未知的，我们用与模型相对应的基于某种信息准则的模型权重 ${\omega }_i$ 来估计(2.3)中的 $c_i$，则参数 $\theta$ 的置信区间 $\left[{\theta }_L,{\theta }_U\right]$ 满足下列方程

$\{\begin{array}{l}{\displaystyle {\sum }_{i=1}^R{\omega }_i{G}_i\left({\theta }_L\right)}=\alpha \\ {\displaystyle {\sum }_{i=1}^R{\omega }_i{G}_i\left({\theta }_U\right)}=1-\alpha \end{array}$ (2.4)

如果不知道 $\widehat{{\theta }_i}$ 的真实分布，假设每个模型的学生化的参数估计 $\frac{\widehat{{\theta }_i}-\theta }{\widehat{se}\left(\widehat{{\theta }_i}\right)}$ 服从标准正态分布或t分布，即 $G_i(\cdot )$ 为标准正态分布或者t分布的分布函数，获取的参数 $\theta$ 的 $100\left(1-\alpha \right)\%$ 的置信区间即为z版本或者

t版本的MATA区间，本文我们考虑t版本的MATA区间 [8]。

但是如果相应变量是偏态的，或者正态假设不明显时，使用Wald区间和MATA区间会导致比较差的覆盖率。本文利用参数Bootstrap 方法构造变量的分布，即用参数估计 $\widehat{{\theta }_i}$ 的Bootstrap分布 $\widehat{G_i}(\cdot )$ 来近似(2.4)中的分布 $G_i(\cdot )$，从而计算参数的模型平均置信区间。

设样本 $X=\left(X_1,X_2,\cdots ,X_n\right)$ 来自模型 $M_i$ 的分布 $F_i\left(X,\theta \right)$ 容量为n的样本， $\widehat{{\theta }_i}$ 是参数 $\theta$ 的极大似然估计。下面给出计算参数模型平均估计的计算步骤。

l 对于来自分布为 $F_i\left(X,\theta \right)$ 的样本，选取Bootstrap样本量B，比如 $B=2000$ ；

n 对于 $i=1,\cdots ,B$，从 $F_i\left(X,\widehat{{\theta }_i}\right)$ 产生容量为n的Bootstrap样本；

n 计算出 ${\theta }_i$ 的极大似然估计 ${\hat{\theta }}_i^{\left(j\right)}$ ；

n 由 $\left({\hat{\theta }}_i^{\left(1\right)},\cdots ,{\hat{\theta }}_i^{\left(B\right)}\right)$ 给出 $\widehat{{\theta }_i}$ 的Bootstrap分布 $\widehat{G_i}(\cdot )$。

求 ${\theta }_L$ 和 ${\theta }_U$ 满足下列方程

$\{\begin{array}{l}{\displaystyle {\sum }_{i=1}^R{\omega }_i\widehat{G_i}\left({\theta }_L\right)}=\alpha /2\\ {\displaystyle {\sum }_{i=1}^R{\omega }_i\widehat{G_i}\left({\theta }_U\right)}=1-\alpha /2\end{array}$ (2.5)

其中 ${\omega }_i$ 为某种信息准则下模型 $M_i$ 的权重。那么 $\left[{\theta }_L,{\theta }_U\right]$ 为参数 $\theta$ 的Bootstrap模型平均置信区间。

3. 数据模拟

本节对Bootstrap模型平均置信区间的优良性进行评判，并与已有的与Wald区间，MATA区间进行比较，最终从区间的频率性质和区间长度这两方面进行评价。区间的频率性质用以下两个指标：覆盖率；上下错误率。

进行模型选择时，权重的选择是基于信息准则AIC

${\omega }_i=\frac{\exp \left(-{\text{AIC}}_i/2\right)}{{\displaystyle {\sum }_{j=1}^R\exp \left(-{\text{AIC}}_j/2\right)}}$

其中 ${\text{AIC}}_i$ 为第i个模型的AIC。

3.1. 线性模型

考虑正态线性模型，包括两个预测变量 $x_1$ 和 $x_2$，5个候选模型 $\left\{M_1,M_2,\cdots ,M_5\right\}$。所有模型的响应变量 $y_i~N\left({\mu }_i,{\sigma }^2\right)$，并且 ${\mu }_i$ 指定如下：

$M_1:{\mu }_i={\beta }_0$

$M_2:{\mu }_i={\beta }_0+{\beta }_1{x}_{1,i}$

$M_3:{\mu }_i={\beta }_0+{\beta }_2{x}_{2,i}$

$M_4:{\mu }_i={\beta }_0+{\beta }_1{x}_{1,i}+{\beta }_2{x}_{2,i}$

$M_5:{\mu }_i={\beta }_0+{\beta }_1{x}_{1,i}+{\beta }_2{x}_{2,i}+{\beta }_{12}{x}_{1,i}{x}_{2,i}$

随意设置 $\sigma =1$，选择 ${\beta }_0=1,{\beta }_1={\beta }_2=0.3,{\beta }_{12}=0.1$，这样所有的模型都倾向于依据样本容量的范围分配非平凡权重。 $x_1$ 和 $x_2$ 是在每轮模拟中随机产生的，并且 $X_{1,i}~N\left(0,1\right),X_{2,i}~Gamma\left(\alpha =2,\beta =1\right)$，样本量大小分别取 $n=20,25,\cdots ,100$ 这17种情形(见表1)。对于感兴趣的参数选取为预测空间 $\left(x_1,x_2\right)$ 各个点的 $\mu$ 值，为了简单起见，这里只考虑 $x_1$ 生成分布的50%处分位数和 $x_2$ 生成分布的90%处分位数。这里模拟次数为10⁵次。在计算Bootstrap模型平均置信区间时，Bootstrap样本量取 $B=2000$。

图1展示了响应变量均值的三种区间的上下错误率和覆盖率随样本容量n的变化形势，图1的左图上方为三种区间的上错误率。为了表示方便，下错误率取负值后展示在图1的左图下方。表1展示了3种区间估计的平均区间长度。由图1和表1可以看出：

1) 下错误率方面，Bootstrap区间要好于其它两种区间，更接近名义错误率。上错率和覆盖率方面，Bootstrap区间整体要好于Wald区间，略逊于MATA区间。

2) 平均区间长度方面，Bootstrap区间是最窄的。

3) 随着样本量的增加，三种区间估计方法都逐渐接近名义覆盖率，并且平均区间长度区别越来越小。

Figure 1. Error rate and coverage of Bootstrap (red triangle), MATA (green square) and Wald (blue circle) interval estimates of linear model mean

图1. 线性模型均值的Bootstrap (红色三角)、MATA (绿色方形)、Wald (蓝色圆形)区间估计的上下错误率、覆盖率

Table 1. Average interval lengths of three interval estimates of linear model mean

表1. 线性模型均值的三种区间估计的平均区间长度

Figure 2. Upper and lower error rates and coverage of Bootstrap (red triangle), MATA (green square) and Wald (blue circle) interval estimates of the mean of lognormal model

图2. 对数正态模型均值的Bootstrap (红色三角)、MATA (绿色方形)、Wald (蓝色圆形)区间估计的上下错误率、覆盖率

Table 2. Average interval lengths of three interval estimates of lognormal model when n = 50

表2. n = 50时对数正态模型的三种区间估计的平均区间长度

3.2. 对数正态模型

下面考虑对数正态的偏态数据情形。设观察值 $Y_{ij}$ 满足 $\log \left(Y_{ij}\right)~N\left({\beta }_i,{\sigma }^2\right)$，我们希望去估计 $Y_{ij}$ 期望 $E\left(Y_{ij}\right)\equiv {\mu }_i=\exp \left({\beta }_i+\frac{{\sigma }^2}{2}\right)\left(i=1,2;j=1,\cdots ,n\right)$。模型平均基于以下两个模型：

$M_1:{\beta }_1={\beta }_2$

$M_2:{\beta }_1,{\beta }_2$ 未指定

我们设置 ${\sigma }^2=1$， ${\beta }_1=0$。我们分别考虑了 ${\beta }_2=0,0.1,0.2,\cdots ,3$ 及样本容量 $n=15,25,50$ 的不同情形，模型数据由 $M_1$ 产生。这里模拟次数为5000次。在计算Bootstrap模型平均置信区间时，Bootstrap样本量取 $B=2000$。

图2展示了 $n=15,25,50$ 三种情形下均值的三种区间的上下错误率和覆盖率随 ${\beta }_2$ 的变化形势，图2的左图上方为三种区间的上错误率。为了表示方便，下错误率取负值后展示在图1的左图下方。表2给出了样本量 $n=50$ 时三种区间的平均区间长度。由图2和表2可以看出

1) 当样本容量较小时，比如 $n=15$，Bootstrap区间的上、下错误率和覆盖率要明显优于MATA区间和Wald区间，并且更接近名义上、下错误率2.5%和名义置信水平95%。随着样本容量的增大，三种区间的表现性能越来越接近，Bootstrap区间的上、下错误率和覆盖率总是优于MATA区间和Wald₂区间，并且更接近名义水平。

2) 当 ${\beta }_2$ 逐渐增大时，也就是两个模型 $M_1$ 和 $M_2$ 之间的差异越来越大时，三种模型平均置信区间的表现越来越好。

3) 在样本量 $n=50$ 时，Wald区间长度最短。Bootstrap区间的平均区间长度略低于MATA区间，并且随着 ${\beta }_2$ 的增大，两种区间的平均区间长度的相对差异也很小。

4. 结论

无论是Wald区间还是MATA区间，都是在参数估计或者参数估计的变换近似服从正态分布的假设下构建的，本文提出了Bootstrap模型平均置信区间，它是基于将参数估计的Bootstrap分布来近似其真实分布从而得到的一种模型平均置信区间。

当响应变量是正态数据时，我们研究了在线性正态模型下Wald区间、MATA区间和Bootstrap区间的表现性能。模拟结果表明，Bootstrap数区间在下错误率方面优势明显，一直稳定在名义下错误率附近。上错误率方面，Bootstrap区间整体上也是仅次于MATA区间，并且在平均区间长度方面，Bootstrap区间较其他三种区间是最窄的。可以看出数Bootstrap 区间缩短了区间长度，牺牲了一些覆盖率的表现性能。

当响应变量是偏态数据时，我们研究了对数正态模型下三种区间的表现性能。模拟结果表明，Bootstrap区间要优于其他两种置信区间。特别是在上错误率方面，Bootstrap区间提供了更低的更接近名义水平的上错误率。除此之外，下错误率和覆盖率也优于其他两种区间。在样本量很小时，Bootstrap区间的优势最明显。随着样本容量的增加，三种区间的表现越来越接近。在平均区间长度方面，MATA区间和Bootstrap区间相差不大。

因此，当响应变量是偏态的，甚至正态假设不明显时，都会导致MATA区间或者Wald区间的错误率和覆盖率很差。Bootstrap区间适应性比其它两个区间要强。特别是在响应变量是偏态的情况下，覆盖率是最好的。所以在不知道响应变量是否是正态的情况下，Bootstrap区间是一个不错的选择。