Logistische Regression mit kategorischen Prädiktoren mit JAGS

Question

Ich bin neu in JAGS und ich versuche, ein binäres Ergebnis (0/1) mit 9 nicht-kontinuierliche Prädiktoren vorherzusagen. Prädiktorwerte können 0, 1 oder 2 sein. Dies ist mein erstes Mal, und obwohl ich das Modell zum Laufen bringen kann, bin ich mir 100% sicher, dass es hier eine Reihe von Problemen gibt.

Datendatei Probe (Liste)

$y 
[1] 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 

$N 
[1] 50 

$oAnt 
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 
[29] 1 1 1 1 1 1 2 1 1 0 1 1 1 1 1 2 1 1 1 1 1 1 

$nAnt 
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
[29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

$cAnt 
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 

$oPen 
[1] 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 0 1 1 1 1 1 1 0 
[29] 1 0 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 0 2 1 1 

$nPen 
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 2 1 
[29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 

$cPen 
[1] 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0 
[29] 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 

$oFin 
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 
[29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

$nFin 
[1] 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1 
[29] 1 1 1 1 2 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1 3 

$cFin 
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 
[29] 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 

Modell

model { 
    for(i in 1 : N){ 
     y[i] ~ dbern(mu[i]) 
     mu[i] <- 1/(1+exp(-(b0 + b1*oAnt[i] + b2*nAnt[i] + b3*cAnt[i] + b4*oPen[i] + b5*nPen[i] + b6*cPen[i] + b7*oFin[i] + b8*nFin[i] + b9*cFin[i]))) 
} 
b0 ~ dnorm(0, 1.0e-12) 
b1 ~ dnorm(0, 1.0e-12) 
b2 ~ dnorm(0, 1.0e-12) 
b3 ~ dnorm(0, 1.0e-12) 
b4 ~ dnorm(0, 1.0e-12) 
b5 ~ dnorm(0, 1.0e-12) 
b6 ~ dnorm(0, 1.0e-12) 
b7 ~ dnorm(0, 1.0e-12) 
b8 ~ dnorm(0, 1.0e-12) 
b9 ~ dnorm(0, 1.0e-12) 
} 

Ich benutzte Schätzungen von einem glm() Modell als INiTS (wie von A. Gelman vorschlagen) -aber aus Gründen der Einfachheit, Nehmen wir an, ich lasse JAGS die Anfangswerte für die Ketten auswählen.

Lauf Modell usw.

jagsModel = jags.model(file = "antPenFin.txt", data = dataList, n.chains = 2, n.adapt = 500) 

update(jagsModel, n.iter = 500) 

codaSamples = coda.samples(jagsModel, 
variable.names = names(dataList)[3:11], n.iter = 5000) 

Probleme

Der Ausgang meines Modells sieht völlig aus (was deutlich wird, wenn ich versuche, es zu zeichnen). Ich bin sicher, dass es hier ein sehr grundlegendes Problem gibt. Könnte jemand helfen?

Vielen Dank.

Answer 1

Ihr Problem hier ist, dass Sie eine Reihe von Zahlen (von 0 bis 3 in Ihrem Fall) nicht als kategorische Kovariaten angeben können. Momentan interpretiert Ihr Modell diese Zahlen als fortlaufend. Was Sie tun müssen, ist diese in dummy variables konvertieren. Sie können dies einfach mit der Funktion model.matrix in R tun. Ich werde als Beispiel einige Daten generieren, die Sie dann auf Ihre Daten anwenden können.

# generate y 
y <- sample(0:1, 30, replace = TRUE) 
# add three different categorical covariates. Note here that these are all 
# factors. 
oAnt <- factor(sample(0:2, 30, replace = TRUE)) 
cAnt <- factor(sample(0:1, 30, replace = TRUE)) 
nFin <- factor(sample(0:3, 30, replace = TRUE)) 

# create your model matrix 
my_matrix <- model.matrix(y ~ oAnt + cAnt + nFin) 

head(my_matrix) 

    (Intercept) oAnt1 oAnt2 cAnt1 nFin1 nFin2 nFin3 
1   1  1  0  1  0  1  0 
2   1  1  0  1  1  0  0 
3   1  1  0  0  0  0  1 
4   1  1  0  1  0  0  1 
5   1  0  0  1  0  1  0 
6   1  1  0  1  0  0  1 

Für jede Ebene in einem Faktor für einen Prädiktor, werden Sie n-1 Spalten erstellen müssen (zum Beispiel NFIN im Bereich von 0-3 ist, werden Sie 3 Spalten von Dummy-Variablen erstellen). Daher werden Sie viel mehr Regressionskoeffizienten in Ihrem Modell haben. Sobald Sie Ihre Modellmatrix erstellt haben, können Sie sie in eine Liste als solche konvertieren.

# remove the intercept from the list as you already have it in your model 
matrix_list <- as.list(data.frame(my_matrix[,-1])) 

Von dort müssen Sie nur noch weitere Prädiktoren für das Modell erstellen. Auch, wenn nAnt in Ihrem Modell wirklich alle sind, gehen Sie voran und entfernen Sie so gut wie Sie im Grunde nur zwei Abschnitte enthalten.