Test du rapport de vraisemblance dans R

25

Supposons que je vais faire une régression logistique univariée sur plusieurs variables indépendantes, comme ceci:

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

J'ai fait une comparaison de modèle (test de rapport de vraisemblance) pour voir si le modèle est meilleur que le modèle nul par cette commande

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

Ensuite, j'ai construit un autre modèle avec toutes les variables

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

Afin de voir si la variable est statistiquement significative dans le modèle multivarié, j'ai utilisé la lrtestcommande deepicalc

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

Je me demande si la pchisqméthode et la lrtestméthode sont équivalentes pour faire un test de loglikelihood? Comme je ne sais pas comment utiliser lrtestpour le modèle logistique univate.

lokheart
la source
@Gavin merci de me le rappeler, car en comparant avec stackoverflow, j'ai besoin de passer plus de temps à "digérer" la réponse avant de décider si la réponse est appropriée ou non, de toute façon, merci encore.
lokheart
Je ne recommanderais pas d'utiliser waldtest de lmtest. Utilisez le package aod pour tester les modèles. C'est beaucoup plus simple. cran.r-project.org/web/packages/aod/aod.pdf
Mr. Nobody
epicalca été supprimé ( source ). Une alternative pourrait être lmtest.
Martin Thoma

Réponses:

21

Fondamentalement, oui, à condition d'utiliser la bonne différence de log-vraisemblance:

> 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

et non la déviance pour le modèle nul qui est le même dans les deux cas. Le nombre de df est le nombre de paramètres qui diffèrent entre les deux modèles imbriqués, ici df = 1. BTW, vous pouvez regarder le code source lrtest()en tapant simplement

> lrtest

à l'invite R.

chl
la source
merci, et je viens de découvrir que je peux utiliser glm (output ~ NULL, data = z, family = binomial ("logistic")) pour créer un modèle NULL, et donc je peux utiliser lrtest ensuite. Pour info, merci encore
lokheart
2
@lokheart anova(model1, model0)fonctionnera également.
chl
5
@lokheart glm(output ~ 1, data=z, family=binomial("logistic"))serait un modèle nul plus naturel, qui dit que cela outputest expliqué par un terme constant (l'interception) / L'interception est implicite dans tous vos modèles, vous testez donc l'effet de la prise en acompte de l'interception.
Rétablir Monica - G. Simpson le
Ou vous pouvez le faire "manuellement": valeur de p du test LR = 1-pchisq (déviance, ddl)
Umka
22

Une alternative est le lmtestpackage, qui a une lrtest()fonction qui accepte un seul modèle. Voici l'exemple de ?lrtestdans le lmtestpackage, qui est pour un LM mais il existe des méthodes qui fonctionnent avec les GLM:

> 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 
Réintégrer Monica - G. Simpson
la source
+1 C'est bon à savoir (et il semble que j'aie oublié ce package).
chl
2
@GavinSimpson Cela peut sembler idiot, mais comment interpréteriez-vous les résultats 'lrtest (fm2, fm1)'? Le modèle 2 est significativement différent du modèle 1 et donc l'ajout de la variable con1 était utile? Ou le lrtest (fm2) dit que le modèle 2 est significativement différent du modèle 1? Mais quel modèle est le meilleur?
Kerry
5
@Kerry fm1a une probabilité de log plus faible et donc un ajustement moins bon que fm2. Le LRT nous dit que le degré auquel nous avons fait fm1un modèle plus pauvre fm2qu'inattendu est grand si les termes qui sont différents entre les modèles étaient utiles (a expliqué la réponse). lrtest(fm2)n'est pas par rapport fm1à l' ensemble, le modèle fm2est comparé dans ce cas si, comme indiqué dans la sortie, ceci: con ~ 1. Ce modèle, le modèle nul, dit que le meilleur prédicteur de conest la moyenne de l'échantillon de con(le terme interception / constant).
Rétablir Monica - G. Simpson