lmer（来自R包lme4）如何计算对数似然？

Question

lmer（来自R包lme4）如何计算对数似然？

9

我正在尝试理解函数lmer。我已经找到了有关如何使用该命令的大量信息，但几乎没有关于它实际在做什么的信息（除了一些神秘的评论在这里：http://www.bioconductor.org/help/course-materials/2008/PHSIntro/lme4Intro-handout-6.pdf）。我正在尝试运用以下简单的示例：

library(data.table)
library(lme4)
options(digits=15)

n<-1000
m<-100
data<-data.table(id=sample(1:m,n,replace=T),key="id")
b<-rnorm(m)
data$y<-rand[data$id]+rnorm(n)*0.1
fitted<-lmer(b~(1|id),data=data,verbose=T)
fitted

我理解 lmer 拟合的模型形式为 Y_{ij} = beta + B_i + epsilon_{ij}，其中 epsilon_{ij} 和 B_i 是独立正态分布，方差分别为 sigma^2 和 tau^2。如果 theta = tau/sigma 是固定的，我计算出了正确均值和最小方差的 beta 估计值。

c = sum_{i,j} alpha_i y_{ij}

在哪里

alpha_i = lambda/(1 + theta^2 n_i)
lambda = 1/[\sum_i n_i/(1+theta^2 n_i)]
n_i = number of observations from group i

~~我还计算了关于sigma^2的以下无偏估计:~~

~~s^2 = \sum_{i,j} alpha_i (y_{ij} - c)^2 / (1 + theta^2 - lambda)~~

这些估计结果与lmer的产生的结果相符。但是，我无法确定在这种情况下对数似然函数如何定义。我计算出的概率密度为

pd(Y_{ij}=y_{ij}) = \prod_{i,j}[f_sigma(y_{ij}-ybar_i)]
    * prod_i[f_{sqrt(sigma^2/n_i+tau^2)}(ybar_i-beta) sigma sqrt(2 pi/n_i)]

where

ybar_i = \sum_j y_{ij}/n_i (the mean of observations in group i)
f_sigma(x) = 1/(sqrt{2 pi}sigma) exp(-x^2/(2 sigma)) (normal density with sd sigma)

但是，上述的日志并不是lmer生成的。在这种情况下，对数似然值是如何计算的（如果能解释为什么，将获得额外加分）？

编辑：为了保持一致性，更改了符号，并划掉了错误的标准偏差估计公式。

- stewbasic

3

这个软件包是开源的，所以你有查看源代码来了解它是如何计算的吗？ - Joshua Ulrich

哦，我没有意识到。我会看一下的，谢谢。 - stewbasic

你可能会发现这个问答对你有帮助：如何查看函数的源代码？。 - Joshua Ulrich

1

无论是_what_还是_why_，您都可以看一下Doug Bates的lme4草稿书...http://lme4.r-forge.r-project.org/lMMwR/lrgprt.pdf（特别是第1.4节）。不确定书中的代码是否与lme4的最新大更新相关--但这是必读的。 - Nate Pope

1

这是一个非常大、复杂的问题。Doug的书稿是一个合理的起点（但并不容易）。任何一本关于混合模型的书籍（例如Pinheiro和Bates 2000）都是一个很好的起点。 - Ben Bolker

1

谢谢提供的链接。我最终找到了Doug Bates的一篇论文（http://pages.cs.wisc.edu/~bates/reports/MixedComp.pdf），我认为它会回答我的问题。等我看完后，我会在我的简单示例中更新我的问题所涉及的内容... - stewbasic

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- stewbasic · Accepted Answer

评论中的链接包含了答案。以下是这个简单示例中公式的简化结果，因为结果有些直观。 lmer拟合一个形式为 $Y_{ij} = \beta + B_i + \epsilon_{ij}$ 的模型，其中 $\epsilon_{ij}$ 和 $B_i$ 分别是具有方差 $\sigma^2$ 和 $\tau^2$ 的独立正态分布。因此， $Y_{ij}$ 和 $B_i$ 的联合概率分布为 $\left(\prod_{i,j}f_{\sigma^2}(y_{ij}-\beta-b_i)\right)\left(\prod_i f_{\tau^2}(b_i)\right)$ 其中 $f_{\sigma^2}(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{x^2}{2\sigma^2}}$ 。

该概率是通过对 $b_i$ （未被观察到）进行积分得到的，公式如下： $\left(\prod_{i,j}f_{\sigma^2}(y_{ij}-\bar y_i)\right)\left(\prod_i f_{\sigma^2/n_i+\tau^2}(\bar y_i-\beta)\sqrt{2\pi\sigma^2/n_i}\right)$ 其中， $n_i$ 是第 $i$ 组的观察次数， $\bar y_i$ 是第 $i$ 组观察结果的平均值。这个公式有些直观，因为第一项捕捉了每组内部的差异，应该具有方差 $\sigma^2$ ，而第二项则捕捉了组之间的差异。请注意， $\sigma^2/n_i+\tau^2$ 是 $\bar y_i$ 的方差。

然而，默认情况下（REML=T），lmer最大化的不是概率，而是“REML准则”，该准则是通过对 $\beta$ 进行积分得到的。

$\left(\prod_{i,j}f_{\sigma^2}(y_{ij}-\bar y_i)\right)\left(\prod_i f_{\sigma^2/n_i+\tau^2}(\bar y_i-\hat\beta)\sqrt{2\pi\sigma^2/n_i}\right)\sqrt{\frac{2\pi\sigma^2}{\sum_i\frac{n_i}{1+n_i\theta^2}}}$

其中， $\hat\beta$ 如下所示。

最大化似然（REML=F）

如果固定 $\theta=\tau/\sigma$ ，我们可以明确地找到最大化似然的 $\beta$ 和 $\sigma$ 。它们分别为：

$\hat\beta=\frac{\sum_{i,j}y_{ij}/(1+n_i\theta^2)}{\sum_i n_i/(1+n_i\theta^2)}$

$\hat\sigma^2=\frac{1}{n}\left(\sum_{i,j}(y_{ij}-\bar y_i)^2+\sum_i\frac{n_i}{1+n_i\theta^2}(\bar y_i-\hat\beta)^2\right)$

注意： $\hat\sigma^2$ 有两个项表示组内和组间的变化， $\hat\beta$ 介于 $y_{ij}$ 和 $\bar y_i$ 的平均值之间，具体取决于 $\theta$ 的值。

将这些代入似然函数中，我们可以仅用 $\theta$ 来表达对数似然函数 $l$ ：

$-2l=\sum_i\log(1+n_i\theta^2)+n(1+\log(2\pi\hat\sigma^2))$

lmer迭代以找到最小化此值的 $\theta$ 值。在输出中，“偏差”和“对数似然”分别显示在字段“deviance”和“logLik”中（如果REML=F）。

最大化受限制的似然（REML=T）

由于REML准则不依赖于 $\beta$ ，因此我们使用与上面相同的 $\beta$ 估计值。我们估计 $\sigma$ 以最大化REML准则：

$\hat\beta=\frac{\sum_{i,j}y_{ij}/(1+n_i\theta^2)}{\sum_i n_i/(1+n_i\theta^2)}$

$\hat\sigma^2=\frac{1}{n-1}\left(\sum_{i,j}(y_{ij}-\bar y_i)^2+\sum_i\frac{n_i}{1+n_i\theta^2}(\bar y_i-\hat\beta)^2\right)$

求解方程的答案 $l_R$ 是:

$-2l_R=\sum_i\log(1+n_i\theta^2)+(n-1)(1+\log(2\pi\hat\sigma^2))+\log\left(\sum_i\frac{n_i}{1+n_i\theta^2}\right)$

在lmer的输出中，“REMLdev”和“logLik”（如果REML=T）字段分别显示了 $-2l_R$ 和 $l_R$ .