Prueba de razón de verosimilitudes en R

Question

Supongamos que voy a hacer una regresión logística univariante sobre varias variables independientes, así:

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

Hice una comparación de modelos (prueba de razón de verosimilitud) para ver si el modelo es mejor que el modelo nulo mediante este comando

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

Luego construí otro modelo con todas las variables

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

Para ver si la variable es estadísticamente significativa en el modelo multivariante, he utilizado el lrtest comando de epicalc

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

Me pregunto si el pchisq y el método lrtest son equivalentes para hacer la prueba de loglikelihood? Como no sé cómo usar lrtest para el modelo logístico unívoco.

Answer 1

Básicamente, sí, siempre que se utilice la diferencia correcta en la probabilidad logarítmica:

> 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

y no la desviación para el modelo nulo que es el mismo en ambos casos. El número de df es el número de parámetros que difieren entre los dos modelos anidados, aquí df=1. Por cierto, puede consultar el código fuente de lrtest() con sólo teclear

> lrtest

en el indicador R.

Answer 2

Una alternativa es el lmtest que tiene un lrtest() que acepta un único modelo. Este es el ejemplo de ?lrtest en el lmtest paquete, que es para un LM pero el

> 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