Teste da razão de verossimilhança em R

25

Suponha que eu faça uma regressão logística univariada em várias variáveis ​​independentes, como esta:

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

Fiz uma comparação de modelo (teste de razão de verossimilhança) para ver se o modelo é melhor que o modelo nulo por este comando

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

Então eu construí outro modelo com todas as variáveis ​​nele

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

Para ver se a variável é estatisticamente significativa no modelo multivariado, usei o lrtestcomando 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

Gostaria de saber se o pchisqmétodo e o lrtestmétodo são equivalentes para fazer o teste loglikelihood? Como eu não sei como usar lrtestpara o modelo de logística univada.

lokheart
fonte
@ Gavin obrigado por me lembrar que, comparando com o stackoverflow, preciso gastar mais tempo para "digerir" a resposta antes de decidir se a resposta é apropriada ou não, de qualquer forma, obrigado novamente.
Lokheart
Eu não recomendaria o uso do waldtest no lmtest. Use o pacote aod para teste de modelo. É muito mais direto. cran.r-project.org/web/packages/aod/aod.pdf
Sr. Nobody
epicalcfoi removido ( fonte ). Uma alternativa poderia ser lmtest.
22815 Martin Thoma

Respostas:

21

Basicamente, sim, desde que você use a diferença correta na probabilidade de log:

> 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

e não o desvio para o modelo nulo, que é o mesmo em ambos os casos. O número de df é o número de parâmetros que diferem entre os dois modelos aninhados, aqui df = 1. BTW, você pode procurar o código fonte lrtest()digitando apenas

> lrtest

no prompt R.

chl
fonte
obrigado, e acabei de descobrir que posso usar glm (output ~ NULL, data = z, family = binomial ("logistic")) para criar um modelo NULL e, portanto, usar o lrtest posteriormente. FYI, obrigado novamente
lokheart 25/01
2
@lokheart anova(model1, model0)também funcionará.
chl
5
@lokheart glm(output ~ 1, data=z, family=binomial("logistic"))seria um modelo nulo mais natural, que diz que outputé explicado por um termo constante (a interceptação) / A interceptação está implícita em todos os seus modelos, então você está testando o efeito de adepois de contabilizar a interceptação.
Reintegrar Monica - G. Simpson
Ou você pode fazê-lo "manualmente": valor de p do teste LR = 1-pchisq (desvio, dof)
Umka 24/01
22

Uma alternativa é o lmtestpacote, que possui uma lrtest()função que aceita um único modelo. Aqui está o exemplo de ?lrtestno lmtestpacote, que é para um LM mas existem métodos que trabalham com MLG:

> 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 
Restabelecer Monica - G. Simpson
fonte
+1 É bom saber (e parece que eu esqueci esse pacote).
chl
2
@GavinSimpson Isso pode parecer bobagem, mas como você interpretaria os resultados 'lrtest (fm2, fm1)'? O modelo 2 é significativamente diferente do modelo 1 e, portanto, a adição da variável con1 foi útil? Ou o lrtest (fm2) está dizendo que o modelo 2 é significativamente diferente do modelo 1? Mas qual modelo é melhor?
Krypton
5
O @Kerry fm1tem uma menor probabilidade de log e, portanto, um ajuste pior do que fm2. O LRT está nos dizendo que o grau em que criamos fm1um modelo mais pobre do que fm2é inesperadamente grande se os termos diferentes entre os modelos forem úteis (explicou a resposta). lrtest(fm2)Não é comparada com fm1a de tudo, o modelo fm2é comparado com, nesse caso, se, como indicado na saída, este: con ~ 1. Esse modelo, o modelo nulo, diz que o melhor preditor de coné a média da amostra de con(o termo intercepto / constante).
Reinstate Monica - G. Simpson