6

我想在 Stan 中运行稳健的逻辑回归（robit）。该模型在 Gelman & Hill 的“使用回归和多级方法的数据分析”（2006 年，第 124 页）中提出，但我不确定如何实现它。我检查了Stan 的 Github 存储库和参考手册，但不幸的是我仍然感到困惑。这是我用来模拟常规逻辑回归的一些代码。我应该添加什么，以便错误跟随，例如，在具有 7 个自由度的分布时？无论如何，如果我运行多级分析，它会是相同的程序吗？

library(rstan)

set.seed(1)
x1 <- rnorm(100)  
x2 <- rnorm(100)
z <- 1 + 2*x1 + 3*x2      
pr <- 1/(1+exp(-z))       
y <- rbinom(100,1,pr)  

df <- list(N=100, y=y,x1=x1,x2=x2)

# Stan code
model1 <- '
data {                          
  int<lower=0> N;          
  int<lower=0,upper=1> y[N];  
  vector[N] x1;         
  vector[N] x2;
}
parameters {
  real beta_0;     
  real beta_1;        
  real beta_2; 
}
model {
  y ~ bernoulli_logit(beta_0 + beta_1 * x1 + beta_2 * x2);
}
'
# Run the model
fit <- stan(model_code = model1, data = df, iter = 1000, chains = 4)
print(fit)

谢谢！

4

4 回答 4

5

我一定遗漏了一些东西，但我无法适应 danilofreire 从 Luc 发布的解决方案。所以我只是从 JAGS 翻译了一个模型。

我认为这是正确的，尽管看起来与 Luc 的解决方案有点不同。

library(rstan)

N <- 100
x1 <- rnorm(N)
x2 <- rnorm(N)
beta0 <- 1
beta1 <- 2
beta2 <- 3

eta <- beta0 + beta1*x1 + beta2*x2                         # linear predictor
p <- 1/(1 + exp(-eta))                                     # inv-logit
y <- rbinom(N, 1, p)                   

dlist <- list(y = y, x1 = x1, x2 = x2, N = N, nu = 3)      # adjust nu as desired df

mod_string <- "
  data{
    int<lower=0> N;
    vector[N] x1;
    vector[N] x2;
    int<lower=0, upper=1> y[N];
    real nu;
  }
  parameters{
    real beta0;
    real beta1;
    real beta2;
  }
  model{
    vector[N] pi;

    for(i in 1:N){
      pi[i] <- student_t_cdf(beta0 + beta1*x1[i] + beta2*x2[i], nu, 0, 1);
      y[i] ~ bernoulli(pi[i]);
    }
  }
"
fit1 <- stan(model_code = mod_string, data = dlist, chains = 3, iter = 1000)
print(fit1)
于 2014-10-25T17:20:01.067 回答
4

Luc Coffeng 在Stan 邮件列表上给我发了这个答案，我想我应该把它加在这里。他说：

“采用 GLM 作为您的 robit 回归的基础：只需将标准误差项替换为e ~ student_t(7, 0, sigma_e), wheresigma_e ~ cauchy(0, 2)或您认为可以的任何规模（我不会超过 5，因为 (-5,5) 的逆 logit 涵盖了大多数[0,1] 区间的值。除了 t-error 的标度外，您还可以指定 t-error 的 df 作为参数。请参阅下面的建议代码。

但是，我确实希望您的数据包含比您提供的玩具示例更多的信息，即每个人的多个观察结果（如下所示）。每个人/单位只有一次观察，该模型实际上是不可能识别的。”

然后他提供了以下示例：

library(rstan)

set.seed(1)
x1 <- rnorm(100)  
x2 <- rnorm(100)
z <- 1 + 2*x1 + 3*x2 + 0.1 * rt(100, 7)
pr <- 1/(1+exp(-z))       
y <- rbinom(100,10,pr)  

df <- list(N=100, y=y, x1=x1, x2=x2, nu = 7)

# Stan code
model1 <- '
data {                          
   int<lower=0> N;          
   int<lower=0,upper=10> y[N];  
   vector[N] x1;         
   vector[N] x2;
   real nu;
}
parameters {
   real beta_0;     
   real beta_1;        
   real beta_2; 
   real<lower=0> sigma_e;
   vector[N] e;
}
model {
   e ~ student_t(nu, 0, sigma_e);
   sigma_e ~ cauchy(0, 1);
   y ~ binomial_logit(10, beta_0 + beta_1 * x1 + beta_2 * x2 + e);
}
'
# Run the model
fit <- stan(model_code = model1, data = df, iter = 4000, chains = 2)
print(fit)

Bob Carpenter 还简要评论了这个问题：

“[...] 是的，你可以在分层设置中做同样的事情，但你必须小心，因为当你接近常态时，自由度建模可能会很棘手，因为规模会跑到无穷大。”

在回答 Bernd 的问题时，Luc 澄清了他y ~ bernoulli_logit(10...在模型代码中写的原因：

“在我为您提供的示例代码中，10 是样本量。您可能已经注意到，玩具数据包含每个个体/单位的多个观察值（即每个单位 10 个观察值）。

Stan 手册还提供了有关函数参数和采样语句的大量信息。”

于 2014-10-21T20:57:09.150 回答
1

更新：我将 johnmyleswhite 示例翻译成 Stan Synthax 不起作用。我不太了解 Stan Synthax 来翻译代码。也许有人可以帮忙？以下是原答案。

如果您查看 jbaums 提到的johnmyleswhite 示例，您可以看到重要的代码片段是：

y[i] ~ dbern(p[i])
p[i] <- pt(z[i], 0, 1, 1)
z[i] <- a * x[i] + b

正如你所看到的，他使用 invlogit 来计算概率，他使用了 t 分布（实际上是累积的 t）。在 stan 中，只需使用：

student_t_cdf

我不太了解 Stan 合成器，但我认为您可以使用以下内容：

   model {
y ~ bernoulli(theta);
theta <- student_t_cdf(df, mu, sigma)
mu <- beta_0 + beta_1 * x1 + beta_2 * x2;
}

请注意，您必须将先验放在 df 和 sigma 上。像：

df_inv ~ uniform(0, 0.5);
df <- 1 / df_inv;
sigma_z <- sqrt((df-2)/df);

我会在这里尝试看看它是否有效。让我知道稍微调整一下我的答案是否可以使它起作用。

于 2014-10-21T01:19:58.097 回答
1

Stan 2.4 参考手册的第 26 页：

y ~ bernoulli(Phi( beta_0 + beta_1 * x1 + beta_2 * x2 ))

一般的解决方案是y ~ bernoulli(link_function(eta))where link_functionis，例如Phi。恰好有一个特殊的函数bernoulli_logit包装了这个功能并且在数值上更稳定。

如果原因不清楚，我建议阅读广义线性模型。维基百科页面并不是那么糟糕的评论。

于 2014-10-21T04:00:29.230 回答