机器算法验证 - 为什么没有随机效应的 gls 模型会产生与混合效应模型类似的拟合？ - 吾爱随笔录

为什么没有随机效应的 gls 模型会产生与混合效应模型类似的拟合？

机器算法验证混合模式 lme4-nlme 广义最小二乘法

2022-04-07 21:52:39

我正在尝试回答 Pinhiero 和 Bates Mixed Effects Models in S and S-Plus的问题，解释随机效应如何无法比具有混合效应的 gls 模型带来任何好处。

该模型包含随机效应

library(nlme)
fm1Ovar.lme <- lme( follicles ~ 1 + sin(2*pi*Time) + cos(2*pi*Time), 
                         data = Ovary, 
                         random = pdDiag(~sin(2*pi*Time)),
                         corr = corARMA(p = 1, q = 1))

输出

fm1Ovar.lme

Linear mixed-effects model fit by REML
  Data: Ovary 
  Log-restricted-likelihood: -771.9471
  Fixed: follicles ~ 1 + sin(2 * pi * Time) + cos(2 * pi * Time) 
       (Intercept) sin(2 * pi * Time) cos(2 * pi * Time) 
        12.1248728         -2.9198264         -0.8487095 

Random effects:
 Formula: ~sin(2 * pi * Time) | Mare
 Structure: Diagonal
        (Intercept) sin(2 * pi * Time) Residual
StdDev:    2.614435           1.004898 3.733423

Correlation Structure: ARMA(1,1)
 Formula: ~1 | Mare 
 Parameter estimate(s):
      Phi1     Theta1 
 0.7868896 -0.2793591 
Number of Observations: 308
Number of Groups: 11

而这个模型没有

fm1Ovar.gls <- gls( follicles ~ 1 + sin(2*pi*Time) + cos(2*pi*Time), 
                    data = Ovary,
                    corr = corARMA(form = ~ 1 | Mare, p = 1, q = 1) )

输出

fm1Ovar.gls

Generalized least squares fit by REML
  Model: follicles ~ 1 + sin(2 * pi * Time) + cos(2 * pi * Time) 
  Data: Ovary 
  Log-restricted-likelihood: -773.3402

Coefficients:
       (Intercept) sin(2 * pi * Time) cos(2 * pi * Time) 
        12.0587080         -2.8832412         -0.8035591 

Correlation Structure: ARMA(1,1)
 Formula: ~1 | Mare 
 Parameter estimate(s):
      Phi1     Theta1 
 0.8908103 -0.3496050 
Degrees of freedom: 308 total; 305 residual
Residual standard error: 4.597197

比较对数似然比

anova(fm1Ovar.lme, fm1Ovar.gls)

            Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fm1Ovar.lme     1  8 1559.894 1589.657 -771.9471                        
fm1Ovar.gls     2  6 1558.680 1581.002 -773.3402 1 vs 2 2.786262  0.2483

表示具有随机效应的lme模型与没有随机效应的模型相比没有显着优势gls。所以我想我们会倾向于gls参数更少的更简约的模型。

我想知道如何解释这一点？据推测，在模型中包含随机效应没有任何好处说明了这些随机效应，但是什么？是不是母马间（即主体间）的变化非常低？

我想一个更广泛的问题是我不明白在这种情况下实际上随机效应是什么。我通过gls模型了解到，我们正在对组内观察值之间的相关性进行建模（在本例中为 Mares），类似于混合效应模型中发生的情况。那么随机效应建模为我们提供了建模相关结构没有提供的哪些额外信息呢？

1个回答

随机效应也模拟相关性。为了更正式地解释这一点，两者都适合的lme()模型gls()如下

y_{i} = X_{i} β + ε_{i}, ε_{i} \sim N (0, V_{i}),

$y_i = X_i\beta + \varepsilon_i, \quad \varepsilon_i \sim \mathcal N(0, V_i),$ 在哪里

y_{i}

$y_i$ 表示结果变量

i

$i$ -th 组（即，follicles在您的示例中，组由Mare变量指示），

X_{i}

$X_i$ 是固定效应的设计矩阵

β

$\beta$ （在你的例子中

X_{i}

$X_i$ 具有三列，截距sin(2 * pi * Time)、和cos(2 * pi * Time))。

在您的问题中，重点是 $V_i$ ，这是误差项的方差-协方差矩阵 $i$ -第组。您有以下选择：

当你使用gls()你假设 $V_i = \Sigma_i$ ，其中的结构 $\Sigma_i$ 由您在correlation参数中选择的选项确定。即，您为它指定了一个 ARMA 模型。
当您使用lme()其中仅指定random参数而不指定correlation参数时，您假设 $V_i = Z_i D Z_i^\top + \sigma^2 \texttt{I}$ ，在哪里 $Z_i$ 表示随机效应的设计矩阵， $D$ 是随机效应的方差-协方差矩阵，并且 $\sigma^2$ 是组内误差项的方差 $\texttt I$ 表示单位矩阵。在您的示例中，您指定了 $Z_i$ 有两列（intercept 和sin(2 * pi * Time)），并且 $2 \times 2$ 随机效应的协方差矩阵是对角线，协方差设置为零。
当您使用lme()其中指定random和correlation参数时，您会说 $V_i = Z_i D Z_i^\top + \Sigma_i$ ，在哪里 $\Sigma_i$ 与中的矩阵相同gls()。也就是说，您对边际协方差矩阵进行建模 $V_i$ 使用随机效应和序列相关结构。在您的示例中，您将第一个项目符号中提到的序列相关结构与第二个项目符号中提到的随机效应结构相结合。

两个模型之间的似然比检验结果表明，包含随机效应，即包含 $Z_i D Z_i^\top$ 与仅使用 ARMA 序列相关项相比，边际协方差矩阵中的项不会提高模型的拟合度 $\Sigma_i$ .

其它你可能感兴趣的问题

上一篇使用高斯过程的超参数优化下一篇在因果推理中使用种族/民族何时有效？