Regressão linear múltipla em Python

Um pequeno guia para realizar a análise e o diagnóstico de uma regressão múltipla.

O modelo de regressão linear é uma técnica que obtém equação que melhor ajuste a variável de interesse que está sendo estudada (dependente) e um conjunto de variáveis (independentes ou regressores). Dessa forma, é desenvolvida uma relação matemática entre essas variáveis.

Quando queremos ajustar uma reta entre uma variável independente e outra variável dependente, chamamos o modelo de regressão linear simples. Quando o modelo de regressão que possui mais de um regressor, o modelo é chamado de regressão linear múltipla. No caso do modelo múltiplo, objetivo é mensurar o impacto sobre a média de y de mais de uma covariável. Esse modelo é dado conforme abaixo

Essa estrutura com K regressores também pode ser chamada de função de regressão real ou populacional. Note que o modelo de regressão simples é obtido como um caso particular quando tomamos k = 2.

O primeiro beta é chamado de constate (ou intercepto) e é a média de y quando os valores das variáveis dependentes são conjuntamente iguais a zero. os demais betas são chamados de coeficientes parciais de regressão ou coeficientes de regressão reais (ou Populacionais).

Suposições do modelo

O modelo de regressão múltipla possui algumas suposições e essas suposições garantem que o modelo tenha parâmetros que tragam uma boa interpretação para e que sejam possível captar, o mais próximo da realidade, o impacto das variáveis independente na variável dependente. Essas suposições são:

1) O modelo está corretamente especificado (Ausência de viés);

2) O modelo é linear nos parâmetros;

3) Os valores de X são independentes do termo de erro. A covariância entre os erros e cada variável de X é igual a zero (Ortogonalidade);

4) A esperança dos erros tem média zero : 𝐸(𝑒)=0, para todo 𝑡;

5) A variância dos erros é constante (homocedasticidade):

6) A covariância entre os erros é zero (Ausência de autocorrelação residual):

7) O número de observações deve ser maior que o número de parâmetros;

8) Cada variável X deve ser linearmente independente uma da outra (Ausência de multicolinearidade);

9) Os erros são normalmente distribuídos.

O ideal é que o modelo de regressão siga todas essas suposições, entretanto podem ter situações (e isso é muito comum) que hajam violações e, para cada, violação há um tratamento.

Importando as bibliotecas

bibliotecas e módulos que serão usados nessa pequeno projeto.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
pd.set_option('display.max_columns', 500)
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings("ignore")
from scipy.stats import kurtosis
import statsmodels.stats.api as sms
from statsmodels.compat import lzip
from statsmodels.graphics.gofplots import qqplot
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.stats.diagnostic import het_breuschpagan, het_goldfeldquandt,het_white
from statsmodels.stats.diagnostic import linear_harvey_collier, linear_reset, spec_white
from statsmodels.stats.diagnostic import linear_rainbow
from statsmodels.graphics.regressionplots import plot_leverage_resid2
from yellowbrick.regressor import CooksDistance
from statsmodels.stats.outliers_influence import OLSInfluence, variance_inflation_factor
from sklearn.linear_model import LinearRegression

Dados

Os dados (amostra) foram coletados em São Paulo — Brasil, em uma área universitária, onde acontecem algumas festas com turmas de alunos de 18 a 28 anos (média). O conjunto de dados utilizado para esta atividade possui 7 variáveis, sendo uma variável dependente, com período de um ano. Os dados podem ser obtidos aqui e aqui.

Importando os dados

cerveja = pd.read_csv('Consumo_cerveja.csv')

Análise exploratória dos dados

Primeiras observações da base de dados.

cerveja.head()

Últimas observações.

cerveja.tail()

Dimensão da base de dados: são 365 observações e 7 variáveis.

cerveja.shape
(365, 7)

Não há valores faltantes na base de dados.

cerveja.isna().sum()

Tipo das variáveis.

cerveja.dtypes

Correlação entre as variáveis.

cerveja.corr()

Tabela descritiva das variáveis.

cerveja.describe()

A maioria dos dias não são finais de semana.

plt.figure(figsize=(10,6))
sns.countplot(x='Final de Semana', data=cerveja)
plt.xlabel('Final de Semana')
plt.ylabel('');

Gráfico das temperaturas média, mínima e máxima em graus Celsius.

cerveja[['Temperatura Media (C)', 
         'Temperatura Minima (C)',
         'Temperatura Maxima (C)']].plot(figsize=(15,7));
plt.title('Séries de temperaturas máximas, Médias e Mínimas', size=15);

Gráfico da precipitação diária.

cerveja['Precipitacao (mm)'].plot(figsize=(15,7))
plt.title('Precipitação em mm',size=15);

Gráfico do consumo de cerveja.

cerveja['Consumo de cerveja (litros)'].plot(figsize=(15,7), color='black');
plt.title('Consumo de cerveja',size=15);

Gráfico de correlograma com a correlação de Pearson.

plt.figure(figsize=(20,8))
sns.heatmap(cerveja.corr(), annot = True, cmap= "RdYlGn");
plt.title('Correlação de Pearson',size=15);

Correlograma com a correlação de Spearman.

plt.figure(figsize=(20,8))
plt.title('Correlação de Spearman',size=15)
sns.heatmap(cerveja.corr('spearman'), annot = True, cmap= "RdYlGn");

Boxplots: Não há presença de valores extremos (outliers).

cerveja[['Temperatura Media (C)', 'Temperatura Minima (C)',
       'Temperatura Maxima (C)', 'Precipitacao (mm)',
       'Consumo de cerveja (litros)']].plot.box(figsize=(20,6));

Histograma das variáveis.

cerveja[['Temperatura Media (C)', 'Temperatura Minima (C)',
       'Temperatura Maxima (C)', 'Precipitacao (mm)',
       'Consumo de cerveja (litros)']].hist(figsize=(20,8), bins=50);

Gráfico de dispersão entre as variáveis.

fig,ax = plt.subplots(2,2, figsize=(20,10))
sns.scatterplot(x='Consumo de cerveja (litros)',y='Temperatura Media (C)',data = cerveja,ax=ax[0][0]);
sns.scatterplot(x='Consumo de cerveja (litros)',y='Temperatura Minima (C)',data = cerveja,ax=ax[0][1]);
sns.scatterplot(x='Consumo de cerveja (litros)',y='Temperatura Maxima (C)',data = cerveja,ax=ax[1][0]);
sns.scatterplot(x='Consumo de cerveja (litros)',y='Precipitacao (mm)',data = cerveja,ax=ax[1][1]);

Regressão linear múltipla em Python

Existem três formas de estimar o modelo de regressão múltipla em Python:

• Biblioteca Scikit Learn : que é mais voltada para resolução de problemas de machine learning e por esse motivo é limitada, principalmente para gerar modelos inferenciais;
• Biblioteca Pingouin : é mais sofisticada que a Scikit Learn gerando mais estatísticas e resultados do modelo;
• Biblioteca Statsmodels : é a biblioteca que considero mais completa para gerar o modelo, permitindo a realização de testes para análise e diagnóstico.

Para fins de demonstração será realizado, primeiramente, o procedimento com a biblioteca Scikit Learn e depois com a Statsmodels.

Antes será realizada a separação entre a variável dependente e as variáveis independentes.

#Variáveis independentes
X = cerveja.drop(['Consumo de cerveja (litros)','Data'],axis=1)
#Variável dependentes
y = cerveja['Consumo de cerveja (litros)']

Criando modelo com a Scikit Learn

modelo = LinearRegression()
modelo.fit(X,y)
LinearRegression() #esse é o objeto criado

Coeficiente de Determinação (R²)

modelo.score(X,y)
0.7226497614758338

Intercepto (ou constante) do modelo.

modelo.intercept_
6.444696360572017

Coeficientes (ou parâmetros) do modelo.

modelo.coef_
array([ 0.03079559, -0.01903491,  0.65600076, -0.05746938,  5.18318073])

Não há como gerar outras informações além dessas.

Gerando os modelos com a Pingouin

Irei gerar dois modelos: um com intercepto e outro sem o intercepto. Porém, antes, terei que transformar as variáveis para o formato numpy.

Xp = X.to_numpy()
yp = y.to_numpy()

Modelo de regressão com intercepto:

lm1=pg.linear_regression(
X,y,add_intercept=True,relimp=True).round(4)
lm1

Gerou-se os coeficientes, erros-padrão, teste-t, p-valor, R² e R² ajustado, intervalos de confiança e as contribuições de cada variável na variação da variável dependente. Todas essas estatísticas serão explicadas mais a frente.

Para acessar informações que se encontram oculta, basta colocar o parâmetro as_dataframe como False

lm1 = pg.linear_regression(X, y, add_intercept=True, relimp=True, as_dataframe=False)
print(lm1['df_model']) #graus de liberdade do modelo
print(lm1['df_resid']) #graus de liberdade dos resíduos
#Saídas
5
359

Acessando os valores preditos pelo modelo

#Usando lm1['pred'] é possível acessar as previsões do modelo
x = lm1['pred'].tolist()
Y = y.tolist()

Plotando os valores reais com os valores preditos.

%matplotlib inline
plt.figure(figsize=(20,5))
plt.plot(x, linewidth=2, color='r')
plt.plot(Y, linewidth=0.5,color='b')
plt.title('Valores preditos e os valores reais',size=15)
plt.legend(['Predições','Real'],fontsize=15)
plt.show()

Acessando os resíduos do modelo com intercepto.

%matplotlib inline
plt.figure(figsize=(20,5))
plt.plot(lm1['residuals'].tolist(), linewidth=2, color='g')
plt.legend(['resíduos'],fontsize=15)
plt.show()

Modelos de regressão sem intercepto:

Para o modelo sem o intercepto basta add_intercept=False. Gera-se as mesmas estatísticas, porém sem a constante.

lm2 = pg.linear_regression(X, y, add_intercept=False, relimp=True).round(4)
lm2

Para acessar outras informações basta inserir as_dataframe=False.

lm2 = pg.linear_regression(X, y, add_intercept=False, relimp=True, as_dataframe=False)
print(lm2['df_model'])
print(lm2['df_resid'])
#Saída
5
360

Valores preditos vs Valores reais.

x2 = lm2['pred'].tolist()
Y2 = y.tolist()
plt.figure(figsize=(20,5))
plt.plot(x2, linewidth=2, color='r')
plt.plot(Y2, linewidth=0.5,color='b')
plt.title('Valores preditos e os valores reais',size=15)
plt.legend(['Predições','Real'],fontsize=15)
plt.show()

Resíduos do modelo sem intercepto:

plt.figure(figsize=(20,5))
plt.plot(lm2['residuals'].tolist(), linewidth=2, color='orange')
plt.legend(['resíduos'],fontsize=15)
plt.show()

Mais uma vez ressalto que mesmo sendo melhor, em termos de análise, do que a Scikit-Learn ainda há limitações, como a falta de testes estatísticos nativos dessa biblioteca para o diagnóstico dos modelos.

Gerando os modelos com a Statsmodels

Nessa etapa são gerado dois modelos novamente: um com intercepto (ou constante) e ou sem. A presença ou não do intercepto gera mudanças consideráveis nas estatísticas geradas.

modelo1 = (sm.OLS(y,sm.add_constant(X)).fit())
modelo1.summary(title='Sumário do modelo com intercepto')

modelo2 = sm.OLS(y,X).fit()
modelo2.summary(title='Sumário do modelo sem intercepto')

Agora as próximas etapas serão a análise de cada parte desses relatórios gerados pela Statsmodels.

Analisando a primeira parte dos relatórios

Informações: a variável dependente, o tipo de modelo, a data e hora que ele foi gerado, o número de observações, os graus de liberdade (degrees freedom) dos resíduos e dos modelos e a covariância do modelo é não robusta, ou seja, não está calculada para minimizar ou eliminar variáveis.

𝑅² e 𝑅² ajustado — Coeficiente de determinação

O que é o 𝑅² ?

O coeficiente de determinação é uma proporção que ajuda a entender o quanto as variáveis explicativas explicam a variação da média do consumo de cerveja. No sumário do modelo com intercepto, o coeficiente de determinação foi de 72.3%; já no modelo sem intercepto, o valor foi de 99.1%.

O 𝑅² é um métrica que varia de 0 a 1, se o modelo tiver intercepto; caso contrário usa-se o 𝑅² não centrado (caso do modelo 2).

O 𝑅² varia entre 0 e 1, então quanto maior o 𝑅² melhor é o modelo de regressão, pois teria uma maior a capacidade de explicação.

Uma limitação dessa medida é que com a inserção de regressores ao modelo o 𝑅² tende a aumentar.

𝑅² ajustado

Diferente do 𝑅², o 𝑅² ajustado não sofre a limitação de nunca decair. Caso seja inserido um modelo de regressão uma variável que não seja importante o 𝑅² ajustado irá diminuir.

Uma característica do 𝑅² ajustado é que ele pode ser negativo e por isso ele não pode ser interpretado como uma proporção. Além disso essa medida serve para fazer a comparação entre modelos diferentes.

No nosso exemplo o 𝑅² ajustado caiu um pouco no primeiro modelo, porém manteve-se no segundo modelo.

Fórmula do 𝑅²

SST = Soma dos quadrados totais (sum of squares total).

SSR = Soma dos quadrados da regressão (sum of squares due to regression).

SSE = Soma dos quadrados dos resíduos (sum of squares error).

Fórmula do 𝑅² ajustado

T = Número de observações.

K = Número de parâmetros.

F-statistic e Prob(F-statistic)

Teste de significância conjunta dos parâmetros do modelo.

Hipóteses do teste

Testa se os coeficientes são conjuntamente nulos e para esse teste deve-se rejeitar a hipótese nula e para rejeitar a hipótese nula precisamos que o valor da estatística gerado no modelo esteja fora da região de aceitação da hipótese nula. Esse teste é conhecido como o teste F.

A regra de rejeição é que

O valor foi de 187.1 e observando a tabela da distribuição F de Snedecor a 5% aqui, veremos que esse valor fica muito acima dos valores limites seja qual for o grau de liberdade (os graus de liberdade são definidos como 𝑛−𝑘).

Conclusão rejeitamos a hipótse nula e os parâmetros/coeficientes são conjuntamente significativos.

O Prob(F-statistic) diz a mesma coisa que o F-statistic : probabilidade da hipótese nula ser verdadeira.

Log-Likelihood

1) Usada para fazer comparação entre modelos;

2) Varia de −∞ a +∞;

3) Quanto maior o valor, melhor o modelo.

No caso, entre os dois modelos, o modelo 1 é melhor que o modelo 2.

Critérios de Informação

AIC (Akaike information criterion) e BIC (Bayesian information criterion) são usado para comparar modelos e possuem uma fundamentação estatística e matemática mais rigorosa do que o 𝑅² ajustado.

Existem ainda o AICc (versão alternativa do AIC), o HQIC (Critério de informação Hannan–Quinn), FIC (Critério de informação focada) e o DIC (Critério de informação de desvio). Qualquer um desses pode ser usado, porém o AIC é o mais frequentemente.

O modelo que possui o menor valor do critério de informação é o melhor e no nosso caso, assim como na estatística Log-Likelihood, o melhor modelo é o 1.

Mais informações sobre critério de informação aqui.

Fórmulas dos critérios de informação

Analisando a segunda parte do relatório

Nessa segunda parte constam os parâmetros, o erro padrão de cada coeficiente, a estatística t, seus respectivos p-valores e os intervalos de confiança.

Coeficientes ou parâmetros

Interpretando os coeficientes

A constante é a média da variável dependente, quando todos os valores das outras variáveis forem zero. Os parâmetros das variáveis independentes medem o impacto na variação média da variável dependente.

• Modelo 1 : o aumento de uma unidade na temperatura máxima, aumenta o consumo de cerveja em 0.65 litros;
• Modelo 2 : o aumento de uma unidade na temperatura máxima, aumenta o consumo de cerveja em 0.73 litros.

Para a variável dummy Final de Semana:

• Modelo 1 : Se o dia for em um fim de semana o consumo de cerveja aumenta em 5.1832 litros;
• Modelo 2 : Se o dia for em um fim de semana o consumo de cerveja aumenta em 5.4816 litros.

(Os valores da variável final de semana fazem total sentido, visto que se trata de universitário).

Erro padrão

Erro padrão é uma estimativa do desvio padrão do coeficiente, uma medida da quantidade de variação no coeficiente ao longo de seus pontos de dados.

Teste de significância individual dos parâmetros

Esse resultado faz parte do teste de hipótese do parâmetro. Nesse teste de hipótese testamos se o parâmetro é estatisticamente igual a um determinado valor, ou seja,

O ideal para esse teste é que rejeitemos a hipótese nula. A regra de rejeição de é

Para a variável temperatura máxima, conforme a tabela t-student aqui, os nossos limites da região de aceitação são -2.009 e +2.009 (mais de 50 graus de liberdade e nível de significância de 0.025%, por ser um teste bilateral). Como o resultado excede esses limites da região de aceitação da hipótese nula, rejeitamos a hipótese nula.

Apenas as variáveis Temperatura média e Temperatura Mínima os valores t encontram-se dentro desses intervalos, sendo assim estatisticamente nulos.

P-valor

O que é um p-valor ?

É uma estatística que relata os resultados de um teste de hipótese, essa estatística satisfaz 0≤𝑝(𝑥)≤1. É conhecido como nível de significância exato ou observado ou probabilidade exata de cometer um erro de Tipo I (rejeitar a hipótese nula quando ela é verdadeira).

Quanto maior o valor dessa estatística, maior a evidência a favor da hipótese alternativa do teste. Aqui o teste é da significância individual dos parâmetros.

Para constante, Temperatura Média, Temperatura Mínima, Precipitação e Final de semana o p-valor é menor que 0.05, então rejeitamos a hipótese nula, logo esses parâmetros são estatisticamente diferentes de zero.

Intervalos de confiança

Diagnóstico do modelo

Nem tudo são flores, pois como foi falado no início desse artigo, podem existir violações das suposições do modelo. Por esse motivo é necessário realizar algumas análises por meio de testes estatísticos para avaliar o modelo.

Para os diagnósticos dos modelos serão gerados os resíduos, que são a diferença entre o real e o predito pelo modelo, conforme código abaixo

modelo1.resid
modelo2.resid

Observar o ajuste dos modelos

Predicoes = pd.DataFrame(modelo1.predict(), columns=['Predições 1'])
Predicoes['Predições 2'] = modelo2.predict()
Predicoes['Consumo de cerveja (litros)']=cerveja['Consumo de cerveja (litros)']

Ajuste do modelo 2.

plt.figure()
Predicoes[['Predições 2','Consumo de cerveja (litros)']].plot(figsize=(20,10), color=['r','g']);

Breve análise gráfica dos resíduos

residuos1 = modelo1.resid
fig, ax = plt.subplots(2,2,figsize=(15,6))
residuos1.plot(title="Resíduos do modelo 1", ax=ax[0][0])
sns.distplot(residuos1,ax=ax[0][1])
plot_acf(residuos1,lags=40, ax=ax[1][0])
qqplot(residuos1,line='s', ax=ax[1][1]);

residuos2 = modelo2.resid
fig, ax = plt.subplots(2,2,figsize=(15,6))
residuos2.plot(title="Resíduos do modelo 2", ax=ax[0][0])
sns.distplot(residuos2,ax=ax[0][1])
plot_acf(residuos2,lags=40, ax=ax[1][0])
qqplot(residuos2,line='s', ax=ax[1][1]);

Analisando a terceira parte do relatório

A terceira parte trata da análise residual e identificação de multicolinearidade.

Teste Omnibus

Descreve a normalidade da distribuição de nossos resíduos usando inclinação e curtose como medidas. Um 0 indicaria normalidade perfeita.

Já a Prob(Omnibus) é um teste estatístico que mede a probabilidade de os resíduos serem normalmente distribuídos. Um 1 indicaria uma distribuição perfeitamente normal.

Calculando o teste Omnibus para os modelos.

nome = ['Estatística', 'Probabilidade']
teste = sms.omni_normtest(modelo1.resid)
lzip(nome, teste)
[('Estatística', 39.36215025856847), ('Probabilidade', 2.8354217956923733e-09)] #SAÍDA

Modelo 2.

nome2 = ['Estatística', 'Probabilidade']
teste2 = sms.omni_normtest(modelo2.resid)
lzip(nome2, teste2)
[('Estatística', 20.752436672582746),
 ('Probabilidade', 3.1164892524943026e-05)] #SAÍDA

Assimetria e Curtose

A assimetria é uma medida da falta de simetria de uma distribuição de probabilidade e é obtida a partir do terceiro momento central; em uma distribuição normal a assimetria é igual a zero.

A curtose é uma medida que caracteriza o achatamento de uma distribuição de probabilidade e é obtida a partir do quarto momento central da distribuição; em uma distribuição normal a curtose é igual ao valor 3.

Para os modelos é ideal que a distribuição dos dados tenham uma distribuição normal (assimetria 0 e curtose 3).

Para os modelos a assimetria está próxima de zero, mas a curtose não está próxima de 3. Pode-se concluir que as distribuições dos erros dos modelos não são normais.

Para cálculo dessas estatísticas a parte usar as funções skew() e kurtosis() da biblioteca scipy.

Teste de autocorrelação Durbin-Watson

No teste Durbin-Watson

• se o valor do teste estiver próximo de 4, então há evidência para autocorrelação negativa;
• se estiver próximo de 2, evidência para ausência de autocorrelação;
• e se estiver perto de 0, há evidência para autocorrelação positiva.

Para o primeiro modelo há evidência para ausência de autocorrelação. Para o segundo modelo, a evidência de ausência de autocorrelação é um pouco mais fraca.

Teste de normalidade Jarque-Bera

Mais um teste de normalidade dos resíduos.

O Jarque-Bera testa se a distribuição dos dados é uma distribuição normal (𝐻0) em comparação com uma hipótese alternativa (𝐻1) em que os dados seguem alguma outra distribuição. A estatística do teste é baseada em dois momentos dos dados, a assimetria e a curtose, e possui uma distribuição Chi-quadrado com dois graus de liberdade.

A estatística do teste é

onde alpha_1 estimado é o coeficiente de assimetria e o alpha_2 estimado o coeficiente de curtose.

Regra de rejeição do teste

Para o primeiro modelo o p-valor teve um valor de 0.00155, bem abaixo do nível de significância de 5%, então rejeitamos 𝐻0. Para o segundo modelo, 0.00771.

Multicolinearidade

Regressores com alta correlação entre si. Há um aumento da variância do estimador tornando-o menos eficientes.

Tipos multicolinearidade:

• Exata (perfeita) : duas variáveis, onde uma é função da outra;

Consequências: Não é posssível calcular o estimador MQO.

• Quase exata (imperfeita)

Consequências 1: atrapalha a estimação dos efeitos dos coeficientes.

Consequências 2: parâmetros individualmente significativos, mas não conjuntamente.

Consequências 3: sensibilidade nas estimações.

Formas de detectar multicolinearidade

• Coeficiente de correlação entre as variáveis acima de 90%;
• Determinante de 𝑋′𝑐𝑋𝑐, onde

Se o determinante for igual a zero, há presença de multicolinearidade.

• Fatores de inflação da variância.
• Número Condição : Se o valor for maior que 900, então há evidência para multicolinearidade.

print('Número condição do modelo 1 :',np.linalg.cond(modelo1.model.exog))
Número condição do modelo 1 : 270.83428302907566
print('Número condição do modelo 2 :',np.linalg.cond(modelo2.model.exog))
Número condição do modelo 2 : 85.793834807165

O que fazer em caso de multicolinearidade ?

1) Não fazer nada (não é o ideal);

2) Excluir um dos regressores que está altamente correlacionado com outro (é o ideal a se fazer antes de gerar o modelo);

3) Usar um método de regularização chamado Regressão Ridge (para mais informações acessar aqui).

Outros diagnóstico do modelo

• Heterocedasticidade
• Análise de influência de outliers
• Testes de linearidade

O que é heterocedasticidade ?

Violação da suposição de que a variância dos resíduos são constantes. Ao invés de termos

temos

Formas de identificação

• Gráfica : gerando um gráfico de dispersão/regressão entre os resíduos e a variável dependente;
• Testes estatísticos: Goldfeld-Quandt, Breusch-Pagan e White.

Antes irei fazer uma cópia da base de dados original apenas para auxiliar na análise gráfica.

cerveja2 = cerveja
cerveja2['residuos1'] = modelo1.resid
cerveja2['residuos2'] = modelo2.resid

Forma gráfica

fig, ax = plt.subplots(1,2,figsize=(20,7))
sns.regplot(x='Consumo de cerveja (litros)',y='residuos1',data=cerveja2, ax=ax[0])
sns.regplot(x='Consumo de cerveja (litros)',y='residuos2',data=cerveja2, ax=ax[1]);

Conforme o gráfico acima, a relação entre os resíduos dos modelos e a variável dependente possui uma tendência, logo há uma evidência para presença de heterocedasticidade nos modelos.

Teste estatísticos

Teste Goldfeld-Quandt : testa se a variância entre duas amostras são não constantes. O p-valor da hipótese é de que a variância em uma sub amostra é maior do que na outra sub amostra.

nomes = ['Estatística F','p-valor','Situação da variância']
testeh = het_goldfeldquandt(modelo1.resid, modelo1.model.exog)
lzip(nomes, testeh)
[('Estatística F', 1.0446927666799872),
 ('p-valor', 0.38584666972430953),
 ('Situação da variância', 'increasing')]
nomes = ['Estatística F','p-valor','Situação da variância']
testeh = het_goldfeldquandt(modelo2.resid, modelo2.model.exog)
lzip(nomes, testeh)
[('Estatística F', 0.953250073927243),
 ('p-valor', 0.6248802785957279),
 ('Situação da variância', 'increasing')]

Conforme resultado acima, para ambos os modelos aceita-se a hipótese nula de diferença de variância entre as sub amostras.

Teste Breusch-Pagan : Testa se os erros quadrados tem relação com os regressores .

Estatística F da hipótese é a de que a variância do erro não depende dos regressores. Como o f p-valor é bem pequeno, podemos considerar que a variância do erro depende de 𝑋, para ambos os modelos.

nomes = ['Estatística Multiplicador de Lagrange','p-valor','Estatística F','f p-valor']
for i,j in zip(nomes,het_breuschpagan(modelo1.resid, modelo1.model.exog)):
    print(i,':',j)
Estatística Multiplicador de Lagrange : 31.352601855902844
p-valor : 7.979156260466449e-06
Estatística F : 6.746993460088667
f p-valor : 5.00426367827051e-06
nomes = ['Estatística Multiplicador de Lagrange','p-valor','Estatística F','f p-valor']
for i,j in zip(nomes,het_breuschpagan(modelo2.resid, modelo2.model.exog)):
    print(i,':',j)
Estatística Multiplicador de Lagrange : 172.56571154054927
p-valor : 2.9428003798144765e-36
f-valor : 64.56609853881449
f p-valor : 5.46898062929374e-48

Teste White : testa a mesma hipóteses do Breusch-Pagan e obtivemos, abaixo, os mesmos resultados.

nomes = ['Estatística Multiplicador de Lagrange','p-valor','Estatística F','f p-valor']
for i,j in zip(nomes,het_white(modelo1.resid, modelo1.model.exog)):
    print(i,':',j)
Estatística Multiplicador de Lagrange : 39.890819533082166
p-valor : 0.003382014584983984
Estatística F : 2.2279693886459695
f p-valor : 0.0024897294986551207

Análise de influência

Análise de pontos de alavanca: observações cujo os regressores apresentam padrão atípico.

fig, ax = plt.subplots(1,2,figsize=(20,5))
plot_leverage_resid2(modelo1, ax = ax[0])
plot_leverage_resid2(modelo2, ax = ax[1]);

Distância de Cook: é uma medida da influência de uma observação ou instâncias em uma regressão linear. Muitos pontos altamente influentes podem não ser adequados para regressão linear.

plt.figure(figsize=(20,5))
CooksDistance().fit(X, y).show();

Realizar outras análises de influência

Existem várias estatísticas a serem geradas para análise de influência e podem ser feitas pelos métodos e propriedades abaixo

OLSInfluence(modelo1)
#.summary_table()
# Métodos: get_resid_studentized_external, plot_index, summary_frame, summary_table
# Propriedades : cooks_distance, cov_ratio, dfbeta, dfbetas, dffits, dffits_internal, ess_press,
#hat_diag_factor, hat_matrix_diag, influence, resid_press, resid_std, #resid_studentized, resid_studentized_external,
#resid_studentized_internal,resid_var, sigma2_not_obsi

Exemplo : identificando graficamente observações influentes.

fig, ax = plt.subplots(1,2,figsize=(20,5))
OLSInfluence(modelo1).plot_influence(ax=ax[0])
OLSInfluence(modelo2).plot_influence(ax=ax[1]);

Teste de linearidade e especificação

Teste RESET: O teste RESET (Ramsey Regression Equation Specification Error Test) usa uma regressão aumentada

Hipóteses do teste:

Se 𝛾 foi igual a zero, então o modelo é linear; caso contrário, o modelo não é linear.

linear_reset(modelo1, power = 3)
<class 'statsmodels.stats.contrast.ContrastResults'>
<Wald test (chi2): statistic=[[1.12330099]], p-value=0.5702670650325437, df_denom=2>
linear_reset(modelo2, power = 3)
<class 'statsmodels.stats.contrast.ContrastResults'>
<Wald test (chi2): statistic=[[57.34772848]], p-value=3.52451194177165e-13, df_denom=2>

Para o modelo 1 houve a aceitação da hipótese nula, então 𝛾 é não significativo e o modelo é linear. Para o modelo 2 houve rejeição da hipótese nula, então o modelo é não linear.

Teste de Especificação de White

𝐻0: O modelo é homocedástico e está bem especificado

𝐻1: O modelo não é homocedástico e não está bem especificado

nome = ['Estatística do teste', 'p-valor','Graus de liberdade']
lzip(nome,spec_white(modelo1.resid, modelo1.model.exog))
#spec_white(modelo2.resid, modelo2.model.exog)
[('Estatística do teste', 35.951279826646534),
 ('p-valor', 0.010701864272732915),
 ('Graus de liberdade', 19)]

Para o modelo 1 o p-valor ficou abaixo de 5%, sendo assim, a hipótese nula foi rejeitada. Esse teste só gera resultados caso o modelo tenha intercepto.

Teste Rainbow : A hipótese nula é que o ajuste do modelo usando amostra completa é o mesmo que usar um subconjunto central. A alternativa é que os ajustes são diferentes.

nome = ['Estatística do teste', 'p-valor']
lzip(nome,linear_rainbow(modelo1))
[('Estatística do teste', 1.2892667239720166),
 ('p-valor', 0.04506824253380087)]
nome = ['Estatística do teste', 'p-valor']
lzip(nome,linear_rainbow(modelo2))
[('Estatística do teste', 1.3851336872387137),
 ('p-valor', 0.014830045462060166)]

Para ambos os modelos a hipótese nula é rejeitada, os ajustes dos modelos são diferentes.

Mais testes

Para mais testes e outras análises acessar aqui.

Conclusão

A biblioteca Statsmodels possui ainda uma variedade de testes a serem explorados. Os modelos gerados e analisados passaram em algumas etapas, mas em outras não foram bem sucedidos, demonstrando que há a necessidade de tratamento desses problemas para que o(s) modelo(s) estejam adequado para inferência.

Esse artigo faz parte do meu aprendizado em estatística e econometria na linguagem Python e visa contribuir para divulgação de conhecimento. Críticas construtivas e sugestões de melhoria serão sempre bem vindos e estou a disposição para assimilá-las.