Test współczynnika wiarygodności w R.

25

Załóżmy, że zamierzam wykonać jednoczynnikową regresję logistyczną dla kilku niezależnych zmiennych, takich jak to:

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

Zrobiłem porównanie modelu (test współczynnika prawdopodobieństwa), aby sprawdzić, czy model jest lepszy niż model zerowy za pomocą tego polecenia

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

Następnie zbudowałem inny model ze wszystkimi zmiennymi

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

Aby sprawdzić, czy zmienna jest istotna statystycznie w modelu wielowymiarowym, użyłem lrtestpolecenia zepicalc

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

Zastanawiam się, czy pchisqmetoda i lrtestmetoda są równoważne z wykonaniem testu wiarygodności? Jak nie wiem, jak korzystać lrtestz univate modelu logistycznego.

lokheart
źródło
@Gavin dziękuję za przypomnienie, ponieważ w porównaniu z przepełnieniem stosu muszę poświęcić więcej czasu na „przetrawienie” odpowiedzi, zanim zdecyduję, czy odpowiedź jest odpowiednia, czy nie, dziękuję jeszcze raz.
lokheart
Nie polecałbym używania waldtest z lmtest. Użyj pakietu aod do testowania modelu. Jest o wiele prostszy. cran.r-project.org/web/packages/aod/aod.pdf
Pan Nikt
epicalczostał usunięty ( źródło ). Alternatywą może być lmtest.
Martin Thoma,

Odpowiedzi:

21

Zasadniczo tak, pod warunkiem, że użyjesz prawidłowej różnicy w prawdopodobieństwie dziennika:

> 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

a nie odchylenie dla modelu zerowego, który jest taki sam w obu przypadkach. Liczba df to liczba parametrów, które różnią się między dwoma modelami zagnieżdżonymi, tutaj df = 1. BTW, możesz spojrzeć na kod źródłowy, lrtest()po prostu wpisując

> lrtest

po znaku zachęty R.

chl
źródło
dzięki i właśnie odkryłem, że mogę używać glm (wyjście ~ NULL, dane = z, rodzina = dwumianowy („logistyczny”)) do tworzenia modelu NULL, dzięki czemu mogę później użyć lrtest. FYI, jeszcze raz dziękuję
lokheart
2
@lokheart też anova(model1, model0)będzie działać.
chl
5
@lokheart glm(output ~ 1, data=z, family=binomial("logistic"))byłby bardziej naturalnym modelem zerowym, który outputtłumaczy to stałym terminem (przechwytywanie) / Przechwytywanie jest implikowane we wszystkich twoich modelach, więc testujesz efekt apo uwzględnieniu przechwytywania.
Przywróć Monikę - G. Simpson
Lub możesz to zrobić „ręcznie”: wartość p testu LR = 1-pchisq (dewiacja, dof)
Umka
22

Alternatywą jest lmtestpakiet, który ma lrtest()funkcję, która akceptuje jeden model. Oto przykład ze ?lrtestw lmtestopakowaniu, które jest na LM, ale istnieją sposoby, że praca z GLMs:

> 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 
Przywróć Monikę - G. Simpson
źródło
+1 Dobrze wiedzieć (i wygląda na to, że zapomniałem o tym pakiecie).
chl
2
@GavinSimpson To może wydawać się głupie, ale jak interpretowałbyś wyniki „lrtest (fm2, fm1)”? Model 2 różni się znacznie od modelu 1 i dlatego dodanie zmiennej con1 było przydatne? Czy też lrtest (fm2) mówi, że model 2 różni się znacznie od modelu 1? Ale który model jest lepszy?
Kerry
5
@Kerry fm1ma niższe prawdopodobieństwo dziennika, a zatem gorsze dopasowanie niż fm2. LRT mówi nam, że stopień, w jakim stworzyliśmy fm1gorszy model, fm2jest nieoczekiwanie duży, jeśli użyteczne byłyby terminy różniące się między modelami (wyjaśnił odpowiedź). lrtest(fm2)nie jest w porównaniu z fm1wcale, model fm2jest w porównaniu z tym w przypadku, gdy, jak stwierdzono na wyjściu, to: con ~ 1. Ten model, model zerowy, mówi, że najlepszym predyktorem conjest średnia próbkowa con(punkt przecięcia / stały).
Przywróć Monikę - G. Simpson