最小范数二次无偏估计法

在统计学中，最小范数二次不偏估计法（minimum norm quadratic unbiased estimation，MINQUE）理论 由卡利安普迪·拉达克里希纳·拉奥提出。最小范数二次不偏估计法与矩估计法或最大似然估计等其他估计方法一起构成估计理论。与最佳线性不偏估计理论类似，最小范数二次不偏估计法是专门针对线性回归的模型。该方法最初是为了估计多元线性回归中的异方差误差变异数而设计的。对于混合效应模型中的变异数分量，最小范数二次不偏估计法估计量也为最大似然估计量或限制极大似然估计量提供了一种替代方案。最小范数二次不偏估计法估计量是响应变量的二次型，用于估计变异数的线性函数。

## 原理

在混合效应模型中，随机向量 $\mathbf{Y} \in \mathbb{R}^n$具有以下线性结构：

$$\mathbf{Y} = \mathbf{X}\boldsymbol\beta + \mathbf{U}_1 \boldsymbol\xi_1 + \cdots + \mathbf{U}_k \boldsymbol\xi_k$$

此处， $\mathbf{X} \in \mathbb{R}^{n\times m}$是固定效应的设计矩阵， $\boldsymbol\beta \in \mathbb{R}^m$代表未知的固定效应参数， $\mathbf{U}_i \in \mathbb{R}^{n\times c_i}$是设计矩阵第$i$个随机效应分量，而$\boldsymbol\xi_i\in\mathbb{R}^{c_i}$是第$i$ 个随机效应分量的随机向量。随机效应的均值被假设为零（$\mathbb{E}[\boldsymbol\xi_i]=\mathbf{0}$）且独立（$\mathbb{V}[\boldsymbol\xi_i]=\sigma^2_i\mathbf{I}_{c_i}$），而任意两个随机效应向量也是独立的（$\mathbb{V}[\boldsymbol\xi_i, \boldsymbol\xi_j]=\mathbf{0}\, \forall i\neq j$）。其中未知变异数$\sigma^2_1,\cdots,\sigma^2_k$代表模型的变异数分量。

这是一个涵盖常用线性回归模型的通用模型。

- 高斯-马尔可夫模型：如果我们考虑一个单分量模型，其中$\mathbf{U}_1=\mathbf{I}_n$。那么，该模型就等价于高斯-马尔可夫模型中 $\mathbf{Y}=\mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon$和$\mathbb{E}[\boldsymbol\epsilon]=\mathbf{0}$和$\mathbb{V}[\boldsymbol\epsilon]=\sigma^2_1 \mathbf{I}_n$。
- 异方差模型：每组随机变量$\mathbf{Y}$具有共同变异数的变量，并可以建模为具有适当系数的个体变异数分量： $\mathbf{U}_i$。

该模型的紧凑表示如下，其中$\mathbf{U} = \left[\begin{array}{c|c|c}\mathbf{U}_1&\cdots&\mathbf{U}_k\end{array}\right]$和$\boldsymbol\xi^\top = \left[\begin{array}{c|c|c} \boldsymbol\xi_1^\top&\cdots&\boldsymbol\xi_k^\top\end{array}\right]$：

$$\mathbf{Y}=\mathbf{X}\boldsymbol\beta+\mathbf{U}\boldsymbol\xi$$

以下模型没有对除了第一动差和第二动差之外的$\mathbf{Y}$分布做任何假设：

$$\mathbb{E}[\mathbf{Y}] = \mathbf{X}\boldsymbol\beta$$

$$\mathbb{V}[\mathbf{Y}]=\sigma^2_1\mathbf{U}_1\mathbf{U}_1^\top + \cdots + \sigma^2_k \mathbf{U}_k \mathbf{U}_k^\top \equiv \sigma^2_1\mathbf{V}_1 + \cdots + \sigma^2_k \mathbf{V}_k$$

最小范数二次不偏估计法的目标透过是使用二次型$\hat{\theta}=\mathbf{Y}^\top \mathbf{A} \mathbf{Y}$估算$\theta = \sum_{i=1}^k p_i \sigma^2_i$。最小范数二次不偏估计法估计量是通过找出一个使得估计量具有一些理想的性质的矩阵$\mathbf{A}$而得到的。下面描述矩阵$\mathbf{A}$寻找方式：

### 约束最小范数二次不偏估计法的最佳估计量性质

#### 固定效应的平移不变性

考虑对原始固定效应的转换，产生一个新的固定效应参数$\boldsymbol\gamma=\boldsymbol\beta - \boldsymbol\beta_0$，则新的等效模型如下所示：

$$\mathbf{Y} - \mathbf{X}\boldsymbol\beta_0 = \mathbf{X}\boldsymbol\gamma + \mathbf{U}\boldsymbol\xi$$

在这个等效模型下，最小范数二次不偏估计法估计量现在是$(\mathbf{Y} - \mathbf{X}\boldsymbol\beta_0)^\top \mathbf{A} (\mathbf{Y} - \mathbf{X}\boldsymbol\beta_0)$。

卡利安普迪·拉达克里希纳·拉奥认为，由于底层模型是等价的，因此该估计量应该等于$\mathbf{Y}^\top \mathbf{A} \mathbf{Y}$ 此估计量可以通过寻找符合条件的$\mathbf{A}$（能使得$\mathbf{A}\mathbf{X} = \mathbf{0}$以确保二次型的展开式中除$\mathbf{Y}^\top \mathbf{A} \mathbf{Y}$系数之外的所有系数均为零）来寻得。

#### 不偏估计

假设如上文所述限制$\mathbf{A}\mathbf{X} = \mathbf{0}$，最小范数二次不偏估计法估计量具有以下形式：

$$\begin{align} \hat{\theta} &= \mathbf{Y}^\top \mathbf{A} \mathbf{Y}\\ &= (\mathbf{X}\boldsymbol\beta + \mathbf{U}\boldsymbol\xi)^\top \mathbf{A} (\mathbf{X}\boldsymbol\beta + \mathbf{U}\boldsymbol\xi)\\ &= \boldsymbol\xi^\top\mathbf{U}^\top\mathbf{A}\mathbf{U}\boldsymbol\xi \end{align}$$

为了确保该估计量是不偏的，需要确保估计量期望值$\mathbb{E}[\hat{\theta}]$必须等于所关注的参数$\theta$，在下方公式中，由于各分量彼此独立，估计量的期望值可以分解为每个分量的期望值之和。此外，迹的循环性质被用于评估与 $\boldsymbol\xi_i$有关的期望值。

$$\begin{align} \mathbb{E}[\hat{\theta}] &= \mathbb{E}[\boldsymbol\xi^\top \mathbf{U}^\top \mathbf{A} \mathbf{U} \boldsymbol\xi]\\ &= \sum_{i=1}^k \mathbb{E}[\boldsymbol\xi_i^\top\mathbf{U}_i^\top\mathbf{A}\mathbf{U}_i\boldsymbol\xi_i]\\ &= \sum_{i=1}^k \sigma_i^2 \mathrm{Tr}[\mathbf{U}_i^\top \mathbf{A} \mathbf{U}_i] \end{align}$$

为了确保该估计量不偏，卡利安普迪·拉达克里希纳·拉奥建议设定$\sum_{i=1}^k \sigma_i^2 \mathrm{Tr}[\mathbf{U}_i^\top \mathbf{A} \mathbf{U}_i] = \sum_{i=1}^k p_i \sigma_i^2$，以藉约束$\mathbf{A}$以对所有组件使得$\mathrm{Tr}[\mathbf{U}_i^\top \mathbf{A} \mathbf{U}_i] = \mathrm{Tr}[\mathbf{A}\mathbf{V}_i] = p_i$。

#### 最小范数

拉奥认为，如果$\boldsymbol\xi$被观测到，因为$\mathbb{E}[\boldsymbol\xi_i^\top\boldsymbol\xi_i]=c_i \sigma_i^2$，一个对$\theta$的“自然”估计量将如以下公式所示。其中， $\boldsymbol\Delta$被定义为对角矩阵：

$$\frac{p_1}{c_1}\boldsymbol\xi_1^\top\boldsymbol\xi_1 + \cdots + \frac{p_k}{c_k}\boldsymbol\xi_k^\top\boldsymbol\xi_k = \boldsymbol\xi^\top\left[\mathrm{diag}\left(\frac{p_1}{c_i},\cdots,\frac{p_k}{c_k}\right)\right]\boldsymbol\xi \equiv \boldsymbol\xi^\top\boldsymbol\Delta\boldsymbol\xi$$

所提出的估计量与自然估计量之间的差是$\boldsymbol\xi^\top (\mathbf{U}^\top \mathbf{A} \mathbf{U} - \boldsymbol\Delta)\boldsymbol\xi$。此差可以通过最小化矩阵范数$\lVert \mathbf{U}^\top\mathbf{A}\mathbf{U}-\boldsymbol\Delta \rVert$被最小化。

### 运算流程

鉴于上述最佳性质推导出的约束条件和优化策略，最小范数二次不偏估计法估计量$\hat{\theta}$为了$\theta=\sum_{i=1}^k p_i\sigma_i^2$是通过选择一个在符合$\mathbf{A}\mathbf{X}=\mathbf{0}$且$\mathrm{Tr}[\mathbf{A}\mathbf{V}_i]=p_i$条件下能最小化$\lVert \mathbf{U}^\top\mathbf{A}\mathbf{U}-\boldsymbol\Delta \rVert$矩阵$\mathbf{A}$推导出来的。

## 估计量范例

### 同方差性误差的标准估计量

在高斯-马尔可夫模型中, 变异数 $\sigma^2$ 由以下公式推导：

$$s^2 = \frac{1}{n-m}(\mathbf{Y}-\mathbf{X}\hat{\boldsymbol\beta})^\top(\mathbf{Y}-\mathbf{X}\hat{\boldsymbol\beta})$$

该估计量是不偏的，并且可以证明它能最小化如下形式的欧几里得范数： $\lVert \mathbf{U}^\top\mathbf{A}\mathbf{U}-\boldsymbol\Delta \rVert$ 因此，高斯-马尔可夫模型中误差变异数的标准估计量是最小范数二次不偏估计法估计量。

### 具有共同均值和异方差性误差的随机变量

对于随机变量$Y_1,\cdots,Y_n$均值相同但变异数不同的$\sigma^2_1,\cdots,\sigma^2_n$，最小范数二次不偏估计法估计量$\sigma^2_i$是$\frac{n}{n-2}(Y_i - \overline{Y})^2 - \frac{s^2}{n - 2}$，其中$\overline{Y} = \frac{1}{n} \sum_{i=1}^n Y_i$且$s^2 = \frac{1}{n-1} \sum_{i=1}^n (Y_i - \overline{Y})^2$。

### 变异数分量估计量

卡利安普迪·拉达克里希纳·拉奥提出了一种基于最小化欧几里得范数的变异数分量模型之最小范数二次不偏估计法估计量。欧氏范数$\lVert \cdot \rVert_2$是矩阵中所有元素平方和的平方根。在评估以下范数时， $\mathbf{V}=\mathbf{V}_1+\cdots+\mathbf{V}_k = \mathbf{U} \mathbf{U}^\top$此外，利用迹的循环性质，可得： $\mathrm{Tr}[\mathbf{U}^\top\mathbf{A}\mathbf{U}\boldsymbol\Delta] = \mathrm{Tr}[\mathbf{A}\mathbf{U}\boldsymbol\Delta\mathbf{U}^\top] = \mathrm{Tr}\left[\sum_{i=1}^k \frac{p_i}{c_i} \mathbf{A}\mathbf{V}_i \right] = \mathrm{Tr}[\boldsymbol\Delta\boldsymbol\Delta]$。

$$\begin{align} \lVert \mathbf{U}^\top\mathbf{A}\mathbf{U} - \boldsymbol\Delta \rVert^2_2 &= (\mathbf{U}^\top\mathbf{A}\mathbf{U} - \boldsymbol\Delta)^\top (\mathbf{U}^\top\mathbf{A}\mathbf{U} - \boldsymbol\Delta)\\ &= \mathrm{Tr}[\mathbf{U}^\top\mathbf{A}\mathbf{U}\mathbf{U}\mathbf{A}\mathbf{U}^\top] - \mathrm{Tr}[2\mathbf{U}^\top\mathbf{A}\mathbf{U}\boldsymbol\Delta] + \mathrm{Tr}[\boldsymbol\Delta\boldsymbol\Delta]\\ &= \mathrm{Tr}[\mathbf{A}\mathbf{V}\mathbf{A}\mathbf{V}] - \mathrm{Tr}[\boldsymbol\Delta\boldsymbol\Delta] \end{align}$$

由于$\mathrm{Tr}[\boldsymbol\Delta\boldsymbol\Delta]$不取决于$\mathbf{A}$，具有欧几里得范数的最小范数二次不偏估计法矩阵可通过能最大限度地减少$\mathrm{Tr}[\mathbf{A}\mathbf{V}\mathbf{A}\mathbf{V}]$且受前面所提最小范数二次不偏估计法限制的识别矩阵$\mathbf{A}$取得。

拉奥指出满足此最佳化问题的矩阵$\mathbf{A}$是：

$\mathbf{A}_\star=\sum_{i=1}^k \lambda_i \mathbf{R}\mathbf{V}_i\mathbf{R}$,

其中$\mathbf{R} = \mathbf{V}^{-1}(\mathbf{I}-\mathbf{P})$， $\mathbf{P}=\mathbf{X}(\mathbf{X}^\top\mathbf{V}^{-1}\mathbf{X})^{-}\mathbf{X}^\top\mathbf{V}^{-1}$两者是投影到$\mathbf{X}$列空间的矩阵，而$(\cdot)^{-}$表示一矩阵的广义逆阵。

因此，最小范数二次不偏估计法估计量如下所示，其中向量 $\boldsymbol\lambda$和$\mathbf{Q}$是根据总和定义的：

$$\begin{align} \hat{\theta} &= \mathbf{Y}^\top \mathbf{A}_\star\mathbf{Y}\\ &= \sum_{i=1}^k \lambda_i \mathbf{Y}^\top\mathbf{R}\mathbf{V}_i\mathbf{R}\mathbf{Y}\\ &\equiv\sum_{i=1}^k \lambda_i Q_i\\ &\equiv \boldsymbol\lambda^\top \mathbf{Q} \end{align}$$

向量$\boldsymbol\lambda$是通过使用约束条件$\mathrm{Tr}[\mathbf{A}_\star\mathbf{V}_i]=p_i$获得的。也就是说，该向量$\boldsymbol\lambda$表示以下方程组的解。$\forall j\in\{1,\cdots,k\}$：

$$\begin{align} \mathrm{Tr}[\mathbf{A}_\star\mathbf{V}_j] &= p_j\\ \mathrm{Tr}\left[ \sum_{i=1}^k \lambda_i \mathbf{R}\mathbf{V}_i\mathbf{R}\mathbf{V}_j \right] &= p_j\\ \sum_{i=1}^k \lambda_i \mathrm{Tr}[\mathbf{R}\mathbf{V}_i\mathbf{R}\mathbf{V}_j] &= p_j \end{align}$$

这可以写成矩阵乘积的形式 $\mathbf{S}\boldsymbol\lambda=\mathbf{p}$，其中$\mathbf{p}=[p_1\,\cdots\,p_k]^\top$而$\mathbf{S}$如下所示：

$$\mathbf{S}=\begin{bmatrix} \mathrm{Tr}[\mathbf{R}\mathbf{V}_1\mathbf{R}\mathbf{V}_1] & \cdots & \mathrm{Tr}[\mathbf{R}\mathbf{V}_k\mathbf{R}\mathbf{V}_1]\\ \vdots & \ddots & \vdots\\ \mathrm{Tr}[\mathbf{R}\mathbf{V}_1\mathbf{R}\mathbf{V}_k] & \cdots & \mathrm{Tr}[\mathbf{R}\mathbf{V}_k\mathbf{R}\mathbf{V}_k] \end{bmatrix}$$

因此， $\boldsymbol\lambda=\mathbf{S}^{-}\mathbf{p}$，代表 最小范数二次不偏估计法中$\hat{\theta}=\boldsymbol\lambda^\top\mathbf{Q}=\mathbf{p}^\top(\mathbf{S}^{-})^\top\mathbf{Q}=\mathbf{p}^\top\mathbf{S}^{-}\mathbf{Q}$。此外，$\theta=\sum_{i=1}^k p_i \sigma_i^2 = \mathbf{p}^\top\boldsymbol\sigma$，其中$\boldsymbol\sigma = [\sigma^2_1\,\cdots\,\sigma^2_k]^\top$。因此，变异数分量的估计量为$\hat{\boldsymbol\sigma}=\mathbf{S}^{-}\mathbf{Q}$。

## 扩展

最小范数二次不偏估计法估计量可以在不满足不变性准则的情况下获得，在这种情况下，估计量只是不偏的，并且最小化范数。这类估计量对最小化问题有略微不同的限制。

此模型可以扩展到估计协变异数分量。在这种模型中，分量的随机效应被假设具有共同的协方差结构$\mathbb{V}[\boldsymbol\xi_i]=\boldsymbol\Sigma$。此外，一种用于混合方差和协方差分量的最小范数二次不偏估计法估计量之模型也被提出。在该模型中，在$i\in \{1,\cdots,s\}$下 $\mathbb{V}[\boldsymbol\xi_i]=\boldsymbol\Sigma$，而$i\in\{s+1,\cdots,k\}$下 $\mathbb{V}[\boldsymbol\xi_i]= \sigma_i^2\mathbf{I}_{c_i}$。