在 rstan 中指定矩阵的先验分布

Question

我在使用贝叶斯混合效应模型来生成固定且混合良好的链时遇到了问题。我已经创建了自己的数据，所以我知道模型应该检索哪些参数。不幸的是，因为参数的有效数量太少而 Rhat 太高，参数估计完全是无稽之谈。

数据设计为有 60 名受试者，分为三组（g1、g2、g3），每组 20 名受试者。每个受试者都暴露于 3 个条件（条件 1、条件 2、条件 3）。我把数据设计成组间没有差异，但是条件有差异，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 中构建模型，使用 cat() 函数将字符串保存为 .stan 文件

###### 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 中的矩阵效率低下 (see here)。最好使用向量的向量：

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]];
}

在 rstan 中指定矩阵的先验分布

Specifying prior distribution on matrix in rstan

r

stan

rstan