Wie führe ich ein Robit-Modell in Stan?

Question

Ich würde gerne eine robuste logistische Regression (robit) in Stan ausführen. Das Modell wird in Gelman & Hill's "Datenanalyse mit Regression und Multilevel-Methoden" (2006, S. 124) vorgeschlagen, aber ich bin mir nicht sicher, wie es zu implementieren ist. Ich habe Stan's Github repository und überprüft, aber leider bin ich immer noch verwirrt. Hier ist ein Code, den ich verwendet habe, um eine reguläre logistische Regression zu modellieren. Was soll ich hinzufügen, damit die Fehler etwa einer t-Verteilung mit 7 Freiheitsgraden folgen? Könnte es bei einer mehrstufigen Analyse zufällig genauso sein?

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) 

Vielen Dank!

Answer 1

ich etwas fehlen muss, aber ich hatte die Lösung Probleme der Anpassung, die von Luc geschrieben danilofreire. Also von I übersetzte nur ein Modell JAGS.

ich denke, das ist richtig, obwohl es ein bisschen anders als Luc-Lösung suchen.

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

Update: Meine Übersetzung von johnmyleswhite Beispiel zu Stan Synthax funktioniert nicht. Ich verstehe Stan Stan Synth nicht sehr gut, um den Code zu übersetzen. Vielleicht kann jemand helfen? Unten ist die ursprüngliche Antwort.

Wenn Sie die johnmyleswhite example überprüfen, von jbaums erwähnt, können Sie sehen, dass der wichtige Teil des Codes ist:

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

Wie Sie sehen können, insted invlogit der Verwendung von Wahrscheinlichkeiten zu berechnen, nutzt er die t-Verteilung (eigentlich, das kumulative t). In stan, benutzen Sie einfach:

student_t_cdf 

Ich weiß nicht, Vey gut Stan Synthax, aber ich nehme an, Sie so etwas wie die Folge verwenden können:

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

Beachten Sie, dass priors setzen müssen auf df und sigma. Something like:

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

Ich werde hier versuchen, um zu sehen, ob es funktioniert. Lass es mich wissen, wenn ich meine Antwort ein wenig verbessern kann.

Answer 3

Seite 26 des Stan 2.4 Referenzhandbuch:

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

Die allgemeine Lösung ist y ~ bernoulli(link_function(eta)) wo link_function ist, z.B. Phi. Es gibt nur eine spezielle Funktion bernoulli_logit, die diese Funktionalität umschließt und numerisch stabiler ist.

Ich empfehle, generalisierte lineare Modelle zu lesen, wenn der Grund dafür nicht klar ist. Die Wikipedia-Seite ist keine so schlechte Bewertung.

Answer 4

Luc Coffeng schickte mir diese Antwort auf die Stan mailing list und ich dachte, ich sollte es hier hinzufügen. Er sagte:

„ein GLM Adopt als Basis für Ihre robit Regression: einfach ersetzt den Standard-Fehlerterm mit e ~ student_t(7, 0, sigma_e), wo sigma_e ~ cauchy(0, 2) oder was auch immer Sie denken, Skala ok ist (ich würde nicht mehr als 5 ohnehin als inversen logit gehen (-5,5) deckt den größten Teil des [0,1] Intervalls ab.Zusätzlich zum Maßstab des t-Fehlers können Sie auch den df des t-Fehlers als Parameter angeben.Siehe unten für den vorgeschlagenen Code.

Ich hoffe jedoch, dass Ihre Daten mehr Informationen enthalten als das von Ihnen bereitgestellte Spielzeugbeispiel, dh mehrere Beobachtungen pro Individuum (wie unten): Mit nur einer Beobachtung pro Individuum/Einheit ist das Modell praktisch nicht zu identifizieren."

Er hat dann das folgende Beispiel:

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 auch auf die Frage kurz kommentiert:

" [...] Und ja, können Sie die gleiche Sache in einer hierarchischen tun Einstellung, aber Sie müssen vorsichtig sein, weil das Modellieren der Freiheitsgrade schwierig sein kann, da die Skala in Richtung Unendlichkeit abläuft, wenn Sie sich der Normalität nähern. "

Auf Bernds Frage antwortete Luc, warum er y ~ bernoulli_logit(10... in dem Modellcode schrieb:

"In dem Beispielcode, den ich Ihnen zur Verfügung gestellt habe, ist 10 die Stichprobengröße. Sie haben vielleicht bemerkt, dass die Spielzeugdaten mehrere Beobachtungen pro Individuum/Einheit enthalten (d. H. 10 Beobachtungen pro Einheit).

Das Stan Handbuch enthält auch umfangreiche Informationen über Argumente an Funktionen und Sampling-Aussagen“.