Я пытаюсь понять функцию 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 соответственно. Если тета = тау / сигма фиксировано, я вычислил оценку бета с правильным средним и минимальной дисперсией, чтобы
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 + theta ^ 2 - lambda)
Эти оценки, кажется, согласуются с тем, что дает Имер. Однако я не могу понять, как в этом контексте определяется вероятность журнала. Я рассчитал, что плотность вероятности
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 производит не лог вышеперечисленного. Как в этом случае рассчитывается логарифмическая вероятность (и почему для бонусных оценок)?
Изменить: изменено обозначение для единообразия, вычеркнута неправильная формула для оценки стандартного отклонения.
