Análise de Regressão

computacional

Author

Affiliation

Universidade Federal de Mato Grosso

Bacharelado em Estatística - UFMT

Modelo de Regressão Linear Simples

Experimento

O artigo “Ecofriendly Dyeing of Silk with Extract of Yerba Mate” (Textile Res. J. 2017: 829–837) descreve um experimento para estudar os efeitos da concentração de corante (mg/L), da temperatura (°C) e do pH na adsorção de corante (mg de corante por grama de tecido). Adsorção de corante é um indicador de cor.

\[ Y = \beta_0 + \beta_1 X + \varepsilon \]

\(Y\): adsorção
\(X\): concentração
\(\varepsilon \sim N(0, \sigma^2)\)

dados <- data.frame(
  concentracao = c(10, 20, 20, 20, 10, 20, 10, 20, 15, 15),
  adsorcao = c(250, 520, 387, 593, 157, 377, 225, 451, 382, 373)
)
dados

   concentracao adsorcao
1            10      250
2            20      520
3            20      387
4            20      593
5            10      157
6            20      377
7            10      225
8            20      451
9            15      382
10           15      373

Estatísticas descritivas

summary(dados)

  concentracao      adsorcao    
 Min.   :10.00   Min.   :157.0  
 1st Qu.:11.25   1st Qu.:280.8  
 Median :17.50   Median :379.5  
 Mean   :16.00   Mean   :371.5  
 3rd Qu.:20.00   3rd Qu.:435.0  
 Max.   :20.00   Max.   :593.0

sd(dados$concentracao)

[1] 4.594683

sd(dados$adsorcao)

[1] 133.3552

Correlação entre as variáveis

plot(dados$concentracao, dados$adsorcao)

ggplot(dados, aes(x = concentracao, y = adsorcao)) +
  geom_point() +
  labs(
    title = "Relação entre concentração e adsorção",
    x = "Concentração (mg/L)",
    y = "Adsorção (mg/g)"
  )

cor(dados$concentracao, dados$adsorcao)

[1] 0.8640818

Teste de hipóteses:

\(H_0: rho = 0\)

cor.test(dados$concentracao, dados$adsorcao)


    Pearson's product-moment correlation

data:  dados$concentracao and dados$adsorcao
t = 4.8554, df = 8, p-value = 0.001263
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5142117 0.9673971
sample estimates:
      cor 
0.8640818

Observe o IC para \(\rho\).

Ajuste do modelo

modelo <- lm(adsorcao ~ concentracao, data = dados)
modelo


Call:
lm(formula = adsorcao ~ concentracao, data = dados)

Coefficients:
 (Intercept)  concentracao  
      -29.76         25.08

Logo, o modelo é dado por:

\[ Y = -29.76 + 25.08 \text{ concentração} \]

anova(modelo)

Analysis of Variance Table

Response: adsorcao
             Df Sum Sq Mean Sq F value   Pr(>F)   
concentracao  1 119501  119501  23.575 0.001263 **
Residuals     8  40551    5069                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimativas, erros-padrões, estatísticas \(t\), valores \(p\), \(\sigma\), \(R^2\), \(R^2_{aj}\), \(F_c\) para regressão (e valor \(p\)):

summary(modelo)


Call:
lm(formula = adsorcao ~ concentracao, data = dados)

Residuals:
   Min     1Q Median     3Q    Max 
-94.82 -53.22  15.28  33.93 121.18 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)   -29.763     85.654  -0.347  0.73719   
concentracao   25.079      5.165   4.855  0.00126 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 71.2 on 8 degrees of freedom
Multiple R-squared:  0.7466,    Adjusted R-squared:  0.715 
F-statistic: 23.58 on 1 and 8 DF,  p-value: 0.001263

Testes de hipóteses:

\[ H_0: \beta_0 = 0 \quad vs \quad H_1: \beta_0 \neq 0 \] \[ H_0: \beta_1 = 0 \quad vs \quad H_1: \beta_1 \neq 0 \] Para o teste \(t\) dos coeficientes, comparamos os valores \(t_c\) com o valor \(t_{tab} = t_{\alpha/2, gl(erro)}\)

qt(0.05/2, df = 8, lower.tail = FALSE) # t tabulado

[1] 2.306004

Para o teste de hipóteses, \(i = 0\) ou \(1\): \[ H_0: \beta_i = \beta_{i,0} \quad vs \quad H_1: \beta_i \neq \beta_{i,0} \] onde \(\beta_{i,0}\) é o valor hipotetizado para o parâmetro (geralmente 0).

A estatística de teste é dada por: \[ t_c = \frac{\hat{\beta}_i - \beta_{i,0}}{se(\hat{\beta}_i)} \]

Para o seguinte exemplo: \[ H_0: \beta_1 = 25 \quad vs \quad H_1: \beta_1 \neq 25 \]

A estatística de teste é: \[ t_c = \frac{\hat{\beta}_1 - 25}{se(\hat{\beta}_1)} \]

t_calc <- (25.079 - 25)/5.165
t_calc

[1] 0.01529526

t_tab <- qt(0.05/2, df = 8, lower.tail = FALSE) # t tabulado
t_tab

[1] 2.306004

t_calc > t_tab

[1] FALSE

# portanto, não rejeitamos H_0!

Conclusões: 1. \(\beta_0 = 0\): Não rejeita \(H_0\) (p = 0.737) 2. \(\beta_1 = 0\): Rejeita \(H_0\) (p = 0.00126)

Coeficientes (betas estimados)

coef(modelo) # coeficientes

 (Intercept) concentracao 
   -29.76316     25.07895

confint(modelo) # IC para os coef.

                  2.5 %    97.5 %
(Intercept)  -227.28135 167.75504
concentracao   13.16815  36.98974

\(\hat{\beta}_0\): adsorção média quando a concentração é zero
\(\hat{\beta}_1\): variação média da adsorção para cada unidade de aumento na concentração

Gráfico com a reta de regressão

dados$ajustado <- fitted(modelo)

ggplot(dados, aes(x = concentracao, y = adsorcao)) +
  geom_point() +
  geom_line(aes(y = ajustado), color = "blue", linewidth = 1) +
  labs(
    title = "Relação entre concentração e adsorção",
    x = "Concentração (mg/L)",
    y = "Adsorção (mg/g)"
  )

plot(dados$concentracao, dados$adsorcao,
     main = "Concentração x Adsorção",
     xlab = "Concentração (mg/L)",
     ylab = "Adsorção (mg/g)")
abline(modelo, col = "blue", lwd = 2)

SQ e QM dos resíduos:

residuos <- residuals(modelo)
SQRes <- sum(residuos^2)
SQRes

[1] 40551.32

GL_Res <- length(dados$adsorcao) - length(coef(modelo))
GL_Res

[1] 8

QMRes <- SQRes / GL_Res
QMRes

[1] 5068.914

sqrt(QMRes)

[1] 71.19631

IC para resposta média

preditos <- predict(modelo, interval = "confidence")
head(preditos)

       fit      lwr      upr
1 221.0263 132.6935 309.3591
2 471.8158 401.3506 542.2810
3 471.8158 401.3506 542.2810
4 471.8158 401.3506 542.2810
5 221.0263 132.6935 309.3591
6 471.8158 401.3506 542.2810

colnames(preditos) <- c("adsorcao_estimada", "lim_inf", "lim_sup") 
head(preditos)

  adsorcao_estimada  lim_inf  lim_sup
1          221.0263 132.6935 309.3591
2          471.8158 401.3506 542.2810
3          471.8158 401.3506 542.2810
4          471.8158 401.3506 542.2810
5          221.0263 132.6935 309.3591
6          471.8158 401.3506 542.2810

dados <- cbind(dados, preditos)
dados

   concentracao adsorcao ajustado adsorcao_estimada  lim_inf  lim_sup
1            10      250 221.0263          221.0263 132.6935 309.3591
2            20      520 471.8158          471.8158 401.3506 542.2810
3            20      387 471.8158          471.8158 401.3506 542.2810
4            20      593 471.8158          471.8158 401.3506 542.2810
5            10      157 221.0263          221.0263 132.6935 309.3591
6            20      377 471.8158          471.8158 401.3506 542.2810
7            10      225 221.0263          221.0263 132.6935 309.3591
8            20      451 471.8158          471.8158 401.3506 542.2810
9            15      382 346.4211          346.4211 293.1544 399.6877
10           15      373 346.4211          346.4211 293.1544 399.6877

Ordernar o banco de dados pela concentração:

dados <- dados[order(dados$concentracao),]
dados

   concentracao adsorcao ajustado adsorcao_estimada  lim_inf  lim_sup
1            10      250 221.0263          221.0263 132.6935 309.3591
5            10      157 221.0263          221.0263 132.6935 309.3591
7            10      225 221.0263          221.0263 132.6935 309.3591
9            15      382 346.4211          346.4211 293.1544 399.6877
10           15      373 346.4211          346.4211 293.1544 399.6877
2            20      520 471.8158          471.8158 401.3506 542.2810
3            20      387 471.8158          471.8158 401.3506 542.2810
4            20      593 471.8158          471.8158 401.3506 542.2810
6            20      377 471.8158          471.8158 401.3506 542.2810
8            20      451 471.8158          471.8158 401.3506 542.2810

Gráfico com IC para a resposta média

grid <- data.frame(concentracao = seq(min(dados$concentracao),
                                      max(dados$concentracao),
                                      length.out = 200))
pred  <- as.data.frame(predict(modelo, newdata = grid, interval = "confidence"))
pred  <- cbind(grid, pred)


ggplot() +
  geom_ribbon(data = pred,
              aes(x = concentracao, ymin = lwr, ymax = upr),
              alpha = 0.18) +
  geom_line(data = pred,
            aes(x = concentracao, y = fit),
            color = "blue",
            linewidth = 1) +
  geom_point(data = dados,
             aes(x = concentracao, y = adsorcao),
             size = 3, alpha = 0.8) +
  labs(
    title = "Concentração x Adsorção",
    x = "Concentração (mg/L)",
    y = "Adsorção (mg/g)"
  )

plot(dados$concentracao, dados$adsorcao,
     main = "Concentração x Adsorção",
     xlab = "Concentração (mg/L)",
     ylab = "Adsorção (mg/g)",
     xlim = c(10,20),
     ylim = c(0, 600))

# valores estimados
lines(dados$concentracao, dados$adsorcao_estimada,
      col = "blue", lty = 1, lwd = 2) 
# IC
lines(dados$concentracao, dados$lim_inf,
      col = "blue", lty = 2, lwd = 2)
lines(dados$concentracao, dados$lim_sup,
      col = "blue", lty = 2, lwd = 2)

IC para novas observações (IP)

pred_new  <- as.data.frame(predict(modelo, newdata = grid, interval = "prediction"))
colnames(pred_new) <- c("fit", "lwr_new", "upr_new")
pred  <- cbind(pred, pred_new[,2:3])


ggplot() +
  geom_ribbon(data = pred,
              aes(x = concentracao, ymin = lwr, ymax = upr),
              alpha = 0.18) +
  geom_ribbon(data = pred,
              aes(x = concentracao, ymin = lwr_new, ymax = upr_new),
              alpha = 0.1) +
  geom_line(data = pred,
            aes(x = concentracao, y = fit),
            color = "blue",
            linewidth = 1) +
  geom_point(data = dados,
             aes(x = concentracao, y = adsorcao),
             size = 3, alpha = 0.8) +
  labs(
    title = "Concentração x Adsorção",
    x = "Concentração (mg/L)",
    y = "Adsorção (mg/g)"
  )

new_data <- data.frame(concentracao = seq(0,30, by = 1))
preditos_new <- predict(modelo, newdata = new_data,
                        interval = "prediction")
preditos_new <- cbind(new_data, preditos_new)
head(preditos_new)

  concentracao        fit       lwr      upr
1            0 -29.763158 -286.6059 227.0795
2            1  -4.684211 -252.8178 243.4494
3            2  20.394737 -219.3054 260.0949
4            3  45.473684 -186.0986 277.0460
5            4  70.552632 -153.2309 294.3362
6            5  95.631579 -120.7388 312.0019

Gráfico com IP para novas observações

# dispersão
plot(dados$concentracao, dados$adsorcao,
     main = "Concentração x Adsorção",
     xlab = "Concentração (mg/L)",
     ylab = "Adsorção (mg/g)",
     xlim = c(0,30),
     ylim = c(-300, 1000))
# valores estimados
lines(preditos_new$concentracao, preditos_new$fit,
      col = "blue", lty = 1, lwd = 2) 
# IP
lines(preditos_new$concentracao, preditos_new$lwr,
      col = "brown", lty = 2, lwd = 2)
lines(preditos_new$concentracao, preditos_new$upr,
      col = "brown", lty = 2, lwd = 2)

Gráfico com IC e IP

# dispersão
plot(dados$concentracao, dados$adsorcao,
     main = "Concentração x Adsorção",
     xlab = "Concentração (mg/L)",
     ylab = "Adsorção (mg/g)",
     xlim = c(0,30),
     ylim = c(-300, 1000))

# valores estimados
lines(preditos_new$concentracao, preditos_new$fit,
      col = "blue", lty = 1, lwd = 2) 

# IC
lines(dados$concentracao, dados$lim_inf,
      col = "blue", lty = 2, lwd = 2)
lines(dados$concentracao, dados$lim_sup,
      col = "blue", lty = 2, lwd = 2)  

# IP
lines(preditos_new$concentracao, preditos_new$lwr,
      col = "brown", lty = 2, lwd = 2)
lines(preditos_new$concentracao, preditos_new$upr,
      col = "brown", lty = 2, lwd = 2)

Regressão pela origem

modelo2 <- lm(adsorcao ~ -1 + concentracao, data = dados)
modelo2


Call:
lm(formula = adsorcao ~ -1 + concentracao, data = dados)

Coefficients:
concentracao  
       23.35

anova(modelo2)

Analysis of Variance Table

Response: adsorcao
             Df  Sum Sq Mean Sq F value    Pr(>F)    
concentracao  1 1499012 1499012  327.75 2.182e-08 ***
Residuals     9   41163    4574                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo2)


Call:
lm(formula = adsorcao ~ -1 + concentracao, data = dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-89.945 -61.341   4.027  29.541 126.055 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
concentracao    23.35       1.29    18.1 2.18e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 67.63 on 9 degrees of freedom
Multiple R-squared:  0.9733,    Adjusted R-squared:  0.9703 
F-statistic: 327.7 on 1 and 9 DF,  p-value: 2.182e-08

Comparação entre modelo completo e modelo reduzido (sem intercepto)

anova(modelo2, modelo)

Analysis of Variance Table

Model 1: adsorcao ~ -1 + concentracao
Model 2: adsorcao ~ concentracao
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1      9 41163                           
2      8 40551  1    612.04 0.1207 0.7372

Se não há diferença entre os modelos, qual deve ser escolhido?

Comparação entre modelo completo e modelo reduzido (sem intercepto)

anova(modelo2, modelo)

Analysis of Variance Table

Model 1: adsorcao ~ -1 + concentracao
Model 2: adsorcao ~ concentracao
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1      9 41163                           
2      8 40551  1    612.04 0.1207 0.7372

Se não há diferença entre os modelos, qual deve ser escolhido?

O modelo mais simples: com menos parâmetros.

Comparação de R2 e R2ajustado

summary(modelo)


Call:
lm(formula = adsorcao ~ concentracao, data = dados)

Residuals:
   Min     1Q Median     3Q    Max 
-94.82 -53.22  15.28  33.93 121.18 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)   -29.763     85.654  -0.347  0.73719   
concentracao   25.079      5.165   4.855  0.00126 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 71.2 on 8 degrees of freedom
Multiple R-squared:  0.7466,    Adjusted R-squared:  0.715 
F-statistic: 23.58 on 1 and 8 DF,  p-value: 0.001263

summary(modelo2)


Call:
lm(formula = adsorcao ~ -1 + concentracao, data = dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-89.945 -61.341   4.027  29.541 126.055 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
concentracao    23.35       1.29    18.1 2.18e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 67.63 on 9 degrees of freedom
Multiple R-squared:  0.9733,    Adjusted R-squared:  0.9703 
F-statistic: 327.7 on 1 and 9 DF,  p-value: 2.182e-08

Modelo de Regressão Linear Múltipla

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 + \varepsilon \]

\(Y\): adsorção
\(X_1\): concentração
\(X_2\): temperatura
\(X_3\): pH
\(\varepsilon \sim N(0, \sigma^2)\)

dados <- data.frame(
  conc = c(10, 20, 20, 20, 10, 20, 10, 20, 15, 15),
  temp = c(70, 70, 80, 90, 70, 70, 90, 90, 80, 80),
  ph = c(3.0, 3.0, 3.0, 3.0, 4.0, 4.0, 4.0, 4.0, 3.5, 3.5),
  adsorcao = c(250, 520, 387, 593, 157, 377, 225, 451, 382, 373)
)
dados

   conc temp  ph adsorcao
1    10   70 3.0      250
2    20   70 3.0      520
3    20   80 3.0      387
4    20   90 3.0      593
5    10   70 4.0      157
6    20   70 4.0      377
7    10   90 4.0      225
8    20   90 4.0      451
9    15   80 3.5      382
10   15   80 3.5      373

Ajuste do modelo (regressão múltipla)

# regressão múltipla
modelo3 <- lm(adsorcao ~ conc + temp + ph, data = dados)
modelo3


Call:
lm(formula = adsorcao ~ conc + temp + ph, data = dados)

Coefficients:
(Intercept)         conc         temp           ph  
     53.759       21.548        3.657      -90.272

anova(modelo3)

Analysis of Variance Table

Response: adsorcao
          Df Sum Sq Mean Sq F value   Pr(>F)   
conc       1 119501  119501 34.6477 0.001066 **
temp       1   5157    5157  1.4952 0.267261   
ph         1  14700   14700  4.2621 0.084530 . 
Residuals  6  20694    3449                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coef(modelo3) # coeficientes

(Intercept)        conc        temp          ph 
  53.758696   21.548478    3.657174  -90.271739

confint(modelo3) # IC para os coef.

                  2.5 %     97.5 %
(Intercept) -507.581522 615.098914
conc          10.518537  32.578419
temp          -1.988525   9.302873
ph          -197.265543  16.722065

Estimativas, erros-padrões, estatísticas \(t\), valores \(p\), \(\sigma\), \(R^2\), \(R^2_{aj}\), \(F_c\) para regressão (e valor \(p\)):

summary(modelo3)


Call:
lm(formula = adsorcao ~ conc + temp + ph, data = dados)

Residuals:
     Min       1Q   Median       3Q      Max 
-119.487   -6.477   -2.215   26.141   50.085 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   53.759    229.408   0.234  0.82252   
conc          21.548      4.508   4.780  0.00306 **
temp           3.657      2.307   1.585  0.16405   
ph           -90.272     43.726  -2.064  0.08453 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 58.73 on 6 degrees of freedom
Multiple R-squared:  0.8707,    Adjusted R-squared:  0.8061 
F-statistic: 13.47 on 3 and 6 DF,  p-value: 0.004493

IC para resposta média

preditos_mod3 <- predict(modelo3, interval = "confidence")
head(preditos_mod3)

       fit      lwr      upr
1 254.4304 147.2270 361.6338
2 469.9152 380.2725 559.5579
3 506.4870 431.8766 581.0974
4 543.0587 445.7325 640.3849
5 164.1587  66.8325 261.4849
6 379.6435 266.7295 492.5575

colnames(preditos_mod3) <- c("adsorcao_estimada", "lim_inf", "lim_sup") 
head(preditos_mod3)

  adsorcao_estimada  lim_inf  lim_sup
1          254.4304 147.2270 361.6338
2          469.9152 380.2725 559.5579
3          506.4870 431.8766 581.0974
4          543.0587 445.7325 640.3849
5          164.1587  66.8325 261.4849
6          379.6435 266.7295 492.5575

dados <- cbind(dados, preditos_mod3)
dados

   conc temp  ph adsorcao adsorcao_estimada  lim_inf  lim_sup
1    10   70 3.0      250          254.4304 147.2270 361.6338
2    20   70 3.0      520          469.9152 380.2725 559.5579
3    20   80 3.0      387          506.4870 431.8766 581.0974
4    20   90 3.0      593          543.0587 445.7325 640.3849
5    10   70 4.0      157          164.1587  66.8325 261.4849
6    20   70 4.0      377          379.6435 266.7295 492.5575
7    10   90 4.0      225          237.3022 129.4725 345.1319
8    20   90 4.0      451          452.7870 354.3144 551.2596
9    15   80 3.5      382          353.6087 306.2310 400.9863
10   15   80 3.5      373          353.6087 306.2310 400.9863

IC para novas observações (IP)

new_data <- data.frame(conc = c(12, 12 , 25, 25),
                       temp = c(50, 80, 50, 80),
                       ph = c(6.0, 6.0, 6.0, 6.0))
preditos_new_mod3 <- predict(modelo3, newdata = new_data,
                        interval = "prediction")
preditos_new_mod3 <- cbind(new_data, preditos_new_mod3)
preditos_new_mod3

  conc temp ph       fit          lwr      upr
1   12   50  6 -46.43130 -406.3418415 313.4792
2   12   80  6  63.28391 -235.1306176 361.6984
3   25   50  6 233.69891 -177.5734979 644.9713
4   25   80  6 343.41413   -0.8326895 687.6610

Análise de resíduos

Example 3.1 The Delivery Time Data (Montgomery, Peck e Vining, 2021, p. 76)

A soft drink bottler is analyzing the vending machine service routes in his distribution system. He is interested in predicting the amount of time required by the route driver to service the vending machines in an outlet. This service activity includes stocking the machine with beverage products and minor maintenance or housekeeping. The industrial engineer responsible for the study has suggested that the two most important variables affecting the delivery time (y) are the number of cases of product stocked (x1) and the distance walked by the route driver (x2). The engineer has collected 25 observations on delivery time:

ID <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
        11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 
        21, 22, 23, 24, 25)

y <- c(16.68, 11.50, 12.03, 14.88, 13.75, 18.11, 8.00, 17.83, 79.24, 21.50,
       40.33, 21.00, 13.50, 19.75, 24.00, 29.00, 15.35, 19.00, 9.50, 35.10,
       17.90, 52.32, 18.75, 19.83, 10.75)

x1 <- c(7, 3, 3, 4, 6, 7, 2, 7, 30, 5,
        16, 10, 4, 6, 9, 10, 6, 7, 3, 17,
        10, 26, 9, 8, 4)

x2 <- c(560, 220, 340, 80, 150, 330, 110, 210, 1460, 605,
        688, 215, 255, 462, 448, 776, 200, 132, 36, 770,
        140, 810, 450, 635, 150)

dados <- data.frame(ID, y, x1, x2)

# Ajuste do modelo linear múltiplo
mod <- lm(y ~ x1 + x2, data = dados)

# estimativas dos parâmetros/coeficientes
summary(mod)


Call:
lm(formula = y ~ x1 + x2, data = dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.7880 -0.6629  0.4364  1.1566  7.4197 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.341231   1.096730   2.135 0.044170 *  
x1          1.615907   0.170735   9.464 3.25e-09 ***
x2          0.014385   0.003613   3.981 0.000631 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.259 on 22 degrees of freedom
Multiple R-squared:  0.9596,    Adjusted R-squared:  0.9559 
F-statistic: 261.2 on 2 and 22 DF,  p-value: 4.687e-16

# Resíduos

# resíduos brutos
res_1 <- residuals(mod)

# resíduos padronizados
res_2 <- res_1/summary(mod)$sigma

# resíduos studentizados
# res_3 <- rstandard(mod, type = "sd.1")
res_3 <- rstandard(mod)

# resíduos PRESS
res_4 <- rstandard(mod, type = "pred")

# resíduos R-Student
res_5 <- rstudent(mod)

residuos <- data.frame(
  "brutos" = res_1,
  "padronizados" = res_2,
  "studentizado" = res_3,
  "PRESS" = res_4,
  "RStudent" = res_5
)

residuos

       brutos padronizados studentizado       PRESS    RStudent
1  -5.0280843  -1.54260631  -1.62767993 -5.59796734 -1.69562881
2   1.1463854   0.35170879   0.36484267  1.23360321  0.35753764
3  -0.0497937  -0.01527661  -0.01609165 -0.05524867 -0.01572177
4   4.9243539   1.51078203   1.57972040  5.38401290  1.63916491
5  -0.4443983  -0.13634053  -0.14176094 -0.48043610 -0.13856493
6  -0.2895743  -0.08884082  -0.09080847 -0.30254339 -0.08873728
7   0.8446235   0.25912883   0.27042496  0.91986749  0.26464769
8   1.1566049   0.35484408   0.36672118  1.23532680  0.35938983
9   7.4197062   2.27635117   3.21376278 14.78889824  4.31078012
10  2.3764129   0.72907878   0.81325432  2.95682585  0.80677584
11  2.2374930   0.68645843   0.71807970  2.44837821  0.70993906
12 -0.5930409  -0.18194377  -0.19325733 -0.66908638 -0.18897451
13  1.0270093   0.31508443   0.32517935  1.09387183  0.31846924
14  1.0675359   0.32751789   0.34113547  1.15815364  0.33417725
15  0.6712018   0.20592338   0.21029137  0.69997845  0.20566324
16 -0.6629284  -0.20338513  -0.22270023 -0.79482144 -0.21782566
17  0.4363603   0.13387449   0.13803929  0.46393280  0.13492400
18  3.4486213   1.05803019   1.11295196  3.81594602  1.11933065
19  1.7931935   0.55014821   0.57876634  1.98460588  0.56981420
20 -5.7879699  -1.77573772  -1.87354643 -6.44313971 -1.99667657
21 -2.6141789  -0.80202492  -0.87784258 -3.13179171 -0.87308697
22 -3.6865279  -1.13101946  -1.44999541 -6.05913500 -1.48962473
23 -4.6075679  -1.41359270  -1.44368977 -4.80585779 -1.48246718
24 -4.5728535  -1.40294240  -1.49605875 -5.20001871 -1.54221512
25 -0.2125839  -0.06522033  -0.06750861 -0.22776283 -0.06596332

Nas análises de resíduos é indicado usar o resíduo studentizado ou o resíduo R-Student. É esperado que os valores dos resíduos sejam uma amostra da distribuição \(t\) com \(n-p-1\) graus de liberdade. Para tamanho amostral suficientemente grande, pode-se comparar com a distribuição normal padrão.

Gráficos de resíduos

# Gráfico qqnorm
qqnorm(res_3)
qqline(res_3, col = "blue", lwd = 2)

# Gráfico de dispersão y_est x resíduos
y_est <- fitted(mod)
plot(y_est, res_3, xlim = c(0,80), ylim = c(-2,5),
     main = "Gráfico de dispersão y_est x resíduos",
     xlab = "y estimado",
     ylab = "Resíduos studentizados")

# Gráfico dos resíduos seguindo a ordem de coleta
# ou execução dos dados
# atenção: supondo que a coleta tenha sido na ordem da variável ID, o que nem sempre é verdade!
plot(res_3 ~ dados$ID)

Testes para os resíduos

Normalidade dos resíduos

Alguns testes são: Qui-quadrado, Kolmogorov-Smirnov, Jarque-Bera, Shapiro-Wilk, Anderson-Darling, Cramér-von Mises, D’Agostino-Pearson, Lilliefors, e Shapiro-Francia.

# Teste de Shapiro-Wilk
shapiro.test(res_3)


    Shapiro-Wilk normality test

data:  res_3
W = 0.92285, p-value = 0.05952

# Teste de Kolmogorov-Smirnov
ks.test(res_3, "pnorm", mean = 0, sd = 1)


    Exact one-sample Kolmogorov-Smirnov test

data:  res_3
D = 0.17188, p-value = 0.4052
alternative hypothesis: two-sided

Independência dos Resíduos

Alguns testes são: Durbin-Watson, Breusch-Godfrey e Ljung-Box.

Para testar a independência dos resíduos é necessário conhecer a ordem de coleta (ou execução) dos dados.

Supondo que a coleta dos dados tenha sido conforme a variável \(ID\) (no exemplo delivery time), então, o teste DW é:

library(lmtest)
dwtest(res_3 ~ dados$ID, alternative = "two.sided")


    Durbin-Watson test

data:  res_3 ~ dados$ID
DW = 1.3994, p-value = 0.07306
alternative hypothesis: true autocorrelation is not 0

Homocedasticidade dos Resíduos

Alguns testes são: Bartlett, Breusch-Pagan, Levene, Samiuddin, O’Neill e Mathews, Layard, Park, White, Cochran, Hartley e Goldfeld-Quandt.

# Teste de Breusch-Pagan
# library(lmtest)
bptest(mod)


    studentized Breusch-Pagan test

data:  mod
BP = 11.988, df = 2, p-value = 0.002493

Valores discrepantes ou outliers

# Análise de outliers
cbind(dados, res_3)

   ID     y x1   x2       res_3
1   1 16.68  7  560 -1.62767993
2   2 11.50  3  220  0.36484267
3   3 12.03  3  340 -0.01609165
4   4 14.88  4   80  1.57972040
5   5 13.75  6  150 -0.14176094
6   6 18.11  7  330 -0.09080847
7   7  8.00  2  110  0.27042496
8   8 17.83  7  210  0.36672118
9   9 79.24 30 1460  3.21376278
10 10 21.50  5  605  0.81325432
11 11 40.33 16  688  0.71807970
12 12 21.00 10  215 -0.19325733
13 13 13.50  4  255  0.32517935
14 14 19.75  6  462  0.34113547
15 15 24.00  9  448  0.21029137
16 16 29.00 10  776 -0.22270023
17 17 15.35  6  200  0.13803929
18 18 19.00  7  132  1.11295196
19 19  9.50  3   36  0.57876634
20 20 35.10 17  770 -1.87354643
21 21 17.90 10  140 -0.87784258
22 22 52.32 26  810 -1.44999541
23 23 18.75  9  450 -1.44368977
24 24 19.83  8  635 -1.49605875
25 25 10.75  4  150 -0.06750861

# é necessário verificar se a observação 9 é um outlier.


# gráfico dos resíduos
plot(res_3, main = "Resíduos studentizados",
     ylab = "Resíduos studentizados",
     xlab = "Observações")
abline(h = 2, col = "red", lty = 2)
abline(h = -2, col = "red", lty = 2)

# Se for, de fato, um outlier, o que fazer?
# Estimar novamente o modelo sem a observação 9.

mod_2 <- lm(y ~ x1 + x2, data = dados[-9,])
summary(mod_2)


Call:
lm(formula = y ~ x1 + x2, data = dados[-9, ])

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0325 -1.2331  0.0199  1.4730  4.8167 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.447238   0.952469   4.669 0.000131 ***
x1          1.497691   0.130207  11.502 1.58e-10 ***
x2          0.010324   0.002854   3.618 0.001614 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.43 on 21 degrees of freedom
Multiple R-squared:  0.9487,    Adjusted R-squared:  0.9438 
F-statistic: 194.2 on 2 and 21 DF,  p-value: 2.859e-14

# comparar os resultados dos dois modelos:
summary(mod)


Call:
lm(formula = y ~ x1 + x2, data = dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.7880 -0.6629  0.4364  1.1566  7.4197 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.341231   1.096730   2.135 0.044170 *  
x1          1.615907   0.170735   9.464 3.25e-09 ***
x2          0.014385   0.003613   3.981 0.000631 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.259 on 22 degrees of freedom
Multiple R-squared:  0.9596,    Adjusted R-squared:  0.9559 
F-statistic: 261.2 on 2 and 22 DF,  p-value: 4.687e-16

summary(mod_2)


Call:
lm(formula = y ~ x1 + x2, data = dados[-9, ])

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0325 -1.2331  0.0199  1.4730  4.8167 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.447238   0.952469   4.669 0.000131 ***
x1          1.497691   0.130207  11.502 1.58e-10 ***
x2          0.010324   0.002854   3.618 0.001614 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.43 on 21 degrees of freedom
Multiple R-squared:  0.9487,    Adjusted R-squared:  0.9438 
F-statistic: 194.2 on 2 and 21 DF,  p-value: 2.859e-14

# Estatística PRESS
PRESS_mod <- sum(res_4^2)
PRESS_mod

[1] 459.0393

res_4_mod_2 <- rstandard(mod_2, type = "pred")
PRESS_mod_2 <- sum(res_4_mod_2^2)
PRESS_mod_2

[1] 165.2484

# Comparar as estatísticas PRESS
PRESS_mod

[1] 459.0393

PRESS_mod_2

[1] 165.2484

Pontos de alavancagem (leverage points)

# valores matriz H
p <- 3 # qde de parâmetros (b0, b1 e b2)
#ou
p <- length(coef(mod))


n <- 25 # qde de observaçoes
# ou
n <- length(dados$y)

2*p/n  # limite para valores da matriz hat

[1] 0.24

h <- hatvalues(mod)
cbind(dados, h, 2*p/n, h > 2*p/n)

   ID     y x1   x2          h 2 * p/n h > 2 * p/n
1   1 16.68  7  560 0.10180178    0.24       FALSE
2   2 11.50  3  220 0.07070164    0.24       FALSE
3   3 12.03  3  340 0.09873476    0.24       FALSE
4   4 14.88  4   80 0.08537479    0.24       FALSE
5   5 13.75  6  150 0.07501050    0.24       FALSE
6   6 18.11  7  330 0.04286693    0.24       FALSE
7   7  8.00  2  110 0.08179867    0.24       FALSE
8   8 17.83  7  210 0.06372559    0.24       FALSE
9   9 79.24 30 1460 0.49829216    0.24        TRUE
10 10 21.50  5  605 0.19629595    0.24       FALSE
11 11 40.33 16  688 0.08613260    0.24       FALSE
12 12 21.00 10  215 0.11365570    0.24       FALSE
13 13 13.50  4  255 0.06112463    0.24       FALSE
14 14 19.75  6  462 0.07824332    0.24       FALSE
15 15 24.00  9  448 0.04111077    0.24       FALSE
16 16 29.00 10  776 0.16594043    0.24       FALSE
17 17 15.35  6  200 0.05943202    0.24       FALSE
18 18 19.00  7  132 0.09626046    0.24       FALSE
19 19  9.50  3   36 0.09644857    0.24       FALSE
20 20 35.10 17  770 0.10168486    0.24       FALSE
21 21 17.90 10  140 0.16527689    0.24       FALSE
22 22 52.32 26  810 0.39157522    0.24        TRUE
23 23 18.75  9  450 0.04126005    0.24       FALSE
24 24 19.83  8  635 0.12060826    0.24       FALSE
25 25 10.75  4  150 0.06664345    0.24       FALSE

# observações 9 e 22 são possíveis pontos de alavancagem.

Distância de Cook

# Distância de Cook
D <- cooks.distance(mod)

# comparar com o valor 1
# Se D > 1, a observação é considerada um ponto influente.

cbind(dados, D, 1, D > 1)

   ID     y x1   x2            D 1 D > 1
1   1 16.68  7  560 1.000921e-01 1 FALSE
2   2 11.50  3  220 3.375704e-03 1 FALSE
3   3 12.03  3  340 9.455785e-06 1 FALSE
4   4 14.88  4   80 7.764718e-02 1 FALSE
5   5 13.75  6  150 5.432217e-04 1 FALSE
6   6 18.11  7  330 1.231067e-04 1 FALSE
7   7  8.00  2  110 2.171604e-03 1 FALSE
8   8 17.83  7  210 3.051135e-03 1 FALSE
9   9 79.24 30 1460 3.419318e+00 1  TRUE
10 10 21.50  5  605 5.384516e-02 1 FALSE
11 11 40.33 16  688 1.619975e-02 1 FALSE
12 12 21.00 10  215 1.596392e-03 1 FALSE
13 13 13.50  4  255 2.294737e-03 1 FALSE
14 14 19.75  6  462 3.292786e-03 1 FALSE
15 15 24.00  9  448 6.319880e-04 1 FALSE
16 16 29.00 10  776 3.289086e-03 1 FALSE
17 17 15.35  6  200 4.013419e-04 1 FALSE
18 18 19.00  7  132 4.397807e-02 1 FALSE
19 19  9.50  3   36 1.191868e-02 1 FALSE
20 20 35.10 17  770 1.324449e-01 1 FALSE
21 21 17.90 10  140 5.086063e-02 1 FALSE
22 22 52.32 26  810 4.510455e-01 1 FALSE
23 23 18.75  9  450 2.989892e-02 1 FALSE
24 24 19.83  8  635 1.023224e-01 1 FALSE
25 25 10.75  4  150 1.084694e-04 1 FALSE

DFBETAS

# DFBETAS
# a influência de cada observação em cada coeficiente
# comparar com o valor 2/sqrt(n)

dfbetas(mod)

    (Intercept)            x1           x2
1  -0.187267279  0.4113118750 -0.434862094
2   0.089793299 -0.0477642427  0.014414155
3  -0.003515177  0.0039483525 -0.002846468
4   0.451964743  0.0882802920 -0.273373097
5  -0.031674102 -0.0133001129  0.024240457
6  -0.014681480  0.0017921068  0.001078986
7   0.078071285 -0.0222783194 -0.011018802
8   0.071202807  0.0333823324 -0.053823961
9  -2.575739806  0.9287433421  1.507550618
10  0.107919369 -0.3381628707  0.341326746
11 -0.034274535  0.0925271540 -0.002686252
12 -0.030268935 -0.0486664488  0.053973390
13  0.072366473 -0.0356212226  0.011335105
14  0.049516699 -0.0670868604  0.061816778
15  0.022279094 -0.0047895025  0.006838236
16 -0.002693186  0.0644208340 -0.084187552
17  0.028855555  0.0064876499 -0.015696507
18  0.248558020  0.1897331043 -0.272430555
19  0.172558506  0.0235737344 -0.098968842
20  0.168036548 -0.2149950233 -0.092915080
21 -0.161928685 -0.2971750929  0.336406248
22  0.398566309 -1.0254140704  0.573140240
23 -0.159852248  0.0372930389 -0.052651959
24 -0.119720216  0.4046225960 -0.465446949
25 -0.016816024  0.0008498979  0.005592192

2/sqrt(n)

[1] 0.4

round(dfbetas(mod), 2)

   (Intercept)    x1    x2
1        -0.19  0.41 -0.43
2         0.09 -0.05  0.01
3         0.00  0.00  0.00
4         0.45  0.09 -0.27
5        -0.03 -0.01  0.02
6        -0.01  0.00  0.00
7         0.08 -0.02 -0.01
8         0.07  0.03 -0.05
9        -2.58  0.93  1.51
10        0.11 -0.34  0.34
11       -0.03  0.09  0.00
12       -0.03 -0.05  0.05
13        0.07 -0.04  0.01
14        0.05 -0.07  0.06
15        0.02  0.00  0.01
16        0.00  0.06 -0.08
17        0.03  0.01 -0.02
18        0.25  0.19 -0.27
19        0.17  0.02 -0.10
20        0.17 -0.21 -0.09
21       -0.16 -0.30  0.34
22        0.40 -1.03  0.57
23       -0.16  0.04 -0.05
24       -0.12  0.40 -0.47
25       -0.02  0.00  0.01

abs(dfbetas(mod)) >= 2/sqrt(n)

   (Intercept)    x1    x2
1        FALSE  TRUE  TRUE
2        FALSE FALSE FALSE
3        FALSE FALSE FALSE
4         TRUE FALSE FALSE
5        FALSE FALSE FALSE
6        FALSE FALSE FALSE
7        FALSE FALSE FALSE
8        FALSE FALSE FALSE
9         TRUE  TRUE  TRUE
10       FALSE FALSE FALSE
11       FALSE FALSE FALSE
12       FALSE FALSE FALSE
13       FALSE FALSE FALSE
14       FALSE FALSE FALSE
15       FALSE FALSE FALSE
16       FALSE FALSE FALSE
17       FALSE FALSE FALSE
18       FALSE FALSE FALSE
19       FALSE FALSE FALSE
20       FALSE FALSE FALSE
21       FALSE FALSE FALSE
22       FALSE  TRUE  TRUE
23       FALSE FALSE FALSE
24       FALSE  TRUE  TRUE
25       FALSE FALSE FALSE

DFFITS

# DFFITS
# a influência de cada observação em cada predição
# comparar com o valor 2*sqrt(p/n)

dffits(mod)

           1            2            3            4            5            6 
-0.570850478  0.098618619 -0.005203676  0.500801817 -0.039458989 -0.018779374 
           7            8            9           10           11           12 
 0.078990030  0.093760764  4.296080927  0.398713071  0.217953207 -0.067670223 
          13           14           15           16           17           18 
 0.081259033  0.097362643  0.042584374 -0.097159801  0.033915978  0.365309285 
          19           20           21           22           23           24 
 0.186167873 -0.671771402 -0.388501185 -1.195036104 -0.307538544 -0.571139627 
          25 
-0.017626149

2*sqrt(p/n)

[1] 0.6928203

round(dffits(mod), 2)

    1     2     3     4     5     6     7     8     9    10    11    12    13 
-0.57  0.10 -0.01  0.50 -0.04 -0.02  0.08  0.09  4.30  0.40  0.22 -0.07  0.08 
   14    15    16    17    18    19    20    21    22    23    24    25 
 0.10  0.04 -0.10  0.03  0.37  0.19 -0.67 -0.39 -1.20 -0.31 -0.57 -0.02

abs(dffits(mod)) >= 2*sqrt(p/n)

    1     2     3     4     5     6     7     8     9    10    11    12    13 
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE 
   14    15    16    17    18    19    20    21    22    23    24    25 
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE

Exercício

Considere os dados da Liga Nacional de Futebol Americano (NFL) na Tabela 1.

equipe <- c("Washington", "Minnesota", "New England", "Oakland", "Pittsburgh", 
            "Baltimore", "Los Angeles", "Dallas", "Atlanta", "Buffalo",
            "Chicago", "Cincinnati", "Cleveland", "Denver", "Detroit",
            "Green Bay", "Houston", "Kansas City", "Miami", "New Orleans",
            "New York Giants", "New York Jets", "Philadelphia", "St. Louis", "San Diego",
            "San Francisco", "Seattle", "Tampa Bay")

y <- c(10, 11, 11, 13, 10,
       11, 10, 11, 4, 2,
       7, 10, 9, 9, 6,
       5, 5, 5, 6, 4, 
       3, 3, 4, 10, 6, 
       8, 2, 0)

x1 <- c(2113, 2003, 2957, 2285, 2971,
        2309, 2528, 2147, 1689, 2566,
        2363, 2109, 2295, 1932, 2128,
        1722, 1498, 1873, 2118, 1775,
        1904, 1929, 2080, 2301, 2040,
        2447, 1416, 1503)

x2 <- c(1985, 2855, 1737, 2905, 1666,
        2927, 2341, 2737, 1414, 1838,
        1480, 2191, 2229, 2204, 2438,
        1730, 2072, 2929, 2268, 1983,
        1792, 1606, 1492, 2835, 2416,
        1638, 2649, 1503)

x3 <- c(38.9, 38.8, 40.1, 41.6, 39.2,
        39.7, 38.1, 37.0, 42.1, 42.3,
        37.3, 39.5, 37.4, 35.1, 38.8,
        36.6, 35.3, 41.1, 38.2, 39.3,
        39.7, 39.7, 35.5, 35.3, 38.7,
        39.9, 37.4, 39.3)


x4 <- c(64.7, 61.3, 60.0, 45.3, 53.8,
        74.1, 65.4, 78.3, 47.6, 54.2,
        48.0, 51.9, 53.6, 71.4, 58.3,
        52.6, 59.3, 55.3, 69.6, 78.3,
        38.1, 68.8, 68.8, 74.1, 50.0,
        57.1, 56.3, 47.0)

x5 <- c(4, 3, 14, -4, 15,
        8, 12, -1, -3, -1,
        19, 6, -5, 3, 6,
        -19, -5, 10, 6, 7,
        -9, -21, -8, 2, 0,
        -8, -22, -9)

x6 <- c(868, 615, 914, 957, 836,
        786, 754, 761, 714, 797,
        984, 700, 1037, 986, 819,
        791, 776, 789, 582, 901,
        734, 627, 722, 683, 576,
        848, 684, 875)

x7 <- c(59.7, 55.0, 65.6, 61.4, 66.1,
        61.0, 66.1, 58.0, 57.0, 58.9,
        67.5, 57.2, 57.8, 58.6, 59.2,
        54.4, 49.6, 54.3, 58.7, 51.7,
        61.9, 52.7, 57.8, 59.7, 54.9,
        65.3, 43.8, 53.5)

x8 <- c(2205, 2096, 1847, 1903, 1457,
        1848, 1564, 1821, 2577, 2476,
        1984, 1917, 1761, 1790, 1901,
        2288, 2072, 2861, 2411, 2289,
        2203, 2592, 2053, 1979, 2048,
        1786, 2876, 2560)

x9 <- c(1917, 1575, 2175, 2476, 1866,
        2339, 2092, 1909, 2001, 2254,
        2217, 1758, 2032, 2025, 1686,
        1835, 1914, 2496, 2670, 2202,
        1988, 2324, 2550, 2110, 2628,
        1776, 2524, 2241)




tabela1 <- data.frame(equipe, y, x1, x2, x3, x4, x5, x6, x7, x8, x9)

knitr::kable(tabela1, caption = "Tabela 1: Desempenho das Equipes da National Football League de 1976")

Tabela 1: Desempenho das Equipes da National Football League de 1976
equipe	y	x1	x2	x3	x4	x5	x6	x7	x8	x9
Washington	10	2113	1985	38.9	64.7	4	868	59.7	2205	1917
Minnesota	11	2003	2855	38.8	61.3	3	615	55.0	2096	1575
New England	11	2957	1737	40.1	60.0	14	914	65.6	1847	2175
Oakland	13	2285	2905	41.6	45.3	-4	957	61.4	1903	2476
Pittsburgh	10	2971	1666	39.2	53.8	15	836	66.1	1457	1866
Baltimore	11	2309	2927	39.7	74.1	8	786	61.0	1848	2339
Los Angeles	10	2528	2341	38.1	65.4	12	754	66.1	1564	2092
Dallas	11	2147	2737	37.0	78.3	-1	761	58.0	1821	1909
Atlanta	4	1689	1414	42.1	47.6	-3	714	57.0	2577	2001
Buffalo	2	2566	1838	42.3	54.2	-1	797	58.9	2476	2254
Chicago	7	2363	1480	37.3	48.0	19	984	67.5	1984	2217
Cincinnati	10	2109	2191	39.5	51.9	6	700	57.2	1917	1758
Cleveland	9	2295	2229	37.4	53.6	-5	1037	57.8	1761	2032
Denver	9	1932	2204	35.1	71.4	3	986	58.6	1790	2025
Detroit	6	2128	2438	38.8	58.3	6	819	59.2	1901	1686
Green Bay	5	1722	1730	36.6	52.6	-19	791	54.4	2288	1835
Houston	5	1498	2072	35.3	59.3	-5	776	49.6	2072	1914
Kansas City	5	1873	2929	41.1	55.3	10	789	54.3	2861	2496
Miami	6	2118	2268	38.2	69.6	6	582	58.7	2411	2670
New Orleans	4	1775	1983	39.3	78.3	7	901	51.7	2289	2202
New York Giants	3	1904	1792	39.7	38.1	-9	734	61.9	2203	1988
New York Jets	3	1929	1606	39.7	68.8	-21	627	52.7	2592	2324
Philadelphia	4	2080	1492	35.5	68.8	-8	722	57.8	2053	2550
St. Louis	10	2301	2835	35.3	74.1	2	683	59.7	1979	2110
San Diego	6	2040	2416	38.7	50.0	0	576	54.9	2048	2628
San Francisco	8	2447	1638	39.9	57.1	-8	848	65.3	1786	1776
Seattle	2	1416	2649	37.4	56.3	-22	684	43.8	2876	2524
Tampa Bay	0	1503	1503	39.3	47.0	-9	875	53.5	2560	2241

Definições das variáveis:

\(y\): Jogos vencidos (por temporada de 14 jogos)
\(x_1\): Jardas conquistadas em corrida (temporada)
\(x_2\): Jardas de passe (temporada)
\(x_3\): Média de punt (jardas/punt)
\(x_4\): Percentual de acerto em field goal (field goals convertidos/field goals tentados na temporada)
\(x_5\): Diferença de turnovers (turnovers ganhos - turnovers perdidos)
\(x_6\): Jardas de penalidade (temporada)
\(x_7\): Percentual de corrida (corridas/jogadas totais)
\(x_8\): Jardas conquistadas em corrida pelos adversários (temporada)
\(x_9\): Jardas de passe dos adversários (temporada)

Ajuste um modelo de regressão linear múltipla relacionando o número de jogos ganhos com as jardas aéreas do time (\(x_2\)), a porcentagem de jogadas terrestres (\(x_7\)) e as jardas terrestres dos adversários (\(x_8\)).

m1 <- lm(y  ~ x2 + x7 + x8, data = tabela1)
m1


Call:
lm(formula = y ~ x2 + x7 + x8, data = tabela1)

Coefficients:
(Intercept)           x2           x7           x8  
  -0.792280     0.003439     0.185537    -0.004911

Construa a tabela de análise de variância (ANOVA) e teste a significância da regressão.

anova(m1)

Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq F value    Pr(>F)    
x2         1  73.256  73.256  23.271 6.497e-05 ***
x7         1 135.885 135.885  43.166 8.534e-07 ***
x8         1  42.273  42.273  13.429  0.001224 ** 
Residuals 24  75.550   3.148                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Calcule as estatísticas \(t\) para testar as hipóteses: \(H_0: \beta_2 = 0\), \(H_0: \beta_7 = 0\) e \(H_0: \beta_8 = 0\). Quais conclusões você pode tirar sobre os papéis das variáveis \(x_2\), \(x_7\) e \(x_8\) no modelo?

summary(m1)


Call:
lm(formula = y ~ x2 + x7 + x8, data = tabela1)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2379 -0.6299 -0.1298  1.0467  3.7201 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.7922801  8.2481313  -0.096  0.92427    
x2           0.0034385  0.0007194   4.780 7.27e-05 ***
x7           0.1855367  0.0917292   2.023  0.05439 .  
x8          -0.0049115  0.0013403  -3.665  0.00122 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.774 on 24 degrees of freedom
Multiple R-squared:  0.7689,    Adjusted R-squared:  0.7401 
F-statistic: 26.62 on 3 and 24 DF,  p-value: 8.279e-08

Calcule o \(R^2\) e o \(R^2_{\text{aj}}\) (ajustado) para este modelo.

summary(m1)


Call:
lm(formula = y ~ x2 + x7 + x8, data = tabela1)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2379 -0.6299 -0.1298  1.0467  3.7201 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.7922801  8.2481313  -0.096  0.92427    
x2           0.0034385  0.0007194   4.780 7.27e-05 ***
x7           0.1855367  0.0917292   2.023  0.05439 .  
x8          -0.0049115  0.0013403  -3.665  0.00122 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.774 on 24 degrees of freedom
Multiple R-squared:  0.7689,    Adjusted R-squared:  0.7401 
F-statistic: 26.62 on 3 and 24 DF,  p-value: 8.279e-08

names(summary(m1))

 [1] "call"          "terms"         "residuals"     "coefficients" 
 [5] "aliased"       "sigma"         "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"

summary(m1)$r.squared

[1] 0.7689342

summary(m1)$adj.r.squared

[1] 0.740051

Usando o teste \(F\) parcial, determine a contribuição de \(x_7\) para o modelo. Como essa estatística \(F\) parcial está relacionada ao teste \(t\) para \(\beta_7\) calculado no item c?

m2 <- lm(y ~ x2 + x8, data = tabela1)
m2


Call:
lm(formula = y ~ x2 + x8, data = tabela1)

Coefficients:
(Intercept)           x2           x8  
  15.064942     0.002968    -0.006837

anova(m2, m1)

Analysis of Variance Table

Model 1: y ~ x2 + x8
Model 2: y ~ x2 + x7 + x8
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     25 88.429                              
2     24 75.550  1    12.879 4.0911 0.05439 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Construa um gráfico de probabilidade normal dos resíduos. Há indícios de violação da suposição de normalidade? Realize os testes de normalidade de Shapiro-Wilk e Kolmogorov-Smirnov. O que você conclui?

res_m1 <- rstandard(m1) # resíduos studentizados
qqnorm(res_m1)
qqline(res_m1)

shapiro.test(res_m1)


    Shapiro-Wilk normality test

data:  res_m1
W = 0.9683, p-value = 0.5358

ks.test(res_m1, "pnorm", mean = 0, sd = 1)


    Exact one-sample Kolmogorov-Smirnov test

data:  res_m1
D = 0.14736, p-value = 0.5294
alternative hypothesis: two-sided

Calcule os resíduos studentizados e os resíduos Rstudent para este modelo. Que tipo de informação é fornecida por esses resíduos?

res_m1 # resíduos studentizados já calculados

          1           2           3           4           5           6 
 2.15516954  1.23879394  1.63300701  1.07819228 -0.02783659 -0.31470640 
          7           8           9          10          11          12 
-1.12142891  0.33781927  1.24773451 -1.38022492 -0.04837882  1.20332312 
         13          14          15          16          17          18 
 0.03160084  0.07830376 -1.88268058  0.58375011 -0.23461680 -0.21896915 
         19          20          21          22          23          24 
-0.03321221 -0.22520769 -1.79683917  0.74508590 -0.58702034 -0.18784834 
         25          26          27          28 
-0.96252871 -0.11057780 -0.21672035 -1.05742134

# resíduos R-Student
res_m1_rstudent <- rstudent(m1)
cbind(res_m1, res_m1_rstudent)

        res_m1 res_m1_rstudent
1   2.15516954      2.34934086
2   1.23879394      1.25344711
3   1.63300701      1.69559868
4   1.07819228      1.08202138
5  -0.02783659     -0.02725093
6  -0.31470640     -0.30871791
7  -1.12142891     -1.12776223
8   0.33781927      0.33149557
9   1.24773451      1.26311855
10 -1.38022492     -1.40820857
11 -0.04837882     -0.04736251
12  1.20332312      1.21521610
13  0.03160084      0.03093613
14  0.07830376      0.07666487
15 -1.88268058     -1.99634394
16  0.58375011      0.57555991
17 -0.23461680     -0.22994078
18 -0.21896915     -0.21457320
19 -0.03321221     -0.03251368
20 -0.22520769     -0.22069926
21 -1.79683917     -1.89077811
22  0.74508590      0.73798335
23 -0.58702034     -0.57883109
24 -0.18784834     -0.18402853
25 -0.96252871     -0.96099363
26 -0.11057780     -0.10827717
27 -0.21672035     -0.21236521
28 -1.05742134     -1.06014758

Construa e interprete o gráfico dos resíduos studentizados em função da resposta predita. Realize o teste de homocedasticidade de Breusch-Pagan. O que você conclui?

y_estimado <- fitted(m1) # valores ajustados
plot(y_estimado, res_m1,
     xlab = "Valores ajustados",
     ylab = "Resíduos studentizados"
     )
abline(h = -2, col = "red", lty = 2)
abline(h = 2, col = "red", lty = 2)

library(lmtest)
bptest(m1)


    studentized Breusch-Pagan test

data:  m1
BP = 2.2021, df = 3, p-value = 0.5315

Calcule as correlações entre as variáveis independentes \(x_2\), \(x_7\) e \(x_8\). Existe multicolinearidade entre essas variáveis? Justifique sua resposta.

cor(tabela1[, c("x2", "x7", "x8")], use = "pairwise.complete.obs")

            x2         x7          x8
x2  1.00000000 -0.1926007 -0.06269288
x7 -0.19260071  1.0000000 -0.68347083
x8 -0.06269288 -0.6834708  1.00000000

Calcule o VIF (Variance Inflation Factor) para as variáveis \(x_2\), \(x_7\) e \(x_8\). O que você conclui? Os valores do VIF são úteis para indicar multicolinearidade. Valores de VIF superior a 5 ou 10 sugerem multicolinearidade. Use a função vif() do pacote car para calcular o VIF.

library(car)

Loading required package: carData

vif(m1)

      x2       x7       x8 
1.121032 2.095503 2.025732

Calcule a estatística PRESS (Press Residual Sum of Squares) para os modelos completo (com \(x_2\), \(x_7\) e \(x_8\)) e reduzido (apenas com \(x_2\) e \(x_8\)). O que você conclui?

# Cálculo PRESS m1
res_press_m1 <- rstandard(m1, type = "pred")
PRESS_m1 <- sum(res_press_m1^2)
PRESS_m1

[1] 93.7735

# Cálculo PRESS m2
res_press_m2 <- rstandard(m2, type = "pred")
PRESS_m2 <- sum(res_press_m2^2)
PRESS_m2

[1] 106.1485

Calcule os valores dos pontos de alavancagem (leverage). Quais observações são consideradas influentes? Lembre-se que valores superiores a \(2p/n\) (onde \(p\) é o número de parâmetros do modelo e \(n\) é o número de observações) são considerados influentes.

# valores matriz H
p <- length(coef(m1)) # qde de parâmetros (b0, b1, b2 e b3)


n <- nrow(tabela1) # qde de observaçoes

2*p/n  # limite para valores da matriz hat

[1] 0.2857143

h <- hatvalues(m1)
cbind(tabela1[, c("y", "x2", "x7", "x8")], y_estimado, h, 2*p/n, h > 2*p/n)

    y   x2   x7   x8 y_estimado          h   2 * p/n h > 2 * p/n
1  10 1985 59.7 2205   6.279885 0.05349006 0.2857143       FALSE
2  11 2855 55.0 2096   8.934726 0.11705751 0.2857143       FALSE
3  11 1737 65.6 1847   8.280112 0.11874564 0.2857143       FALSE
4  13 2905 61.4 1903  11.242004 0.15546351 0.2857143       FALSE
5  10 1666 66.1 1457  10.044226 0.19812848 0.2857143       FALSE
6  11 2927 61.0 1848  11.513569 0.15401656 0.2857143       FALSE
7  10 2341 66.1 1564  11.839697 0.14508199 0.2857143       FALSE
8  11 2737 58.0 1821  10.436251 0.11533567 0.2857143       FALSE
9   4 1414 57.0 2577   1.988467 0.17437071 0.2857143       FALSE
10  2 1838 58.9 2476   4.294979 0.12171963 0.2857143       FALSE
11  7 1480 67.5 1984   7.076058 0.21484443 0.2857143       FALSE
12 10 2191 57.2 1917   7.938887 0.06800333 0.2857143       FALSE
13  9 2229 57.8 1761   8.947065 0.10862272 0.2857143       FALSE
14  9 2204 58.6 1790   8.867098 0.08489476 0.2857143       FALSE
15  6 2438 59.2 1901   9.237858 0.06041149 0.2857143       FALSE
16  5 1730 54.4 2288   4.012064 0.09013164 0.2857143       FALSE
17  5 2072 49.6 2072   5.358344 0.25893477 0.2857143       FALSE
18  5 2929 54.3 2861   5.302011 0.39569550 0.2857143        TRUE
19  6 2268 58.7 2411   6.055682 0.10708379 0.2857143       FALSE
20  4 1983 51.7 2289   4.376149 0.11380529 0.2857143       FALSE
21  3 1792 61.9 2203   6.034254 0.09414052 0.2857143       FALSE
22  3 1606 52.7 2592   1.777183 0.14437079 0.2857143       FALSE
23  4 1492 57.8 2053   4.978722 0.11694558 0.2857143       FALSE
24 10 2835 59.7 1979  10.312623 0.12016535 0.2857143       FALSE
25  6 2416 54.9 2048   7.642414 0.07506052 0.2857143       FALSE
26  8 1638 65.3 1786   8.183638 0.12387496 0.2857143       FALSE
27  2 2649 43.8 2876   2.317419 0.31853977 0.2857143        TRUE
28  0 1503 53.5 2560   1.728612 0.15106502 0.2857143       FALSE

which(h > 2*p/n)

18 27 
18 27

Calcule os valores da distância de Cook. Quais observações são consideradas influentes? Lembre-se que valores superiores a 1 são considerados influentes.

# Distância de Cook
D <- cooks.distance(m1)

# comparar com o valor 1
# Se D > 1, a observação é considerada um ponto influente.

cbind(tabela1[, c("y", "x2", "x7", "x8")], y_estimado, D, 1, D > 1)

    y   x2   x7   x8 y_estimado            D 1 D > 1
1  10 1985 59.7 2205   6.279885 6.562220e-02 1 FALSE
2  11 2855 55.0 2096   8.934726 5.086336e-02 1 FALSE
3  11 1737 65.6 1847   8.280112 8.983230e-02 1 FALSE
4  13 2905 61.4 1903  11.242004 5.349861e-02 1 FALSE
5  10 1666 66.1 1457  10.044226 4.786457e-05 1 FALSE
6  11 2927 61.0 1848  11.513569 4.507718e-03 1 FALSE
7  10 2341 66.1 1564  11.839697 5.335468e-02 1 FALSE
8  11 2737 58.0 1821  10.436251 3.719580e-03 1 FALSE
9   4 1414 57.0 2577   1.988467 8.220019e-02 1 FALSE
10  2 1838 58.9 2476   4.294979 6.600354e-02 1 FALSE
11  7 1480 67.5 1984   7.076058 1.601102e-04 1 FALSE
12 10 2191 57.2 1917   7.938887 2.641316e-02 1 FALSE
13  9 2229 57.8 1761   8.947065 3.042261e-05 1 FALSE
14  9 2204 58.6 1790   8.867098 1.422051e-04 1 FALSE
15  6 2438 59.2 1901   9.237858 5.697380e-02 1 FALSE
16  5 1730 54.4 2288   4.012064 8.439032e-03 1 FALSE
17  5 2072 49.6 2072   5.358344 4.808307e-03 1 FALSE
18  5 2929 54.3 2861   5.302011 7.848943e-03 1 FALSE
19  6 2268 58.7 2411   6.055682 3.307110e-05 1 FALSE
20  4 1983 51.7 2289   4.376149 1.628320e-03 1 FALSE
21  3 1792 61.9 2203   6.034254 8.388304e-02 1 FALSE
22  3 1606 52.7 2592   1.777183 2.341782e-02 1 FALSE
23  4 1492 57.8 2053   4.978722 1.140887e-02 1 FALSE
24 10 2835 59.7 1979  10.312623 1.204850e-03 1 FALSE
25  6 2416 54.9 2048   7.642414 1.879601e-02 1 FALSE
26  8 1638 65.3 1786   8.183638 4.322085e-04 1 FALSE
27  2 2649 43.8 2876   2.317419 5.488612e-03 1 FALSE
28  0 1503 53.5 2560   1.728612 4.974227e-02 1 FALSE

Calcule os valores de DFBETAS. Quais observações são consideradas influentes? Lembre-se que valores superiores a \(2/\sqrt{n}\) merecem atenção.

# DFBETAS
# a influência de cada observação em cada coeficiente
# comparar com o valor 2/sqrt(n)

dfbetas(m1)

    (Intercept)           x2            x7           x8
1  -0.238761486 -0.040456057  2.678571e-01  0.267125084
2   0.040033901  0.313636747 -1.021143e-01 -0.064830777
3  -0.196885826 -0.166337545  3.456996e-01  0.053364028
4  -0.200196651  0.380035349  1.852540e-01  0.054283262
5  -0.004010632  0.005856577  3.938946e-05  0.007745337
6   0.040040046 -0.105304797 -3.693300e-02  0.002588051
7   0.086414024 -0.125821399 -1.763182e-01  0.121286090
8   0.029309162  0.063807642 -2.976293e-02 -0.056095124
9  -0.129007939 -0.301851886  1.335750e-01  0.319286268
10  0.308141383  0.051364852 -2.940507e-01 -0.405393028
11  0.012814808  0.007294419 -1.760943e-02 -0.009328888
12  0.212424295 -0.042828098 -1.835889e-01 -0.221636976
13  0.007610617 -0.001329077 -6.341772e-03 -0.008771910
14  0.014362879 -0.002515488 -1.121180e-02 -0.017565145
15 -0.065620490 -0.198249783  3.844244e-02  0.180321590
16  0.095696698 -0.115663577 -9.170044e-02 -0.029770584
17 -0.124451269  0.047409617  1.258636e-01  0.093791450
18  0.122330111 -0.113225722 -8.484365e-02 -0.141247451
19  0.008207135 -0.004249292 -7.116382e-03 -0.008968230
20 -0.056389087  0.031419481  6.045758e-02  0.028108274
21  0.332694544  0.116525776 -3.939454e-01 -0.335059866
22  0.056493195 -0.166000228 -7.421115e-02  0.083167543
23 -0.117084527  0.172185093  8.220421e-02  0.078963885
24  0.023857025 -0.055240703 -1.982562e-02 -0.006551652
25 -0.152685965 -0.045583775  1.648678e-01  0.136713875
26  0.004088923  0.017322230 -1.505958e-02  0.004446726
27 -0.030507073 -0.031979480  6.794611e-02 -0.027451532
28 -0.083438767  0.279622005  9.472017e-02 -0.112893575

2/sqrt(n)

[1] 0.3779645

round(dfbetas(m1), 2)

   (Intercept)    x2    x7    x8
1        -0.24 -0.04  0.27  0.27
2         0.04  0.31 -0.10 -0.06
3        -0.20 -0.17  0.35  0.05
4        -0.20  0.38  0.19  0.05
5         0.00  0.01  0.00  0.01
6         0.04 -0.11 -0.04  0.00
7         0.09 -0.13 -0.18  0.12
8         0.03  0.06 -0.03 -0.06
9        -0.13 -0.30  0.13  0.32
10        0.31  0.05 -0.29 -0.41
11        0.01  0.01 -0.02 -0.01
12        0.21 -0.04 -0.18 -0.22
13        0.01  0.00 -0.01 -0.01
14        0.01  0.00 -0.01 -0.02
15       -0.07 -0.20  0.04  0.18
16        0.10 -0.12 -0.09 -0.03
17       -0.12  0.05  0.13  0.09
18        0.12 -0.11 -0.08 -0.14
19        0.01  0.00 -0.01 -0.01
20       -0.06  0.03  0.06  0.03
21        0.33  0.12 -0.39 -0.34
22        0.06 -0.17 -0.07  0.08
23       -0.12  0.17  0.08  0.08
24        0.02 -0.06 -0.02 -0.01
25       -0.15 -0.05  0.16  0.14
26        0.00  0.02 -0.02  0.00
27       -0.03 -0.03  0.07 -0.03
28       -0.08  0.28  0.09 -0.11

abs(dfbetas(m1)) >= 2/sqrt(n)

   (Intercept)    x2    x7    x8
1        FALSE FALSE FALSE FALSE
2        FALSE FALSE FALSE FALSE
3        FALSE FALSE FALSE FALSE
4        FALSE  TRUE FALSE FALSE
5        FALSE FALSE FALSE FALSE
6        FALSE FALSE FALSE FALSE
7        FALSE FALSE FALSE FALSE
8        FALSE FALSE FALSE FALSE
9        FALSE FALSE FALSE FALSE
10       FALSE FALSE FALSE  TRUE
11       FALSE FALSE FALSE FALSE
12       FALSE FALSE FALSE FALSE
13       FALSE FALSE FALSE FALSE
14       FALSE FALSE FALSE FALSE
15       FALSE FALSE FALSE FALSE
16       FALSE FALSE FALSE FALSE
17       FALSE FALSE FALSE FALSE
18       FALSE FALSE FALSE FALSE
19       FALSE FALSE FALSE FALSE
20       FALSE FALSE FALSE FALSE
21       FALSE FALSE  TRUE FALSE
22       FALSE FALSE FALSE FALSE
23       FALSE FALSE FALSE FALSE
24       FALSE FALSE FALSE FALSE
25       FALSE FALSE FALSE FALSE
26       FALSE FALSE FALSE FALSE
27       FALSE FALSE FALSE FALSE
28       FALSE FALSE FALSE FALSE

Calcule os valores de DFFITS. Quais observações são consideradas influentes? Lembre-se que valores superiores a \(2\sqrt{p/n}\) merecem atenção.

# DFFITS
# a influência de cada observação em cada predição
# comparar com o valor 2*sqrt(p/n)

dffits(m1)

          1           2           3           4           5           6 
 0.55849589  0.45639350  0.62241674  0.46423820 -0.01354573 -0.13172391 
          7           8           9          10          11          12 
-0.46458167  0.11969342  0.58048144 -0.52424068 -0.02477530  0.32825505 
         13          14          15          16          17          18 
 0.01079930  0.02335078 -0.50620485  0.18115064 -0.13591988 -0.17363134 
         19          20          21          22          23          24 
-0.01125959 -0.07908926 -0.60953470  0.30314013 -0.21064445 -0.06801023 
         25          26          27          28 
-0.27375978 -0.04071417 -0.14519275 -0.44720955

2*sqrt(p/n)

[1] 0.7559289

round(dffits(m1), 2)

    1     2     3     4     5     6     7     8     9    10    11    12    13 
 0.56  0.46  0.62  0.46 -0.01 -0.13 -0.46  0.12  0.58 -0.52 -0.02  0.33  0.01 
   14    15    16    17    18    19    20    21    22    23    24    25    26 
 0.02 -0.51  0.18 -0.14 -0.17 -0.01 -0.08 -0.61  0.30 -0.21 -0.07 -0.27 -0.04 
   27    28 
-0.15 -0.45

abs(dffits(m1)) >= 2*sqrt(p/n)

    1     2     3     4     5     6     7     8     9    10    11    12    13 
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
   14    15    16    17    18    19    20    21    22    23    24    25    26 
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
   27    28 
FALSE FALSE

Construa um modelo de regressão linear relacionando o número de jogos vencidos com:

jardas terrestres dos adversários (\(x_8\)),
porcentagem de jogadas terrestres (\(x_7\)) e
diferença de turnovers (turnovers ganhos - turnovers perdidos) (\(x_5\)).

Especificamente, considere a diferença de turnovers como uma variável indicadora, cuja codificação depende de o valor real do diferencial ser positivo, negativo ou igual a zero. Ou seja, considere \(x_5\) como uma variável indicadora que assume os seguintes valores:

\[\begin{align*} x_5 = 1, & \text{ se } x_5 > 0 \\ x_5 = 0, & \text{ se } x_5 = 0 \\ x_5 = -1, & \text{ se } x_5 < 0. \end{align*}\]

Quais conclusões podem ser obtidas sobre o efeito dos turnovers no número de jogos vencidos?

I_x5 <- ifelse(x5 > 0, 1, ifelse(x5 == 0, 0, -1))
I_x5 <- as.factor(I_x5)

m3 <- lm(y ~ x7 + x8 + I_x5)
summary(m3)


Call:
lm(formula = y ~ x7 + x8 + I_x5)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2975 -1.6308  0.0043  1.2600  5.6613 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 20.708747   9.650005   2.146 0.042658 *  
x7          -0.016338   0.118371  -0.138 0.891426    
x8          -0.006499   0.001719  -3.781 0.000968 ***
I_x50       -0.502574   2.452150  -0.205 0.839412    
I_x51        1.909872   0.958422   1.993 0.058288 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.325 on 23 degrees of freedom
Multiple R-squared:  0.6198,    Adjusted R-squared:  0.5537 
F-statistic: 9.374 on 4 and 23 DF,  p-value: 0.0001202

Considerando o enunciado da letra anterior, e se a variável indicadora for definida como:

\[\begin{align*} x_5 = 1, & \text{ se } x_5 > 0 \\ x_5 = 0, & \text{ se } x_5 \leq 0. \end{align*}\]

Qual é a interpretação do coeficiente estimado para \(x_5\)? O que você conclui sobre o efeito dos turnovers no número de jogos vencidos?

I2_x5 <- ifelse(x5 > 0, 1, 0)
I2_x5 <- as.factor(I2_x5)
m4 <- lm(y ~ x7 + x8 + I2_x5)
summary(m4)


Call:
lm(formula = y ~ x7 + x8 + I2_x5)

Residuals:
   Min     1Q Median     3Q    Max 
-3.288 -1.646  0.008  1.264  5.697 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 20.332528   9.282800   2.190 0.038451 *  
x7          -0.012662   0.114646  -0.110 0.912973    
x8          -0.006438   0.001659  -3.880 0.000713 ***
I2_x51       1.944113   0.924721   2.102 0.046194 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.278 on 24 degrees of freedom
Multiple R-squared:  0.6191,    Adjusted R-squared:  0.5715 
F-statistic:    13 on 3 and 24 DF,  p-value: 3.025e-05

Considerando todas as variáveis disponíveis, utilize o método de seleção de variáveis forward (seleção progressiva) para selecionar um modelo de regressão com subconjunto de variáveis. Comente sobre o modelo final. O método de seleção de variáveis forward (seleção progressiva) é um método de seleção de variáveis que começa com um modelo nulo (sem variáveis preditoras) e adiciona variáveis preditoras uma a uma, com base em critérios estatísticos, até que não haja mais variáveis significativas a serem adicionadas.

# Modelo reduzido
# Modelo mínimo a ser considerado.
# O modelo mínimo pode ser um modelo nulo (sem variáveis preditoras)
modelo_reduzido <- lm(y ~ 1, data = tabela1)

# Forward selection
modelo_forward <- step(modelo_reduzido,
    scope = as.formula(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9),
    direction = "forward")

Start:  AIC=70.81
y ~ 1

       Df Sum of Sq    RSS    AIC
+ x8    1   178.688 148.28 50.673
+ x1    1   116.251 210.71 60.512
+ x7    1    95.849 231.12 63.100
+ x5    1    86.116 240.85 64.255
+ x2    1    73.256 253.71 65.711
+ x9    1    30.167 296.80 70.104
<none>              326.96 70.814
+ x4    1    21.844 305.12 70.878
+ x6    1    16.411 310.55 71.372
+ x3    1     2.135 324.83 72.631

Step:  AIC=50.67
y ~ x8

       Df Sum of Sq     RSS    AIC
+ x2    1    59.848  88.429 38.200
+ x5    1    11.394 136.883 50.434
<none>              148.277 50.673
+ x4    1     6.867 141.410 51.345
+ x1    1     6.507 141.770 51.416
+ x3    1     5.515 142.762 51.611
+ x7    1     0.803 147.474 52.521
+ x6    1     0.309 147.967 52.614
+ x9    1     0.011 148.266 52.670

Step:  AIC=38.2
y ~ x8 + x2

       Df Sum of Sq    RSS    AIC
+ x7    1   12.8786 75.550 35.793
+ x1    1   11.3495 77.079 36.354
+ x3    1    7.5539 80.875 37.700
<none>              88.429 38.200
+ x5    1    5.6023 82.827 38.367
+ x6    1    1.4090 87.020 39.750
+ x9    1    1.1860 87.243 39.822
+ x4    1    0.1458 88.283 40.154

Step:  AIC=35.79
y ~ x8 + x2 + x7

       Df Sum of Sq    RSS    AIC
<none>              75.550 35.793
+ x9    1    3.2411 72.309 36.565
+ x1    1    1.4090 74.141 37.266
+ x4    1    1.3332 74.217 37.294
+ x3    1    1.0597 74.491 37.397
+ x6    1    0.6291 74.921 37.558
+ x5    1    0.3118 75.238 37.677

summary(modelo_forward)


Call:
lm(formula = y ~ x8 + x2 + x7, data = tabela1)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2379 -0.6299 -0.1298  1.0467  3.7201 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.7922801  8.2481313  -0.096  0.92427    
x8          -0.0049115  0.0013403  -3.665  0.00122 ** 
x2           0.0034385  0.0007194   4.780 7.27e-05 ***
x7           0.1855367  0.0917292   2.023  0.05439 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.774 on 24 degrees of freedom
Multiple R-squared:  0.7689,    Adjusted R-squared:  0.7401 
F-statistic: 26.62 on 3 and 24 DF,  p-value: 8.279e-08

Utilize o método de seleção de variáveis backward (eliminacão progressiva) para selecionar um modelo de regressão com subconjunto de variáveis. Comente sobre o modelo final. O método de seleção de variáveis backward (eliminacão progressiva) é um método de seleção de variáveis que começa com um modelo completo (com todas as variáveis preditoras) e remove variáveis preditoras uma a uma, com base em critérios estatísticos, até que não haja mais variáveis insignificativas a serem removidas.

# Modelo completo
# Modelo máximo a ser considerado.
# O modelo máximo pode ser um modelo completo (com todas as variáveis preditoras)

modelo_completo <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9,
                      data = tabela1)

# Backward elimination
modelo_backward <- step(modelo_completo,
                        direction = "backward")

Start:  AIC=44.29
y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9

       Df Sum of Sq     RSS    AIC
- x5    1     0.000  66.676 42.294
- x3    1     0.562  67.238 42.529
- x6    1     0.999  67.675 42.710
- x1    1     1.006  67.681 42.713
- x4    1     2.105  68.780 43.164
- x7    1     2.620  69.296 43.373
- x9    1     4.131  70.806 43.977
<none>               66.676 44.294
- x8    1    12.316  78.992 47.040
- x2    1    55.721 122.397 59.302

Step:  AIC=42.29
y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 + x9

       Df Sum of Sq     RSS    AIC
- x3    1     0.563  67.239 40.529
- x1    1     1.008  67.684 40.714
- x6    1     1.015  67.691 40.717
- x4    1     2.183  68.859 41.196
- x7    1     2.954  69.630 41.508
- x9    1     4.130  70.807 41.977
<none>               66.676 42.294
- x8    1    12.340  79.016 45.049
- x2    1    60.340 127.017 58.339

Step:  AIC=40.53
y ~ x1 + x2 + x4 + x6 + x7 + x8 + x9

       Df Sum of Sq     RSS    AIC
- x6    1     1.069  68.308 38.971
- x4    1     1.646  68.884 39.206
- x1    1     2.162  69.401 39.415
- x7    1     3.300  70.538 39.871
<none>               67.239 40.529
- x9    1     4.980  72.219 40.530
- x8    1    14.416  81.655 43.968
- x2    1    65.897 133.135 57.657

Step:  AIC=38.97
y ~ x1 + x2 + x4 + x7 + x8 + x9

       Df Sum of Sq     RSS    AIC
- x4    1     1.393  69.701 37.536
- x1    1     1.956  70.264 37.762
- x7    1     3.775  72.083 38.477
- x9    1     4.938  73.246 38.925
<none>               68.308 38.971
- x8    1    16.722  85.030 43.103
- x2    1    64.835 133.143 55.658

Step:  AIC=37.54
y ~ x1 + x2 + x7 + x8 + x9

       Df Sum of Sq     RSS    AIC
- x1    1     2.608  72.309 36.565
- x7    1     2.824  72.525 36.648
- x9    1     4.440  74.141 37.266
<none>               69.701 37.536
- x8    1    20.291  89.992 42.690
- x2    1    70.097 139.798 55.024

Step:  AIC=36.56
y ~ x2 + x7 + x8 + x9

       Df Sum of Sq     RSS    AIC
- x9    1     3.241  75.550 35.793
<none>               72.309 36.565
- x7    1    14.934  87.243 39.822
- x8    1    25.945  98.254 43.150
- x2    1    75.165 147.474 54.520

Step:  AIC=35.79
y ~ x2 + x7 + x8

       Df Sum of Sq     RSS    AIC
<none>               75.550 35.793
- x7    1    12.879  88.429 38.200
- x8    1    42.273 117.824 46.236
- x2    1    71.924 147.474 52.521

summary(modelo_backward)


Call:
lm(formula = y ~ x2 + x7 + x8, data = tabela1)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2379 -0.6299 -0.1298  1.0467  3.7201 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.7922801  8.2481313  -0.096  0.92427    
x2           0.0034385  0.0007194   4.780 7.27e-05 ***
x7           0.1855367  0.0917292   2.023  0.05439 .  
x8          -0.0049115  0.0013403  -3.665  0.00122 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.774 on 24 degrees of freedom
Multiple R-squared:  0.7689,    Adjusted R-squared:  0.7401 
F-statistic: 26.62 on 3 and 24 DF,  p-value: 8.279e-08

Utilize o método de seleção de variáveis stepwise (passo a passo) para selecionar um modelo de regressão com subconjunto de variáveis. Comente sobre o modelo final. O método de seleção de variáveis stepwise (passo a passo) é um método de seleção de variáveis que combina os métodos forward e backward. Ele começa com um modelo nulo (sem variáveis preditoras) e adiciona variáveis preditoras uma a uma, com base em critérios estatísticos, mas também pode remover variáveis preditoras que se tornam não significativas à medida que novas variáveis são adicionadas.

# Stepwise regression
modelo_stepwise <- step(modelo_reduzido,
        scope = formula(modelo_completo),
        direction = "both")

Start:  AIC=70.81
y ~ 1

       Df Sum of Sq    RSS    AIC
+ x8    1   178.688 148.28 50.673
+ x1    1   116.251 210.71 60.512
+ x7    1    95.849 231.12 63.100
+ x5    1    86.116 240.85 64.255
+ x2    1    73.256 253.71 65.711
+ x9    1    30.167 296.80 70.104
<none>              326.96 70.814
+ x4    1    21.844 305.12 70.878
+ x6    1    16.411 310.55 71.372
+ x3    1     2.135 324.83 72.631

Step:  AIC=50.67
y ~ x8

       Df Sum of Sq    RSS    AIC
+ x2    1    59.848  88.43 38.200
+ x5    1    11.394 136.88 50.434
<none>              148.28 50.673
+ x4    1     6.867 141.41 51.345
+ x1    1     6.507 141.77 51.416
+ x3    1     5.515 142.76 51.611
+ x7    1     0.803 147.47 52.521
+ x6    1     0.309 147.97 52.614
+ x9    1     0.011 148.27 52.670
- x8    1   178.688 326.96 70.814

Step:  AIC=38.2
y ~ x8 + x2

       Df Sum of Sq     RSS    AIC
+ x7    1    12.879  75.550 35.793
+ x1    1    11.349  77.079 36.354
+ x3    1     7.554  80.875 37.700
<none>               88.429 38.200
+ x5    1     5.602  82.827 38.367
+ x6    1     1.409  87.020 39.750
+ x9    1     1.186  87.243 39.822
+ x4    1     0.146  88.283 40.154
- x2    1    59.848 148.277 50.673
- x8    1   165.280 253.709 65.711

Step:  AIC=35.79
y ~ x8 + x2 + x7

       Df Sum of Sq     RSS    AIC
<none>               75.550 35.793
+ x9    1     3.241  72.309 36.565
+ x1    1     1.409  74.141 37.266
+ x4    1     1.333  74.217 37.294
+ x3    1     1.060  74.491 37.397
+ x6    1     0.629  74.921 37.558
+ x5    1     0.312  75.238 37.677
- x7    1    12.879  88.429 38.200
- x8    1    42.273 117.824 46.236
- x2    1    71.924 147.474 52.521

summary(modelo_stepwise)


Call:
lm(formula = y ~ x8 + x2 + x7, data = tabela1)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2379 -0.6299 -0.1298  1.0467  3.7201 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.7922801  8.2481313  -0.096  0.92427    
x8          -0.0049115  0.0013403  -3.665  0.00122 ** 
x2           0.0034385  0.0007194   4.780 7.27e-05 ***
x7           0.1855367  0.0917292   2.023  0.05439 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.774 on 24 degrees of freedom
Multiple R-squared:  0.7689,    Adjusted R-squared:  0.7401 
F-statistic: 26.62 on 3 and 24 DF,  p-value: 8.279e-08

Usando o método de validação cruzada leave-one-out (LOOCV), estime o erro de previsão para os modelos (1) com as variáveis regressoras \(x_2\), \(x_7\) e \(x_8\) e (2) com as variáveis regressoras \(x_2\) e \(x_8\). O que você conclui? O método de validação cruzada leave-one-out (LOOCV) é uma técnica de validação cruzada onde cada observação do conjunto de dados é usada como um ponto de teste, enquanto todas as outras observações são usadas como pontos de treinamento. Isso é repetido para cada observação, e o erro médio de previsão é calculado. Existe também o método leave-p-out.

library(boot)


Attaching package: 'boot'

The following object is masked from 'package:car':

    logit

# LOOCV
m1 <- glm(y ~ x2 + x7 + x8,
              data = tabela1,
              family = gaussian(link = "identity"))
K <- nrow(tabela1) # k subconjuntos (folds)
cv_m1 <- cv.glm(data = tabela1, glmfit = m1, K = K)
cv_m1$delta  # delta[1] = erro de validação cruzada

[1] 3.349054 3.336679

m2 <- glm(y ~ x2 + x8,
              data = tabela1,
              family = gaussian(link = "identity"))
cv_m2 <- cv.glm(data = tabela1, glmfit = m2, K = K)
cv_m2$delta  # delta[1] = erro de validação cruzada

[1] 3.791018 3.778988

Usando o método de validação cruzada k-fold, com \(k = 3\), estime o erro de previsão para os modelos (1) com as variáveis regressoras \(x_2\), \(x_7\) e \(x_8\) e (2) com as variáveis regressoras \(x_2\) e \(x_8\). O que você conclui? O método de validação cruzada k-fold é uma técnica de validação cruzada onde o conjunto de dados é dividido em k subconjuntos (folds). O modelo é treinado em k-1 folds e testado no fold restante. Isso é repetido k vezes, e o erro médio de previsão é calculado.

# k-fold
m1 <- glm(y ~ x2 + x7 + x8,
              data = tabela1,
              family = gaussian(link = "identity"))

K <- 3 # k subconjuntos (folds)
cv_m1 <- cv.glm(data = tabela1, glmfit = m1, K = K)
cv_m1$delta  # delta[1] = erro de validação cruzada

[1] 3.963470 3.692716

m2 <- glm(y ~ x2 + x8,
              data = tabela1,
              family = gaussian(link = "identity"))
cv_m2 <- cv.glm(data = tabela1, glmfit = m2, K = K)
cv_m2$delta  # delta[1] = erro de validação cruzada

[1] 4.224191 4.000779

@profjulianobortolini www.julianobortolini.com.br linkedin github lattes