Jak testować równoczesność równości wybranych współczynników w modelu logit lub probit?

Jak przetestować równość wybranych współczynników w modelu logit lub probit? Jakie jest standardowe podejście i jakie jest najnowsze podejście?

hypothesis-testing logit probit

— Qbik
źródło

Test Walda

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: $\beta_{gre}=\beta_{gpa}$ $\beta_{gre}\neq \beta_{gpa}$ $\beta_{gre} - \beta_{gpa} = 0$

W. = \frac{(\hat{β} - β_{0})}{\hat{se} (\hat{β})} \sim N. (0, 1)

$W=\frac{(\hat{\beta}-\beta_{0})}{\widehat{\operatorname{se}}(\hat{\beta})}\sim \mathcal{N}(0,1)$

lub

W^{2} = \frac{(\hat{θ} - θ_{0})^{2}}{Var (\hat{θ})} \sim χ_{1}^{2}

$W^2 = \frac{(\hat{\theta}-\theta_{0})^2}{\operatorname{Var}(\hat{\theta})}\sim \chi_{1}^2$

Nasz tutaj to i . Potrzebujemy więc tylko standardowego błędu . Możemy obliczyć błąd standardowy za pomocą metody Delta : $\widehat{\theta}$ $\beta_{gre} - \beta_{gpa}$ $\theta_{0}=0$ $\beta_{gre} - \beta_{gpa}$

\hat{s e} (β_{g r e} - β_{g p a}) \approx \sqrt{Var (β_{g r e}) + Var (β_{g p a}) - 2 \cdot Cov (β_{g r e}, β_{g p a})}

$\hat{se}(\beta_{gre} - \beta_{gpa})\approx \sqrt{\operatorname{Var}(\beta_{gre}) + \operatorname{Var}(\beta_{gpa}) - 2\cdot \operatorname{Cov}(\beta_{gre},\beta_{gpa})}$

Potrzebujemy także kowariancji i . Macierz wariancji-kowariancji można wyodrębnić za pomocą polecenia po uruchomieniu regresji logistycznej: $\beta_{gre}$ $\beta_{gpa}$ vcov

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 $z$

wald.z <- (gre-gpa)/se
wald.z
[1] -2.413564

Aby uzyskać wartość , wystarczy użyć standardowego rozkładu normalnego: $p$

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 Rużyciu multcomppakietu. 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 .

Test wskaźnika wiarygodności (LRT)

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: $LL_{a}$ $LL_{0}$

D = 2 \cdot (L L_{a} - L L_{0}) \sim χ_{d f 1 - d f 2}^{2}

$D=2\cdot (LL_{a} - LL_{0})\sim \chi_{df1-df2}^{2}$

Statystyka testu prawdopodobieństwa jest zgodna z przy czym stopnie swobody są różnicą w liczbie zmiennych. W naszym przypadku jest to 2. $\chi^{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ć: . $\beta_{gre}=\beta_{gpa}$

\log (\frac{p_{i}}{1 - p_{i}}) = β_{0} + β_{1} \cdot g r e + β_{2} \cdot g p a + β_{3} \cdot r a n k_{2} + β_{4} \cdot r a n k_{3} + β_{5} \cdot r a n k_{4}

$\log\left(\frac{p_{i}}{1-p_{i}}\right)=\beta_{0}+\beta_{1}\cdot \mathrm{gre} + \beta_{2}\cdot \mathrm{gpa}+\beta_{3}\cdot \mathrm{rank_{2}} + \beta_{4}\cdot \mathrm{rank_{3}}+\beta_{5}\cdot \mathrm{rank_{4}}$

\log (\frac{p_{i}}{1 - p_{i}}) = β_{0} + β_{1} \cdot (g r e + g p a) + β_{2} \cdot r a n k_{2} + β_{3} \cdot r a n k_{3} + β_{4} \cdot r a n k_{4}

$\log\left(\frac{p_{i}}{1-p_{i}}\right)=\beta_{0}+\beta_{1}\cdot (\mathrm{gre} + \mathrm{gpa})+\beta_{2}\cdot \mathrm{rank_{2}} + \beta_{3}\cdot \mathrm{rank_{3}}+\beta_{4}\cdot \mathrm{rank_{4}}$

mylogit2 <- glm(admit ~ I(gre + gpa) + rank, data = mydata, family = "binomial")

W naszym przypadku możemy użyć logLikdo 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 grei gpama 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ść : $\chi^{2}_{2}$ $p$

1-pchisq(D, df=1)
[1] 0.01458625

Wartość jest bardzo mała, co wskazuje, że współczynniki są różne. $p$

R ma wbudowany test współczynnika prawdopodobieństwa; możemy użyć anovafunkcji 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 grei gpasą znacznie różne od siebie.

Test punktowy (aka Rao's Score test, czyli test mnożnika Lagrange'a)

Funkcja score jest pochodną funkcji logarytmu wiarygodności ( ), gdzie to parametry, a dane (tutaj pokazano ilustrację wielkości pojedynczej cele): $U(\theta)$ $\text{log} L(\theta|x)$ $\theta$ $x$

U (θ) = \frac{\partial log L. (θ | x)}{\partial θ}

$U(\theta) = \frac{\partial \text{log} L(\theta|x)}{\partial \theta}$

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: $I(\theta)$ $\theta$

S. (θ_{0}) = \frac{U (θ_{0}^{2)})}{ja (θ_{0})} \sim χ_{1}^{2)}

$S(\theta_{0})=\frac{U(\theta_{0}^{2})}{I(\theta_{0})}\sim\chi^{2}_{1}$

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. $\geq$ $\geq$

— COOLSerdash
źródło

Zastanawiam się, dlaczego zredukowany model po prostu wyklucza grei gpa? Czy to nie testowanie , a nie ? Dla mnie, aby poprawnie przetestować , musimy zachować, a tymczasem narzucić .

β_{1} = β_{2} = 0

$\beta_1=\beta_2=0$

β_{1} = β_{2}

$\beta_1=\beta_2$

β_{1} = β_{2}

$\beta_1=\beta_2$ gregpa

β_{gre} = β_{gpa}

$\beta_{\text{gre}}=\beta_{\text{gpa}}$

— Sibbs Gambling

@SibbsGambling Dobry połów! Zaktualizowałem odpowiednio swoją odpowiedź.

— COOLSerdash

Czy ogranicza się to tylko do ciągłych predyktorów, czy też - na przykład - mogę również sprawdzić, czy dwa poziomy zmiennej kategorialnej są znacząco różne? Powiedzmy, czy różnica między rangą 3 a rangą 4 jest znacząca?

— Daniel

@Daniel Tak, to podejście może być również zastosowane do poziomów zmiennej jakościowej. Te multcomppakiety 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 podniesienie rank3poziomu odniesienia (użycie mydata$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 dla rank4byłaby pożądanym porównaniem.

— COOLSerdash

@Daniel Wartości p z danych wyjściowych modelu (zmieniony poziom odniesienia) glhtsą 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.

0.591

$0.591$ linfct = c("rank3 - rank4= 0")mcp(rank="Tukey")rank

— COOLSerdash,

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.

— Jen Bohold
ź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).

clear all
set more off

// data preparation
sysuse nlsw88, clear

gen byte edcat = cond(grade <  12, 1,     ///
                 cond(grade == 12, 2, 3)) ///
                 if grade < .
label define edcat 1 "less than high school" ///
                   2 "high school"           ///
                   3 "more than high school"
label value edcat edcat
label variable edcat "education in categories"

// create cascading dummies, i.e.
// edcat2 compares high school with less than high school
// edcat3 compares more than high school with high school
gen byte edcat2 = (edcat >= 2) if edcat < .
gen byte edcat3 = (edcat >= 3) if edcat < .

keep union edcat2 edcat3 race south
bsample 150 if !missing(union, edcat2, edcat3, race, south)

// constraining edcat2 = edcat3 is equivalent to adding 
// a linear effect (in the log odds) of edcat
constraint define 1 edcat2 = edcat3

// estimate the constrained model
logit union edcat2 edcat3 i.race i.south, constraint(1)

// predict the probabilities
predict pr
gen byte ysim = .
gen w = .

program define sim, rclass
    // create a dependent variable such that the null hypothesis is true
    replace ysim = runiform() < pr

    // estimate the constrained model
    logit ysim edcat2 edcat3 i.race i.south, constraint(1)
    est store constr

    // score test
    tempname b0
    matrix `b0' = e(b)
    logit ysim edcat2 edcat3 i.race i.south, from(`b0') iter(0)
    matrix chi = e(gradient)*e(V)*e(gradient)'
    return scalar p_score = chi2tail(1,chi[1,1])

    // estimate unconstrained model
    logit ysim edcat2 edcat3 i.race i.south 
    est store full

    // Wald test
    test edcat2 = edcat3
    return scalar p_Wald = r(p)

    // likelihood ratio test
    lrtest full constr
    return scalar p_lr = r(p)
end

simulate p_score=r(p_score) p_Wald=r(p_Wald) p_lr=r(p_lr), reps(2000) : sim
simpplot p*, overall reps(20000) scheme(s2color) ylab(,angle(horizontal))

wprowadź opis zdjęcia tutaj

— Maarten Buis
źródło

test punktowy to inna nazwa tego, co @ jen-bohold nazwał testem mnożnika Lagrange'a (LM).

— Maarten Buis

Dobra odpowiedź (+1). Szczególnie podoba mi się wysiłek symulacji. Nie wiedziałem, jak obliczyć wynik testu w Stacie. Dzięki.

— COOLSerdash