Diagnosticando a sua Regressão Linear

Verificando a robustez do seu modelo com statsmodels em Python

Published in

Turing Talks

8 min readJun 5, 2022

Olá, querido leitor, seja muito bem-vindo a mais um Turing Talks. Imagino que você, que não perde nenhum dos nossos textos, já deve estar familiarizado com um dos modelos estatísticos mais clássicos de todos os tempos, a regressão linear. No entanto, caso o conceito esteja um pouco nebuloso, você pode sempre pode reler o nosso texto aqui no Turing Talks sobre o tema.

A proposta de hoje é se aprofundar um pouco em alguns métodos de diagnóstico de uma regressão linear e, assim, verificar a robustez do modelo. Podemos entender robustez de diversas formas, mas aqui, especificamente, damos foco na capacidade do modelo de explicar os dados e de prever de forma segura observações que estejam no seu escopo. É importante pontuar que a regressão linear é muito utilizada não só para prever, mas também para explicar os nossos dados. Nesse sentido, ainda que o nosso modelo produza métricas de negócio satisfatórias, precisamos ficar atentos aos pressupostos e às anormalidades, para que não cheguemos a conclusões erradas sobre os nossos dados.

Começando…

Bom, começando, quando ajustamos um modelo de regressão linear, assumimos:

A nossa variável dependente é aleatória, portanto, haverá um erro e_i na previsão de um valor Y_i.

2. Tais erros são normalmente distribuídos.

3. Os erros das observações são independentes entre si, não havendo nenhuma estrutura temporal, espacial ou outra correlação entre eles.

4. Os erros tem variância constante, ou seja, se estou prevendo preços de casas, a variância do erro numa faixa de 100k-200k BRL não deve ser diferente da variância do erro da previsão de uma casa na faixa de 1M-1.5M BRL.

Lembrando que esses pressupostos permitem a construção de procedimentos de inferência, mas não são necessários para uma convergência numérica estimativas de mínimos quadrados (OLS).

Por último, antes de falarmos em detalhe sobre esses pontos, é importante relembrar que não temos acesso ao erro em si, mas ao chamado resíduo, que é uma aproximação do erro calculada pela diferença do valor previsto e o valor real:

O exemplo

Para ilustrar esse Turing Talks, vamos utilizar um dataset clássico, o California Housing Prices. Utilizaremos apenas uma covariável, a medium income, de forma a tornar mais simples as análises.

Erros não normalmente distribuídos

Quando o erro segue a distribuição normal, podemos garantir que nossos parâmetros beta estimados com mínimos quadrados são iguais a uma estimativa gerada via máxima verossimilhança. Entretanto, a recíproca não é verdadeira, ou seja, a não normalidade dos resíduos não nos permite afirmar que, com certeza, nossos parâmetros estão errados, inclusive, com um número grande de dados, as estimativas podem estar adequadas.

Ainda assim, é importante estar atento a esse ponto, já que uma distribuição muito distante do esperado pode significar que um modelo linear não é suficiente para explicar os dados e, nesse caso, vale a pena inserir transformações nos dados (como log) para tentar resolver o problema. Se isso ainda assim não resolver a questão, talvez seja melhor utilizar outros métodos que podem lidar com a distribuição real.

Para verificar a normalidade dos resíduos, nós podemos utilizar alguns métodos, tanto gráficos como testes estatísticos. Vamos ilustrar com o nosso exemplo. Primeiro, vamos começar ajustando um modelo de regressão linear, com a ajuda da lib statsmodels:

import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy import statsY = data[‘median_house_value’]X = data[‘median_income’]X = sm.add_constant(X)model = sm.OLS(Y, X).fit()

QQ Plot

O gráfico QQ (quantil-quantil) é um gráfico de probabilidade, que é um método gráfico para comparar duas distribuições de probabilidade plotando seus quantis um contra o outro. Nesse caso, a linha é a distribuição teórica (normal) e os pontos são referentes a distribuição encontrada dos resíduos.

res = model.resid
fig = sm.qqplot(res, stats.norm, fit=True, line="45")
plt.show()

Olhando para o gráfico, vemos que nossos resíduos se aproximam de uma distribuição normal nos quantis mais centrais, entretanto, parecem ter uma cauda longa, o que os distanciam do caso ideal.

Histograma

Outra possível estratégia é plotar o histograma dos resíduos e verificar “no olho” se a sua distribuição se aproxima da normal.

plt.hist(res)
plt.plot()

Nesse histograma, fica clara a existência de uma cauda longa na nossa distribuição.

Testes de normalidade

Por último, uma opção mais rigorosa é utilizar um teste estatístico. Você pode utilizar qualquer teste de normalidade que quiser. Aqui vamos utilizar um que convenientemente já está implementado na lib que estamos utilizando.

name = ['Jarque-Bera', 'p-value', 'Skew', 'Kurtosis']
test = sms.jarque_bera(res)
print(dict(zip(name, test))){'Jarque-Bera': 1273.3059054411217, 'p-value': 3.199895134467452e-277, 'Skew': 1.2073436978383232, 'Kurtosis': 5.087039056819328}

O teste rejeita a hipótese nula, o que indica que os dados violam a suposição de normalidade dos resíduos. Ainda assim, vamos continuar com esse modelo, já que o gráfico QQ e o histograma tem um aspecto considerado satisfatório.

Heterocedasticidade

Caso seja observada heterocedasticidade nos resíduos, isto é, caso a sua variância não seja constante, isso pode significar que observações incertas têm muita influência sobre as estimativas. Além disso, os intervalos de previsão estarão, com certeza, errados, já que a eles são calculados utilizando diretamente com a variância dos resíduos.

De novo, podemos verificar se esse é o caso de uma forma gráfica e também com um teste estatístico. Aqui começa a ficar mais claro que a “linha” de valores no topo do gráfico “atrapalha” o nosso modelo.

Plotando valores previstos vs resíduos

Essa é uma forma imediata e visual de verificar se existe alguma tendência nos resíduos em relação ao valor previsto.

Aqui claramente vemos que há uma tendência sistemática, em forma de linha.

Teste estatístico

A lib que utilizamos, novamente, contém um teste estatístico para que verifiquemos se a variância dos nossos resíduos são constantes.

name = ['Lagrange multiplier statistic', 'p-value','f-value', 'f p-value']
test = sms.het_breuschpagan(model.resid, model.model.exog)
print(dict(zip(name, test))){'Lagrange multiplier statistic': 17.044323033060248, 'p-value': 3.651738255523092e-05, 'f-value': 17.13028485393767, 'f p-value': 3.586259560159815e-05}

O teste é significativo, o que significa que os dados violam a suposição de homocedasticidade, ou seja, a heterocedasticidade está presente nos resíduos.

Erros correlacionados

A última coisa que queremos verificar em relação aos resíduos é se eles tem alguma correlação. Um jeito de verificar isso é com um gráfico de auto correlação. Basicamente, interpretamos os resíduos como uma série temporal. Então calculamos a correlação entre as observações em função do intervalo de “tempo” entre elas, ou seja, trasladamos a série em um delta T e calculamos a correlação da série trasladada com a original. Obviamente, quando o lag é zero, a correlação será 1, já que é a correlação da série com ela mesma, isto é, do valor de cada resíduo com ele mesmo.

sm.graphics.tsa.plot_acf(res, lags=10)
plt.plot()

Vemos aqui que os resíduos não parecem estar correlacionados. A faixa azul representa um intervalo de confiança de 95% ao redor do zero.

Algumas outras técnicas úteis para detectar observações estranhas

Agora vamos falar de algumas técnicas mais “avançadas” que podem te ajudar a identificar pontos estranhos na sua regressão linear.

Leverage

De forma simplificada, podemos entender que os valores Y_hat previstos em uma regressão linear são uma combinação linear dos valores Y utilizados para ajustar esse modelo. O leverage é uma medida que nos indica quanto que uma observação influi na previsão do seu próprio valor. Nesse sentido, as observações com alto leverage não conseguem ser explicadas de forma eficiente pelas demais. Um ponto com leverage elevado não necessariamente é um outlier, mas um dado que tem um valor fora do escopo dos dados de treinamento e, por isso, os demais pontos “falham” em descrevê-lo.

influence = model.get_influence()
leverage = influence.hat_matrix_diagplt.scatter(predictions, leverage)
plt.xlabel('predictions')
plt.ylabel('Leverage')
plt.show()

Aqui vemos que as observações com previsão mais alta tem leverage mais alto.

Resíduo Studentizado

Todos os erros estão incluídos em na estimativa da variância dos resíduos, de modo que um grande resíduo contribuirá para uma grande variância, afetando todos os outros resíduos padronizados. Isso pode ser mitigado usando os resíduos studentizado. Um resíduo studentizado é o quociente resultante da divisão de um resíduo por uma estimativa de seu desvio padrão. Quando este valor elevado, isso pode gerar suspeitas em relação ao comportamento daquela observação em específico. De forma geral, é uma forma de “superar” o problema da heterocedasticidade dos resíduos na hora de verificar os pontos de atenção, ou seja, aqueles que tem um resíduo “normalizado” pela variância elevado.

Em código, temos:

student_resid = influence.resid_studentized_externalplt.scatter(predictions, student_resid)
plt.xlabel('predictions')
plt.ylabel('Studentized Residuals')
plt.show()

Distância de Cook

A distância de Cook é uma expressão que avalia o impacto que a retirada de uma observação na estimação dos parâmetros beta da seguinte forma:

Ou equivalentemente:

Sem entrar muitos nos detalhes, intuitivamente, a distância de Cook mede globalmente quanto que os parâmetros se alteram com a retirada de um ponto.

Em código, temos:

(cooks, p) = influence.cooks_distance
plt.scatter(Y, cooks)
plt.xlabel('Y')
plt.ylabel('Cooks Distance')
plt.show()

Aqui verificamos quais pontos estão causando maior mudança nos parâmetros estimados. Poderíamos tentar entender o que difere esses pontos dos demais.

DFBETAS

DFBETAS é uma medida que, assim como a distância de Cook, mede a influência dos pontos nos parâmetros, mas de forma individualizada. Assim, numa regressão múltipla, podemos ver a influência de cada observação em betas específicos, e não como uma medida global.

No nosso exemplo, vamos verificar a influência das observações no beta correspondente a renda média.

dfbetas = influence.dfbetasplt.scatter( data['median_income'], dfbetas[:,1]) 
plt.xlabel('Median income')
plt.ylabel('DFBETAS (beta-median income)')
plt.show()

Conclusão

Nesse texto, você pôde descobrir algumas formas de diagnosticar sua regressão linear e identificar pontos que podem ser outliers ou devem ser melhor estudados.

Esperamos que você tenha gostado e agradecemos por chegar até aqui. Se quiserem conhecer um pouco mais sobre o que fazemos no Turing-USP, não deixem de nos seguir nas redes sociais: Facebook, Instagram, LinkedIn e, claro, acompanhar nossos posts no Medium. Para acompanhar também mais de perto e participar de nossas discussões e eventos, entre no nosso servidor do Discord.

Por hoje é só, até mais!