Regressão linear múltipla em Python

Prévia do material em texto

Ivanildo Batista Follow
Jul 16, 2021 · 20 min read
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.
Q d j iá l i d d
Get unlimited access Open in app
https://ivanildo-batista13.medium.com/?source=post_page-----eb4b6603a3-----------------------------------
https://ivanildo-batista13.medium.com/?source=post_page-----eb4b6603a3-----------------------------------
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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:
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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);
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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_whitewhite 
from statsmodels.stats.diagnostic import linear_harvey_collier, 
linear_reset, spec_whitewhite 
from statsmodels.stats.diagnostic import linear_rainbow 
from statsmodels.graphics.regressionplots import plot_leverage_resid2 
from yellowbrick.regressor import CooksDistance 
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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()
Get unlimited access Open in app
https://www.kaggle.com/dongeorge/beer-consumption-sao-paulo
https://edisciplinas.usp.br/mod/resource/view.php?id=3239600&forceview=1
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
Ú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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
cerveja.isna().sum()
Tipo das variáveis.
cerveja.dtypes
Correlação entre as variáveis.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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) 
i
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
plt.xlabel('Final de Semana') 
plt.ylabel('');
Gráfico das temperaturas média, mínima e máxima em graus Celsius.
cerveja[['TemperaturaMedia (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);
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
Gráfico da precipitação diária.
cerveja['Precipitacao (mm)'].plot(figsize=(15,7)) 
plt.title('Precipitação em mm',size=15);
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
Gráfico do consumo de cerveja.
cerveja['Consumo de cerveja (litros)'].plot(figsize=(15,7), 
color='black'); 
plt.title('Consumo de cerveja',size=15);
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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);
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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");
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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));
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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]); 
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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²)
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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,
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/listsç q p
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()
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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:
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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()
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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')
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
modelo2 = sm.OLS(y,X).fit() 
modelo2.summary(title='Sumário do modelo sem intercepto')
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
p p ç j
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://edisciplinas.usp.br/pluginfile.php/3333945/mod_resource/content/1/Distribuicao%20F%205%25.pdf
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://en.wikipedia.org/wiki/Akaike_information_criterion#Comparison_with_BIC
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
(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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://edisciplinas.usp.br/pluginfile.php/1786946/mod_resource/content/1/Tabelat-student.pdf
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
Diagnósticodo 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.
di ó i d d l ã d íd ã dif
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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)']
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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]);
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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]);
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
Analisando a terceira parte do relatório
A terceira parte trata da análise residual e identificação de multicolinearidade.
Teste Omnibus
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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).
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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 aestimação dos efeitos dos coeficientes.
Consequências 2: parâmetros individualmente significativos, mas não
conjuntamente.
Consequências 3: sensibilidade nas estimações.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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 
i
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
:',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 ?
Get unlimited access Open in app
https://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.OLS.fit_regularized.html
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
q
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 WhiteWhite.
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
22
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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 WhiteWhite : 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)):white 
 print(i,':',j)
Estatística Multiplicador de Lagrange : 39 890819533082166
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBarhttps://medium.com/
https://medium.com/search
https://medium.com/me/lists
q p g
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
Teste de Especificação de WhiteWhite
𝐻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))white 
#spec_white(modelo2.resid, modelo2.model.exog)white
[('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)]
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
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.
Get unlimited access Open in app
https://www.statsmodels.org/stable/stats.html
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists
Get unlimited access Open in app
https://medium.com/
https://medium.com/plans?source=upgrade_membership---nav_full-------------------------------------
https://rsci.app.link/?$canonical_url=https%3A%2F%2Fmedium.com/p/eb4b6603a3&~feature=LiOpenInAppButton&~channel=ShowPostUnderUser&~stage=mobileNavBar
https://medium.com/
https://medium.com/search
https://medium.com/me/lists