r - 如何在Stan中运行robit模型？

Question

我想在斯坦(Stan)中运行健壮的逻辑回归(robit)。该模型在Gelman＆Hill的“使用回归和多级方法进行数据分析”(2006年，第124页)中提出，但是我不确定如何实现。我检查了Stan's Github repository和the reference manual，但不幸的是我仍然感到困惑。这是我用来对常规逻辑回归建模的一些代码。我应该添加什么以使误差遵循具有7个自由度的t分布？如果我进行多级分析，是否有相同的过程？

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)

谢谢!

Answer 1

我一定想念一些东西，但是我无法适应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)