<jupyter_start><jupyter_text>Regressão Múltipla II  Imports<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
import seaborn as sns
import scipy.stats as stats

from scipy.stats import f
from statsmodels.graphics.gofplots import ProbPlot
from statsmodels.formula.api import ols<jupyter_output><empty_output><jupyter_text>1. A loja "GT Auto" decidiu aprimorar os treinamentos de vendedores inexperientes de carros, criando uma base de dados dos veículos disponíveis que apresente algumas informações importantes para treinar os novos vendedores. A base de dados foi criada em um documento de texto e contém as seguintes informações:*   Marca/Modelo/Ano do carro: brand/model/year*   Capacidade volumétrica (cc) ou volume de deslocamento do motor: cap_vol.*   Consumo de combustível (km/l): consumo.*   Potência (cv): power.*   Peso (kg): weight.*   0-100 (s), tempo que o carro demora para atingir a velocidade de 100km/h: cemm. *   Número de cilindros: nu_cy.*   Tipo de motor, aspirado (0) ou turbo (1): Etype. Após ter realizado a análise dos dados apresentados, obter o modelo de regressão múltipla, escolha as variáveis adequadas.<jupyter_code>#Observar os nomes das variáveis
cars.columns
#Observar o tipo de dado de cada coluna
cars.dtypes
#Observar os índices
cars.index
# Reduzir a base de dados para usar apenas 
# o motor naturalmente aspirado
dfcars = cars[cars['Etype'] == 0] 
dfcars
#número de instâncias 
na = dfcars.shape[0]

# número de variáveis
pa = dfcars.shape[1]

na, pa
#Regressão Primeira opção
#Regressão com a fórmula import statsmodels.formula.api as smf 
regmul = smf.ols('consumo ~ cap_vol + power + weight', data = dfcars) 

#Realizar o processo de modelagem
res = regmul.fit()

#Resultado detalhado
print(res.summary())
#cálculo do Fcrit (tabela)
import scipy.stats
F=res.fvalue
k=res.df_model # grau do modelo
n=res.nobs # num. amostras
dfn=k
dfd=n-(k+1)
alpha = 0.1 #nível de confiança.
F_critico=scipy.stats.f.ppf(1-alpha, dfn, dfd)
print("F_crit=",F_critico) #tabela F-dist

#como f-statistic > f_crit, pelo menso um doa parâmetros é diferente de 0
#Fazer um teste para cada parâmetro
#cálculo t_crit (tabela t-student)
from scipy.stats import t
alpha = 0.1  # significia = 5% 
df = 11 # graus de liberdade                                      
v = t.ppf(1 - alpha/2, df) 
tt=v
print(f't_crit=: {v}')

#controntar com os valores t da tabela OLS
#Regressão Segunda opção
regmul2 = smf.ols('consumo ~ cap_vol + power + weight + nu_cy', data = dfcars) 
#Realizar o processo de modelagem
res2 = regmul2.fit()
#Resultado detalhado
print(res2.summary())
#regressão Terceira opção
regmul3 = smf.ols('consumo ~ power + weight + nu_cy', data = dfcars) 
#Realizar o processo de modelagem
res3 = regmul3.fit()
#Resultado detalhado
print(res3.summary())
res3.params
#gráficos dispersão e resíduos
#Valores previstos e residuos
y_pred=list(res3.predict())
resi=list(res3.resid)

#crear instancia influence
influence = res3.get_influence()

#obter residuos standardizados
stdresid = list(influence.resid_studentized_internal)
prop=np.divide(resi,stdresid) #e/e*
y=list(dfcars['consumo'])
dftab = pd.DataFrame(list(zip(y,y_pred,resi,stdresid,prop)), columns =["y","y_p","e","e*","e/e*"])
dftab.head(10)
#Gráficos de diagnóstico 2 (y vs. yc)
sns.lmplot(x='y', y='y_p', data=dftab);plt.grid(True)
plt.xlabel('Valores observados: y= Consumo')
plt.ylabel('Valores esperados:  y_pred= Consumo')
plt.axhline(y=0, color='black', linestyle='--', linewidth=1)
plt.show()
#Gráficos de diagnóstico 3 (resíduos padronizados vs. x)
sns.scatterplot(x='y', y='e*', data=dftab);plt.grid(True)
plt.xlabel('Valores observados: y= Consumo')
plt.ylabel('e*= Resíduos padronizados')
plt.axhline(y=0, color='black', linestyle='--', linewidth=1)
plt.show()
#Gráficos de diagnóstico 4
#Verificar a normalidade dos resíduos
#qqpolot vs. normal distribution 
QQ = ProbPlot(influence.resid_studentized_internal)
plot_lm_2 = QQ.qqplot(line='45', alpha=0.5, color='#4C72B0', lw=1)
plot_lm_2.axes[0].set_xlabel('Percentil')
plot_lm_2.axes[0].set_ylabel('Resíduos padronizados')
plt.grid(True)
#Jarque-Bera normality test import statsmodels.api as sm 
JB, JBpv,skw,kurt = sm.stats.stattools.jarque_bera(res3.resid) 
print(JB,JBpv,skw,kurt)
#Analisar amostras influentes
infl = res3.get_influence()
#leverage 
print(infl.hat_matrix_diag)
#Valores dos resíduos (influencias internas)
#residus = res3.resid.as_matrix() #residuals
leviers = infl.hat_matrix_diag #leverage
sigma_err = np.sqrt(res3.scale) #regression standard error
res_stds = stdresid/(sigma_err*np.sqrt(1.0-leviers))
print(res_stds)
#Limiar 
#threshold leverage
seuil_levier = 2*(pa+1)/na
print(seuil_levier)
#identification
atyp_levier = leviers > seuil_levier
print(atyp_levier)
#Quais carros ultrapassam esse limiar
print(dfcars.index[atyp_levier],leviers[atyp_levier])
#dimensões 
print(dfcars.shape) 
#número de instâncias 
na = dfcars.shape[0] 
# número de variáveis
pa = dfcars.shape[1] - 4
na,pa
#valores fornecidos pelo objeto
print(infl.resid_studentized_external)
#verificando com a fórmula
res_studs = res_stds*np.sqrt((na-pa-2)/(na-pa-1-res_stds**2))
print(res_studs)
#Representação dos pontos influentes
sm.graphics.influence_plot(res3)
#limiares dos resíduos (influencias externas)
import scipy
seuil_stud = scipy.stats.t.ppf(0.975,df=na-pa-2)
print(seuil_stud)
#detection - absolute value > threshold
atyp_stud = np.abs(res_studs) > seuil_stud
#which ones?
print(dfcars.index[atyp_stud],res_studs[atyp_stud])
#Observações suspeitas considerando ambos critérios
pbm_infl = np.logical_or(atyp_levier,atyp_stud)
print(dfcars.index[pbm_infl])
#Outros critérios para determinar amostras influentes
print(infl.summary_frame().filter(["hat_diag","student_resid","dffits","cooks_d"]))
# Calcula a correlação entre os atributos numéricos
corr = dfcars.corr()
# Resultado
plt.figure(figsize=(8, 4))
sns.heatmap(corr, cmap='coolwarm', annot=True, fmt=".2f", vmin=-1, vmax=1)
plt.show()
infl.exog[:,1:4]
#correlation matrix 
import scipy 
mc = scipy.corrcoef(infl.exog[:,1:4], rowvar=0) 
print(mc)
#Klein’s rule of thumb 
mc2 = mc**2 
print(mc2)
#VIF criterion 
vif = np.linalg.inv(mc) 
print(vif)<jupyter_output>[[ 34.35864218  -6.87792632 -30.22491181]
 [ -6.87792632   2.62471127   5.49425794]
 [-30.22491181   5.49425794  27.83640092]]<jupyter_text>Extras<jupyter_code>#ANOVA
table = sm.stats.anova_lm(res, typ=2) # 
print(table)
#estimated coefficients 
print(res.params) 
#R2
print(res.rsquared) 
#calculating the F-statistic 
F = res.mse_model/res.mse_resid 
print(F) 
#F provided by the res object 
print(res.fvalue)<jupyter_output>Intercept    13.897284
cap_vol      -0.000544
power        -0.006088
weight       -0.002053
dtype: float64
0.8583177713910948
18.174144629536404
18.174144629536404<jupyter_text>Analisando todas as variáveis do modelo final<jupyter_code>dfex = dfcars.drop(columns=["Etype",       # nome
                      "cemm",       # nome
                      "cap_vol"])     # nome
dfex.head()
# Calcula a correlação entre os atributos numéricos
corr = dfex.corr()
# Resultado
plt.figure(figsize=(8, 4))
sns.heatmap(corr, cmap='coolwarm', annot=True, fmt=".2f", vmin=-1, vmax=1)
plt.show()
print(corr)
#Klein’s rule of thumb 
mc2 = corr**2 
print(mc2)
sns.heatmap(mc2, cmap='coolwarm', annot=True, fmt=".2f", vmin=-1, vmax=1)
#VIF criterion 
vif = np.linalg.inv(mc2) 
print(vif)
sns.heatmap(vif, cmap='coolwarm', annot=True, fmt=".2f", vmin=-4, vmax=4)<jupyter_output><empty_output><jupyter_text>Analisando apenas as variáveis preditoras do modelo final<jupyter_code>dfex1 = dfcars.drop(columns=["Etype",       # nome
                      "cemm",       # nome
                      "consumo",       # nome
                      "cap_vol"])     # nome
dfex1.head()
# Calcula a correlação entre os atributos numéricos
corr1 = dfex1.corr()
# Resultado
plt.figure(figsize=(8, 4))
sns.heatmap(corr1, cmap='coolwarm', annot=True, fmt=".2f", vmin=-1, vmax=1)
plt.show()
print(corr1)
#Klein’s rule of thumb 
mc21 = corr1**2 
print(mc21)
sns.heatmap(mc21, cmap='coolwarm', annot=True, fmt=".2f", vmin=-1, vmax=1)
#VIF criterion 
vif1 = np.linalg.inv(mc21) 
print(vif1)
sns.heatmap(vif1, cmap='coolwarm', annot=True, fmt=".2f", vmin=-4, vmax=4)<jupyter_output><empty_output><jupyter_text>Transformação de variáveis<jupyter_code>#@title Bibliotecas
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.formula.api import ols
import statsmodels.formula.api as smf #adicionada na semana 6
import seaborn as sns
import scipy.stats as stats #adicionada na semana 6
from scipy.stats import f #adicionada na semana 6
from statsmodels.graphics.gofplots import ProbPlot
#@title Conjunto de dados original
dx1 = (0.93, 1.11, 0.93, 1.11, 0.93, 1.11, 0.93, 1.11,
       1.02, 1.02, 1.02, 1.02)
dx2 = (1.00, 1.00, 1.00, 1.00, 1.40, 1.40, 1.40, 1.40, 
       1.18, 1.18, 1.18, 1.18)
dx3= (0.20, 0.20, 0.50, 0.50, 0.20, 0.20, 0.50, 0.50,
      0.31, 0.31, 0.31, 0.31)
dy= (32.95, 38.72, 35.20, 38.72, 32.27, 39.71, 33.67, 
     38.72, 35.20, 33.67, 36.02, 32.27)
# Construir o DataFrame e nomear as colunas
dfo = pd.DataFrame(list(zip(dx1, dx2, dx3, dy)),
               columns =["x1","x2","x3","y"])
dfo.head(3)
#@title Conjunto de dados transformado 
lnx1 = np.log(dx1)
lnx2 = np.log(dx2)
lnx3= np.log(dx3)
lny= np.log(dy)
# Construir o DataFrame e nomear as colunas
dft = pd.DataFrame(list(zip(lnx1, lnx2, lnx3, lny)),
               columns =["lnx1","lnx2","lnx3","lny"])
dft.head(4)
#@title Regressão modelo de múltiplas variáveis
#regressão com a fórmula import statsmodels.formula.api as smf 
regmul = smf.ols('lny ~ lnx1 + lnx2 + lnx3', data = dft) 
#Realizar o processo de modelagem
res = regmul.fit()
#Resultado detalhado
print(res.summary())
#@title Cálculo do Fcrit (tabela)
import scipy.stats
F=res.fvalue
k=res.df_model # grau do modelo
n=res.nobs # num. amostras
dfn=k
dfd=n-(k+1)
alpha = 0.05 #nível de confiança.
F_critico=scipy.stats.f.ppf(1-alpha, dfn, dfd)
print("F_crit=",F_critico) #tabela F-dist
#@title Alguns comandos importantes
#parâmetros estimados
print(res.params) 
#R2
print(res.rsquared) 
#calculando F-statistic 
F = res.mse_model/res.mse_resid 
print(F) 
#F provided by the res object 
print(res.fvalue)
#@title Regressão linear simples após aceitar hipótese nula para b2 e b3
#Apenas para verificar refazer a regresão com uma úica variável
#Devido a que os preditores b2 e b3 foram cancelados
#regressão com a fórmula import statsmodels.formula.api as smf 
regmulin = smf.ols('lny ~ lnx1', data = dft) 
#Realizar o processo de modelagem
reslin = regmulin.fit()
#Resultado detalhado
print(reslin.summary())
#@title Gráficos dispersão e resíduos
#Valores previstos e residuos
y_pred=list(reslin.predict())
resi=reslin.resid
#crear instancia influence
influence = reslin.get_influence()
#obter residuos standardizados
stdresid = list(influence.resid_studentized_internal)
prop=np.divide(resi,stdresid) #e/e*
dftab = pd.DataFrame(list(zip(lnx1, lny, y_pred, resi, stdresid, prop)),
               columns =["lnx1","lny","lnyc","e","e*","e/e*"])
dftab
sns.scatterplot(x='lnx1', y='e*', data=dftab);plt.grid(True)
plt.xlabel('lnx= Velocidade de corte')
plt.ylabel('e*= Resíduos padronizados')
plt.title('Resíduos padronizados vs x1`')
plt.axhline(y=0, color='black', linestyle='--', linewidth=1)
plt.show()
#Verificar a normalidade dos resíduos
#qqpolot vs. normal distribution 
QQ = ProbPlot(influence.resid_studentized_internal)
plot_lm_2 = QQ.qqplot(line='45', alpha=0.5, color='#4C72B0', lw=1)
plot_lm_2.axes[0].set_xlabel('Percentil')
plot_lm_2.axes[0].set_ylabel('Resíduos padronizados')
plt.grid(True)

#@title extras 
#ANOVA
table = sm.stats.anova_lm(res, typ=1) # 
print(table)
# Calcula a correlação entre as variáveis do modelo
corr = dft.corr()
# Resultado
plt.figure(figsize=(8, 4))
sns.heatmap(corr, cmap='coolwarm', annot=True, fmt=".2f", vmin=-1, vmax=1)
plt.show()
import scipy 
mc = scipy.corrcoef(dft, rowvar=0) 
print(mc)
#Klein’s rule of thumb 
mc2 = mc**2 
print(mc2)
#VIF criterion 
vif = np.linalg.inv(mc) 
print(vif)
#Gráfico do VIF, um critério indicador de multicolinearidade é considerar VIF>4
sns.heatmap(vif, cmap='coolwarm', annot=True, fmt=".2f", vmin=-10, vmax=4)<jupyter_output><empty_output><jupyter_text>Modelo logístico<jupyter_code>#@title Definir os dados
lstx1 = (1.80, 1.65, 2.70, 3.67, 1.41, 1.76, 2.10, 2.10, 4.57, 3.59,
         8.33, 2.86, 2.58, 2.90, 3.89, 0.80, 0.60, 1.30, 0.83, 0.57, 
         1.44, 2.08, 1.50, 1.38, 0.94, 1.58, 1.67, 3.00, 2.21)
lstx2 = (2.40, 2.54, 0.84, 1.68, 2.41, 1.93, 1.77, 1.50, 2.43, 5.55,
         5.58, 2.00, 3.68, 1.13, 2.49, 1.37, 1.27, 0.87, 0.97, 0.94, 
         1.00, 0.78, 1.03, 0.82, 1.30, 0.83, 1.05, 1.19, 0.86)
lsty= (1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
# Construir o DataFrame e nomear as colunas
df = pd.DataFrame(list(zip(lstx1, lstx2,lsty)),
               columns =["x1","x2","y"])
x=df[['x1','x2']]
y=df['y']
df.head(3)
#@title Regressão logística
#adicionar uma constante preditora
x = sm.add_constant(x)
# Construir o modelo e ajustar os dados 
model = sm.Logit(y, x).fit()
print(model.summary())
model.params
#@title Razão das chances
# ... Define and fit model
odds_ratios = pd.DataFrame(
    {
        "OR": model.params,
        "Lower CI": model.conf_int()[0],
        "Upper CI": model.conf_int()[1],
    }
)
print(odds_ratios)
odds_ratios = np.exp(odds_ratios)
print(odds_ratios)
#@title Valores previstos vs. observados
# performing predictions on the test datdaset
yhat = model.predict()
prediction = list(map(round, yhat))
 
# comparing original and predicted values of y
print('Valores observados:', list(y))
print('Valores previstos:', prediction)
df2 = pd.DataFrame(list(zip(y, prediction)),
               columns =["x","y"])
#df2.to_csv(index=False)
#@title Acurácia do modelo
from sklearn.metrics import (confusion_matrix,accuracy_score)
 
# confusion matrix
cm = confusion_matrix(y, prediction)
print ("Confusion Matrix : \n", cm)
 
# accuracy score of the model
print('Test accuracy = ', accuracy_score(y, prediction))
#usar a tabela tstudent pata t
from scipy.stats import t
alpha = 0.05  # significia = 5% 
df = 11 # graus de liberdade                                      
v = t.ppf(1 - alpha/2, df) 
tt=v
print(f't_crit=: {v}')<jupyter_output>t_crit=: 2.200985160082949<jupyter_text>Extras<jupyter_code># Calcula a correlação entre as variáveis do modelo
corr = df.corr()
# Resultado
plt.figure(figsize=(8, 4))
sns.heatmap(corr, cmap='coolwarm', annot=True, fmt=".2f", vmin=-1, vmax=1)
plt.show()
import scipy 
mclo = scipy.corrcoef(df, rowvar=0) 
print(mclo)
#Klein’s rule of thumb 
mc2lo = mc**2 
print(mc2lo)
#VIF criterion 
viflo = np.linalg.inv(mc) 
print(viflo)
sns.heatmap(viflo, cmap='coolwarm', annot=True, fmt=".2f", vmin=-4, vmax=4)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13999465669249264, 'fraction_numerical': 0.05283195297889393, 'mean_word_length': 4.752208989627353, 'pattern_counts': {'":': 3, '<': 30, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '430755', 'license': 'MIT', 'n_tokens_mistral': 6820, 'n_tokens_neox': 6212, 'n_words': 1452, 'path': '07_regressao_multipla_II.ipynb', 'repo_name': 'felipeninek/statistical-modeling-inference', 'size': '14972'}