Jak obliczyć standardowe błędy współczynników regresji logistycznej

18

Korzystam ze scikit-learn Pythona do trenowania i testowania regresji logistycznej.

scikit-learn zwraca współczynniki regresji zmiennych niezależnych, ale nie podaje standardowych błędów współczynników. Potrzebuję tych standardowych błędów, aby obliczyć statystykę Walda dla każdego współczynnika i z kolei porównać te współczynniki ze sobą.

Znalazłem jeden opis, w jaki sposób obliczyć standardowe błędy dla współczynników regresji logistycznej ( tutaj ), ale nieco trudniej jest podążać.

Jeśli zdarzy ci się znać proste, dokładne wyjaśnienie, w jaki sposób obliczyć te standardowe błędy i / lub możesz mi je podać, bardzo bym to docenił! Nie mam na myśli konkretnego kodu (choć prosimy o opublikowanie dowolnego kodu, który może być pomocny), ale raczej algorytmiczne wyjaśnienie związanych z tym kroków.

Gyan Veda
źródło
1
Czy pytasz o kod Pythona, aby uzyskać standardowe błędy, lub o sposób obliczania SE (matematycznie / algorytmicznie), abyś mógł to zrobić sam? Jeśli to pierwsze, to pytanie Q byłoby nie na temat CV (zobacz nasze centrum pomocy ), ale może dotyczyć przepełnienia stosu . Jeśli to drugie, będzie to temat na ten temat (ale możesz nie otrzymać żadnych sugestii kodu). Przeprowadź edycję swojego Q, aby to wyjaśnić. Jeśli jest to ten pierwszy, możemy go dla Ciebie przenieść do SO ( proszę jednak nie przesyłać pocztą ).
gung - Przywróć Monikę
1
Dzięki, Gung. Celowo opublikowałem tutaj, ponieważ oczekuję tego drugiego, ale dokonam edycji, aby wyjaśnić. Wspomniałem, że pracuję w Pythonie ze scikit-learn na wypadek, gdyby ktoś, kto korzysta z tego oprogramowania, mógł udzielić mi wskazówek na ten temat.
Gyan Veda,
Cześć @GyanVeda, Mam teraz ten sam problem, jakie jest twoje ostateczne rozwiązanie, proszę?
zyxue

Odpowiedzi:

12

Czy twoje oprogramowanie daje ci macierz kowariancji parametrów (lub wariancji-kowariancji)? Jeśli tak, standardowe błędy to pierwiastek kwadratowy przekątnej tej macierzy. Prawdopodobnie zechcesz zajrzeć do podręcznika (lub google do notatek z wykładów uniwersyteckich), aby dowiedzieć się, jak uzyskać macierz dla liniowych i uogólnionych modeli liniowych.V.β

użytkownik_ogólny
źródło
1
Nie udało mi się znaleźć niczego online dla przypadku uogólnionego modelu liniowego (może nie znam odpowiednich wyszukiwanych haseł?). Wsparcie?
Kevin H. Lin.
3
Oto taki, który znalazłem po kilku minutach google. Radzę najpierw zrozumieć, w jaki sposób obliczana jest wariancja parametru w podstawowym modelu liniowym. Gdy to otrzymasz, rozszerzenie do GLM jest łatwiejsze. Niemniej jednak wiedza o tym, jak to obliczyć i jak zdobyć go w pakiecie oprogramowania, to nie to samo. www.sagepub.com/upm-data/21121_Chapter_15.pdf
generic_user
18

Standardowe błędy współczynników modelu to pierwiastki kwadratowe przekątnych wpisów macierzy kowariancji. Rozważ następujące:

  • Matryca projektowa:

X = [1x1,1x1,p1x2),1x2),p1xn,1xn,p]xja,jotjotja

(UWAGA: Zakłada to model z przechwytywaniem.)

  • V = [π^1(1-π^1)000π^2)(1-π^2))000π^n(1-π^n)]π^jaja

Macierz kowariancji można zapisać jako:

(XT.V.X)-1

Można to zaimplementować za pomocą następującego kodu:

import numpy as np
from sklearn import linear_model

# Initiate logistic regression object
logit = linear_model.LogisticRegression()

# Fit model. Let X_train = matrix of predictors, y_train = matrix of variable.
# NOTE: Do not include a column for the intercept when fitting the model.
resLogit = logit.fit(X_train, y_train)

# Calculate matrix of predicted class probabilities.
# Check resLogit.classes_ to make sure that sklearn ordered your classes as expected
predProbs = resLogit.predict_proba(X_train)

# Design matrix -- add column of 1's at the beginning of your X_train matrix
X_design = np.hstack([np.ones((X_train.shape[0], 1)), X_train])

# Initiate matrix of 0's, fill diagonal with each predicted observation's variance
V = np.diagflat(np.product(predProbs, axis=1))

# Covariance matrix
# Note that the @-operater does matrix multiplication in Python 3.5+, so if you're running
# Python 3.5+, you can replace the covLogit-line below with the more readable:
# covLogit = np.linalg.inv(X_design.T @ V @ X_design)
covLogit = np.linalg.inv(np.dot(np.dot(X_design.T, V), X_design))
print("Covariance matrix: ", covLogit)

# Standard errors
print("Standard errors: ", np.sqrt(np.diag(covLogit)))

# Wald statistic (coefficient / s.e.) ^ 2
logitParams = np.insert(resLogit.coef_, 0, resLogit.intercept_)
print("Wald statistics: ", (logitParams / np.sqrt(np.diag(covLogit))) ** 2)

To powiedziawszy, statsmodelsprawdopodobnie będzie lepszym pakietem do użycia, jeśli chcesz uzyskać dostęp do wielu „gotowych” narzędzi diagnostycznych.

j_sack
źródło
2
Aby uniknąć problemów z pamięcią i uwzględnić przypadek pojedynczej macierzy, możesz zaktualizować kod w następujący sposób -V = np.product(predProbs, axis=1); covLogit = np.linalg.pinv(np.dot(X_design.T * V), X_design)
steadyfish
6

Jeśli jesteś zainteresowany wnioskowaniem, prawdopodobnie będziesz chciał rzucić okiem na statsmodels . Dostępne są standardowe błędy i wspólne testy statystyczne. Oto przykład regresji logistycznej .

Jseabold
źródło
Dzięki za rekomendację! Zajrzę do statystyk. Szkoda, że ​​scikit-learn nie zapewnia tego rodzaju danych wyjściowych.
Gyan Veda
1
Tak. Po prostu zwykle nie jest celem zestawów narzędzi uczenia maszynowego dostarczanie narzędzi do (częstych) testów hipotez. Jeśli natrafisz na ograniczenia wielkości danych, które nie działają dobrze w statsmodels, ale działają w scikit-learn, chciałbym usłyszeć o nich na github.
jseabold
@jseabold Jednak jeśli chcesz uzyskać pojęcie ad hoc znaczenia funkcji w regresji logistycznej, nie możesz po prostu odczytać wielkości efektu (współczynników) bez myślenia o ich standardowych błędach. Więc nawet jeśli nie wykonujesz częstych testów, a chcesz tylko wskazać wielkość efektu i solidność, brak zmienności wyjściowej jest trudny.
ely