関数 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
また、シグマ ^ 2 の次の公平な推定値も計算しました。
s^2 = \sum_{i,j} alpha_i (y_{ij} - c)^2 / (1 + シータ^2 - ラムダ)
これらの推定値は、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 が生成するものではありません。この場合、対数尤度はどのように計算されますか (ボーナス マークの場合はなぜですか)。
編集:一貫性のために表記を変更し、標準偏差推定の誤った式を削除しました。