0

我正在尝试为多级逻辑回归编写 stan 代码。我尝试的模型是具有两个预测变量的混合截距逻辑模型。第一级是儿童级，第二级是妈妈级。当我尝试将我编写的代码的摘要结果与由 function 生成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)

我得到了以下结果： 写模型的结果 但是，如果我使用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)

但结果不匹配，尤其是固定截距的结果。 stan_glmer() 函数生成的结果

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

4

1 回答 1

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

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

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

例如，考虑：

  • 定义：x_m = x + 1
  • 您的型号（米）：yhat_m = alpha_m + x_m1*beta_m1 + x_m2*beta_m2
  • 斯坦格尔默：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将归结为先验、分层参数的参数化、采样质量等。

另外：有多种方法可以将分类协变量编码到 model.matrix中，每种方法都用于对效果参数的特定解释。这些模型通常是等效的，这意味着可以使用上述效果的线性变换从一种参数化转换为另一种参数化。

于 2021-06-19T07:18:42.610 回答