Oprócz odpowiedzi @ gung postaram się podać przykład tego, co anova
funkcja faktycznie testuje. Mam nadzieję, że pozwoli ci to zdecydować, które testy są odpowiednie dla hipotez, które chcesz przetestować.
yx1x2)x3)my.mod <- glm(y~x1+x2+x3, family="binomial")
anova(my.mod, test="Chisq")
glm(y~1, family="binomial")
vs. glm(y~x1, family="binomial")
glm(y~x1, family="binomial")
vs. glm(y~x1+x2, family="binomial")
glm(y~x1+x2, family="binomial")
vs. glm(y~x1+x2+x3, family="binomial")
Dlatego sekwencyjnie porównuje mniejszy model z kolejnym bardziej złożonym modelem, dodając jedną zmienną na każdym etapie. Każde z tych porównań odbywa się za pomocą testu współczynnika wiarygodności (test LR; patrz przykład poniżej). O ile mi wiadomo, te hipotezy rzadko są interesujące, ale musisz o tym zadecydować.
Oto przykład w R
:
mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
mydata$rank <- factor(mydata$rank)
my.mod <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")
summary(my.mod)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.989979 1.139951 -3.500 0.000465 ***
gre 0.002264 0.001094 2.070 0.038465 *
gpa 0.804038 0.331819 2.423 0.015388 *
rank2 -0.675443 0.316490 -2.134 0.032829 *
rank3 -1.340204 0.345306 -3.881 0.000104 ***
rank4 -1.551464 0.417832 -3.713 0.000205 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# The sequential analysis
anova(my.mod, test="Chisq")
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 399 499.98
gre 1 13.9204 398 486.06 0.0001907 ***
gpa 1 5.7122 397 480.34 0.0168478 *
rank 3 21.8265 394 458.52 7.088e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# We can make the comparisons by hand (adding a variable in each step)
# model only the intercept
mod1 <- glm(admit ~ 1, data = mydata, family = "binomial")
# model with intercept + gre
mod2 <- glm(admit ~ gre, data = mydata, family = "binomial")
# model with intercept + gre + gpa
mod3 <- glm(admit ~ gre + gpa, data = mydata, family = "binomial")
# model containing all variables (full model)
mod4 <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")
anova(mod1, mod2, test="LRT")
Model 1: admit ~ 1
Model 2: admit ~ gre
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 399 499.98
2 398 486.06 1 13.92 0.0001907 ***
anova(mod2, mod3, test="LRT")
Model 1: admit ~ gre
Model 2: admit ~ gre + gpa
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 398 486.06
2 397 480.34 1 5.7122 0.01685 *
anova(mod3, mod4, test="LRT")
Model 1: admit ~ gre + gpa
Model 2: admit ~ gre + gpa + rank
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 397 480.34
2 394 458.52 3 21.826 7.088e-05 ***
psummary(my.mod)
- Dla współczynnika
x1
: glm(y~x2+x3, family="binomial")
vs.
glm(y~x1+x2+x3, family="binomial")
- Dla współczynnika
x2
: glm(y~x1+x3, family="binomial")
vs.glm(y~x1+x2+x3, family="binomial")
- Dla współczynnika
x3
: glm(y~x1+x2, family="binomial")
vs.glm(y~x1+x2+x3, family="binomial")
Zatem każdy współczynnik w stosunku do pełnego modelu zawierającego wszystkie współczynniki. Testy Walda są przybliżeniem testu współczynnika wiarygodności. Możemy również wykonać testy współczynnika wiarygodności (test LR). Oto jak:
mod1.2 <- glm(admit ~ gre + gpa, data = mydata, family = "binomial")
mod2.2 <- glm(admit ~ gre + rank, data = mydata, family = "binomial")
mod3.2 <- glm(admit ~ gpa + rank, data = mydata, family = "binomial")
anova(mod1.2, my.mod, test="LRT") # joint LR test for rank
Model 1: admit ~ gre + gpa
Model 2: admit ~ gre + gpa + rank
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 397 480.34
2 394 458.52 3 21.826 7.088e-05 ***
anova(mod2.2, my.mod, test="LRT") # LR test for gpa
Model 1: admit ~ gre + rank
Model 2: admit ~ gre + gpa + rank
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 395 464.53
2 394 458.52 1 6.0143 0.01419 *
anova(mod3.2, my.mod, test="LRT") # LR test for gre
Model 1: admit ~ gpa + rank
Model 2: admit ~ gre + gpa + rank
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 395 462.88
2 394 458.52 1 4.3578 0.03684 *
psummary(my.mod)
rank
anova(my.mod, test="Chisq")
rank
anova(mod1.2, my.mod, test="Chisq")
wartość jest taka sama,7.088 ⋅ 10- 5. Jest to każdorazowe porównanie między modelem bez rank
modelu zawierającego go.