【基础还得练】H-似然

H似然（H-likelihood）是统计学中用于处理具有隐变量模型的似然函数的一种扩展方法。在统计学中，似然函数是一种衡量在给定参数的情况下，观察到的数据出现的概率。当模型中包含不可观测的变量，即隐变量时，标准的最大似然估计方法可能不再适用，因为隐变量的存在使得似然函数变得复杂且难以处理。

H似然是由日本统计学家Hirotugu Akaike提出的，旨在解决含有隐变量的情况下的参数估计问题。H似然方法的核心思想是将隐变量视为缺失数据，并利用期望最大化（EM）算法或者其它相关算法来估计模型参数。
具体来说，H似然的主要特点和步骤如下：

1. 隐变量处理：将模型中的隐变量视为缺失数据，不需要明确地对其进行推断。
2. 积分消除：通过对隐变量进行积分，消除似然函数中的隐变量，从而得到一个关于模型参数的函数。
3. 似然函数的近似：在积分难以计算的情况下，H似然方法通过一些近似手段来简化计算，比如使用牛顿-拉夫森方法或者泰勒展开。
4. 参数估计：通过最大化H似然函数来估计模型参数。
   H似然方法在时间序列分析、结构方程模型以及其它统计模型中有着广泛的应用，特别是在处理具有隐变量的问题时。通过使用H似然，研究者可以更加灵活地处理复杂的统计模型，并得到可靠的参数估计。

---

H似然方法的核心在于如何处理包含隐变量的似然函数。以下是一些基本的公式和概念，用于解释H似然方法：

基本似然函数

首先，我们考虑一个包含观测变量 $Y$ 和隐变量 $Z$ 的模型。其联合似然函数可以表示为：
$$L(\theta ;Y,Z)=P(Y,Z;\theta )$$
其中，$\theta$ 是模型参数。

完全数据的对数似然

如果我们能够观测到 $Z$，那么完全数据的对数似然函数为：
$$\log L(\theta ;Y,Z)=\log P(Y,Z;\theta )$$

H似然函数

然而，由于 $Z$ 是隐变量，我们不能直接观测到它。H似然函数是通过考虑 $Z$ 的期望值来近似对数似然函数。H似然函数定义为：
$$H(\theta ;Y)=\int \log P(Y,z;\theta )P(z|Y;\theta )dz$$
这里，$P(z|Y;\theta )$ 是给定观测数据 $Y$ 下隐变量 $Z$ 的后验概率密度。

EM算法中的应用

在期望最大化（EM）算法中，H似然函数可以用来迭代地估计参数 $\theta$。EM算法的步骤如下：

1. E步（期望步）：计算隐变量 $Z$ 的期望值，即后验概率 $P(z|Y;{\theta }^{(t)})$。
2. M步（最大化步）：最大化关于 $Z$ 的期望对数似然，即最大化 $Q(\theta ,{\theta }^{(t)})$，其中：
   $$Q(\theta ,{\theta }^{(t)})=\int \log P(Y,z;\theta )P(z|Y;{\theta }^{(t)})dz$$
   在H似然方法中，$Q$ 函数实际上就是H似然函数的一个表达式。

H似然的近似

在实际操作中，H似然函数可能很难直接计算，因此需要近似。例如，可以使用泰勒展开来近似对数似然函数：
$$\log P(Y,z;\theta )\approx \log P(Y,z;{\theta }_0)+\frac{\partial \log P(Y,z;\theta )}{\partial \theta }{|}_{\theta ={\theta }_0}(\theta -{\theta }_0)+\text{higher order terms}$$
这里，${\theta }_0$ 是参数的一个初始估计。

参数估计

最终，我们通过最大化H似然函数来估计参数 $\theta$：
$$\hat{\theta }=\arg \underset{\theta }{\max }H(\theta ;Y)$$
通过迭代地进行E步和M步，直到参数估计收敛，我们可以得到模型参数的估计值。
需要注意的是，H似然方法的具体公式和推导过程会根据不同的统计模型而有所不同。上述公式提供了一个通用的框架，但实际应用时需要根据具体模型进行相应的调整。

---

让我们考虑一个简单的线性混合效应模型（也称为随机效应模型）作为H似然的一个具体案例。在这个模型中，我们假设观测数据是由固定效应和随机效应共同作用的结果。

模型设定

假设我们有以下线性混合效应模型：
$$Y_i=\beta X_i+b_i+{ϵ}_i$$
其中：

• $Y_i$ 是第 $i$ 个观测值。
• $X_i$ 是对应的固定效应设计矩阵。
• $\beta$ 是固定效应参数。
• $b_i$ 是第 $i$ 个观测的随机效应，假设 $b_i\sim N(0,{\sigma }_b^2)$。
• ${ϵ}_i$ 是观测误差，假设 ${ϵ}_i\sim N(0,{\sigma }^2)$。

似然函数

对于这个模型，完整的似然函数（包含随机效应 $b_i$）可以写为：
$$L(\beta ,{\sigma }^2,{\sigma }_b^2;Y,b)=\prod _{i=1}^N\left[\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(Y_i-\beta X_i-b_i{)}^2}{2{\sigma }^2}\right)\right]\times \prod _{i=1}^N\left[\frac{1}{\sqrt{2\pi {\sigma }_b^2}}\exp \left(-\frac{b_i^2}{2{\sigma }_b^2}\right)\right]$$

H似然函数

由于随机效应 $b_i$ 是未观测的，我们需要消除它们来得到H似然函数。H似然函数是通过对随机效应 $b_i$ 进行积分得到的：
$$H(\beta ,{\sigma }^2,{\sigma }_b^2;Y)=\int L(\beta ,{\sigma }^2,{\sigma }_b^2;Y,b)db$$
具体来说，对于每个 $i$，我们需要对 $b_i$ 进行积分：
$$H_i(\beta ,{\sigma }^2,{\sigma }_b^2;Y_i)=\int \left[\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(Y_i-\beta X_i-b_i{)}^2}{2{\sigma }^2}\right)\right]\left[\frac{1}{\sqrt{2\pi {\sigma }_b^2}}\exp \left(-\frac{b_i^2}{2{\sigma }_b^2}\right)\right]d{b}_i$$

积分的计算

通过对 $b_i$ 进行积分，我们可以得到 $H_i$ 的表达式。这个积分实际上是两个高斯分布的乘积的积分，可以通过合并两个高斯分布来简化：
$$H_i(\beta ,{\sigma }^2,{\sigma }_b^2;Y_i)=\frac{1}{\sqrt{2\pi ({\sigma }^2+{\sigma }_b^2)}}\exp \left(-\frac{(Y_i-\beta X_i{)}^2}{2({\sigma }^2+{\sigma }_b^2)}\right)$$

最终的H似然函数

将所有观测值的 $H_i$ 相乘，我们得到最终的H似然函数：
$$H(\beta ,{\sigma }^2,{\sigma }_b^2;Y)=\prod _{i=1}^N\frac{1}{\sqrt{2\pi ({\sigma }^2+{\sigma }_b^2)}}\exp \left(-\frac{(Y_i-\beta X_i{)}^2}{2({\sigma }^2+{\sigma }_b^2)}\right)$$

参数估计

最后，我们通过最大化H似然函数来估计参数 $\beta$，${\sigma }^2$，和 ${\sigma }_b^2$：
$$\hat{\beta },{\hat{\sigma }}^2,{\hat{\sigma }}_b^2=\arg \underset{\beta ,{\sigma }^2,{\sigma }_b^2}{\max }H(\beta ,{\sigma }^2,{\sigma }_b^2;Y)$$
这通常通过数值优化方法来实现。
这个案例展示了如何将H似然方法应用于一个具体的统计模型，即线性混合效应模型，以估计模型参数。

---

如果在自变量 $X$ 上直接添加随机效应，这意味着 $X$ 本身是一个随机变量，其包含固定部分和随机部分。这种情况在统计建模中是可行的，通常被称为随机系数模型。以下是在自变量 $X$ 上添加随机效应的情况以及求解思路：

模型设定

假设我们的模型是：
$$Y_i={\beta }_0+{\beta }_1(X_{1i}+b_{1i})+{ϵ}_i$$
这里，$X_{1i}$ 是第 $i$ 个观测的自变量，$b_{1i}$ 是对应的随机效应，假设 $b_{1i}\sim N(0,{\sigma }_{b1}^2)$，而 ${ϵ}_i$ 是误差项，假设 ${ϵ}_i\sim N(0,{\sigma }^2)$。

求解思路

1. 模型表述：首先，我们需要明确模型的形式，包括固定效应和随机效应。
2. 条件分布：给定 $b_{1i}$，$Y_i$ 的条件分布是：
   $$Y_i|b_{1i}\sim N({\beta }_0+{\beta }_1(X_{1i}+b_{1i}),{\sigma }^2)$$
3. 边际分布：我们需要从条件分布中消除随机效应 $b_{1i}$ 来得到 $Y_i$ 的边际分布。这通常涉及到积分操作：
   $$p(Y_i|{\beta }_0,{\beta }_1,{\sigma }^2,{\sigma }_{b1}^2,X_{1i})=\int p(Y_i|b_{1i},{\beta }_0,{\beta }_1,{\sigma }^2,X_{1i})p(b_{1i}|{\sigma }_{b1}^2)d{b}_{1i}$$
4. H似然函数：H似然函数可以通过对 $b_{1i}$ 的条件分布取期望来构建：
   $$H({\beta }_0,{\beta }_1,{\sigma }^2,{\sigma }_{b1}^2;Y,X)=\sum _i\log \left(\int p(Y_i|b_{1i},{\beta }_0,{\beta }_1,{\sigma }^2,X_{1i})p(b_{1i}|{\sigma }_{b1}^2)d{b}_{1i}\right)$$
5. 参数估计：通过最大化H似然函数来估计固定效应参数 ${\beta }_0$ 和 ${\beta }_1$，以及随机效应的方差参数 ${\sigma }_{b1}^2$ 和误差项的方差 ${\sigma }^2$。这通常通过数值优化方法实现，例如使用迭代重加权最小二乘（IRWLS）或者梯度上升法。
6. 迭代过程：在实际操作中，可能需要使用EM算法（期望最大化算法）来迭代求解。E步涉及计算随机效应的期望，M步则是最大化H似然函数。
   在具体实现时，由于涉及到复杂的积分运算，可能需要使用数值积分方法或者近似方法来简化计算。这种方法在处理具有随机系数的统计模型时非常有用，尤其是在生物统计和心理学等领域，其中个体差异可能导致自变量的效应具有随机性。

---

下面我将提供一个简化的例子，说明如何在混合效应模型中求解H似然。

模型设定

考虑以下混合效应模型：
$$Y_{ij}={\beta }_0+{\beta }_1(X_{ij}+b_j)+{ϵ}_{ij}$$
这里：

• $Y_{ij}$ 是第 $j$ 所学校中第 $i$ 个学生的成绩。
• $X_{ij}$ 是第 $j$ 所学校中第 $i$ 个学生的家庭背景指标。
• ${\beta }_0$ 是全局截距项。
• ${\beta }_1$ 是固定效应系数，表示家庭背景对成绩的平均影响。
• $b_j$ 是第 $j$ 所学校的随机效应，它直接作用于家庭背景指标 $X_{ij}$。
• ${ϵ}_{ij}$ 是第 $j$ 所学校中第 $i$ 个学生的随机误差项。

H似然的求解过程

1. 写出完整的似然函数：
   由于 $b_j$ 是随机效应，我们需要对 $b_j$ 进行积分以得到边际似然函数（H似然）。
   $$L(\theta ;Y,X)=\int \prod _{j=1}^J\prod _{i=1}^{n_j}\left[\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{1}{2{\sigma }^2}(Y_{ij}-{\beta }_0-{\beta }_1(X_{ij}+b_j){)}^2\right)\right]\varphi (b_j;0,{\sigma }_b^2)d{b}_j$$
   其中 $\varphi (b_j;0,{\sigma }_b^2)$ 是 $b_j$ 的正态密度函数，${\sigma }^2$ 是误差项的方差，${\sigma }_b^2$ 是随机效应的方差。
2. 最大化H似然：
   由于直接最大化上述积分似然函数非常困难，通常使用数值方法，如限制最大似然（REML）或期望最大化（EM）算法。
3. 使用EM算法：
   • E步：计算随机效应 $b_j$ 的期望值和条件方差。
   • M步：最大化关于固定效应参数和随机效应方差的期望似然函数。
4. 迭代：
   重复E步和M步直到参数估计收敛。

结果

假设我们已经通过EM算法或其他数值优化方法得到了参数的估计值。以下是一些可能的估计结果：

• ${\hat{\beta }}_0$：全局截距项的估计值。
• ${\hat{\beta }}_1$：固定效应系数的估计值，表示家庭背景对成绩的平均影响。
• ${\hat{\sigma }}^2$：误差项方差的估计值。
• ${\hat{\sigma }}_b^2$：随机效应方差的估计值。
  这些结果可以用来解释模型，例如，${\hat{\beta }}_1$ 告诉我们家庭背景对成绩的平均影响，而 ${\hat{\sigma }}_b^2$ 则告诉我们学校之间对家庭背景影响的差异程度。
  实际操作中，这些计算通常使用专门的统计软件来完成，如R的lme4包或nlme包，它们内部实现了EM算法和其他数值优化技术来求解混合效应模型。

---

H似然（H-likelihood）是一种结合固定效应和随机效应的完整模型似然方法，广泛用于广义线性混合模型（GLMMs）等随机效应模型的参数估计。以下是一个H似然案例及其求解过程。

案例背景

我们研究一个作物实验，其中作物产量受到肥料种类和随机的地块效应的影响。数据如下：

• 固定效应：肥料种类（A和B，因子变量）
• 随机效应：地块效应（每个地块的作物生长条件不同）
• 响应变量：作物产量（连续变量）

假设模型为：
$$y_{ij}={\beta }_0+{\beta }_1{x}_{ij}+u_i+{ϵ}_{ij},u_i\sim N(0,{\sigma }_u^2),{ϵ}_{ij}\sim N(0,{\sigma }_{ϵ}^2),$$
其中：

• $y_{ij}$ 是第 $i$ 个地块在第 $j$ 次试验中的产量；
• $x_{ij}$ 是肥料种类的指标变量（0表示A，1表示B）；
• $u_i$ 是地块的随机效应；
• ${\beta }_0,{\beta }_1$ 是固定效应的回归系数；
• ${\sigma }_u^2,{\sigma }_{ϵ}^2$ 是随机效应和误差的方差。

求解过程

1. H似然的定义

H似然由两个部分组成：

• 固定效应的似然（基于数据的观测值）：
  $$L_f(\beta ;y,u)=\prod _{i,j}\varphi \left(y_{ij}\mid {\beta }_0+{\beta }_1{x}_{ij}+u_i,{\sigma }_{ϵ}^2\right),$$
  其中 $\varphi (\cdot )$ 是正态分布的概率密度函数。

• 随机效应的似然（基于随机效应的分布）：
  $$L_r(u;{\sigma }_u^2)=\prod _i\varphi \left(u_i\mid 0,{\sigma }_u^2\right).$$

H似然为两者的乘积：
$$H(\beta ,u;y,{\sigma }_u^2,{\sigma }_{ϵ}^2)=L_f(\beta ;y,u)\cdot L_r(u;{\sigma }_u^2).$$

2. 参数估计的目标

通过最大化H似然，估计${\beta }_0,{\beta }_1,u_i,{\sigma }_u^2,{\sigma }_{ϵ}^2$。

3. 计算步骤

1. 固定效应参数和随机效应的初步估计：
   给定初始值 ${\sigma }_u^2,{\sigma }_{ϵ}^2$，使用联合最大化的方法，估计 ${\beta }_0,{\beta }_1,u_i$：
   $${\hat{u}}_i=\arg \underset{u}{\max }H(\beta ,u;y,{\sigma }_u^2,{\sigma }_{ϵ}^2).$$

2. 方差分量的估计：
   方差分量 ${\sigma }_u^2,{\sigma }_{ϵ}^2$ 可以通过期望最大化（EM）算法或REML（限制最大似然）方法估计。

4. 数值示例

假设：

• 数据：$y=[10,12,15,13,14]$，肥料种类 $x=[0,1,0,1,0]$，地块分组为 $[1,1,2,2,3]$；
• 初始值：${\beta }_0=10,{\beta }_1=2,{\sigma }_u^2=1,{\sigma }_{ϵ}^2=1$。

通过迭代优化，得到：

• 固定效应估计：${\hat{\beta }}_0=11.2,{\hat{\beta }}_1=1.8$；
• 随机效应估计：${\hat{u}}_1=0.5,{\hat{u}}_2=-0.3,{\hat{u}}_3=0.2$；
• 方差分量估计：${\hat{\sigma }}_u^2=0.8,{\hat{\sigma }}_{ϵ}^2=1.2$。

通过这种方法，H似然允许我们同时估计固定效应和随机效应的贡献，并通过联合优化提高估计的效率和精确度。

---

是的，模型 $$Y_{ij}={\beta }_0+{\beta }_1(X_{ij}+b_j)+{ϵ}_{ij}$$ 仍然可以使用 H似然 方法进行求解，因为它包含了固定效应（${\beta }_0,{\beta }_1$）和随机效应（$b_j$）。下面是针对该模型的具体求解过程。

模型说明

• 响应变量：$Y_{ij}$，第 $j$ 个随机组内第 $i$ 个观测值；
• 固定效应：
  • ${\beta }_0$：截距；
  • ${\beta }_1$：固定效应斜率；
• 随机效应：
  • $b_j$：第 $j$ 个随机组的效应，假设 $b_j\sim N(0,{\sigma }_b^2)$；
• 误差项：
  • ${ϵ}_{ij}\sim N(0,{\sigma }_{ϵ}^2)$，为独立同分布的随机误差。

将 $X_{ij}+b_j$ 看作 ${\beta }_1$ 的“混合效应因子”，模型属于 线性混合效应模型 的一种。

H似然的分解

H似然的核心是结合固定效应和随机效应的联合似然。对于该模型，H似然由以下部分构成：

1. 固定效应的似然

$$L_f(\beta ;Y,b)=\prod _j\prod _i\varphi \left(Y_{ij}\mid {\beta }_0+{\beta }_1(X_{ij}+b_j),{\sigma }_{ϵ}^2\right),$$
其中，$\varphi (\cdot )$ 是正态分布的概率密度函数。

2. 随机效应的似然

$$L_r(b;{\sigma }_b^2)=\prod _j\varphi \left(b_j\mid 0,{\sigma }_b^2\right).$$

3. H似然公式

两者结合后，H似然为：
$$H(\beta ,b;Y,{\sigma }_b^2,{\sigma }_{ϵ}^2)=L_f(\beta ;Y,b)\cdot L_r(b;{\sigma }_b^2).$$

求解过程

我们希望通过最大化H似然来估计参数 ${\beta }_0,{\beta }_1,b_j,{\sigma }_b^2,{\sigma }_{ϵ}^2$。常用的优化方法如下：

1. 步骤一：初始化参数

给定初始值：

• 固定效应：${\beta }_0,{\beta }_1$；
• 随机效应：$b_j=0$；
• 方差分量：${\sigma }_b^2,{\sigma }_{ϵ}^2$。

2. 步骤二：条件优化

采用迭代方法优化固定效应和随机效应：

1. 固定效应：在给定随机效应 $b_j$ 的条件下，优化固定效应 ${\beta }_0,{\beta }_1$：
   $$({\hat{\beta }}_0,{\hat{\beta }}_1)=\arg \underset{\beta }{\max }H(\beta ,b;Y,{\sigma }_b^2,{\sigma }_{ϵ}^2).$$

2. 随机效应：在给定固定效应的条件下，优化随机效应 $b_j$：
   $${\hat{b}}_j=\arg \underset{b_j}{\max }H(\beta ,b;Y,{\sigma }_b^2,{\sigma }_{ϵ}^2).$$

3. 方差分量：更新 ${\sigma }_b^2,{\sigma }_{ϵ}^2$，可以通过期望最大化（EM）或限制最大似然（REML）方法估计。

3. 步骤三：迭代收敛

重复步骤二，直到参数估计收敛。

示例

假设数据如下：

• $Y_{ij}=[2.1,3.5,2.8,3.9,4.0]$；
• $X_{ij}=[0.5,1.0,0.7,1.2,1.0]$；
• $j=[1,1,2,2,3]$（分组标识）；
• 初始值：${\beta }_0=2.0,{\beta }_1=1.0,{\sigma }_b^2=0.5,{\sigma }_{ϵ}^2=0.2$。

通过迭代优化，估计结果可能为：

• 固定效应：${\hat{\beta }}_0=2.1,{\hat{\beta }}_1=0.9$；
• 随机效应：${\hat{b}}_1=0.2,{\hat{b}}_2=-0.1,{\hat{b}}_3=0.05$；
• 方差分量：${\hat{\sigma }}_b^2=0.4,{\hat{\sigma }}_{ϵ}^2=0.25$。

总结

这种模型的H似然方法通过联合优化固定效应和随机效应，提供了一种自然的估计方法。与经典的REML方法相比，H似然的最大优点是直接将随机效应视为参数进行估计，适用于复杂随机效应结构的模型。