Jak przetestować równość wybranych współczynników w modelu logit lub probit? Jakie jest standardowe podejście i jakie jest najnowsze podejście?
Jak przetestować równość wybranych współczynników w modelu logit lub probit? Jakie jest standardowe podejście i jakie jest najnowsze podejście?
Jednym standardowym podejściem jest test Walda . To właśnie robi polecenie Stata test
po regresji logit lub probit. Zobaczmy, jak to działa w R, patrząc na przykład:
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv") # Load dataset from the web
mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial") # calculate the logistic regression
summary(mylogit)
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
Powiedzmy, że chcesz przetestować hipotezę vs. . Jest to równoważne testowaniu . Statystyka testu Walda to: β g r e ≠ β g p a β g r e - β g p a = 0
lub
Nasz tutaj to i . Potrzebujemy więc tylko standardowego błędu . Możemy obliczyć błąd standardowy za pomocą metody Delta : βgre-βgpθ0=0βgre-βgp
Potrzebujemy także kowariancji i . Macierz wariancji-kowariancji można wyodrębnić za pomocą polecenia po uruchomieniu regresji logistycznej: β g p avcov
var.mat <- vcov(mylogit)[c("gre", "gpa"),c("gre", "gpa")]
colnames(var.mat) <- rownames(var.mat) <- c("gre", "gpa")
gre gpa
gre 1.196831e-06 -0.0001241775
gpa -1.241775e-04 0.1101040465
Wreszcie możemy obliczyć błąd standardowy:
se <- sqrt(1.196831e-06 + 0.1101040465 -2*-0.0001241775)
se
[1] 0.3321951
Więc twoja wartość Wald jest
wald.z <- (gre-gpa)/se
wald.z
[1] -2.413564
Aby uzyskać wartość , wystarczy użyć standardowego rozkładu normalnego:
2*pnorm(-2.413564)
[1] 0.01579735
W tym przypadku mamy dowody, że współczynniki różnią się od siebie. To podejście można rozszerzyć na więcej niż dwa współczynniki.
Za pomocą multcomp
To dość żmudne obliczenia można wygodnie wykonać przy R
użyciu multcomp
pakietu. Oto ten sam przykład, co powyżej, ale wykonano za pomocą multcomp
:
library(multcomp)
glht.mod <- glht(mylogit, linfct = c("gre - gpa = 0"))
summary(glht.mod)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
gre - gpa == 0 -0.8018 0.3322 -2.414 0.0158 *
confint(glht.mod)
Przedział ufności dla różnicy współczynników można również obliczyć:
Quantile = 1.96
95% family-wise confidence level
Linear Hypotheses:
Estimate lwr upr
gre - gpa == 0 -0.8018 -1.4529 -0.1507
Aby zobaczyć dodatkowe przykłady multcomp
, zobacz tutaj lub tutaj .
Współczynniki regresji logistycznej znajdują się na podstawie największego prawdopodobieństwa. Ponieważ jednak funkcja wiarygodności obejmuje wiele produktów, prawdopodobieństwo logarytmiczne jest zmaksymalizowane, co zamienia produkty w sumy. Model, który lepiej pasuje, ma większe prawdopodobieństwo dziennika. Model obejmujący więcej zmiennych ma co najmniej takie samo prawdopodobieństwo jak model zerowy. Oznacz logarytmiczne prawdopodobieństwo modelu alternatywnego (model zawierający więcej zmiennych) z i logarytmiczne prawdopodobieństwo modelu zerowego z , statystyka testu współczynnika wiarygodności wynosi: L L 0
Statystyka testu prawdopodobieństwa jest zgodna z przy czym stopnie swobody są różnicą w liczbie zmiennych. W naszym przypadku jest to 2.
Aby wykonać test współczynnika wiarygodności, musimy również dopasować model do ograniczenia aby móc porównać dwa prawdopodobieństwa. Pełny model ma postać . Nasz model ograniczeń ma postać: . log ( p i
mylogit2 <- glm(admit ~ I(gre + gpa) + rank, data = mydata, family = "binomial")
W naszym przypadku możemy użyć logLik
do wyodrębnienia prawdopodobieństwa logarytmicznego dwóch modeli po regresji logistycznej:
L1 <- logLik(mylogit)
L1
'log Lik.' -229.2587 (df=6)
L2 <- logLik(mylogit2)
L2
'log Lik.' -232.2416 (df=5)
Model zawierający ograniczenie gre
i gpa
ma nieco większe prawdopodobieństwo logarytmu (-22,24) w porównaniu do pełnego modelu (-229.26). Nasza statystyka testu współczynnika wiarygodności to:
D <- 2*(L1 - L2)
D
[1] 16.44923
Możemy teraz użyć CDF z aby obliczyć wartość :
1-pchisq(D, df=1)
[1] 0.01458625
Wartość jest bardzo mała, co wskazuje, że współczynniki są różne.
R ma wbudowany test współczynnika prawdopodobieństwa; możemy użyć anova
funkcji do obliczenia testu współczynnika wiarygodności:
anova(mylogit2, mylogit, test="LRT")
Analysis of Deviance Table
Model 1: admit ~ I(gre + gpa) + rank
Model 2: admit ~ gre + gpa + rank
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 395 464.48
2 394 458.52 1 5.9658 0.01459 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Ponownie mamy mocne dowody na to, że współczynniki gre
i gpa
są znacznie różne od siebie.
Funkcja score jest pochodną funkcji logarytmu wiarygodności ( ), gdzie to parametry, a dane (tutaj pokazano ilustrację wielkości pojedynczej cele):
Jest to w zasadzie nachylenie funkcji logarytmu wiarygodności. Dalej, niech będzie macierzą informacji Fishera, która jest ujemnym oczekiwaniem na drugą pochodną funkcji logarytmu prawdopodobieństwa w odniesieniu do . Statystyki testu punktowego to:
Test wyników można również obliczyć za pomocą anova
(statystyki testu wyników nazywane są „Rao”):
anova(mylogit2, mylogit, test="Rao")
Analysis of Deviance Table
Model 1: admit ~ I(gre + gpa) + rank
Model 2: admit ~ gre + gpa + rank
Resid. Df Resid. Dev Df Deviance Rao Pr(>Chi)
1 395 464.48
2 394 458.52 1 5.9658 5.9144 0.01502 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Wniosek jest taki sam jak poprzednio.
Uwaga
Interesująca zależność między różnymi statystykami testu, gdy model jest liniowy, to (Johnston i DiNardo (1997): Metody ekonometryczne ): Wald LR Score.
gre
igpa
? Czy to nie testowanie , a nie ? Dla mnie, aby poprawnie przetestować , musimy zachować, a tymczasem narzucić .gre
gpa
multcomp
pakiety czyni go szczególnie łatwe. Na przykład, spróbuj tego:glht.mod <- glht(mylogit, linfct = c("rank3 - rank4= 0"))
. Ale o wiele łatwiejszym sposobem byłoby podniesienierank3
poziomu odniesienia (użyciemydata$rank <- relevel(mydata$rank, ref="3")
), a następnie użycie zwykłego wyniku regresji. Każdy poziom współczynnika jest porównywany z poziomem odniesienia. Wartość p dlarank4
byłaby pożądanym porównaniem.glht
są dla mnie takie same (około ). Odnośnie drugiego pytania: przetestuj tylko jedną hipotezę liniową, podczas gdy przetestujesz wszystkie 6 porównań parami . Dlatego wartości p należy dostosować do wielu porównań. Oznacza to, że wartości p za pomocą testu Tukeya są na ogół wyższe niż pojedyncze porównanie.linfct = c("rank3 - rank4= 0")
mcp(rank="Tukey")
rank
Nie określiłeś swoich zmiennych, jeśli są one binarne lub coś innego. Myślę, że mówisz o zmiennych binarnych. Istnieją również wielomianowe wersje modelu probit i logit.
Zasadniczo można użyć pełnej trójcy podejść testowych, tj
Test ilorazu wiarygodności
Test LM
Wald-Test
Każdy test wykorzystuje inne statystyki testowe. Standardowym podejściem byłoby wykonanie jednego z trzech testów. Wszystkie trzy można wykorzystać do wykonania wspólnych testów.
Test LR wykorzystuje różnicę logarytmu prawdopodobieństwa modelu ograniczonego i nieograniczonego. Tak więc ograniczonym modelem jest model, w którym określone współczynniki są ustawione na zero. Nieograniczony jest model „normalny”. Zaletą testu Wald jest to, że szacowany jest tylko model bezstratny. Zasadniczo pyta, czy ograniczenie jest prawie spełnione, jeśli zostanie oszacowane przy niezmienionym MLE. W przypadku testu mnożnika Lagrange'a należy oszacować tylko model ograniczony. Ograniczony estymator ML służy do obliczania wyniku nieograniczonego modelu. Wynik ten zwykle nie będzie równy zero, więc ta rozbieżność jest podstawą testu LR. Test LM może być w twoim kontekście również używany do testowania heteroscedastyczności.
źródło
Standardowe podejścia to test Walda, test współczynnika prawdopodobieństwa i test punktowy. Asymptotycznie powinny być takie same. Z mojego doświadczenia wynika, że testy współczynnika prawdopodobieństwa wydają się nieco lepiej wykonywać w symulacjach na próbkach skończonych, ale przypadki, w których ma to znaczenie, byłyby w bardzo ekstremalnych (małych próbkach) scenariuszach, w których wszystkie te testy traktowałbym jedynie jako przybliżone przybliżenie. Jednak w zależności od modelu (liczba zmiennych towarzyszących, obecność efektów interakcji) i danych (wielokolaryzacja, rozkład krańcowy zmiennej zależnej), „cudowne królestwo Asymptotii” można dobrze przybliżyć dzięki zaskakująco małej liczbie obserwacji.
Poniżej znajduje się przykład takiej symulacji w Stacie z wykorzystaniem Walda, współczynnika wiarygodności i testu punktowego w próbie jedynie 150 obserwacji. Nawet w tak małej próbce trzy testy dają dość podobne wartości p, a rozkład próbkowania wartości p, gdy hipoteza zerowa jest prawdziwa, wydaje się być zgodna z rozkładem równomiernym, tak jak powinna (lub przynajmniej z odchyleniami od rozkładu równomiernego nie są większe niż można by się spodziewać ze względu na losowość odziedziczoną w eksperymencie Monte Carlo).
źródło