Скажем, у меня есть 2 факторные переменные foo и bar, которые содержат одинаковые уровни "a", "b" и "c".Есть ли способ указать в lme4 (или любом другом пакете) модель со случайными перехватами для foo и bar с корреляцией между перехватами с одинаковым уровнем?Другими словами, я думаю, что эффект "a" в foo должен коррелироваться с "a" в bar (аналогично для "b" и "c").Формально это может выглядеть примерно так:
для каждого уровня k в ["a", "b", "c"].
Вот некоторый кодкоторый оценивает sigma^2_foo и sigma^2_bar:
library(lme4)
levs <- c("a", "b", "c")
n <- 1000
df <- data.frame(y = rpois(n, 3.14),
foo = sample(levs, n, TRUE),
bar = sample(levs, n, TRUE))
mod <- glmer(y ~ (1 | foo) + (1 | bar), df, poisson)
> mod
Formula: y ~ (1 | foo) + (1 | bar)
Random effects:
Groups Name Std.Dev.
foo (Intercept) 0.009668
bar (Intercept) 0.006739
, но, конечно, пропускает корреляционный член rho.Можно ли добавить эту структуру корреляции к этой модели?
ОБНОВЛЕНИЕ
В надежде, что это полезно для людей, которые знакомы с Stan В Stan базовая реализация этой модели случайных эффектов будет выглядеть так:
data {
int<lower = 1> num_data;
int<lower = 1> num_levels;
int<lower = 0> y[num_data];
int<lower = 1, upper = num_levels> foo_ix[num_data];
int<lower = 1, upper = num_levels> bar_ix[num_data];
}
parameters {
real alpha;
vector[num_levels] alpha_foo;
vector[num_levels] alpha_bar;
real<lower = 0.0> sigma_foo;
real<lower = 0.0> sigma_bar;
real<lower = -1.0, upper = 1.0> rho;
}
transformed parameters {
matrix[2, 2] Sigma;
Sigma[1, 1] = square(sigma_foo);
Sigma[2, 1] = rho * sigma_foo * sigma_bar;
Sigma[1, 2] = rho * sigma_foo * sigma_bar;
Sigma[2, 2] = square(sigma_bar);
}
model {
for (i in 1:num_levels) {
[alpha_foo[i], alpha_bar[i]] ~ multi_normal([0.0, 0.0], Sigma);
}
y ~ poisson_log(alpha + alpha_foo[foo_ix] + alpha_bar[bar_ix]);
}