关于Stan规范中多级逻辑模型(混合截取逻辑模型)的问题

Question

我正在尝试编写多层逻辑回归的stan代码。我尝试的模型是具有两个预测器的混合截取逻辑模型。第一层是孩子层，第二层是妈妈层。当我尝试将我编写的代码的汇总结果与函数stan_glmer()生成的汇总结果进行匹配时，固定截取的结果不匹配。首先，我使用的数据如下：

library(rstanarm)
library(rstan)

data(guImmun, package = "mlmRev")
summary(guImmun)
require(dplyr)
guImmun <- guImmun %>%
  mutate(immun = ifelse(immun == "N",0,1))

其次，stan代码编写如下：

data {
    int N; // number of obs 
    int M; // number of groups 
    int K; // number of predictors
    
    int y[N]; // outcome
    row_vector[K] x[N]; // predictors
    int g[N];    // map obs to groups (kids to women)
}
parameters {
    real alpha;
    real a[M]; 
    vector[K] beta;
    real<lower=0,upper=10> sigma;  
}
model {
  alpha ~ normal(0,1);
  a ~ normal(0,sigma);
  beta ~ normal(0,1);
  for(n in 1:N) {
    y[n] ~ bernoulli(inv_logit( alpha + a[g[n]] + x[n]*beta));
  }
}

将数据拟合到模型：

guI_data <- list(g=as.integer(guImmun$mom),
                y=guImmun$immun,
                x=data.frame(guImmun$kid2p, guImmun$mom25p),
                N=nrow(guImmun),
                K=2,
                M=nlevels(guImmun$mom))
ranIntFit <- stan(file = "first_model.stan", data = guI_data,
                  iter = 500, chains = 1)
summary(ranIntFit, pars = c("alpha", "beta", "a[1]", "a[2]", "a[3]", "sigma"),
                          probs = c(0.025, 0.975),
                          digits = 2)

我得到了如下结果：results of written model但是，如果我使用stan_glmer()函数，结果将如下所示。

M1_stanglmer <- stan_glmer(immun ~ kid2p + mom25p + (1 | mom), 
                           family = binomial("logit"), 
                           data = guImmun,
                           iter = 500,
                           chains = 1,
                           seed = 349)
print(M1_stanglmer, digits = 2)

但结果并不匹配，特别是固定截距的结果。Results generated by the stan_glmer() function

有人能帮我找出我的代码出了什么问题吗？谢谢!

Answer 1

因此，我不期望Stan中的模型和stan_glmer中实现的版本之间存在精确的等价性，但对于样本良好的模型，合理地期望估计值相似。

然而，在您的情况下，还有另一个问题会影响您的估计：

您在guI_Data$x对象中使用的协变量具有{1,2}中的值，其中典型的实现将使用{0,1}中的值来表示二进制协变量。这就是在stan_glmer中所做的。

如果您使用glimpse检查数据结构，则此编码非常明显：

> library(tidyverse)
> glimpse(guI_data)
List of 6
 $ g: int [1:2159] 1 2 3 4 5 5 6 7 7 8 ...
 $ y: num [1:2159] 1 0 0 0 0 1 1 1 1 1 ...
 $ x:'data.frame':  2159 obs. of  2 variables:
  ..$ guImmun.kid2p : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 1 2 2 ...
  ..$ guImmun.mom25p: Factor w/ 2 levels "N","Y": 1 1 1 1 2 1 1 2 2 2 ...
 $ N: int 2159
 $ K: num 2
 $ M: int 1595

这对intercept参数的影响最大，因为当所有协变量都为0时，intercept表示预期的线性预测器。当协变量被转换或添加时，该值通常会发生变化。

实际上，如果考虑到这种变换，我希望从拟合和stan_glmer模型得到的估计系数实际上是相似的。

例如，考虑以下内容：

模型定义：x_m = x + 1

• Your (m)：yhat_m = alpha_m + x_m1*beta_m1 + x_m2*beta_m2
• Stan_glmer：yhat = alpha + x_1*beta_1 + x_2*beta_2

并替换为：

yhat_m = alpha_m + (x_1 + 1)*beta_m1 + (x_2 + 1)*beta_m1

yhat_m = alpha_m + x_1*beta_m1 + beta_m1 + x_2*beta_m2 + beta_m2

yhat_m = alpha_m + beta_m1 + beta_m2 + x_1*beta_m1 + x_2*beta_m2

如果我们假设yhat_m ~= yhat、beta_m1 ~= beta_1和beta_m2 ~= beta_2...然后

alpha = alpha_m + beta_m1 + beta_m2

因此，我预计stan_glmer alpha (-1.7)将接近手工编码的Stan alpha +两个betas (-3.2 + 1.7 - 0.1)。

它确实是(-1.6)。

如果您进一步更新Stan数据以将这些协变量缩放为{0,1}而不是{1,2}：

guI_data2 <- list(g=as.integer(guImmun$mom),
                y=guImmun$immun,
                x=data.frame(guImmun$kid2p == "Y", guImmun$mom25p == "Y"),
                N=nrow(guImmun),
                K=2,
                M=nlevels(guImmun$mom))
ranIntFit2 <- stan(file = "first_model.stan", data = guI_data2,
                  iter = 500, chains = 1)

然后看一下输出：

> summary(ranIntFit2, pars = c('alpha', 'beta'))
$summary
              mean     se_mean        sd       2.5%        25%        50%           75%      97.5%     n_eff      Rhat
alpha   -1.5110714 0.022982199 0.1903571 -1.8974997 -1.6318370 -1.5038593 -1.3861628334 -1.1729671  68.60488 1.0505237
beta[1]  1.5224756 0.025017739 0.1737332  1.2260666  1.4058789  1.5118314  1.6492158203  1.8673450  48.22471 1.0592955
beta[2] -0.1206084 0.009410305 0.1640406 -0.4267987 -0.2368855 -0.1267984 -0.0003187197  0.1894375 303.87510 0.9964177

你可以自己确认一下，你在正确的范围内。

在此之后，您的模型和stan_glmer之间的差异将归结为先验、分层参数的参数化、采样质量等。

旁白：categorical covariates can be coded into a model.matrix有多种方式，每种方式都针对特定的效果参数进行解释。这些模型通常是等效的，这意味着可以使用如上所述的效果线性变换将一种参数化转换为另一种参数化。