Regressão Linear com Python

Olá. Hoje vamos revisar como estimar um modelo de regressão linear por MQO no Python. Para isso, vamos usar pandas, scipy e a biblioteca statsmodels. Há algumas outras bibliotecas para estimação de modelos estatísticos em Python mas considero statsmodels a melhor delas pela facilidade e praticidade de uso. Outras bibliotecas como scikit-learn exigem um pouco mais de linhas de programação e, consequentemente, uma base matemática mais consolidada. Aliás, este é um ótimo exercício para quem está estudando modelos lineares: tente programá-lo no Python sem a ajuda dos comandos do pacote, from scratch. Sua compreensão vai melhorar bastante.

Antes de avançarmos para o código, vamos fazer uma breve revisão sobre os modelos lineares por MQO.

Regressão Linear por mínimos quadrados ordinários (MQO)

A regressão linear é um método de análise estatística que nos permite estimar o valor de uma determinada variável resposta (variável dependente) como função de outras variáveis preditoras (variáveis independentes). Essa relação pode ser expressa na seguinte equação:

onde yi representa a variável dependente, xi1, xi2, xip representam as variáveis independentes, os preditores da regressão, b0 representa o intercepto da regressão, ou seja, o valor esperado de y quando todos os preditores tem valor = 0, b1, b2, bp representam os coeficientes estimados da regressão e epsilon representa o termo do erro, ou seja, a diferença entre os valores esperados e os valores observados de yi.

O método de estimação dos mínimos quadrados ordinários - MQO (do inglês ordinary least squares - OLS) é uma técnica de otimização que busca o melhor ajuste do modelo através da minimização dos quadrados do erro da regressão. A regressão linear é um método de análise bastante poderoso. Contudo, ela está ancorada em alguns pressupostos que precisam ser respeitados. Caso contrário, seus resultados tendem a ficar enviesados. São eles:

1. A relação entre a variável resposta e os preditores é linear.

2. A variável resposta possui distribuição normal. Os testes estatísticos desse método são todos baseados no pressuposto da normalidade. Caso tentemos analisar uma variável que não possui distribuição normal com esse método, os resultados não serão bons.

3. O termo do erro é não correlacionado com a variável resposta.

4. O termo do erro possui distribuição normal com média 0.

5. O termo do erro é homoscedástico, ou seja, possui variância constante em toda a sua extensão.

6. Os parâmetros populacionais são constantes, ou seja, valores fixos desconhecidos os quais serão estimados pela equação de regressão.

Estimando MQO no Python

Vamos aos códigos. Primeiro, vamos gerar algumas variáveis aleatórias distribuídas normalmente, uma variável categórica binária e um erro normalmente distribuído para servir de banco de dados para nossa regressão.

%matplotlib inline
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as sm

x1 = sp.random.normal(size=100)
x2 = sp.random.normal(size=100)
x3 = sp.random.normal(size=100)
x4 = sp.random.normal(size=100)
categoria = sp.random.binomial(n=1, p=.5, size=100)
erro = sp.random.normal(size=100)

Em seguida, vamos construir a variável resposta com parâmetros conhecidos para fins de comparação. A equação verdadeira pode ser expressa assim:

y = 3 + 2(x1) + 0.65(x2) + 3.14(x3) - 1.75(x4) + 1.5(dummy) + epsilon

Vamos verificar, no final deste tutorial, o desempenho da análise de regressão em relação aos parâmetros verdadeiros.

y = 3 + 2*x1 + 0.65*x2 + 3.14*x3 - 1.75*x4 + 1.5*categoria + erro

Em seguida, vamos criar um DataFrame com essas variáveis.

dic = {"y":y, "x1":x1, "x2":x2, "x3":x3, "x4":x4, "categoria":categoria}
dados = pd.DataFrame(data=dic)
dados   #Visualiza o DataFrame

Prontinho. Agora, vamos dar uma olhada nas estatísticas descritivas e na distribuição da variável resposta. Lembre-se de que ela deve possuir distribuição normal.

print("Estatísticas descritivas de y:")
dados['y'].describe()

Estatísticas descritivas de y:

count    100.000000
mean       3.709140
std        4.417438
min       -9.847082
25%        0.665959
50%        3.753046
75%        6.709719
max       13.142979
Name: y, dtype: float64

# Elabora um histograma
plt.hist(y, color='red', bins=15)
plt.title('Histograma da variável resposta')
plt.show()

A distribuição da variável parece “bem normal”. Agora vamos estimar o modelo e mostrar os resultados.

reg = sm.ols(formula='y~x1+x2+x3+x4+categoria', data=dados).fit()
print(reg.summary())

Vamos começar pela primeira tabela apresentada. As medidas mais importantes para nós neste momento são o R2 ajustado, a estatística de teste F, o p-valor dessa estatística e, caso queiramos comparar diferentes modelos, o log-likelihood, o Akaike Information Criterion (AIC) e o Bayesian Information Criterion (BIC).

O R2 ajustado nos dá uma medida de porcentagem da variância explicada pelo modelo. Este modelo tem um poder explicativo de 95,7%. Nada mau, hum? Mas não se engane! O R2 não é uma medida de ajuste mas apenas nos diz o quanto de nossos dados esse modelo explica. O restante do que falta explicar é atribuído ao termo do erro, ou seja, todas as outras coisas que tem influência sobre nossa variável resposta mas que não estão no modelo ou não podem de alguma forma ser mensuradas. Em ciências naturais é comum obtermos um R2 próximo de 100%. Já nas ciências sociais, costumamos ficar muito felizes quando conseguimos 30% da variância dos dados explicados. A natureza dos dados estudados é bem diferente. A estatística de teste F e seu p-valor < 0.001 basicamente nos mostram que esse modelo é estatisticamente válido.

Quanto aos coeficientes da regressão, os betas estimados, vejam os resultados. Parecem bem confiáveis, não é mesmo? O modelo pegou muito bem as rlações entre as variáveis e os valores verdadeiros estão todos dentro do intervalo de confiança estimado. Podemos interpretá-los da seguinte forma: quando todos os preditores são iguais a 0, o valor esperado de y é 2,64. Para cada aumento de 1 unidade em x1, y aumenta, em média, 2,04. Para cada aumento de 1 unidade em x2, y aumenta, em média, 0,74. Para cada aumento de 1 unidade em x3, y aumenta, em média, 3,06. Para cada aumento de 1 unidade em x4, y diminui (veja o sinal negativo no coeficiente), em média, 1,83. A variável categórica possui uma análise diferente. Os casos da categoria 1 y têm um valor esperado médio 1,74 maior do que os casos da categoria 0. O p-valor para todas os coeficientes estimados é menor que 0.001 o que mostra que eles são estatisticamente significativos.

Vamos agora plotar alguns gráficos de avaliação. Primeiro, vamos plotar um histograma dos resíduos (epsilon = y - ŷ). Lembre-se de que os resíduos precisam ter distribuição normal com média 0. Em seguida, vamos plotar um gráfico de dispersão (scatterplot) dos dos valores preditos (ŷ) e resíduos. É preciso que os resíduos e a variável resposta sejam não correlacionados.

y_hat = reg.predict()
res = y - y_hat

plt.hist(res, color='orange', bins=15)
plt.title('Histograma dos resíduos da regressão')
plt.show()

plt.scatter(y=res, x=y_hat, color='green', s=50, alpha=.6)
plt.hlines(y=0, xmin=-10, xmax=15, color='orange')
plt.ylabel('$\epsilon = y - \hat{y}$ - Resíduos')
plt.xlabel('$\hat{y}$ ou $E(y)$ - Predito')
plt.show()

Agora, vamos observar apresentar apenas os coeficientes da regressão se quisermos facilitar a análise e traçar a reta de regressão entre as variáveis y e x1.

coefs = pd.DataFrame(reg.params)
coefs.columns = ['Coeficientes']
print(coefs)

           Coeficientes
Intercept      2.964863
x1             2.083678
x2             0.751478
x3             3.229035
x4            -1.880821
categoria      1.378816

plt.scatter(y=dados['y'], x=dados['x1'], color='blue', s=50, alpha=.5)
X_plot = sp.linspace(min(dados['x1']), max(dados['x1']), len(dados['x1']))
plt.plot(X_plot, X_plot*reg.params[1] + reg.params[0], color='r')
plt.ylim(-11,16)
plt.xlim(-2.5,3)
plt.title('Reta de regressão')
plt.ylabel('$y$ - Variável Dependente')
plt.xlabel('$x1$ - Preditor')
plt.show()