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

8

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

logL(β)=i=1nyilogΛ(ηi)+(1yi)log(1Λ(ηi))

Dabei werden Λ(ηi)=11+exp(ηi) und ηi=Xiβ 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

L=j=1JPr(y1j,...,ynjj|x,bj)f(bj)dbj

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 seingrwrr=1,...,RR=10bjRgrbjwr

L=j=1Jr=1RPr(y1j,...,ynjj|x,bj=gr)wr

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.bjN(0,σb2)ησjθjθjN(0,1)θgβ

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.

Steve
quelle
Ich habe jetzt keine Zeit, einen guten Blick darauf zu werfen, aber dies scheint eine gute Verwendung für MCMC zu sein.
Shadowtalker
@ssdecontrol danke, MCMC ist eine gute Alternative. Aber ich möchte den klassischen Ansatz anwenden.
Steve
Was ist nicht klassisch an MCMC, um ein Integral zu bewerten?
Shadowtalker
@ssdecontrol Nun, ich bin mir nicht sicher ... Aber alle Bücher, die ich überprüft habe, und die lme4, ordinale R-Pakete, verwenden kein MCMC. Also möchte ich dabei bleiben. Zumindest am Anfang.
Steve
1
Haben Sie dies in der R-sig-ME-Liste ([email protected]) gefragt? Einige Leute könnten Ihnen dort möglicherweise weiterhelfen. Darüber hinaus: Ich möchte Sie dringend bitten, die Arbeit Efficient Laplacian and Adaptive Gaussian Quadrature Algorithms für mehrstufige generalisierte lineare gemischte Modelle von Pinheiro und Chao zu studieren . Es enthält Ergebnisse zur AGQ-Leistung sowie zur dahinter stehenden linearen Algebra. Wenn Sie sie codieren möchten, machen Sie sich bereit für einige ernsthafte Zerlegungen mit spärlicher Matrix. : D
usεr11852

Antworten:

2

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().

Randel
quelle
Das Buch scheint reichhaltige Informationen darüber zu haben, woran ich arbeite. Der Code selbst sagt nicht viel aus - ich habe gerade das Buch bestellt, also muss ich darauf warten. Die Sache mit den Vektoren ist, dass, wenn wir die Terme, die 2 bdurch die 10 Knoten ersetzen , wie wir die Matrizen multiplizieren können Zund g? Oder habe ich es völlig falsch?
Steve
Ich weiß, dass ich noch weiter gehen sollte, um genaue Schätzungen zu erhalten, aber ich hatte gehofft, zuerst den GQ als ersten Schritt zu verstehen.
Steve
Sie können zuerst eine Vorschau der beiden Editionen in Google Books anzeigen. In Ihrem Code ist eine Matrix? Aber nur einen skalaren? Sie können zuerst mit dem Zufallsschnittmodell beginnen . Wenn die Dimension der zufälligen Effekte zwei ist, benötigen Sie insgesamt 100 zweidimensionale Knoten, dann 10 Knoten in jeder Dimension. Zn×2bZ = rep(1,n)
Randel
Entschuldigung, aber je mehr ich darüber nachdenke, desto mehr verwirrt es mich. Wenn ich das mache, Z=rep(1,n)bekomme ich einen zufälligen Abschnitt für jede Zeile, was bedeutet, dass jede Person eine Gruppe ist? In meinem Beispiel habe ich zwei Gruppen, also haben wir und , um die zu geben, die wir brauchen. Nein? Z%*%b n×22×1n×1
Steve
Oh, ich habe gerade bemerkt, dass dies Zverwendet wird, um den zufälligen Achsenabschnitt für jeden Cluster zu trennen, nicht die Entwurfsmatrix für den zufälligen Effekt. Dann haben Sie Recht, aber Sie sollten das Integral auswerten und die Quadratur für jeden Cluster separat verwenden. Sie brauchen nicht Zmehr, bewerten Sie einfach das Integral für jeden Cluster und addieren Sie sie dann. Die Entwurfsmatrix für den zufälligen Achsenabschnitt ist nur der Vektor von 1s.
Randel