问

lmer(来自R包lme4)如何计算对数似然？

yellow start 发布于 2023-02-05 12:59

我正在努力理解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是固定的,我用正确的均值和最小方差计算β的估计值

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产生的.在这种情况下如何计算对数似然(对于奖励标记,为什么)？

编辑:改变符号的一致性,删除标准偏差估计的错误公式.

1 个回答

评论中的链接包含答案.下面我在这个简单的例子中给出了公式简化的内容,因为结果有些直观.

lmer适合表格的模型 $Y_ {ij} =\beta + B_i +\epsilon_ {ij}$ ,哪里 $\ epsilon_ {IJ}$ 和 $双$ 是具有差异的独立法线 $\西格马^ 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_ {\西格马^ 2}(X)= \压裂{1} {\ SQRT {2\PI \西格马^ 2}}ë^ { - \压裂{X ^ 2} {2 \西格玛^ 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$ 是来自群体观察的平均值 $一世$ .这有点直观,因为第一项捕获在每个组内传播,这应该有差异 $\西格马^ 2$ ,第二个捕获了群体之间的传播.注意 $\西格玛^ 2/n_i个+\tau蛋白^ 2$ 是方差 $\ bar y_i$ .

然而,默认情况下(REML = T)lmer最大化不是可能性而是最大化"REML标准",通过另外整合它来获得 $\公测$ 给

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

哪里 $\帽子\公测$ 如下.

最大化可能性(REML = F)

如果 $\ THETA =\TAU/\西格玛$ 是固定的,我们可以明确地找到 $\公测$ 和 $\西格玛$ 最大化可能性.他们结果是

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

注意 $\帽子\西格马^ 2$ 对于组内和组之间的变化有两个术语,和 $\帽子\公测$ 介于两者之间 $Y_ {} IJ$ 和的意思 $\ bar y_i$ 取决于的价值 $\ THETA$ .

将这些代入可能性,我们可以表示对数似然 $升$ 就......而言 $\ THETA$ 只要:

$-2L =\sum_i \日志(1 + n_i个\ THETA ^ 2)+ N(1 + \日志(2\PI \帽子\西格马^ 2))$

lmer迭代找到的值 $\ THETA$ 这最小化了这一点.在输出中, $-2L$ 和 $升$ 分别显示在"deviance"和"logLik"字段中(如果REML = F).

最大限制可能性(REML = T)

由于REML标准不依赖于 $\公测$ ,我们使用相同的估计 $\公测$ 如上.我们估计 $\西格玛$ 最大化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 \日志(1 + n_i个\ THETA ^ 2)+(N-1)(1 + \日志(2\PI \帽子\西格马^ 2))+\LOG \左(\ sum_i \压裂{ n_i个} {1个+ n_i个\ THETA ^ 2} \右)$

在lmer的输出中, $-2l_R$ 和 $L_R$ 分别显示在"REMLdev"和"logLik"字段中(如果REML = T).

2023-02-05 13:05 回答

阿穆尔湾的阿娜莎

撰写答案

今天，你开发时遇到什么问题呢？

立即提问

热门标签

lmer(来自R包lme4)如何计算对数似然？

最大化可能性(REML = F)

最大限制可能性(REML = T)