computacional
Bacharelado em Estatística - UFMT
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 \]
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 373summary(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.594683sd(dados$adsorcao)[1] 133.3552plot(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\(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 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 ' ' 1Estimativas, 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.001263Testes 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.306004Para 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.01529526t_tab <- qt(0.05/2, df = 8, lower.tail = FALSE) # t tabulado
t_tab[1] 2.306004t_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)
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.98974dados$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)
residuos <- residuals(modelo)
SQRes <- sum(residuos^2)
SQRes[1] 40551.32GL_Res <- length(dados$adsorcao) - length(coef(modelo))
GL_Res[1] 8QMRes <- SQRes / GL_Res
QMRes[1] 5068.914sqrt(QMRes)[1] 71.19631preditos <- 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.2810colnames(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.2810dados <- 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.6877Ordernar 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.2810grid <- 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) 
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# 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) 
# 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) 
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 ' ' 1summary(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-08anova(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.7372Se não há diferença entre os modelos, qual deve ser escolhido?
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.7372Se não há diferença entre os modelos, qual deve ser escolhido?
O modelo mais simples: com menos parâmetros.
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.001263summary(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-08Seja \(Y\) a adsorção de corante, \(X_1\) a concentração, \(X_2\) a temperatura e \(X_3\) o pH. O modelo de regressão linear múltipla, que relaciona a adsorção com as variáveis explicativas é dado por: \[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 + \varepsilon, \]
em que, \(\varepsilon \sim N(0, \sigma^2)\), \(Y\) é a variável resposta (adsorção) e \(X_1\), \(X_2\) e \(X_3\) são as variáveis explicativas (concentração, temperatura e pH, respectivamente).
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 373O modelo de regreão linear múltipla é estimado pela função lm() do R, que é uma função genérica para ajuste de modelos lineares. O modelo é ajustado pela função lm() usando a fórmula adsorcao ~ conc + temp + ph, que indica que a variável resposta é adsorcao e as variáveis explicativas são conc, temp e ph.
# 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 A análise de variância do modelo é realizada pela função anova(), que é uma função genérica para análise de variância de modelos ajustados. O modelo é passado como argumento para a função anova(), que retorna a tabela de análise de variância do modelo, incluindo as somas de quadrados, graus de liberdade, quadrados médios, estatística \(F\) e valor \(p\) para o teste de significância global do modelo.
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 ' ' 1coef(modelo3) # coeficientes(Intercept) conc temp ph
53.758696 21.548478 3.657174 -90.271739 A função que calcula os intervalos de confiança para os coeficientes é a função confint(), que é uma função genérica para calcular intervalos de confiança para os parâmetros de um modelo ajustado. O modelo é passado como argumento para a função confint(), que retorna os intervalos de confiança para os coeficientes do modelo. O padrão do nível de confiança é 95%, mas pode ser alterado usando o argumento level.
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.722065Uma função que apresenta os coeficientes estimados, erros-padrões, estatísticas \(t\) e valores \(p\) é a função summary(), que é uma função genérica para sumarizar os resultados de um modelo ajustado. O modelo é sumarizado pela função summary(modelo3), que apresenta as informações mencionadas.
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.004493Será apresentado a seguir como calcular os intervalos de confiança para a resposta média, ou seja, para a predição da adsorção de corante para cada combinação de concentração, temperatura e pH presentes no banco de dados. Para isso, é utilizada a função predict(), que é uma função genérica para fazer predições a partir de um modelo ajustado. O modelo é passado como argumento para a função predict(), juntamente com o argumento interval = "confidence", que indica que os intervalos de confiança devem ser calculados para os valores médios da resposta. A função retorna uma matriz com as predições da resposta média, os limites inferior e superior dos intervalos de confiança para cada combinação de variáveis explicativas presentes no banco de dados.
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.5575colnames(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.5575dados <- 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.9863Quando queremos calcular os intervalos de predição para novas observações, ou seja, para combinações de concentração, temperatura e pH que não estão presentes no banco de dados, utilizamos a função predict() com o argumento interval = "prediction". A função retorna uma matriz com as predições da resposta para as novas observações, os limites inferior e superior dos intervalos de predição para cada combinação de variáveis explicativas presentes no banco de dados. Para isso, é necessário criar um novo banco de dados com as combinações de concentração, temperatura e pH para as quais se deseja calcular os intervalos de predição. A seguir, é apresentado um exemplo de como criar um novo banco de dados e calcular os intervalos de predição para essas novas observações.
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.6610A análise de resíduos é uma etapa importante na avaliação do ajuste do modelo de regressão. Os resíduos são as diferenças entre os valores observados e os valores ajustados pelo modelo. A análise de resíduos envolve a verificação de suposições do modelo, como normalidade, homocedasticidade e independência dos resíduos, além da identificação de possíveis outliers ou pontos influentes.
Há diversos tipos de resíduos, como resíduos brutos, resíduos padronizados, resíduos studentizados, resíduos PRESS e resíduos R-Student. Cada tipo de resíduo tem suas próprias características e é utilizado para diferentes propósitos na análise de resíduos.
O resíduo bruto é a diferença simples entre o valor observado e o valor ajustado pelo modelo. Ou seja, \(e_i = y_i - \hat{y}_i\), onde \(y_i\) é o valor observado e \(\hat{y}_i\) é o valor ajustado pelo modelo para a observação \(i\).
O resíduo padronizado é o resíduo bruto dividido pelo desvio padrão dos resíduos, o que permite comparar os resíduos em diferentes escalas. Isso é dado por \(d_i^* = \frac{e_i}{\hat{\sigma}}\), onde \(\hat{\sigma}\) é a estimativa do desvio padrão dos resíduos. Em geral, o valor de \(\sigma\) é estimado pelo modelo de regressão, e é dado por \(\hat{\sigma} = \sqrt{QMRes}\), onde \(QMRes\) é o quadrado médio dos resíduos.
Na definição do resíduo studentizado (internamente), considera-se variância do resíduo bruto e a matriz \(H\). O resíduo studentizado é definido por \(r_i = \frac{e_i}{\hat{\sigma} \sqrt{1 - h_{ii}}}\), onde \(h_{ii}\) é o elemento da matriz \(H\) correspondente à observação \(i\). O resíduo studentizado é uma medida de quão grande é o resíduo em relação à variabilidade dos resíduos, levando em consideração a influência da observação \(i\) no ajuste do modelo. A justificativa para esse resíduo é que \(Var(e) = Var(Y - \hat{Y}) = Var((I - H)Y) = (I - H)\sigma^2\), consequentemente, \(Var(e_i) = (1 - h_{ii})\sigma^2\), sendo \(h_{ii}\) o elemento da diagonal da matriz \(H\), que representa a influência da observação \(i\) no ajuste do modelo. Assim, o resíduo studentizado leva em consideração a variabilidade dos resíduos e a influência da observação \(i\) no ajuste do modelo, o que o torna uma medida mais robusta para identificar outliers e pontos influentes.
O resíduo PRESS é o resíduo studentizado calculado para cada observação, considerando a exclusão dessa observação do ajuste do modelo. Matematicamente, tem-se que \(e_(i) = y_i - \hat{y}_{(i)}\), onde \(\hat{y}_{(i)}\) é o valor ajustado pelo modelo de regressão quando as observações do indivíduo \(i\) são excluídas do ajuste. O resíduo PRESS é utilizado para avaliar a capacidade de predição do modelo, pois considera a influência de cada observação no ajuste do modelo e a variabilidade dos resíduos.
O resíduo R-Student é o resíduo studentizado calculado para cada observação, considerando a exclusão dessa observação do ajuste do modelo e a variabilidade dos resíduos. Ou seja, \(t_i = \frac{e_i}{\hat{\sigma}_{(i)} \sqrt{1 - h_{ii}}}\), onde \(\hat{\sigma}_{(i)}\) é a estimativa do desvio padrão dos resíduos quando as observações do indivíduo \(i\) são excluídas do ajuste do modelo. O resíduo R-Student é utilizado para identificar outliers e pontos influentes, pois leva em consideração a influência das observações do indivíduo \(i\) no ajuste do modelo e a variabilidade dos resíduos, mesmo quando essa observação é excluída do ajuste.
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)
head(dados) ID y x1 x2
1 1 16.68 7 560
2 2 11.50 3 220
3 3 12.03 3 340
4 4 14.88 4 80
5 5 13.75 6 150
6 6 18.11 7 330# descritiva dos dados
summary(dados) ID y x1 x2
Min. : 1 Min. : 8.00 Min. : 2.00 Min. : 36.0
1st Qu.: 7 1st Qu.:13.75 1st Qu.: 4.00 1st Qu.: 150.0
Median :13 Median :18.11 Median : 7.00 Median : 330.0
Mean :13 Mean :22.38 Mean : 8.76 Mean : 409.3
3rd Qu.:19 3rd Qu.:21.50 3rd Qu.:10.00 3rd Qu.: 605.0
Max. :25 Max. :79.24 Max. :30.00 Max. :1460.0 # Ajuste do modelo linear múltipla
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# valores preditos
y_hat <- predict(mod)
y_hat 1 2 3 4 5 6 7 8
21.708084 10.353615 12.079794 9.955646 14.194398 18.399574 7.155376 16.673395
9 10 11 12 13 14 15 16
71.820294 19.123587 38.092507 21.593041 12.472991 18.682464 23.328798 29.662928
17 18 19 20 21 22 23 24
14.913640 15.551379 7.706807 40.887970 20.514179 56.006528 23.357568 24.402854
25
10.962584 dados <- cbind(dados, y_hat)
head(dados) ID y x1 x2 y_hat
1 1 16.68 7 560 21.708084
2 2 11.50 3 220 10.353615
3 3 12.03 3 340 12.079794
4 4 14.88 4 80 9.955646
5 5 13.75 6 150 14.194398
6 6 18.11 7 330 18.399574# estimativa de sigma^2
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# Análise de variância do modelo
anova(mod)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x1 1 5382.4 5382.4 506.619 < 2.2e-16 ***
x2 1 168.4 168.4 15.851 0.0006312 ***
Residuals 22 233.7 10.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Testes de hipóteses para parâmetros
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# Intervalo de confiança para os parâmetros
confint(mod) 2.5 % 97.5 %
(Intercept) 0.066751987 4.61571030
x1 1.261824662 1.96998976
x2 0.006891745 0.02187791# alterando o nível de confiança
confint(mod, level = 0.90) 5 % 95 %
(Intercept) 0.457987107 4.22447518
x1 1.322730706 1.90908372
x2 0.008180636 0.02058902# IC para y_hat(x) (x1 = 8 e x2 = 275)
predict(mod,
newdata = data.frame(x1 = 8, x2 = 275),
interval = "confidence") fit lwr upr
1 19.22432 17.6539 20.79474# IP para y_hat(x) (x1 = 8 e x2 = 275)
predict(mod,
newdata = data.frame(x1 = 8, x2 = 275),
interval = "prediction") fit lwr upr
1 19.22432 12.28456 26.16407# IP para diversos valores de x1 e x2:
x_novo <- data.frame(x1 = c(8, 10, 12), x2 = c(275, 300, 325))
predict(mod, newdata = x_novo, interval = "prediction", level = 0.95) fit lwr upr
1 19.22432 12.28456 26.16407
2 22.81575 15.81734 29.81416
3 26.40719 19.30571 33.50866# modelo de regressão polinomial (grau 2) com x1
# $Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i1}^2 + \epsilon_i$
mod2 <- lm(y ~ x1 + I(x1^2), data = dados)
summary(mod2)
Call:
lm(formula = y ~ x1 + I(x1^2), data = dados)
Residuals:
Min 1Q Median 3Q Max
-8.7819 -2.4305 -0.0255 1.5719 6.7452
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.21264 2.30994 3.555 0.00177 **
x1 1.13562 0.42985 2.642 0.01489 *
I(x1^2) 0.03456 0.01379 2.507 0.02005 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.771 on 22 degrees of freedom
Multiple R-squared: 0.9459, Adjusted R-squared: 0.941
F-statistic: 192.4 on 2 and 22 DF, p-value: 1.155e-14# modelo de regressão polinomial (grau 2) com x2
# $Y_i = \beta_0 + \beta_1 X_{i2} + \beta_2 X_{i2}^2 + \epsilon_i$
mod3 <- lm(y ~ x2 + I(x2^2), data = dados)
summary(mod3)
Call:
lm(formula = y ~ x2 + I(x2^2), data = dados)
Residuals:
Min 1Q Median 3Q Max
-8.218 -3.340 -1.617 3.096 16.096
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.161e+01 2.721e+00 4.268 0.000314 ***
x2 8.800e-03 1.045e-02 0.842 0.408772
I(x2^2) 2.665e-05 7.708e-06 3.457 0.002245 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.909 on 22 degrees of freedom
Multiple R-squared: 0.8672, Adjusted R-squared: 0.8551
F-statistic: 71.83 on 2 and 22 DF, p-value: 2.265e-10# modelo linear múltipla com interação entre x1 e x2
# $Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \beta_3 X_{i1}X_{i2} + \epsilon_i$
mod4 <- lm(y ~ x1 + x2 + x1:x2, data = dados)
# ou
# mod4 <- lm(y ~ x1 * x2, data = dados)
summary(mod4)
Call:
lm(formula = y ~ x1 + x2 + x1:x2, data = dados)
Residuals:
Min 1Q Median 3Q Max
-3.7316 -1.5387 0.0606 1.4375 4.7841
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.1390846 1.3997413 5.100 4.73e-05 ***
x1 1.0144063 0.1912517 5.304 2.93e-05 ***
x2 0.0058273 0.0033825 1.723 0.099622 .
x1:x2 0.0007419 0.0001750 4.240 0.000366 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.449 on 21 degrees of freedom
Multiple R-squared: 0.9782, Adjusted R-squared: 0.9751
F-statistic: 314.6 on 3 and 21 DF, p-value: < 2.2e-16# Comparando dois modelos encaixados:
# Definição: O Modelo 1 (M1) é encaixado no Modelo 2 (M2) se todos os parâmetros de M1 estão contidos em M2, e M2 possui parâmetros adicionais.
# mod é encaixado em mod4, pois mod4 tem os mesmos parâmetros de mod mais um parâmetro adicional (interação x1:x2).
anova(mod, mod4)Analysis of Variance Table
Model 1: y ~ x1 + x2
Model 2: y ~ x1 + x2 + x1:x2
Res.Df RSS Df Sum of Sq F Pr(>F)
1 22 233.73
2 21 125.92 1 107.81 17.979 0.0003659 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod5 <- lm(y ~ x1 + x1:x2, data = dados)
anova(mod5, mod4)Analysis of Variance Table
Model 1: y ~ x1 + x1:x2
Model 2: y ~ x1 + x2 + x1:x2
Res.Df RSS Df Sum of Sq F Pr(>F)
1 22 143.72
2 21 125.92 1 17.797 2.968 0.09962 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_nulo <- lm(y ~ 1, data = dados)
anova(mod_nulo, mod)Analysis of Variance Table
Model 1: y ~ 1
Model 2: y ~ x1 + x2
Res.Df RSS Df Sum of Sq F Pr(>F)
1 24 5784.5
2 22 233.7 2 5550.8 261.24 4.687e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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.06596332Nas 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á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)
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-sidedAlguns 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 0Alguns 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# Análise de outliers
cbind(dados, res_3) ID y x1 x2 y_hat res_3
1 1 16.68 7 560 21.708084 -1.62767993
2 2 11.50 3 220 10.353615 0.36484267
3 3 12.03 3 340 12.079794 -0.01609165
4 4 14.88 4 80 9.955646 1.57972040
5 5 13.75 6 150 14.194398 -0.14176094
6 6 18.11 7 330 18.399574 -0.09080847
7 7 8.00 2 110 7.155376 0.27042496
8 8 17.83 7 210 16.673395 0.36672118
9 9 79.24 30 1460 71.820294 3.21376278
10 10 21.50 5 605 19.123587 0.81325432
11 11 40.33 16 688 38.092507 0.71807970
12 12 21.00 10 215 21.593041 -0.19325733
13 13 13.50 4 255 12.472991 0.32517935
14 14 19.75 6 462 18.682464 0.34113547
15 15 24.00 9 448 23.328798 0.21029137
16 16 29.00 10 776 29.662928 -0.22270023
17 17 15.35 6 200 14.913640 0.13803929
18 18 19.00 7 132 15.551379 1.11295196
19 19 9.50 3 36 7.706807 0.57876634
20 20 35.10 17 770 40.887970 -1.87354643
21 21 17.90 10 140 20.514179 -0.87784258
22 22 52.32 26 810 56.006528 -1.44999541
23 23 18.75 9 450 23.357568 -1.44368977
24 24 19.83 8 635 24.402854 -1.49605875
25 25 10.75 4 150 10.962584 -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-16summary(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.0393res_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.0393PRESS_mod_2[1] 165.2484Ponto de alavancagem é uma observação que tem valores incomuns nas variáveis preditoras \(x_i\), \(i=1,\cdots,k\), ou seja, é um ponto que está “longe” do centro de massa das demais observações no espaço das preditoras. Já um ponto influente é uma observação que, se removida do conjunto de dados, provoca uma mudança substancial nos resultados da análise (coeficientes estimados, valores preditos, erros padrão ou estatísticas de ajuste como \(R^2\)).
O valor de alavancagem \(\hat{h}_{ii}\) é o \(i\)-ésimo elemento da diagonal da matriz chapéu (hat matrix) \(H = X(X^TX)^{-1}X^T\). O valor de alavancagem é uma medida de quão longe a observação \(x_i\) está do centro das preditoras. O valor de alavancagem varia entre 0 e 1, e a média dos valores de alavancagem é dada por \(\bar{h} = \frac{p+1}{n}\), onde \(p\) é o número de preditores e \(n\) é o número de observações. Valores de alavancagem maiores que \(\frac{2p}{n}\) são merecem atenção.
# 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.24h <- hatvalues(mod)
cbind(dados, h, 2*p/n, h > 2*p/n) ID y x1 x2 y_hat h 2 * p/n h > 2 * p/n
1 1 16.68 7 560 21.708084 0.10180178 0.24 FALSE
2 2 11.50 3 220 10.353615 0.07070164 0.24 FALSE
3 3 12.03 3 340 12.079794 0.09873476 0.24 FALSE
4 4 14.88 4 80 9.955646 0.08537479 0.24 FALSE
5 5 13.75 6 150 14.194398 0.07501050 0.24 FALSE
6 6 18.11 7 330 18.399574 0.04286693 0.24 FALSE
7 7 8.00 2 110 7.155376 0.08179867 0.24 FALSE
8 8 17.83 7 210 16.673395 0.06372559 0.24 FALSE
9 9 79.24 30 1460 71.820294 0.49829216 0.24 TRUE
10 10 21.50 5 605 19.123587 0.19629595 0.24 FALSE
11 11 40.33 16 688 38.092507 0.08613260 0.24 FALSE
12 12 21.00 10 215 21.593041 0.11365570 0.24 FALSE
13 13 13.50 4 255 12.472991 0.06112463 0.24 FALSE
14 14 19.75 6 462 18.682464 0.07824332 0.24 FALSE
15 15 24.00 9 448 23.328798 0.04111077 0.24 FALSE
16 16 29.00 10 776 29.662928 0.16594043 0.24 FALSE
17 17 15.35 6 200 14.913640 0.05943202 0.24 FALSE
18 18 19.00 7 132 15.551379 0.09626046 0.24 FALSE
19 19 9.50 3 36 7.706807 0.09644857 0.24 FALSE
20 20 35.10 17 770 40.887970 0.10168486 0.24 FALSE
21 21 17.90 10 140 20.514179 0.16527689 0.24 FALSE
22 22 52.32 26 810 56.006528 0.39157522 0.24 TRUE
23 23 18.75 9 450 23.357568 0.04126005 0.24 FALSE
24 24 19.83 8 635 24.402854 0.12060826 0.24 FALSE
25 25 10.75 4 150 10.962584 0.06664345 0.24 FALSE# observações 9 e 22 são possíveis pontos de alavancagem.A distância de Cook é uma medida que combina a alavancagem e o resíduo para avaliar a influência de cada observação no ajuste do modelo. Ela é calculada como: \[D_i = \frac{(\hat{\boldsymbol{\beta}} - \hat{\boldsymbol{\beta}}_{(-i)})^T (X^TX) (\hat{\boldsymbol{\beta}} - \hat{\boldsymbol{\beta}}_{(-i)})}{(p+1)\hat{\sigma}^2}.\]
# 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 y_hat D 1 D > 1
1 1 16.68 7 560 21.708084 1.000921e-01 1 FALSE
2 2 11.50 3 220 10.353615 3.375704e-03 1 FALSE
3 3 12.03 3 340 12.079794 9.455785e-06 1 FALSE
4 4 14.88 4 80 9.955646 7.764718e-02 1 FALSE
5 5 13.75 6 150 14.194398 5.432217e-04 1 FALSE
6 6 18.11 7 330 18.399574 1.231067e-04 1 FALSE
7 7 8.00 2 110 7.155376 2.171604e-03 1 FALSE
8 8 17.83 7 210 16.673395 3.051135e-03 1 FALSE
9 9 79.24 30 1460 71.820294 3.419318e+00 1 TRUE
10 10 21.50 5 605 19.123587 5.384516e-02 1 FALSE
11 11 40.33 16 688 38.092507 1.619975e-02 1 FALSE
12 12 21.00 10 215 21.593041 1.596392e-03 1 FALSE
13 13 13.50 4 255 12.472991 2.294737e-03 1 FALSE
14 14 19.75 6 462 18.682464 3.292786e-03 1 FALSE
15 15 24.00 9 448 23.328798 6.319880e-04 1 FALSE
16 16 29.00 10 776 29.662928 3.289086e-03 1 FALSE
17 17 15.35 6 200 14.913640 4.013419e-04 1 FALSE
18 18 19.00 7 132 15.551379 4.397807e-02 1 FALSE
19 19 9.50 3 36 7.706807 1.191868e-02 1 FALSE
20 20 35.10 17 770 40.887970 1.324449e-01 1 FALSE
21 21 17.90 10 140 20.514179 5.086063e-02 1 FALSE
22 22 52.32 26 810 56.006528 4.510455e-01 1 FALSE
23 23 18.75 9 450 23.357568 2.989892e-02 1 FALSE
24 24 19.83 8 635 24.402854 1.023224e-01 1 FALSE
25 25 10.75 4 150 10.962584 1.084694e-04 1 FALSEA medida DFBETAS (do inglês Difference in Betas) avalia a influência de cada observação em cada coeficiente estimado. Ela é calculada como: \[\text{DFBETAS}_{ij} = \frac{\hat{\beta}_j - \hat{\beta}_{j(-i)}}{\hat{\sigma}_{(-i)} \sqrt{(X^TX)^{-1}_{jj}}},\] em que \(\hat{\beta}_j\) é o coeficiente estimado para o preditor \(j\) usando todos os dados, \(\hat{\beta}_{j(-i)}\) é o coeficiente estimado para o preditor \(j\) quando a observação \(i\) é removida, e \(\hat{\sigma}_{(-i)}\) é a estimativa do desvio padrão dos resíduos quando a observação \(i\) é removida. Valores absolutos de DFBETAS maiores que \(2/\sqrt{n}\) indicam que a observação tem uma influência significativa no coeficiente estimado.
# 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.0055921922/sqrt(n)[1] 0.4round(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.01abs(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 FALSEA medida DFFITS (do inglês Difference in Fits) avalia a influência de cada observação em cada predição. Ela é calculada como: \[\text{DFFITS}_i = \frac{\hat{Y}_i - \hat{Y}_{i(-i)}}{\hat{\sigma}_{(-i)} \sqrt{\hat{h}_{ii}}},\] em que \(\hat{Y}_i\) é o valor predito para a observação \(i\) usando todos os dados, \(\hat{Y}_{i(-i)}\) é o valor predito para a observação \(i\) quando a observação \(i\) é removida, \(\hat{\sigma}_{(-i)}\) é a estimativa do desvio padrão dos resíduos quando a observação \(i\) é removida, e \(\hat{h}_{ii}\) é o valor de alavancagem para a observação \(i\). Valores absolutos de DFFITS maiores que \(2\sqrt{p/n}\) indicam que a observação tem uma influência significativa na predição.
# 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.6928203round(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