Указание предварительного распределения по матрице в rstan

Question

У меня проблемы с получением байесовской модели смешанных эффектов для получения стационарных и хорошо смешанных цепочек. Я создал свои собственные данные, поэтому я знаю, какие параметры должны быть получены моделью. К сожалению, потому что эффективное число параметров настолько мало, а Rhat настолько высоко, что оценки параметров - полная ерунда.

Данные составлены таким образом, что имеется 60 субъектов, разделенных на три группы (g1, g2, g3) по 20 субъектов в каждой. Каждый субъект подвергается 3 условиям (cond1, cond2, cond3). Я разработал данные таким образом, чтобы между группами не было различий, но есть различия между условиями: cond1 набирает в среднем 100 баллов, cond2 набирает 75 баллов в среднем и cond3 набирает 125 баллов.

df <- data.frame(id = factor(rep(1:60, 3)),
                 group = factor(rep(c("g1", "g2", "g3"), each = 20, length.out = 180)),
                 condition = factor(rep(c("cond1", "cond2", "cond3"), each = 60)),
                 score = c(ceiling(rnorm(60, 100, 15)), ceiling(rnorm(60, 75, 15)), ceiling(rnorm(60, 125, 15))))

Вот описание

library(dplyr)
df %>% group_by(group, condition) %>% summarise(m = mean(score), sd = sd(score))

# group condition     m    sd
# <fct> <fct>     <dbl> <dbl>
# 1 g1    cond1     108    12.4
# 2 g1    cond2      79.4  13.1
# 3 g1    cond3     128    11.5
# 4 g2    cond1     105    15.5
# 5 g2    cond2      71.6  10.6
# 6 g2    cond3     127    17.7
# 7 g3    cond1     106    13.3
# 8 g3    cond2      75.8  17.6
# 9 g3    cond3     124    14.5

Все выглядит правильно, различия между условиями хорошо сохраняются между группами.

Теперь для модели. Модель, которую я запускаю, имеет большое среднее значение, параметр для группы, параметр для условия, параметр для взаимодействия условия группы x и параметр субъекта.

Вот список данных

##### Step 1: put data into a list
mixList <- list(N = nrow(df),
                nSubj = nlevels(df$id),
                nGroup = nlevels(df$group),
                nCond = nlevels(df$condition),
                nGxC = nlevels(df$group)*nlevels(df$condition),
                sIndex = as.integer(df$id),
                gIndex = as.integer(df$group),
                cIndex = as.integer(df$condition),
                score = df$score)

Теперь для построения модели в rstan, сохраняя строку в виде файла .stan с помощью функции cat()

###### Step 2: build model
cat("
    data{
    int<lower=1> N;
    int<lower=1> nSubj;
    int<lower=1> nGroup;
    int<lower=1> nCond;
    int<lower=1,upper=nSubj> sIndex[N];
    int<lower=1,upper=nGroup> gIndex[N];
    int<lower=1,upper=nCond> cIndex[N];
    real score[N];
    }

    parameters{
    real a0;
    vector[nGroup] bGroup;
    vector[nCond] bCond;
    vector[nSubj] bSubj;
    matrix[nGroup,nCond] bGxC;
    real<lower=0> sigma_s;
    real<lower=0> sigma_g;
    real<lower=0> sigma_c;
    real<lower=0> sigma_gc;
    real<lower=0> sigma;
    }

    model{
    vector[N] mu;
    bCond ~ normal(100, sigma_c);
    bGroup ~ normal(100, sigma_g);
    bSubj ~ normal(0, sigma_s);
    sigma ~ cauchy(0,2)T[0,];
    for (i in 1:N){
    mu[i] = a0 + bGroup[gIndex[i]] + bCond[cIndex[i]] + bSubj[sIndex[i]] + bGxC[gIndex[i],cIndex[i]];
    }
    score ~ normal(mu, sigma);
    }

    ", file = "mix.stan")

Далее следует сгенерировать цепочки в rstan

##### Step 3: generate the chains
mix <- stan(file = "mix.stan",
            data = mixList,
            iter = 2e3,
            warmup = 1e3,
            cores = 1,
            chains = 1)

А вот и вывод

###### Step 4: Diagnostics
print(mix, pars = c("a0", "bGroup", "bCond", "bGxC", "sigma"), probs = c(.025,.975))

#                mean se_mean      sd      2.5%    97.5% n_eff Rhat
# a0         -1917.21  776.69 2222.64  -5305.69  1918.58     8 1.02
# bGroup[1]   2368.36 2083.48 3819.06  -2784.04  9680.78     3 1.54
# bGroup[2]   7994.87  446.06 1506.31   4511.22 10611.46    11 1.00
# bGroup[3]   7020.78 2464.68 4376.83     81.18 14699.90     3 1.91
# bCond[1]   -3887.06  906.99 1883.45  -7681.24  -247.48     4 1.60
# bCond[2]    4588.50  676.28 1941.92   -594.56  7266.09     8 1.10
# bCond[3]      73.91 1970.28 3584.74  -5386.96  5585.99     3 2.13
# bGxC[1,1]   3544.02  799.91 1819.18  -1067.27  6327.68     5 1.26
# bGxC[1,2]  -4960.08 1942.57 3137.33 -10078.84   317.07     3 2.66
# bGxC[1,3]   -396.35  418.34 1276.44  -2865.39  2543.45     9 1.42
# bGxC[2,1]  -2085.90 1231.36 2439.58  -5769.81  3689.38     4 1.46
# bGxC[2,2] -10594.89 1206.58 2560.42 -14767.50 -5074.33     5 1.02
# bGxC[2,3]  -6024.75 2417.43 4407.09 -12002.87  4651.14     3 1.71
# bGxC[3,1]  -1111.81 1273.66 2853.08  -4843.38  5572.87     5 1.48
# bGxC[3,2]  -9616.85 2314.56 4020.02 -15775.40 -4262.64     3 2.98
# bGxC[3,3]  -5054.27  828.77 2245.68  -8666.01  -321.74     7 1.00
# sigma         13.81    0.14    0.74     12.36    15.17    27 1.00

Низкое количество эффективных сэмплов и высокие Rhats говорят мне, что я здесь делаю что-то ужасно неправильное, но что?

Разве это не указывает априор на bGxC?

Как определяет один априор для матрицы?

Answer 1

Матрицы в Stan неэффективны ( см. Здесь ).Лучше использовать вектор векторов:

vector[nCond] bGxC[nGroup];

И установить приоритет:

for(i in 1:nGroup){
  bGxC[i] ~ normal(0, sigma_gc);
}

И:

for (i in 1:N){
  mu[i] = a0 + bGroup[gIndex[i]] + bCond[cIndex[i]] + bSubj[sIndex[i]] + bGxC[gIndex[i]][cIndex[i]];
}