7.2 Regressão Linear Simples
A palavra “regressão” não constitui escolha feliz sob o ponto de vista etimológico, mas já está tão arraigada na literatura estatística que não tentamos substituí-la por outra que exprimisse mais convenientemente suas propriedades essenciais. Foi introduzida por Galton, em conexão com a hereditariedade da estatura. Galton verificou que os filhos de pais cujas estaturas se afastam \(x\) polegadas em relação à estatura média de todos os pais, apresentam um desvio em relação à estatura média de todos os filhos inferior a \(x\) polegadas, isto é, verifica-se o que Galton chamou de “regressão à mediocridade”. Em geral, a idéia ligada à palavra “regressão” nada tem a ver com êste sentido, devendo ser considerada apenas como expressão convencional. (G. Udny Yule and Kendall 1948, 246)
Para maiores detalhes veja (Galton 1889, 95–99) e (Stigler 1997).
7.2.1 Reta
A equação da reta é uma relação matemática utilizada para descrever uma reta no plano cartesiano. Pode ser apresentada de formas distintas, sendo que na Seção 7.2.2 é utilizada a notação da reta reduzida, fazendo \(a' = \beta_{1}\) e \(b' = \beta_{0}\).
| Tipo | Equação |
|---|---|
| Geral | \(ax + by + c = 0\) |
| Segmentária | \(\frac{x}{-c/a} + \frac{y}{-c/b} = 1\) |
| Reduzida | \(y = -\frac{a}{b} x - \frac{c}{b} \; \therefore \; y = a'x + b'\) |
Exemplo 7.15 (Equação reduzida) Considere a reta que passa pelos pontos \(A = (0,-3)\) e \(B = (1.5,0)\).
Uma maneira de descobrir a equação reduzida é substituir os pontos \(A\) e \(B\) em \(y = a'x + b'\):
| Tipo | Equação |
|---|---|
| Ponto \(A\) | \(-3 = a' \times 0 + b' \; \therefore \; b' = -3\) |
| Ponto \(B + b'\) | \(0 = a' \times 1.5 + (-3) \; \therefore \; a' = \frac{3}{1.5} = 2\) |
Assim, a equação reduzida da reta é \(y = 2x - 3\), onde a inclinação é \(a' = 2\) e o intercepto (constante ou coeficiente linear) é \(b' = -3\). Para cada aumento de 1 unidade em \(x\), \(y\) aumenta 2 unidades.
Exemplo 7.16 (Equação segmentária) Do Exemplo 7.15, pode-se obter a equação segmentária da reta a partir da forma reduzida. \[y = 2x - 3 \; \therefore \; 2x - y = 3 \; \therefore \; \frac{2}{3} x - \frac{1}{3} y = \frac{3}{3} \; \therefore \; \frac{x}{3/2} + \frac{y}{3/-1} = 1 \; \therefore \; \frac{x}{1.5} + \frac{y}{-3} = 1.\] Assim, \(-c/a = 1.5\) e \(-c/b = -3\). Note que \(x_{B} = 1.5\) e \(y_{A} = -3\).
Exemplo 7.17 (Equação geral) Do Exemplo 7.15, pode-se obter a equação geral da reta a partir da forma reduzida. \[y = 2x - 3 \; \therefore \; 2x - y - 3 = 0.\] Assim, \(a = 2\), \(b = -1\) e \(c = -3\). Note que \(a' = - \frac{2}{-1} = 2\) e \(b' = - \frac{-3}{-1} = -3\).
7.2.2 Modelo
O modelo de Regressão Linear Simples (RLS) populacional/universal é construído, na abordagem clássica, com todos os \(N\) pares ordenados da população ou universo, i.e., \((x_i,y_i), \; i \in \{1,\ldots,N\}\). Pode ser descrito pela Eq. (7.13).
\[\begin{equation} Y = \beta_0 + \beta_1 x + \varepsilon, \tag{7.13} \end{equation}\] onde \(\varepsilon \sim \mathcal{N}(0,\sigma_{\varepsilon})\).
Note que por simplicidade considera-se \(Y \equiv Y_i\), \(x \equiv x_i\) e \(\varepsilon \equiv \varepsilon_i\). É importante ainda ter em mente que a variável aleatória \(Y\) representa a esperança condicional aos valores observados de \(x\), denotada por \(E(Y|x)\). A estimativa de \(E(Y|x) = Y\) será anotada por \(\widehat{E(Y|x)} = y\).
Exemplo 7.18 Em Python.
import pandas as pd
# Dados (se ainda não estiverem carregados)
dr = pd.read_csv('https://filipezabala.com/data/drinks.txt', sep='\t')
# Calculando a média de 'gar'
m = np.mean(dr['gar'])
print(f"Média de gar: {m:.4f}") # Output: 130.0000Exemplo 7.19 No gráfico a seguir pode-se notar que a linha azul refere-se à média \(\bar{y}=180\). Este valor é constante – ou incondicional – para qualquer valor de \(x\) (no caso a temperatura máxima do dia na escala Celsius). Já a linha vermelha pontilhada indica que espera-se um valor baixo de \(y\) (no caso o número de garrafas vendidas no dia) quanto menor o \(x\), e alto \(y\) quanto maior o \(x\).
plot(dr$temp, dr$gar, las = 1)
abline(h = mean(dr$gar), col = 'blue')
abline(a = -19.11, b = 5.91, col = 'red', lty = 2)
legend(26, 235, legend=c('Incondicional','Condicional'),
col = c('blue', 'red'), lty = 1:2, cex = 0.8, box.lty = 0)
Exemplo 7.20 Em Python.
import matplotlib.pyplot as plt
import numpy as np
# Dados (se ainda não estiverem carregados)
dr = pd.read_csv('https://filipezabala.com/data/drinks.txt', sep='\t')
# Plotando o gráfico de dispersão
plt.scatter(dr['temp'], dr['gar'])
plt.xlabel('Temperatura')
plt.ylabel('Garrafas')
plt.title('Relação entre Temperatura e Garrafas Vendidas')
# Adicionando as linhas
plt.axhline(y=np.mean(dr['gar']), color='blue', linestyle='-', label='Incondicional')
plt.plot(dr['temp'], -19.11 + 5.91 * dr['temp'], color='red', linestyle='--', label='Condicional')
# Adicionando a legenda
plt.legend(loc='upper left')
plt.show()7.2.2.1 Estimativa dos parâmetros
Na maioria dos casos práticos trabalha-se com amostras, sendo necessário estimar os valores de \(\beta_0\) e \(\beta_1\). O método dos mínimos quadrados pode ser utilizado para calcular estas estimativas. O princípio do método é minimizar a soma de quadrado dos erros, i.e., \[\begin{equation} \text{minimizar} \sum_{i=1}^n \varepsilon_{i}^{2}. \tag{7.14} \end{equation}\]
Basicamente utiliza-se \(\varepsilon = Y - \beta_0 - \beta_1 x\) da Eq. (7.13) e deriva-se (parcialmente) em relação a \(\beta_0\) e \(\beta_1\), fazendo cada uma das derivadas parciais igual a zero. Para maiores detalhes recomenda-se (DeGroot and Schervish 2012, 689). As estimativas por mínimos quadrados são enfim dadas por \[\begin{equation} \hat{\beta}_1 = \frac{n \sum{xy} - \sum{x} \sum{y}}{n \sum{x^2} - (\sum{x})^2} \tag{7.15} \end{equation}\] e \[\begin{equation} \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}. \tag{7.16} \end{equation}\]
Exemplo 7.21 Considere novamente os dados do Exemplo 7.3. Pode-se ajustar um modelo linear simples conforme Equações (7.13), (7.15) e (7.16). Note a estimativa na forma \(y = -19.11 + 5.91 x\).
##
## Call:
## lm(formula = gar ~ temp, data = dr)
##
## Coefficients:
## (Intercept) temp
## -19.110 5.915Exemplo 7.22 Em Python.
import pandas as pd
import statsmodels.formula.api as smf
# Dados (se ainda não estiverem carregados)
dr = pd.read_csv('https://filipezabala.com/data/drinks.txt', sep='\t')
# Ajustando o modelo de regressão linear
fit = smf.ols('gar ~ temp', data=dr).fit()
# Imprimindo o resumo do modelo
print(fit.summary())Modelo RPO
Um modelo na forma da Eq. (7.17) é chamado Regressão Pela Origem (RPO) pelo fato de a reta ajustada passar pelo ponto \((0,0)\), a origem do plano cartesiano. Basta fazer \(\beta_0=0\) na Equação (7.13), e para mais detalhes veja (Eisenhauer 2003). \[\begin{equation} Y = \beta_1^{RPO} x + \varepsilon, \tag{7.17} \end{equation}\] onde \(\varepsilon \sim \mathcal{N}(0,\sigma_{\varepsilon})\). A estimativa do parâmetro \(\beta_1^{RPO}\) é dada por \[\begin{equation} \hat{\beta}_1^{RPO} = \frac{\sum{xy} }{\sum{x^2} }. \tag{7.18} \end{equation}\]
Exemplo 7.23 Considere novamente os dados do Exemplo 7.3. Pode-se ajustar um modelo RPO conforme Equações (7.17) e (7.18), resultando em \(y = 5.36 x\). Na função stats::lm (linear model) do R é possível fazer o ajuste utilizando + 0 no lugar de - 1.
##
## Call:
## lm(formula = gar ~ temp - 1, data = dr)
##
## Coefficients:
## temp
## 5.359Exemplo 7.24 Em Python.
import pandas as pd
import statsmodels.formula.api as smf
# Dados (se ainda não estiverem carregados)
dr = pd.read_csv('https://filipezabala.com/data/drinks.txt', sep='\t')
# Ajustando o modelo de regressão linear sem intercepto
fit0 = smf.ols('gar ~ temp - 1', data=dr).fit()
# Imprimindo o resumo do modelo
print(fit0.summary())Exercício 7.7 Considere novamente os dados dos Exemplos 7.3 e 7.21.
a. Calcule à mão a estimativa dos parâmetros do modelo descrito na Equação (7.13) considerando as Equações (7.15) e (7.16).
b. Calcule à mão a estimativa dos parâmetros do modelo descrito na Equação (7.17) considerando a Equação (7.18).
\(\\\)
7.2.3 Diagnóstico inferencial
A análise inferencial consiste em avaliar a qualidade de um modelo em relação à estimação de seus parâmetros (\(\theta\)). Desta forma pode-se pensar no prisma do entendimento do modelo. Este entendimento usualmente passa pela avaliação da aderência do modelo aos dados, atribuindo graus de evidência a hipóteses associadas aos parâmetros modelados ou minimizando algum critério teórico de informação.
Teste para \(\beta_0\)
As hipóteses usuais do teste para \(\beta_0\) são \(H_0: \beta_0 = 0\) vs \(H_1: \beta_0 \ne 0\). Sob \(H_0\) \[\begin{equation} T_0 = \frac{\hat{\beta}_0}{ep(\hat{\beta}_0)} \sim t_{n-2}, \tag{7.19} \end{equation}\] onde \(\hat{\beta}_0\) é dado pela Eq. (7.16) e \[\begin{equation} ep(\hat{\beta}_0) = \sqrt{\hat{\sigma}^2 \left[ \frac{1}{n} + \frac{\bar{x}^2}{S_{xx}} \right] } = \sqrt{\frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2}{n-2} \left[ \frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=1}^n (x_i - \bar{x})^2} \right] }. \tag{7.20} \end{equation}\]
Teste para \(\beta_1\)
O teste para \(\beta_1\) é fundamental na análise de diagnóstico. É com ele que decide-se a respeito da presença ou ausência de relação linear entre \(X\) e \(Y\). As hipóteses usuais do teste para \(\beta_1\) são \(H_0: \beta_1 = 0\) vs \(H_1: \beta_1 \ne 0\). Sob \(H_0\) \[\begin{equation} T_1 = \frac{\hat{\beta}_1}{ep(\hat{\beta}_1)} \sim t_{n-2}, \tag{7.21} \end{equation}\] onde \(\hat{\beta}_1\) é dado pela Eq. (7.15) e \[\begin{equation} ep(\hat{\beta}_1) = \sqrt{\frac{\hat{\sigma}^2}{S_{xx}}} = \sqrt{\frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2 / (n-2)}{\sum_{i=1}^n (x_i - \bar{x})^2} }. \tag{7.22} \end{equation}\]
Teste para \(\beta_1^{RPO}\)
\[\begin{equation} T_1^{RPO} = \frac{\hat{\beta}_1^{RPO}}{ep(\hat{\beta}_1^{RPO})} \sim t_{n-1}, \tag{7.23} \end{equation}\] onde \(\hat{\beta}_1^{RPO}\) é dado pela Eq. (7.18) e \[\begin{equation} ep(\hat{\beta}_1^{RPO}) = \sqrt{\frac{\hat{\sigma}^2}{S_{xx}}} = \sqrt{\frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2 / (n-1)}{\sum_{i=1}^n x_{i}^{2} } }. \tag{7.24} \end{equation}\]
Critérios de informação
AIC - Critério de Informação de Akaike
(Akaike 1974) propõe um critério teórico de informação conforme Eq. (7.25). Quanto menor o valor da medida, melhor (no sentido de mais parcimonioso) é o modelo.
\[\begin{equation} \text{AIC} = -2 l(\hat{\beta}) + 2n_{par}, \tag{7.25} \end{equation}\]
onde \(n_{par}\) indica o número de parâmetros ajustados independentemente no modelo e \(l(\hat{\beta})\) é o logaritmo do valor de máxima verossimilhança do modelo. Para mais detalhes veja ?stats::AIC e ?stats::logLik.
BIC - Critério de Informação Bayesiano
(Schwarz 1978) propõe um critério teórico de informação conforme Eq. (7.26). O autor argumenta36 que tanto o AIC quanto o BIC fornecem uma formulação matemática do princípio da parcimônia na construção de modelos. Porém, para um grande número de observações, os dois critérios podem diferir bastante entre si, sendo mais indicado dar preferência ao BIC.
\[\begin{equation} \text{BIC} = -2 l(\hat{\beta}) + \log(n) n_{par}, \tag{7.26} \end{equation}\]
Exemplo 7.25 Considere o Exemplo 7.21. Pode-se utilizar a função genérica base::summary do R para realizar o diagnóstico dos modelos considerados. O modelo RLS na forma \(y=-19.110+5.915x\) foi atribuído ao objeto fit.
##
## Call:
## lm(formula = gar ~ temp, data = dr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.204 -8.261 3.518 10.796 33.540
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -19.1096 22.9195 -0.834 0.411
## temp 5.9147 0.6736 8.780 1.57e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.24 on 28 degrees of freedom
## Multiple R-squared: 0.7336, Adjusted R-squared: 0.7241
## F-statistic: 77.1 on 1 and 28 DF, p-value: 1.565e-09Exemplo 7.26 Em Python.
import pandas as pd
import statsmodels.formula.api as smf
# Dados (se ainda não estiverem carregados)
dr = pd.read_csv('https://filipezabala.com/data/drinks.txt', sep='\t')
# Ajustando o modelo de regressão linear (se ainda não estiver ajustado)
fit = smf.ols('gar ~ temp', data=dr).fit()
# Imprimindo o resumo do modelo
print(fit.summary())O modelo RPO na forma \(y=5.3590x\) foi atribuído ao objeto fit0.
##
## Call:
## lm(formula = gar ~ temp - 1, data = dr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.810 -10.537 3.503 11.979 30.267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## temp 5.35904 0.09734 55.05 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.14 on 29 degrees of freedom
## Multiple R-squared: 0.9905, Adjusted R-squared: 0.9902
## F-statistic: 3031 on 1 and 29 DF, p-value: < 2.2e-16Exemplo 7.27 Em Python.
import pandas as pd
import statsmodels.formula.api as smf
# Dados (se ainda não estiverem carregados)
dr = pd.read_csv('https://filipezabala.com/data/drinks.txt', sep='\t')
# Ajustando o modelo de regressão linear sem intercepto (se ainda não estiver ajustado)
fit0 = smf.ols('gar ~ temp - 1', data=dr).fit()
# Imprimindo o resumo do modelo
print(fit0.summary())Os critérios de informação de Akaike (AIC) e Schwartz (BIC) podem ser calculados respectivamente via stats::AIC e stats::BIC. Os dois critérios sugerem que o modelo RPO é mais parcimonioso do que o modelo RLS.
## [1] 263.2736## [1] 267.4772## [1] 262.0094## [1] 264.8118Exemplo 7.28 Em Python.
import pandas as pd
import statsmodels.formula.api as smf
# Dados (se ainda não estiverem carregados)
dr = pd.read_csv('https://filipezabala.com/data/drinks.txt', sep='\t')
# Ajustando os modelos (se ainda não estiverem ajustados)
fit = smf.ols('gar ~ temp', data=dr).fit()
fit0 = smf.ols('gar ~ temp - 1', data=dr).fit()
# Calculando AIC e BIC
print(f"AIC (RLS): {fit.aic:.4f}")
print(f"BIC (RLS): {fit.bic:.4f}")
print(f"AIC (RPO): {fit0.aic:.4f}")
print(f"BIC (RPO): {fit0.bic:.4f}")Pode-se ainda realizar uma inspeção visual como complemento da análise inferencial. Para um modelo bem ajustado espera-se que os gráficos de resíduos por valores ajustados (primeira coluna) não indiquem qualquer tipo de tendência ou comportamento modelável. É também esperado para um bom ajuste que o gráfico ‘Normal Q-Q’ apresente os pontos o mais próximo possível da linha pontilhada, que indica normalidade perfeita dos resíduos. A distância de Cook do último gráfico indica possíveis pontos influentes, que podem estar contribuindo de maneira desproporcional ao modelo. Para mais detalhes veja a documentação de ?stats::lm.influence, (R. D. Cook 1977) e (Paula 2013).
Exemplo 7.29 Em Python.
import matplotlib.pyplot as plt
import statsmodels.api as sm
# Ajustando os modelos (se ainda não estiverem ajustados)
# ... (código para ajustar fit e fit0)
# Criando os gráficos de diagnóstico
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
# Gráficos para o modelo 'fit' (RLS)
sm.graphics.plot_regress_exog(fit, 'temp', fig=fig)
# Ajustando os títulos dos subplots
axs[0, 0].set_title("Resíduos vs. Valores Ajustados (RLS)")
axs[0, 1].set_title("Normal Q-Q (RLS)")
axs[1, 0].set_title("Scale-Location (RLS)")
axs[1, 1].set_title("Resíduos vs. Leverage (RLS)")
plt.tight_layout()
plt.show()
# Criando os gráficos de diagnóstico para o modelo 'fit0' (RPO)
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
sm.graphics.plot_regress_exog(fit0, 'temp', fig=fig)
# Ajustando os títulos dos subplots
axs[0, 0].set_title("Resíduos vs. Valores Ajustados (RPO)")
axs[0, 1].set_title("Normal Q-Q (RPO)")
axs[1, 0].set_title("Scale-Location (RPO)")
axs[1, 1].set_title("Resíduos vs. Leverage (RPO)")
plt.tight_layout()
plt.show()Exercício 7.9 Considerando o Exemplo 7.25, verifique:
a. \(\hat{\beta}_0 = -19.110\).
b. \(\hat{\beta}_1 = 5.915\).
c. \(ep(\hat{\beta}_0) = 22.919\).
d. \(ep(\hat{\beta}_1) = 0.674\).
e. \(T_0=-0.83\).
f. \(T_1=8.78\).
g. O p-value associado à hipótese \(H_0:\beta_0=0\) é \(0.41\).
h. O p-value associado à hipótese \(H_0:\beta_1=0\) é \(1.6e-09\).
i. \(\hat{\beta}_1^{RPO} = 5.3590\).
j. \(ep(\hat{\beta}_1^{RPO}) = 0.0973\).
k. \(T_1^{RPO} = 55\).
l. O p-value associado à hipótese \(H_0:\beta_1^{RPO}=0\) é \(<2e-16\).
Exemplo 7.30 Considere novamente os dados do Exemplo 7.1.
a. Faça a regressão de M\(f0 por M\)f0_mhs e vice-versa.
b. O que você percebe após avaliar o summary dos modelos?
Exercício 7.10 Considere novamente os dados do Exemplo 7.3.
a. Realize à mão o teste dos parâmetros do modelo descrito na Equação (7.13) considerando as Equações (7.19), (7.20), (7.21) e (7.22).
b. Realize à mão o teste dos parâmetros do modelo descrito na Equação (7.17) considerando as Equações (7.23) e (7.24).
7.2.3.1 Coeficiente de determinação
Original \[\begin{equation} r^2 = 1 - \frac{SQ_{res}}{SQ_{tot}} = 1 - \frac{\sum_{i=1}^n (y_i-\hat{y}_i)^2}{\sum_{i=1}^n (y_i-\bar{y})^2} \tag{7.27} \end{equation}\]
Ajustado \[\begin{equation} r_{a}^2 = 1 - \frac{SQ_{res}/gl_{res}}{SQ_{tot}/gl_{tot}} = 1 - \frac{(\sum_{i=1}^n (y_i-\hat{y}_i)^2)/(n-p)}{(\sum_{i=1}^n (y_i-\bar{y})^2)/(n-1)} = 1 - (1-r^2) \frac{n-1}{n-p-1} \tag{7.28} \end{equation}\]
7.2.4 Diagnóstico preditivo
A análise preditiva consiste em avaliar a qualidade do modelo em relação à capacidade de predição de novos valores (\(X\)). Desta forma pode-se pensar no prisma da capacidade preditiva do modelo. Tal capacidade pode ser avaliada por diversas medidas, a depender do tipo de variável envolvida.
Variáveis categóricas tipicamente requerem ser agrupadas ou classificadas, e para mensurar a capacidade preditiva utilizam-se medidas tais como acurácia, sensibilidade, sensitividade, etc. Para mais detalhes veja a documentação de ?caret::confusionMatrix e este artigo da Wikipedia.
Modelos para variáveis numéricas tipicamente fazem uso do erro quadrático médio (de previsão), ou EQM, dado pela Eq. (7.29).
\[\begin{equation} \text{EQM} = \frac{\sum_i^{n} (y_i - \hat{y}_i)^2}{n} \tag{7.29} \end{equation}\]
Pode-se escrever uma função para o cálculo do EQM.
Exemplo 7.31 Em Python.
import numpy as np
def eqm(y, yh):
"""
Calcula o Erro Quadrático Médio (EQM).
Args:
y: Valores observados.
yh: Valores previstos.
Returns:
O EQM.
"""
return np.mean((y - yh)**2)Exemplo 7.32 Considere novamente os dados do Exemplo 7.3. Utilizando a ideia de treino-teste, pode-se comparar os modelos quanto aos seus EQMs. Para isso são realizadas \(M=100\) replicações ajustando os modelos RLS e RPO com \(ntrain=20\) observações de treino, confrontadas com as \(ntest=10\) observações de teste restantes. Note que o sorteio da amostra é realizado através da função base::sample associada a sed.seed indexada por i, permitindo a replicação exata dos resultados.
M <- 100 # Número de repetições
ntrain <- 20 # Tamanho da amostra de treino
fit_eqm <- vector(length = M) # Vetor para guardar os M EQMs de fit
fit0_eqm <- vector(length = M) # Vetor para guardar os M EQMs de fit0
for(i in 1:M){
set.seed(i); itrain <- sort(sample(1:nrow(dr), ntrain))
train <- dr[itrain,]
test <- dr[-itrain,]
xtest <- data.frame(temp = dr[-itrain,'temp'])
fit_temp <- lm(gar ~ temp, data = train)
fit0_temp <- lm(gar ~ temp-1, data = train)
fit_eqm[i] <- eqm(test$gar, predict(fit_temp, xtest))
fit0_eqm[i] <- eqm(test$gar, predict(fit0_temp, xtest))
}Exemplo 7.33 Em Python.
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
import statsmodels.formula.api as smf
# Dados (se ainda não estiverem carregados)
dr = pd.read_csv('https://filipezabala.com/data/drinks.txt', sep='\t')
# Parâmetros
M = 100 # Número de repetições
ntrain = 20 # Tamanho da amostra de treino
# Vetores para guardar os EQMs
fit_eqm = np.zeros(M)
fit0_eqm = np.zeros(M)
# Função EQM (definida anteriormente)
# def eqm(y, yh):
# return np.mean((y - yh)**2)
# Loop para calcular os EQMs
for i in range(M):
# Dividindo os dados em treino e teste
train, test = train_test_split(dr, train_size=ntrain, random_state=i)
# Ajustando os modelos
fit_temp = smf.ols('gar ~ temp', data=train).fit()
fit0_temp = smf.ols('gar ~ temp - 1', data=train).fit()
# Calculando as previsões
fit_pred = fit_temp.predict(test)
fit0_pred = fit0_temp.predict(test)
# Calculando os EQMs
fit_eqm[i] = eqm(test['gar'], fit_pred)
fit0_eqm[i] = eqm(test['gar'], fit0_pred)Pode-se comparar numericamente os EQMs dos modelos através da função genérica base::summary e visualmente por graphics::boxplot.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 145.9 301.4 362.1 386.7 455.0 815.2## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 137.8 271.6 349.6 365.4 453.6 779.3
Exemplo 7.34 Em Python.
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
# Dados (se ainda não estiverem carregados)
# ... (código para carregar os dados)
# Parâmetros
# ... (código para definir M e ntrain)
# Vetores para guardar os EQMs
# ... (código para criar fit_eqm e fit0_eqm)
# Função EQM (definida anteriormente)
# ... (código da função eqm)
# Loop para calcular os EQMs
# ... (código do loop for)
# Resumo dos EQMs
print("Resumo de fit_eqm:")
print(pd.Series(fit_eqm).describe())
print("\nResumo de fit0_eqm:")
print(pd.Series(fit0_eqm).describe())
# Boxplot dos EQMs
plt.boxplot([fit_eqm, fit0_eqm], labels=['fit_eqm', 'fit0_eqm'])
plt.title('Comparação dos EQMs')
plt.ylabel('EQM')
plt.show()Nota-se que os valores de EQM associados ao modelo RPO fit0 são ligeriamente inferiores aos do modelo RLS fit. Desta maneira sugere-se o uso da RPO em detrimento da RLS, ainda que os testes de Student e Mann-Whitney (com seus pressupostos básicos devidamente admitidos) sugiram que não haja diferença significativa entre as médias e medianas dos EQMs dos modelos avaliados.
##
## Welch Two Sample t-test
##
## data: fit_eqm and fit0_eqm
## t = 1.1401, df = 197.99, p-value = 0.2556
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -15.56988 58.24303
## sample estimates:
## mean of x mean of y
## 386.7044 365.3678##
## Wilcoxon rank sum test with continuity correction
##
## data: fit_eqm and fit0_eqm
## W = 5445, p-value = 0.2774
## alternative hypothesis: true location shift is not equal to 0Exemplo 7.35 Em Python.
from scipy.stats import ttest_ind, wilcoxon
# ... (código para calcular fit_eqm e fit0_eqm)
# Teste t para amostras independentes
t_stat, p_value = ttest_ind(fit_eqm, fit0_eqm)
print(f"Teste t: t = {t_stat:.4f}, p = {p_value:.4f}")
# Teste de Wilcoxon para amostras independentes
w_stat, p_value = wilcoxon(fit_eqm, fit0_eqm)
print(f"Teste de Wilcoxon: w = {w_stat:.4f}, p = {p_value:.4f}")Referências
Qualitatively both our procedure and Akaike’s give “a mathematical formulation of the principle of parsimony in model building.” Quantitatively, since our procedure differs from Akaike’s only in that the dimension is multiplied by \(\frac{1}{2}\log n\), our procedure leans more than Akaike’s towards lower-dimensional models (when there are 8 or more observations). For large numbers of observations the procedures differ markedly from each other. If the assumptions we made in Section 2 are accepted, Akaike’s criterion cannot be asymptotically optimal. This would contradict any proof of its optimality, but no such proof seems to have been published, and the heuristics of Akaike [1] and of Tong [4] do not seem to lead to any such proof.↩︎