x <- c(408, 408, 554, 554, 680, 680, 812, 812)
y <- c(1100, 1300, 1600, 2500, 3000, 4300, 4200, 4700)Análise de Regressão
Bacharelado em Estatística - UFMT
Atividade prática no R
O aquecimento global é uma questão importante, e as emissões de \(CO_2\) desempenham um papel fundamental nessa discussão. O artigo “Effects of Atmospheric \(CO_2\) Enrichment on Biomass Accumulation and Distribution in Eldarica Pine Trees” (J. Exp. Bot. 1994: 345–349) descreve os resultados do crescimento de pinheiros sob diferentes níveis de \(CO_2\) no ar.
Os dados abaixo correspondem às observações onde:
- \(x\) = concentração atmosférica de \(CO_2\) (partes por milhão)
- \(y\) = massa em gramas após 11 meses de experimento.
Tabela. Concentração de \(CO_2\) (\(x\)) e Massa (\(y\))
| \(x\) | 408 | 408 | 554 | 554 | 680 | 680 | 812 | 812 |
|---|---|---|---|---|---|---|---|---|
| \(y\) | 1100 | 1300 | 1600 | 2500 | 3000 | 4300 | 4200 | 4700 |
- Calcule e interprete o valor do coeficiente de correlação amostral.
cor(x,y)[1] 0.9392405- Realize o teste de hipóteses para decidir se a concentração de \(CO_2\) e a massa dos pinheiros estão linearmente relacionadas e interprete os resultados, considere \(\alpha = 0,05\).
cor.test(x,y)
Pearson's product-moment correlation
data: x and y
t = 6.7024, df = 6, p-value = 0.0005355
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6936940 0.9892025
sample estimates:
cor
0.9392405 - Calcule o intervalo de confiança de \(95\%\) para o coeficiente de correlação amostral e interprete-o.
cor.test(x,y)
Pearson's product-moment correlation
data: x and y
t = 6.7024, df = 6, p-value = 0.0005355
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6936940 0.9892025
sample estimates:
cor
0.9392405 - Ajuste um modelo de regressão linear simples relacionando \(x\) e \(y\) e interprete o coeficiente linear (\(\beta_1\)).
plot(x,y)
modelo <- lm(y ~ x)
modelo
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
-2349.295 8.454 - Realize a análise de variância (ANOVA) para o modelo ajustado e interprete os resultados, considere \(\alpha = 0,05\).
anova(modelo) Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 12808045 12808045 44.922 0.0005355 ***
Residuals 6 1710705 285118
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# a diferença entre os valores do cálculo no R e "na mão" é por causa das aproximações.- Realize o teste de hipóteses \(H_0: \beta_0 = 0\) versus \(H_1: \beta_0 \neq 0\), considere \(\alpha = 0,05\).
summary(modelo)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-734.46 -336.71 82.71 188.19 900.28
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2349.295 796.567 -2.949 0.025637 *
x 8.454 1.261 6.702 0.000536 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 534 on 6 degrees of freedom
Multiple R-squared: 0.8822, Adjusted R-squared: 0.8625
F-statistic: 44.92 on 1 and 6 DF, p-value: 0.0005355- Realize o teste de hipóteses \(H_0: \beta_1 = 0\) versus \(H_1: \beta_1 \neq 0\), considere \(\alpha = 0,05\).
summary(modelo)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-734.46 -336.71 82.71 188.19 900.28
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2349.295 796.567 -2.949 0.025637 *
x 8.454 1.261 6.702 0.000536 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 534 on 6 degrees of freedom
Multiple R-squared: 0.8822, Adjusted R-squared: 0.8625
F-statistic: 44.92 on 1 and 6 DF, p-value: 0.0005355- Realize o teste de hipóteses \(H_0: \beta_1 = 8\) versus \(H_1: \beta_1 \neq 8\), considere \(\alpha = 0,05\).
beta1 <- 8.454
beta1_h0 <- 8
ep_beta1 <- 1.261
tc <- (beta1 - beta1_h0) / ep_beta1
tc # t calculado[1] 0.3600317ttab <- qt(0.975, df = 6)
ttab # t tabelado ou valor crítico.[1] 2.446912abs(tc) > ttab # qual é a decisão? Rej ou nao H0?[1] FALSE- Calcule o intervalo de confiança de \(95\%\) para o coeficiente linear (\(\beta_1\)) e interprete-o.
confint(modelo, level = 0.95) 2.5 % 97.5 %
(Intercept) -4298.424521 -400.16580
x 5.367883 11.54098- Apresente os valores estimados de \(y\) para todos os valores de \(x\).
y_hat <- predict(modelo)
y_hat 1 2 3 4 5 6 7 8
1100.114 1100.114 2334.461 2334.461 3399.720 3399.720 4515.705 4515.705 data.frame(x, y, y_hat) x y y_hat
1 408 1100 1100.114
2 408 1300 1100.114
3 554 1600 2334.461
4 554 2500 2334.461
5 680 3000 3399.720
6 680 4300 3399.720
7 812 4200 4515.705
8 812 4700 4515.705- Calcule o intervalo de confiança para \(y\), quando \(x = 680\), ao nível de \(95\%\) de confiança. Ou seja, \(IC_{95\%}(y(680))\), sendo \(y(680) = E(Y| x = 680)\).
x_conf <- data.frame(x = 680)
predict(modelo, newdata = x_conf, interval = "confidence", level = 0.95) fit lwr upr
1 3399.72 2894.232 3905.208- Calcule os intervalos de confiança para \(y\), para todos os valores de \(x\), ao nível de \(95\%\) de confiança.
x_conf <- data.frame(x = x)
y_hat_conf <- predict(modelo, newdata = x_conf,
interval = "confidence", level = 0.95)
# data.frame com valores de x, y, predito e IC:
tabela <- data.frame(x, y, y_hat_conf)
tabela x y fit lwr upr
1 408 1100 1100.114 315.4437 1884.784
2 408 1300 1100.114 315.4437 1884.784
3 554 1600 2334.461 1837.3540 2831.568
4 554 2500 2334.461 1837.3540 2831.568
5 680 3000 3399.720 2894.2318 3905.208
6 680 4300 3399.720 2894.2318 3905.208
7 812 4200 4515.705 3748.3945 5283.016
8 812 4700 4515.705 3748.3945 5283.016# alterar o nome dos colunas:
colnames(tabela) <- c("x", "y", "y_hat", "LI", "LS")
tabela x y y_hat LI LS
1 408 1100 1100.114 315.4437 1884.784
2 408 1300 1100.114 315.4437 1884.784
3 554 1600 2334.461 1837.3540 2831.568
4 554 2500 2334.461 1837.3540 2831.568
5 680 3000 3399.720 2894.2318 3905.208
6 680 4300 3399.720 2894.2318 3905.208
7 812 4200 4515.705 3748.3945 5283.016
8 812 4700 4515.705 3748.3945 5283.016- Construa um gráfico com os pontos observados, a reta de regressão ajustada e os intervalos de confiança para \(y\).
plot(x,y)
abline(modelo)
lines(tabela$x, tabela$LI, lty = 2)
lines(tabela$x, tabela$LS, lty = 2)
plot(x, y,
pch = 19, col = "black",
xlab = "Concentração de CO2 (ppm)",
ylab = "Massa (g)",
main = "Regressão Linear Simples com IC para y")
# adicionar a reta de regressão:
abline(modelo, col = "red", lwd = 2)
# adicionar os intervalos de confiança:
lines(tabela$x, tabela$LI, col = "brown", lty = 2)
lines(tabela$x, tabela$LS, col = "brown", lty = 2)
legend("topleft",
legend = c("Observações", "Reta de Regressão", "IC para y"),
col = c("black", "red", "brown"),
pch = c(19, NA, NA),
lty = c(NA, 1, 2),
lwd = c(NA, 2, 1),
bty = "n")
# Usando o pacote ggplot2
library(ggplot2)
ggplot(tabela, aes(x = x, y = y)) +
geom_point(color = "black") +
geom_line(aes(y = y_hat), color = "red", linewidth = 1) +
geom_line(aes(y = LI), color = "brown", linewidth = 1, lty = 2) +
geom_line(aes(y = LS), color = "brown", linewidth = 1, lty = 2) +
labs(title = "Regressão Linear Simples com IC para y",
x = "Concentração de CO2 (ppm)",
y = "Massa (g)")
- Calcule a estimativa de \(y\) quando \(x = 600\). Ou seja, \(y(600)\).
x0 <- data.frame(x = 600)
y_hat_x0 <- predict(modelo, newdata = x0)
y_hat_x0 1
2723.365 - Calcule as estimativas de \(y\) quando \(x = 500, 600\) e \(700\).
x0 <- data.frame(x = c(500, 600, 700))
y_hat_x0 <- predict(modelo, newdata = x0)
y_hat_x0 1 2 3
1877.922 2723.365 3568.809 tabela_x0 <- data.frame(x = x0, y_hat = y_hat_x0)
tabela_x0 x y_hat
1 500 1877.922
2 600 2723.365
3 700 3568.809- Calcule as estimativas de \(y\) quando \(x = 500, 600\) e \(700\) e os intervalos de predição para cada um desses valores de \(x\), ao nível de \(95\%\) de confiança.
x0 <- data.frame(x = c(500, 600, 700))
y_hat_x0 <- predict(modelo, newdata = x0,
interval = "prediction", level = 0.95)
y_hat_x0 fit lwr upr
1 1877.922 448.5088 3307.335
2 2723.365 1336.9196 4109.811
3 3568.809 2157.5053 4980.112tabela_x0 <- data.frame(x = x0, y_hat = y_hat_x0)
tabela_x0 x y_hat.fit y_hat.lwr y_hat.upr
1 500 1877.922 448.5088 3307.335
2 600 2723.365 1336.9196 4109.811
3 700 3568.809 2157.5053 4980.112# alterando os nomes das colunas:
colnames(tabela_x0) <- c("x", "y_hat", "LI_pred", "LS_pred")
tabela_x0 x y_hat LI_pred LS_pred
1 500 1877.922 448.5088 3307.335
2 600 2723.365 1336.9196 4109.811
3 700 3568.809 2157.5053 4980.112- Construa o gráfico com os intervalos de confiança e predição, ao nível de \(95\%\) de confiança, para os valores de \(x\) entre \(400\) e \(900\).
x_seq <- seq(400, 900, by = 10)
x_seq_df <- data.frame(x = x_seq)
y_hat_seq <- predict(modelo, newdata = x_seq_df,
interval = "confidence", level = 0.95)
y_pred_seq <- predict(modelo, newdata = x_seq_df,
interval = "prediction", level = 0.95)
tabela_seq <- data.frame(x = x_seq, y_hat = y_hat_seq[, "fit"],
LI_conf = y_hat_seq[, "lwr"], LS_conf = y_hat_seq[, "upr"],
LI_pred = y_pred_seq[, "lwr"], LS_pred = y_pred_seq[, "upr"])
head(tabela_seq) x y_hat LI_conf LS_conf LI_pred LS_pred
1 400 1032.478 227.7168 1837.240 -502.04035 2566.997
2 410 1117.023 337.3341 1896.711 -404.49658 2638.542
3 420 1201.567 446.5220 1956.612 -307.47224 2710.606
4 430 1286.111 555.2373 2016.986 -210.98031 2783.203
5 440 1370.656 663.4314 2077.880 -115.03364 2856.345
6 450 1455.200 771.0502 2139.350 -19.64488 2930.045ggplot(tabela_seq, aes(x = x)) +
geom_line(aes(y = y_hat), color = "red", linewidth = 1) +
geom_line(aes(y = LI_conf), color = "blue", lty = 2) +
geom_line(aes(y = LS_conf), color = "blue", lty = 2) +
geom_line(aes(y = LI_pred), color = "green", lty = 2) +
geom_line(aes(y = LS_pred), color = "green", lty = 2) +
labs(title = "Reta de Regressão com IC e IP",
x = "Concentração de CO2 (ppm)",
y = "Massa (g)") +
theme_minimal()
# sugestão de outro gráfico (que também pode ser aprimorado):
ggplot(tabela_seq, aes(x = x)) +
geom_line(aes(y = y_hat, color = "Y = -2349.295 + 8.454x"), linewidth = 1) +
geom_line(aes(y = LI_conf, color = "IC 95%"), linetype = 2) +
geom_line(aes(y = LS_conf, color = "IC 95%"), linetype = 2) +
geom_line(aes(y = LI_pred, color = "IP 95%"), linetype = 2) +
geom_line(aes(y = LS_pred, color = "IP 95%"), linetype = 2) +
scale_color_manual(
values = c(
"Y = -2349.295 + 8.454x" = "black",
"IC 95%" = "black",
"IP 95%" = "grey"
)
) +
labs(
title = "",
x = "Concentração de CO_2 (ppm)",
y = "Massa (g)",
color = ""
) +
theme_minimal() +
theme(
legend.position = c(0.05, 0.95),
legend.justification = c(0, 1),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(color = "black")
)