1

我使用 rethinking 包和 rstan() 复制了分层模型的结果，我只是好奇为什么 n_eff 不接近。

这是使用 rethinking 包对 2 个组 (intercept_x2) 进行随机截距的模型：

代码：

response = c(rnorm(500,0,1),rnorm(500,200,10))
predicotr1_continuous = rnorm(1000)
predictor2_categorical = factor(c(rep("A",500),rep("B",500) ))
data = data.frame(y = response, x1 = predicotr1_continuous, x2 = predictor2_categorical)
head(data)
library(rethinking)
m22 <- map2stan( 
  alist(
    y ~ dnorm( mu , sigma ) ,
    mu  <- intercept +  intercept_x2[x2] + beta*x1  ,
    intercept ~ dnorm(0,10),
    intercept_x2[x2] ~ dnorm(0, sigma_2),
    beta ~ dnorm(0,10),
    sigma ~ dnorm(0, 10),
    sigma_2 ~ dnorm(0,10)
  ) ,
  data=data , chains=1 , iter=5000 , warmup=500 )

precis(m22, depth = 2)

                  Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
intercept         9.96   9.59      -5.14      25.84  1368    1
intercept_x2[1]  -9.94   9.59     -25.55       5.43  1371    1
intercept_x2[2] 189.68   9.59     173.28     204.26  1368    1
beta              0.06   0.22      -0.27       0.42  3458    1
sigma             6.94   0.16       6.70       7.20  2927    1
sigma_2          43.16   5.01      35.33      51.19  2757    1

现在这是 rstan() 中的相同模型：

  # create a numeric vector to indicate the categorical groups
data$GROUP_ID  = match(  data$x2, levels( data$x2   )  ) 
library(rstan)
standat <- list(
  N = nrow(data),
  y = data$y,
  x1 = data$x1,
  GROUP_ID = data$GROUP_ID,
  nGROUPS = 2
)


stanmodelcode = '
data {
int<lower=1> N; 
int nGROUPS;
real y[N];                                 
real x1[N];   
int<lower=1, upper=nGROUPS> GROUP_ID[N];
}


transformed data{

}

parameters {
real intercept;  
vector[nGROUPS] intercept_x2;
real beta;                   
real<lower=0> sigma;    
real<lower=0>  sigma_2;
}

transformed parameters {  // none needed
}

model {
real mu;

// priors
intercept~ normal(0,10);
intercept_x2 ~ normal(0,sigma_2);
beta ~ normal(0,10);
sigma ~ normal(0,10);
sigma_2 ~ normal(0,10);

// likelihood

for(i in 1:N){
mu = intercept + intercept_x2[ GROUP_ID[i] ] +  beta*x1[i];
y[i] ~ normal(mu, sigma);
}

}
'

fit22 = stan(model_code=stanmodelcode, data=standat, iter=5000, warmup=500, chains = 1)
fit22


    Inference for Stan model: b212ebc67c08c77926c59693aa719288.
    1 chains, each with iter=5000; warmup=500; thin=1; 
    post-warmup draws per chain=4500, total post-warmup draws=4500.

                        mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
    intercept          10.14    0.30 9.72    -8.42     3.56    10.21    16.71    29.19  1060    1
    intercept_x2[1]   -10.12    0.30 9.73   -29.09   -16.70   -10.25    -3.50     8.36  1059    1
    intercept_x2[2]   189.50    0.30 9.72   170.40   182.98   189.42   196.09   208.05  1063    1
    beta                0.05    0.00 0.21    -0.37    -0.10     0.05     0.20     0.47  3114    1
    sigma               6.94    0.00 0.15     6.65     6.84     6.94     7.05     7.25  3432    1
    sigma_2            43.14    0.09 4.88    34.38    39.71    42.84    46.36    53.26  3248    1
    lp__            -2459.75    0.05 1.71 -2463.99 -2460.68 -2459.45 -2458.49 -2457.40  1334    1

    Samples were drawn using NUTS(diag_e) at Thu Aug 31 15:53:09 2017.
    For each parameter, n_eff is a crude measure of effective sample size,
    and Rhat is the potential scale reduction factor on split chains (at 
    convergence, Rhat=1).

我的问题：

  • 使用 rethinking() 时 n_eff 更大。存在模拟差异，但您认为这里发生了其他事情吗？
  • 除了 n_eff 不同之外，后验分布的百分位数也不同。我在想 rethinking() 和 rstan() 应该通过 5000 次迭代返回相似的结果，因为 rethinking 只是调用 rstan。两种实现之间的差异是正常的还是不同的？
  • 我创建了 data$GROUP_ID 来指示分类分组。这是将分类变量合并到 rstan() 中的层次模型中的正确方法吗？我有 2 个组，如果我有 50 个组，我使用相同的 data$GROUP_ID 向量，但这是标准方式吗？

谢谢你。

4

0 回答 0