Como calcular os erros padrão dos coeficientes de uma regressão logística

18

Estou usando o scikit-learn do Python para treinar e testar uma regressão logística.

O scikit-learn retorna os coeficientes da regressão das variáveis ​​independentes, mas não fornece os erros padrão dos coeficientes. Eu preciso desses erros padrão para calcular uma estatística de Wald para cada coeficiente e, por sua vez, comparar esses coeficientes entre si.

Eu encontrei uma descrição de como calcular erros padrão para os coeficientes de uma regressão logística ( aqui ), mas é um pouco difícil de seguir.

Se você conhece uma explicação simples e sucinta de como calcular esses erros padrão e / ou pode me fornecer um, eu realmente aprecio isso! Não quero dizer código específico (embora fique à vontade para publicar qualquer código que possa ser útil), mas uma explicação algorítmica das etapas envolvidas.

Gyan Veda
fonte
1
Você está solicitando código Python para obter os erros padrão ou como os SEs são computados (matematicamente / algoritmicamente) para que você possa fazer isso sozinho? Se o primeiro, esse Q seria fora de tópico para o CV (consulte nossa Central de Ajuda ), mas pode estar no tópico sobre estouro de pilha . Neste último caso, o tópico será abordado aqui (mas você pode não receber sugestões de código). Edite seu Q para esclarecer isso. Se for o primeiro, podemos migrar para SO para você ( por favor, não faça postagens cruzadas ).
gung - Restabelece Monica
1
Obrigado, Gung. Eu propositalmente postei aqui porque estou esperando o último, mas vou editar para esclarecer. Mencionei que estou trabalhando em Python com o scikit-learn, caso alguém que use este software possa me dar dicas específicas.
Gyan Veda # 03:
Olá @GyanVeda, estou enfrentando o mesmo problema agora, qual é a sua solução final, por favor?
zyxue

Respostas:

12

Seu software fornece uma matriz de covariância a parâmetros (ou covariância a variações)? Nesse caso, os erros padrão são a raiz quadrada da diagonal dessa matriz. Você provavelmente deseja consultar um livro didático (ou o google para obter notas de aula universitárias) sobre como obter a matriz para modelos lineares e lineares generalizados.Vβ

generic_user
fonte
1
Não consegui encontrar nada on-line para o caso do modelo linear generalizado (talvez eu não conheça os termos de pesquisa corretos?). Socorro?
Kevin H. Lin
3
Aqui está um que eu encontrei depois de alguns minutos pesquisando. Meu conselho é primeiro entender como a variação de parâmetro é calculada em um modelo linear básico. Depois de conseguir isso, a extensão para GLM é mais fácil. Mesmo assim, saber calculá-lo e como obtê-lo em um pacote de software não é a mesma coisa. www.sagepub.com/upm-data/21121_Chapter_15.pdf
generic_user
18

Os erros padrão dos coeficientes do modelo são as raízes quadradas das entradas diagonais da matriz de covariância. Considere o seguinte:

  • Matriz de projeto:

X = [1x1,1...x1,p1x2,1...x2,p1xn,1...xn,p]xEu,jjEu

(NOTA: Isso pressupõe um modelo com interceptação.)

  • V = [π^1(1-π^1)0 0...0 00 0π^2(1-π^2)...0 00 00 0...π^n(1-π^n)]π^EuEu

A matriz de covariância pode ser escrita como:

(XTVX)-1

Isso pode ser implementado com o seguinte código:

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)

Tudo isso dito, statsmodelsprovavelmente será um pacote melhor para usar se você desejar acessar MUITOS diagnósticos "prontos para uso".

j_sack
fonte
2
Para evitar problemas de memória e levar em conta caso de matriz singular, você pode atualizar seu código da seguinte maneira -V = np.product(predProbs, axis=1); covLogit = np.linalg.pinv(np.dot(X_design.T * V), X_design)
stablefish
6

Se você estiver interessado em fazer inferência, provavelmente desejará dar uma olhada nos modelos de estatísticas . Erros padrão e testes estatísticos comuns estão disponíveis. Aqui está um exemplo de regressão logística .

jseabold
fonte
Obrigado pela recomendação! Vou dar uma olhada nos modelos estatísticos. Pena que o scikit-learn não fornece esse tipo de saída.
Gyan Veda
1
Sim. Geralmente, não é o objetivo das caixas de ferramentas do tipo aprendizado de máquina fornecer ferramentas para testes de hipóteses (freqüentistas). Se você encontrar restrições de tamanho de dados que não funcionam bem nos modelos de estatísticas, mas funcionam no scikit-learn, eu estaria interessado em ouvi-las no github.
jseabold
@jseabold No entanto, se você deseja obter alguma noção ad hoc de importância do recurso na regressão logística, não pode simplesmente ler os tamanhos dos efeitos (os coeficientes) sem pensar nos erros padrão. Portanto, mesmo que você não esteja realizando um teste freqüentista e deseje apenas uma indicação do tamanho e da robustez do efeito, a falta de variação sklearn é um desafio.
Ely