lmer(来自R包lme4)如何计算对数似然率? [英] How does lmer (from the R package lme4) compute log likelihood?
问题描述
我正在尝试了解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适合以下形式的模型:，其中和是具有方差和. 和是
其中
.
通过相对于(未观察到)给予
其中是来自和是来自组的观察值的平均值.这有点直观，因为第一项捕获在每个组中的散布，每个组应具有方差，第二个捕获组之间的传播.请注意，是.
但是，默认情况下(REML = T)，lmer不会使可能性最大化，而会使"REML准则"最大化，这是通过相对于给出
"> p
下面给出了.
最大可能性(REML = F)
如果是固定的，我们可以明确找到和可将可能性最大化.他们原来是
请注意对于组内和组间有两个变化术语，而介于和取决于的值.
将这些可能性替换为可能性，我们可以用来表达对数可能性 ://://latex.codecogs.com/gif.latex?％5ctheta"alt =" \ theta>仅:
lmer反复查找的值，以将其最小化.在输出中，和分别显示在字段"deviance"和"logLik"(如果REML = F)中.
最大化受限可能性(REML = T)
由于REML标准不取决于，因此对如上.我们估计以最大化REML标准:
受限制的对数可能性由
给出
在lmer的输出中，和分别显示在"REMLdev"和"logLik"字段中(如果REML = T).
I'm trying to understand the function lmer. I've found plenty of information about how to use the command, but not much about what it's actually doing (save for some cryptic comments here: http://www.bioconductor.org/help/course-materials/2008/PHSIntro/lme4Intro-handout-6.pdf). I'm playing with the following simple example:
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
I understand that lmer is fitting a model of the form Y_{ij} = beta + B_i + epsilon_{ij}, where epsilon_{ij} and B_i are independent normals with variances sigma^2 and tau^2 respectively. If theta = tau/sigma is fixed, I computed the estimate for beta with the correct mean and minimum variance to be
c = sum_{i,j} alpha_i y_{ij}
where
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
I also computed the following unbiased estimate for sigma^2:
s^2 = \sum_{i,j} alpha_i (y_{ij} - c)^2 / (1 + theta^2 - lambda)
These estimates seem to agree with what lmer produces. However, I can't figure out how log likelihood is defined in this context. I calculated the probability density to be
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)
But log of the above is not what lmer produces. How is log likelihood computed in this case (and for bonus marks, why)?
Edit: Changed notation for consistency, striked out incorrect formula for standard deviation estimate.
The links in the comments contained the answer. Below I've put what the formulae simplify to in this simple example, since the results are somewhat intuitive.
lmer fits a model of the form , where and are independent normals with variances and respectively. The joint probability distribution of and is therefore
where
.
The likelihood is obtained by integrating this with respect to (which isn't observed) to give
where is the number of observations from group , and is the mean of observations from group . This is somewhat intuitive since the first term captures spread within each group, which should have variance , and the second captures the spread between groups. Note that is the variance of .
However, by default (REML=T) lmer maximises not the likelihood but the "REML criterion", obtained by additionally integrating this with respect to to give
where is given below.
Maximising likelihood (REML=F)
If is fixed, we can explicitly find the and which maximise likelihood. They turn out to be
Note has two terms for variation within and between groups, and is somewhere between the mean of and the mean of depending on the value of .
Substituting these into likelihood, we can express the log likelihood in terms of only:
lmer iterates to find the value of which minimises this. In the output, and are shown in the fields "deviance" and "logLik" (if REML=F) respectively.
Maximising restricted likelihood (REML=T)
Since the REML criterion doesn't depend on , we use the same estimate for as above. We estimate to maximise the REML criterion:
The restricted log likelihood is given by
In the output of lmer, and are shown in the fields "REMLdev" and "logLik" (if REML=T) respectively.
这篇关于lmer(来自R包lme4)如何计算对数似然率?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
