Likelihood-Ratio-Test in R

25

Angenommen, ich werde eine univariate logistische Regression für mehrere unabhängige Variablen durchführen:

mod.a <- glm(x ~ a, data=z, family=binominal("logistic"))
mod.b <- glm(x ~ b, data=z, family=binominal("logistic"))

Ich habe einen Modellvergleich (Likelihood Ratio Test) durchgeführt, um festzustellen, ob das Modell mit diesem Befehl besser ist als das Nullmodell

1-pchisq(mod.a$null.deviance-mod.a$deviance, mod.a$df.null-mod.a$df.residual)

Dann baute ich ein anderes Modell mit allen Variablen darin

mod.c <- glm(x ~ a+b, data=z, family=binomial("logistic"))

Um festzustellen, ob die Variable im multivariaten Modell statistisch signifikant ist, habe ich den lrtestBefehl von verwendetepicalc

lrtest(mod.c,mod.a) ### see if variable b is statistically significant after adjustment of a
lrtest(mod.c,mod.b) ### see if variable a is statistically significant after adjustment of b

Ich frage mich, ob die pchisqMethode und die lrtestMethode für die Durchführung des Loglikelihood-Tests gleichwertig sind. Ich weiß nicht, wie ich es lrtestfür das Univate-Logistikmodell verwenden soll.

lokheart
quelle
@Gavin, danke, dass du mich daran erinnert hast, dass ich im Vergleich zu Stackoverflow mehr Zeit aufwenden muss, um die Antwort zu "verdauen", bevor ich mich entscheide, ob die Antwort angemessen ist oder nicht. Trotzdem, danke noch einmal.
Lokheart
Ich würde nicht empfehlen, Waldtest von lmtest zu verwenden. Verwenden Sie das aod-Paket zum Testen von Modellen. Es ist viel einfacher. cran.r-project.org/web/packages/aod/aod.pdf
Mr. Nobody
epicalcwurde entfernt ( Quelle ). Eine Alternative könnte sein lmtest.
Martin Thoma

Antworten:

21

Grundsätzlich ja, vorausgesetzt, Sie verwenden den richtigen Unterschied in der Log-Wahrscheinlichkeit:

> library(epicalc)
> model0 <- glm(case ~ induced + spontaneous, family=binomial, data=infert)
> model1 <- glm(case ~ induced, family=binomial, data=infert)
> lrtest (model0, model1)
Likelihood ratio test for MLE method 
Chi-squared 1 d.f. =  36.48675 , P value =  0 
> model1$deviance-model0$deviance
[1] 36.48675

und nicht die Abweichung für das Nullmodell, die in beiden Fällen gleich ist. Die Anzahl von df ist die Anzahl von Parametern, die sich zwischen den beiden verschachtelten Modellen unterscheiden, hier ist df = 1. Übrigens, Sie können sich den Quellcode ansehen, lrtest()indem Sie einfach etwas eingeben

> lrtest

an der R-Eingabeaufforderung.

chl
quelle
Danke, und ich habe gerade herausgefunden, dass ich glm (output ~ NULL, data = z, family = binomial ("logistic")) zum Erstellen eines NULL-Modells verwenden kann, um anschließend den lrtest zu verwenden. Zu
Ihrer Information
2
@lokheart anova(model1, model0)wird auch funktionieren.
chl
5
@lokheart glm(output ~ 1, data=z, family=binomial("logistic"))wäre ein natürlicheres Nullmodell, das besagt, dass dies outputdurch einen konstanten Term (den Achsenabschnitt) erklärt wird. / Der Achsenabschnitt ist in all Ihren Modellen impliziert a.
Setzen Sie Monica - G. Simpson am
Oder Sie können es "manuell" machen: p-Wert des LR-Tests = 1-pchisq (Abweichung, dof)
Umka
22

Eine Alternative ist das lmtestPaket, das eine lrtest()Funktion hat, die ein einzelnes Modell akzeptiert. Hier ist das Beispiel aus ?lrtestdem lmtestPaket, das sich auf einen LM bezieht, aber es gibt Methoden, die mit GLMs funktionieren:

> require(lmtest)
Loading required package: lmtest
Loading required package: zoo
> ## with data from Greene (1993):
> ## load data and compute lags
> data("USDistLag")
> usdl <- na.contiguous(cbind(USDistLag, lag(USDistLag, k = -1)))
> colnames(usdl) <- c("con", "gnp", "con1", "gnp1")
> fm1 <- lm(con ~ gnp + gnp1, data = usdl)
> fm2 <- lm(con ~ gnp + con1 + gnp1, data = usdl)
> ## various equivalent specifications of the LR test
>
> ## Compare two nested models
> lrtest(fm2, fm1)
Likelihood ratio test

Model 1: con ~ gnp + con1 + gnp1
Model 2: con ~ gnp + gnp1
  #Df  LogLik Df  Chisq Pr(>Chisq)    
1   5 -56.069                         
2   4 -65.871 -1 19.605  9.524e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 
>
> ## with just one model provided, compare this model to a null one
> lrtest(fm2)
Likelihood ratio test

Model 1: con ~ gnp + con1 + gnp1
Model 2: con ~ 1
  #Df   LogLik Df  Chisq Pr(>Chisq)    
1   5  -56.069                         
2   2 -119.091 -3 126.04  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 
Setzen Sie Monica - G. Simpson wieder ein
quelle
+1 Gut zu wissen (und anscheinend habe ich dieses Paket vergessen).
chl
2
@GavinSimpson Das mag albern erscheinen, aber wie würden Sie die Ergebnisse von 'lrtest (fm2, fm1)' interpretieren? Modell 2 unterscheidet sich erheblich von Modell 1 und daher war die Hinzufügung der Variablen con1 sinnvoll? Oder der lrtest (fm2) besagt, dass Modell 2 sich erheblich von Modell 1 unterscheidet? Aber welches Modell ist besser?
Kerry
5
@ Kerry fm1hat eine geringere Log-Wahrscheinlichkeit und damit eine schlechtere Passform als fm2. Das LRT sagt uns, dass der Grad, in dem wir fm1ein ärmeres Modell erstellt haben, fm2unerwartet groß ist, wenn die Begriffe, die sich zwischen den Modellen unterscheiden, nützlich waren (erklärte die Reaktion). lrtest(fm2)nicht gegenüber fm1überhaupt, das Modell fm2ist mit in diesem Fall verglichen , wenn, wie dies in der Ausgabe angegeben: con ~ 1. Dieses Modell, das Nullmodell, sagt aus, dass der beste Prädiktor conder Stichprobenmittelwert con(der Achsenabschnitt / konstante Term) ist.
Setzen Sie Monica - G. Simpson