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

stewbasic 发表于 Dev

炖的

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

编辑：更改了表示法的一致性，删除了用于标准偏差估计的错误公式。

炖的

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

lmer拟合形式的模型 $Y_ {ij} = \ beta + B_i + \ epsilon_ {ij}$ ，其中 $\ epsilon_ {ij}$ 和 $双$ 是分别具有方差 $\ sigma ^ 2$ 和的独立法线 $\你的^ 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}（\ y_i- \ beta）\ sqrt {2 \ pi \ sigma ^ 2 / n_i} \右）$

其中， $你$ 是来自组的观察次数 $一世$ ， $\栏y_i$ 是来自组的观察平均值 $一世$ 。这有点直观，因为第一个术语捕获了每个组内的差异 $\ sigma ^ 2$ ，第二个术语捕获了组之间的差异。请注意， $\ sigma ^ 2 / n_i + \ tau ^ 2$ 是的方差 $\栏y_i$ 。

但是，默认情况下（REML = T），lmer不会使可能性最大化，而是使“ REML准则”最大化，而该“ 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}（\ 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}}}$

在 $\帽子\ 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）$

注意 $\帽子\西格玛^ 2$ 在组内和组之间有两个变化项，并且 $\帽子\ Beta$ 在的平均值 $y_ {ij}$ 和的平均值之间， $\栏y_i$ 取决于的值 $\ theta$ 。

将这些替换为可能性，我们可以仅用以下 $升$ 形式表达对数可能性 $\ theta$ ：

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

lmer反复查找以将其值 $\ theta$ 最小化。在输出中， $-2升$ 并且 $升$ 分别在字段“ 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} \右）$

11聚物中的输出， $-2l_R$ 并且 $l_R$ 被示出在字段“REMLdev”和（REML = T IF）分别为“logLik”。

本文收集自互联网，转载请注明来源。

如有侵权，请联系[email protected] 删除。

编辑于2020-11-29

我来说两句

0条评论

登录后参与评论

上一篇：在Swift中从数组中获取随机元素

来自分类Dev

Related 相关文章

文章