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

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

Answer 1

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

lmer 适合 形式的模型，其中和 是具有方差的独立正态分布 和 分别。 的联合概率分布和 因此是

哪里

.

通过将其与 积分来获得可能性。 (未观察到)给予

哪里是组 的观察数，和是组观察值的平均值。这有点直观，因为第一项捕获了每个组内的分布，应该有方差 ，第二个捕获组之间的传播。请注意是 的方差.

但是，默认情况下 (REML=T) lmer 不是最大化可能性，而是最大化“REML 标准”，这是通过将其与 进行额外积分而获得的。给予

哪里下面给出。

最大化可能性 (REML=F)

如果固定后，我们可以显式找到 和 最大化可能性。事实证明他们是

注意有两个术语表示组内和组间的变化，并且 介于 的平均值之间和 的平均值取决于 的值.

将这些代入似然，我们可以表示对数似然 就而言仅:

lmer 迭代查找 的值这最大限度地减少了这种情况。在输出中，和 分别显示在“偏差”和“logLik”字段中(如果 REML=F)。

最大化限制似然 (REML=T)

因为 REML 标准不依赖于 ，我们对 使用相同的估计如上。我们估计最大化 REML 标准:

受限对数似然 由下式给出:

在 lmer 的输出中，和 分别显示在“REMLdev”和“logLik”字段中(如果 REML=T)。