Wahrscheinlichkeit und Schätzungen für gemischte Effekte Logistische Regression

Lassen Sie uns zunächst einige Daten für eine logistische Regression mit festen und zufälligen Teilen simulieren:

set.seed(1)
n    <- 100
x    <- runif(n)
z    <- sample(c(0,1), n, replace=TRUE)
b    <- rnorm(2)
beta <- c(0.4, 0.8)
X    <- model.matrix(~x)
Z    <- cbind(z, 1-z)
eta  <- X%*%beta + Z%*%b
pr   <- 1/(1+exp(-eta))
y    <- rbinom(n, 1, pr)

Wenn wir nur eine logistische Regression ohne zufällige Teile anpassen wollten, könnten wir die folgende glmFunktion verwenden:

glm(y~x, family="binomial")

glm(y~x, family="binomial")$coefficients
# (Intercept)           x 
#  -0.2992785   2.1429825 

Oder unsere eigene Funktion der Log-Wahrscheinlichkeit konstruieren

Dabei werden $\Lambda(\eta_i)=\frac{1}{1+\exp(-\eta_i)}$ und $\eta_i=\boldsymbol{X_i' \beta}$ verwendet optim(), um die Parameter zu schätzen, die es maximieren, wie im Folgenden Beispielcode:

ll.no.random <- function(theta,X,y){
  beta <- theta[1:ncol(X)]
  eta  <- X%*%beta
  p    <- 1/(1+exp(-eta))
  ll   <- sum( y*log(p) + (1-y)*log(1-p) )
  -ll
}

optim(c(0,1), ll.no.random, X=X, y=y)

optim(c(0,1), ll.no.random, X=X, y=y)$par
# -0.2992456  2.1427484

Dies liefert natürlich die gleichen Schätzungen und maximiert die Log-Wahrscheinlichkeit für den gleichen Wert. Für gemischte Effekte möchten wir so etwas wie

library(lme4)
glmer(y~x + (1|z), family="binomial")

Aber wie können wir dasselbe mit unserer eigenen Funktion tun? Da ist die Wahrscheinlichkeit

und das Integral hat keinen Ausdruck in geschlossener Form, wir müssen eine numerische Integration wie die Gaußsche Quadratur verwenden. Wir können das Paket verwenden statmod, um einige Quadraturen zu erhalten, sagen wir 10

library(statmod)    
gq <- gauss.quad(10)
w  <- gq$weights
g  <- gq$nodes

UPDATE: Unter Verwendung dieser Quadraturpositionen und der Gewichte für (hier ) können wir das Integral über durch eine Summe der Terme und das Ganze durch ersetzt werden Term multipliziert mit den jeweiligen Gewichten . Daher sollte unsere Wahrscheinlichkeitsfunktion jetzt sein$g_r$$w_r$$r=1,...,R$$R=10$$b_j$$R$$g_r$$b_j$$w_r$

Außerdem müssen wir die Varianz des zufälligen Teils berücksichtigen. Ich habe gelesen, dass dies erreicht werden kann, indem das in unserer Funktion durch wobei , also ersetzen wir in der obigen Wahrscheinlichkeitsfunktion tatsächlich 's durch ' s und nicht 's.$b_j \sim N(0,\sigma_b^2)$$\eta$$\sigma_j \theta_j$$\theta_j \sim N(0,1)$$\theta$$g$$\beta$

Ein Rechenproblem, das ich nicht bekomme, ist das Ersetzen der Begriffe, da die Vektoren nicht die gleiche Länge haben. Aber wahrscheinlich verstehe ich das nicht, weil mir hier etwas Entscheidendes fehlt oder ich falsch verstanden habe, wie diese Methode funktioniert.

Ich habe nicht gesehen, wie "die Vektoren nicht gleich lang sein werden", bitte klären Sie Ihre Frage.

Erstens sind für das Integral mit einer Dimension von weniger als 4 die direkten numerischen Methoden wie Quadratur effizienter als MCMC. Ich habe diese Fragen eine Weile studiert und würde dieses Problem gerne mit Ihnen besprechen.

Für die logistische Regression mit gemischten Effekten ist der einzige explizite RCode, den ich gefunden habe, aus Prof. Demidenkos Buch " Gemischte Modelle: Theorie und Anwendungen" . Sie können den Code über die Spalte "SOFTWARE UND DATEN" auf der Webseite herunterladen. Die logMLEgh()finden Sie in \mixed_models_data.zip\MixedModels\Chapter07. Er benutzte das statmodPaket nicht, um die Quadraturen zu erhalten, sondern schrieb seine eigene Funktion gauher(). Es gibt einige kleinere Fehler im Code und ich habe sie mit dem Autor besprochen, aber es ist immer noch sehr hilfreich, von seinem Code und seinem Buch auszugehen. Bei Bedarf kann ich die korrigierte Version bereitstellen.

Ein weiteres Problem ist, dass, wenn Sie genaue Schätzungen erhalten möchten, dies optim()nicht ausreicht, Sie möglicherweise Methoden wie Fisher Scoring verwenden müssen, wie in glm().