스탠에서 로빗 모델을 실행하는 방법은 무엇입니까?

Question

Stan에서 강력한 로지스틱 회귀 (robit)를 실행하고 싶습니다. 모델은 Gelman & Hill의 "회귀 및 다단계 방법을 사용한 데이터 분석"(2006, pp. 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 솔루션을 채택했다

. 그래서 난 그냥에서 모델을 번역 장애가 있었던 경우.

나는

이

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) 

Answer 2

업데이트 : johnmyleswhite 예제를 Stan Synthax로 번역 한 것이 작동하지 않습니다. 나는 코드를 번역하기 위해 Stan Synthax를 잘 이해하지 못한다. 어쩌면 누군가가 도울 수 있습니까? 아래는 원래의 대답입니다.

는 jbaums에 의해 언급 된 johnmyleswhite example을 선택하면, 당신은 코드의 중요한 부분이라고 볼 수 있습니다

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

당신이 볼 수 있듯이이 probabilites을 계산하는 invlogit를 사용 insted, 그는 t 분포를 사용 (실제로 누적 t).

이

student_t_cdf 

내가 참 잘 스탠 synthax 모르지만, 난 당신이 다음과 같은 것을 사용할 수 있습니다 가정 : 스탠에서 바로 사용

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

참고가 전과를 넣어해야합니다 df 및 시그마. Something like :

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

나는 그것이 작동하는지 여기에서 시험해 볼 것이다. 내 대답을 약간 조정하면 효과가 있음을 알려주세요. 스탠 2.4 레퍼런스 매뉴얼

Answer 3

페이지 26 : link_function이고

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

일반적인 해결책은, 예를 y ~ bernoulli(link_function(eta)) 인 Phi. 이 기능을 랩핑하고 수치 적으로 더 안정적인 특별한 함수 bernoulli_logit이 있습니다.

일반 선형 모델을 읽는 것이 좋습니다. Wikipedia 페이지는 그렇게 나쁜 리뷰는 아닙니다.

Answer 4

Luc Coffeng이 Stan mailing list에서이 대답을 보내 주었고 여기에 추가해야한다고 생각했습니다. 그는 말했습니다 :

"로 비트 회귀의 기본으로 GLM을 채택하십시오 : e ~ student_t(7, 0, sigma_e)으로 대체하십시오. sigma_e ~ cauchy(0, 2) 또는 어떤 스케일로도 괜찮습니까? (어쨌든 나는 inverse logit of (-5,5)는 [0,1] 간격의 대부분을 커버합니다 .t- 오차의 척도 이외에, t- 오차의 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) 

밥 카펜터는 간단히 질문에 댓글을 달았습니다 :

그는 다음 예제를 제공

"[...] 그리고 그래, 당신은 계층에서 같은 일을 할 수 설정을해야만하지만 자유도를 모델링하는 것이 정상에 다다를 때 무한대로 도망 가기 때문에 자유도를 모델링하는 것이 까다로울 수 있으므로주의해야합니다. "라고 Bernd의 질문에 대답하면서 Luc은 모델 코드에 y ~ bernoulli_logit(10...을 쓴 이유를 분명히했습니다.

"내가 제공 한 예제 코드에서 10은 샘플 크기입니다. 장난감 데이터에 개인/단위 (즉, 단위당 10 회의 관찰)마다 여러 개의 관찰이 포함되어 있음을 알 수 있습니다.

스탠 설명서는 함수의 인수 및 샘플링 문에 대한 광범위한 정보를 제공합니다. "