Erros padrão pontuais para um ajuste de regressão logística com modelos de estatísticas

Aug 17 2020

Fonte

Uma introdução ao aprendizado estatístico com aplicativos em R , encontrada aqui:https://faculty.marshall.usc.edu/gareth-james/ISL/ISLR%20Seventh%20Printing.pdf

Tarefa

Estou tentando replicar o exemplo de uma regressão logística polinomial no conjunto de dados "Salário" na página 267/8.

Esboço da teoria

De acordo com o livro, uma vez feitas as previsões, os intervalos de confiança podem ser calculados dessa forma. Para um modelo do formulário$$\hat{f}(x_0)=\hat{\beta_0}+\hat{\beta_1}x_0+\hat{\beta_2}x_0^2+\hat{\beta_3}x_0^3+\hat{\beta_4}x_0^4,$$ com um $5\times 5$ matriz de covariância $C$ e vetor $l_0^T=(1, x_0, x_0^2, x_0^3, x_0^4)$, o erro padrão pontual é a raiz quadrada de $\text{Var}[\hat{f}(x_0)]=l_0^TCl_0$. Então, para cada$x_0$ em nosso conjunto de dados, temos um gráfico de previsões $\hat{f}(x_0)$ e um gráfico de intervalos de confiança superior e inferior $\hat{f}(x_0)\pm(2\times \text{Var}[\hat{f}(x_0)])$.

Para uma regressão logística, o mesmo princípio pode ser aplicado, mas a confiança gira em torno da função logit de probabilidade condicional, em oposição às previsões que vêm diretamente da fórmula acima.

Dados e abordagem / código reutilizável

Em primeiro lugar, este é o código para gerar o modelo de regressão logística e traçar os resultados. Esta parte está bem e eu reproduzi com sucesso o que está no livro:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from statsmodels.datasets import get_rdataset
from statsmodels.discrete import discrete_model

from sklearn.preprocessing import PolynomialFeatures
polynomial_feat = PolynomialFeatures(degree=4)

# Get dataset from the R package
data = get_rdataset("Wage", package="ISLR") 
df = data.data.reset_index() 

# Split data into wage (response, y) and age (predictor, X_orig)
y = df.wage
X_orig = df.filter(['age'], axis=1)
# Get the polynomial features from the predictor variable
X = polynomial_feat.fit_transform(X_orig) 

# Set up the test ages for a smooth results plot
X_test = np.linspace(18, 80, 1000)
X_test = X_test[:,np.newaxis] 
X_test_poly = polynomial_feat.fit_transform(X_test) 

# Create a dummy response variable, 1 if wage > 250k and 0 otherwise
y_dummy = pd.DataFrame({'wage': y[:]}) 
y_dummy['wage_split'] = np.where(y_dummy['wage'] > 250, 1, 0) 
y_dummy = y_dummy.drop(['wage'], axis=1)

# Fit a logistic regression model with statsmodels
logit_model = discrete_model.Logit(y_dummy, X).fit() 
# Get predictions, i.e. Pr(Wage > 250 | Age)
y_preds = logit_model.predict(X_test_poly)

# Plot the results
plt.figure(figsize=(8, 8)) 
plt.plot(X_test, y_preds, 'b-') 
plt.ylim(top=0.2) 
plt.xlabel("Age")
plt.ylabel("P(Wage > 250 | Age)")
plt.title("Probability of Earning > 250k with Logistic Regression")

Portanto, agora tento traçar os intervalos de confiança. Não acho que haja um método para fazer isso diretamente nos modelos de estatísticas (corrija-me se estiver errado).

Meu problema

Meu problema aqui é o cálculo dos erros padrão pontuais e os intervalos de confiança. Sabemos que os valores de resposta para o modelo de regressão logística devem ser$y\in [0, 1]$, uma vez que é uma probabilidade condicional.

O problema é que para cada $x_0$, o valor de $$\sqrt{l_0^TCl_0}$$vai ser relativamente grande. Posso demonstrar isso usando o primeiro valor de idade,$x_0=18$:

# Get the covariance matrix from the model class
C = logit_model.normalized_cov_params
x = 18.
L_T = np.array([1, x, x**2, x**3, x**4])

# Compute the pointwise standard error, as outlined above
L_T = np.matrix(L_T)
L = np.transpose(L_T)
C = np.matrix(C)

var_f = np.matmul(np.matmul(L_T, C), L)
var_f = np.asarray(var_f)[0][0]
pointwise_se = np.sqrt(var_f) 
print(pointwise_se)

A saída disso é pointwise_se = 6.14.

No gráfico acima, posso ver que a previsão de $\text{Pr}(\text{Wage} > 250 | x=18)$ está próximo de zero e, a partir do exemplo fornecido no livro, posso ver que o intervalo de confiança em torno desse valor não é amplo e, definitivamente, não fica negativo ou maior que 1.

Se eu fosse obter um intervalo de confiança de um erro padrão pontual de $6.14$, o enredo seria bobo, e não uma réplica disso no livro.

Minha pergunta

O que estou fazendo de errado no cálculo do erro padrão pontual?

Respostas

4 PedroSebe Aug 17 2020 at 23:35

Uma vez que você está fazendo regressão logística e não regressão linear simples, a equação $\hat f(x_0)=\hat\beta_0+\hat\beta_1x_0+\hat\beta_2x_0^2+\hat\beta_3x_0^3+\hat\beta_4x_0^4$não se refere à probabilidade de ganhar> 250 mil, mas ao logit dessa probabilidade. Isso é o mesmo que dizer que a regressão logística é um modelo linear que usa logit como uma função de ligação.

Então, você tem que definir funções para converter entre probabilidades e logits (talvez eles já estejam implementados no Numpy ou algo assim, mas são simples de digitar):

def logit(p):
    return np.log(p/(1-p))

def invlogit(x):
    # inverse function of logit
    return 1/(1+np.exp(-x))

Agora, temos que aplicar o SE pontual que você calculou ao logit das estimativas pontuais e, em seguida, converter de volta para probabilidades:

upper_limit = invlogit(logit(y_pred)+1.96*std_err)
lower_limit = invlogit(logit(y_pred)-1.96*std_err)

Onde std_errestá uma matriz com os erros padrão de$\hat f(x)$que você calculou corretamente. Então, upper_limite lower_limitdará um intervalo em torno da probabilidade estimada.