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

Question

我试图理解函数 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)]

在哪里

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 产生的。在这种情况下如何计算对数似然（对于奖励分数，为什么）？

编辑：更改了一致性符号，删除了标准偏差估计的错误公式。

score 14 · Accepted Answer

评论中的链接包含答案。下面我将公式简化的内容放在这个简单的示例中，因为结果有点直观。

lmer 拟合形式为的模型 $Y_{ij} = \beta + B_i + \epsilon_{ij}$ ，其中 $\epsilon_{ij}$ 和 $双$ 分别是具有方差 $\sigma^2$ 和的独立法线 $\tau^2$ 。和的联合概率分布 $Y_{ij}$ 因此 $双$ 是

$\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}}$ .

可能性是通过将其与 $双$ （未观察到的）积分来获得的

$\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)$

其中 $你$ 是来自组的观察数 $一世$ ，并且 $\bar y_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}}}$