关于Stan code中多级logistic模型混合截距logistic模型的问题

如何解决关于Stan code中多级logistic模型混合截距logistic模型的问题

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

library(rstanarm)
library(rstan)

data(guImmun,package = "mlmRev")
summary(guImmun)
require(dplyr)
guImmun <- guImmun %>%
  mutate(immun = ifelse(immun == "N",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,chains = 1,seed = 349)
print(M1_stanglmer,digits = 2)

但是结果不匹配，尤其是固定截取的结果。 Results generated by the stan_glmer() function

谁能帮我弄清楚我的代码有什么问题？谢谢！

解决方法

因此，我不希望您在 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

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

实际上，一旦考虑到这种转换，我希望您的拟合和 stan_glmer 模型的估计系数实际上相似。

例如，考虑：

定义：x_m = x + 1
您的模型 (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 + 两个 beta (-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有多种方式，每种方式都用于对效果参数的特定解释。这些模型通常是等效的，这意味着可以使用上述效果的线性变换从一种参数化转换为另一种参数化。