Prova - Não Paramétrica

---

This question uses the variables dis (the weighted mean of distances to five Boston employment centers) and nox (nitrogen oxides concen- tration in parts per 10 million) from the Boston data. We will treat dis as the predictor and nox as the response.

	crim	zn	indus	chas	nox	rm	age	dis	rad	tax	ptratio	lstat	medv
1	0.00632	18.0	2.31	0	0.538	6.575	65.2	4.0900	1	296	15.3	4.98	24.0
2	0.02731	0.0	7.07	0	0.469	6.421	78.9	4.9671	2	242	17.8	9.14	21.6
3	0.02729	0.0	7.07	0	0.469	7.185	61.1	4.9671	2	242	17.8	4.03	34.7
4	0.03237	0.0	2.18	0	0.458	6.998	45.8	6.0622	3	222	18.7	2.94	33.4
5	0.06905	0.0	2.18	0	0.458	7.147	54.2	6.0622	3	222	18.7	5.33	36.2
6	0.02985	0.0	2.18	0	0.458	6.430	58.7	6.0622	3	222	18.7	5.21	28.7
7	0.08829	12.5	7.87	0	0.524	6.012	66.6	5.5605	5	311	15.2	12.43	22.9
8	0.14455	12.5	7.87	0	0.524	6.172	96.1	5.9505	5	311	15.2	19.15	27.1
9	0.21124	12.5	7.87	0	0.524	5.631	100.0	6.0821	5	311	15.2	29.93	16.5
10	0.17004	12.5	7.87	0	0.524	6.004	85.9	6.5921	5	311	15.2	17.10	18.9

(a) Use the poly() function to fit a cubic polynomial regression to predict nox using dis. Report the regression output, and plot the resulting data and polynomial fits.

modelo <- lm(nox ~ poly(dis, 3))
summary <- summary(modelo)

predicoes <- data.frame(dis = seq(min(dis), max(dis), length.out = 100))
predicoes$pred <- predict(modelo, newdata = predicoes, se.fit = TRUE)$fit
predicoes$intervalo <- 2 * predict(modelo, newdata = predicoes, se.fit = TRUE)$se.fit


ggplot() +
  geom_point(data = data.frame(dis = dis, nox = nox), aes(x = dis, y = nox), color = "gray15", size = 2, shape = 16) +
  geom_line(data = predicoes, aes(x = dis, y = pred), color = "red", size = 1.5) +
  geom_ribbon(data = predicoes, aes(x = dis, ymin = pred - intervalo, ymax = pred + intervalo), fill = "red", alpha = 0.4) +
  labs(x = "Distances", y = "NOX", title = "Regressão Polinomial Cúbica") + theme_minimal() + theme(plot.title = element_text(hjust = 0.5))

O modelo de regressão polinomial tem como objetivo substituir um modelo linear \[y_i = \beta_o+\beta_1x_i+\epsilon_i\] por um modelo que permita produzir uma curva não linear, através de uma função polinomial: \[y_i = \beta_o+\beta_1x_i+\beta_2x_i^2+...+\beta_kx_i^k+\epsilon_i\].

summary

## 
## Call:
## lm(formula = nox ~ poly(dis, 3))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.121130 -0.040619 -0.009738  0.023385  0.194904 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.554695   0.002759 201.021  < 2e-16 ***
## poly(dis, 3)1 -2.003096   0.062071 -32.271  < 2e-16 ***
## poly(dis, 3)2  0.856330   0.062071  13.796  < 2e-16 ***
## poly(dis, 3)3 -0.318049   0.062071  -5.124 4.27e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06207 on 502 degrees of freedom
## Multiple R-squared:  0.7148, Adjusted R-squared:  0.7131 
## F-statistic: 419.3 on 3 and 502 DF,  p-value: < 2.2e-16

Nesse sentido, produzindo o modelo com a variável resposta nox e a variável preditora dis elevada na terceira potência, temos como output da função summary que todas as potências da variável dis são significativas. Além disso, temos um \(R^2\) de 0.7147737, que podemos interpretar como a variação de nox que é explicada por este modelo.

(b) Plot the polynomial fits for a range of different polynomial degrees (say, from 1 to 10), and report the associated residual sum of squares.

graus <- 1:10

rss_valores <- numeric(length(graus))

lista_graficos <- list()


for (i in seq_along(graus)) {
  
  modelo_graus <- lm(nox ~ poly(dis, graus[i]))
  predicoes$pred <- predict(modelo_graus, newdata = predicoes, se.fit = TRUE)$fit
  predicoes$intervalo <- 2 * predict(modelo_graus, newdata = predicoes, se.fit = TRUE)$se.fit

  lista_graficos[[i]] <- ggplot() +
    geom_point(data = data.frame(dis = dis, nox = nox), aes(x = dis, y = nox), color = "gray15", size = 2, shape = 16) +
    geom_line(data = predicoes, aes(x = dis, y = pred), color = "red", size = 1.5) +
    geom_ribbon(data = predicoes, aes(x = dis, ymin = pred - intervalo, ymax = pred + intervalo), fill = "red", alpha = 0.4) +
    labs(x = "Distances", y = "NOX", title = paste("Regressão Polinomial", graus[i])) +
    theme_minimal() + 
    theme(plot.title = element_text(hjust = 0.5))
  
  rss_valores[i] <- sum(modelo_graus$residuals^2)
  
  
}

for (j in seq(1, length(lista_graficos), 2)){
  
  grid.arrange(grobs = lista_graficos[j:(j+1)], ncol = 2,  bottom = "_________________________________________________________________________", top = "_________________________________________________________________________")
  
}

Fazendo o mesmo procedimento da questão anterior, obtemos os gráficos das curvas para cada potência de dis. Vale reparar que, conforme aumentamos o grau do polinômio, o comportamento da curva fica um pouco esquisito, principalmente nas caudas, o que também acontece com o seus respectivos intervalos de confiança.

dados_rss <- data.frame(Grau_polinômio = graus, Soma_quadrados = rss_valores)

kable(dados_rss)

Grau_polinômio	Soma_quadrados
1	2.768563
2	2.035262
3	1.934107
4	1.932981
5	1.915290
6	1.878257
7	1.849484
8	1.835630
9	1.833331
10	1.832171

Para a tabela de soma residual dos quadrados, temos que ela diminui conforme aumentamos o grau do polinômio mas, apesar de termos essa medida, é válido utilizar técnicas de validação cruzada para avaliarmos de fato o desempenho deste modelo.

(c) Perform cross-validation or another approach to select the optimal degree for the polynomial, and explain your results.

cv.error <- rep(0, 10)

for (i in 1:10) {
  modelo2 <- glm(nox ~ poly(dis, graus[i]))
  cv.error[i] <- cv.glm(dados, modelo2)$delta[1]
}

df <- data.frame(graus = graus, cv.error = cv.error)

ggplot(data = df, aes(x = graus, y = cv.error)) +
  geom_line(color = "blue") +
  geom_point(shape = 15, size = 3, col = "steelblue") +
  geom_text(aes(label = round(cv.error, 5)), vjust = 1.4, size = 3, color = "black") +
  labs(x = "Grau do Polinômio", y = "Estimativa do erro",
       title = "Erro da Validação Cruzada vs. Grau Polinomial da Regressão") +
  scale_x_continuous(breaks = 1:10, labels = 1:10)+
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5))

Considerando o método de Leave-One-Out Cross-Validation (LOOCV), podemos observar pelo gráfico que a regressão polinomial cúbica é a que possui a menor estimativa do erro de predição, seguida pelas regressões de grau 4 e 2.

O método usado acima funciona de forma a separar todos pontos dos dados e treinar o modelo incluindo todos estes pontos (exceto um). Comparando com os outros métodos, ele possui uma certa vantagem pois reduz o viés, mas, por outro lado, é um procedimento com a variância maior e mais caro computacionalmente.

(d) Use the bs() function to fit a regression spline to predict nox using dis. Report the output for the fit using four degrees of freedom. How did you choose the knots? Plot the resulting fit.

modelo3 <- lm(nox ~ bs(dis, df = 4))
pred <- data.frame(dis = dis, fit = predict(modelo3, newdata = dados, se = TRUE)$fit, se = predict(modelo3, newdata = dados, se = TRUE)$se.fit)
summary3 <- summary(modelo3)

ggplot(data = pred, aes(x = dis, y = nox)) +
  geom_point(color = "gray15", size = 2, shape = 16) +
  geom_line(aes(y = fit), color = "red", size = 1.5) +
  geom_ribbon(aes(ymin = fit - 2 * se, ymax = fit + 2 * se), fill = "red", alpha = 0.4) +
  labs(x = "Distance", y = "NOX", title = "Regressão com bases B-Spline")  + theme_minimal() + theme(plot.title = element_text(hjust = 0.5))

A motivação de fazer um modelo de Basis Functions é substituir o modelo linear por um modelo que contenha uma família de funções da variável X. Assim, temos: \[y_i = \beta_o+\beta_1b_1(x_i)+\beta_2b_2(x_i)+...+\beta_kb_k(x_i)+\epsilon_i\]

summary3

## 
## Call:
## lm(formula = nox ~ bs(dis, df = 4))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.124622 -0.039259 -0.008514  0.020850  0.193891 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.73447    0.01460  50.306  < 2e-16 ***
## bs(dis, df = 4)1 -0.05810    0.02186  -2.658  0.00812 ** 
## bs(dis, df = 4)2 -0.46356    0.02366 -19.596  < 2e-16 ***
## bs(dis, df = 4)3 -0.19979    0.04311  -4.634 4.58e-06 ***
## bs(dis, df = 4)4 -0.38881    0.04551  -8.544  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06195 on 501 degrees of freedom
## Multiple R-squared:  0.7164, Adjusted R-squared:  0.7142 
## F-statistic: 316.5 on 4 and 501 DF,  p-value: < 2.2e-16

Nesse sentido, produzindo o modelo com a variável resposta nox e a variável preditora dis com quatro graus de liberdade, temos como output da função summary que todos os fatores do modelo são significativos. Além disso, temos um \(R^2\) de \(0.7164449\), que podemos interpretar como a variação de nox que é explicada por este modelo.

Quanto aos nós, ao definirmos os graus de liberdade como 4, o R automaticamente escolhe 3 nós, que são localizados no primeiro, segundo e terceiro quartil de dis, que correspondem aos valores de \(2.1002, 3.2074, 5.1884\).

(e) Now fit a regression spline for a range of degrees of freedom, and plot the resulting fits and report the resulting RSS. Describe the results obtained.

graus_lib <- 1:10

rss_valores2 <- numeric(length(graus_lib))

lista_graficos2 <- list()

for (i in seq_along(graus_lib)) {
  modelo_spline <- lm(nox ~ bs(dis, df = graus_lib[i]))
  predicoes_spline <- data.frame(dis = dis, fit = predict(modelo_spline, newdata = dados, se = TRUE)$fit, se = predict(modelo_spline, newdata = dados, se = TRUE)$se.fit)
  
  rss_valores2[i] <- sum(modelo_spline$residuals^2)
  
  lista_graficos2[[i]] <- ggplot() +
  geom_point(data = data.frame(dis = dis, nox = nox), aes(x = dis, y = nox), color = "gray15", size = 2, shape = 16) +
  geom_line(data = predicoes_spline, aes(y = fit, x = dis), color = "red", size = 1.5) +
  geom_ribbon(data = predicoes_spline, aes(ymin = fit - 2 * se, ymax = fit + 2 * se, x = dis), fill = "red", alpha = 0.4) +
  labs(x = "Distances", y = "NOX", title = paste("Regression Spline (df =", graus_lib[i], ")")) + 
    theme_minimal() + 
    theme(plot.title = element_text(hjust = 0.5))

}

for (j in seq(1, length(lista_graficos2), 2)){
  
  grid.arrange(grobs = lista_graficos2[j:(j+1)], ncol = 2,  bottom = "_________________________________________________________________________", top = "_________________________________________________________________________")
  
}

Assim como na sequência de plots para os modelos de regressão polinomial da questão b) percebemos que conforme aumentamos o grau de liberdade a curva assume um formato esquisito. Vale ressaltar que o R considerou os graus 1 e 2 muito pequenos e os substituiu pelo grau 3, por isso os 3 primeiros modelos produziram resultados iguais.

dados_rss2 <- data.frame(Grau_de_liberdade = graus_lib, Soma_quadrados = rss_valores2)

kable(dados_rss2)

Grau_de_liberdade	Soma_quadrados
1	1.934107
2	1.934107
3	1.934107
4	1.922775
5	1.840173
6	1.833966
7	1.829884
8	1.816995
9	1.825652
10	1.792535

Novamente, conforme diminuímos o grau de liberdade, obtemos uma soma residual dos quadrados menor.

(f) Perform cross-validation or another approach in order to select the best degrees of freedom for a regression spline on this data. Describe your results.

cv.error2 <- rep(0, 10)

for (i in 1:10) {
  modelo4 <- glm(nox ~ bs(dis, df = graus_lib[i]))
  cv.error2[i] <- cv.glm(dados, modelo4, K = 10)$delta[1]
}

df2 <- data.frame(graus = graus, cv.error = cv.error2)

ggplot(data = df2, aes(x = graus, y = cv.error)) +
  geom_line(color = "blue") +
  geom_point(shape = 15, size = 3, col = "steelblue") +
  geom_text(aes(label = round(cv.error, 5)), vjust = 1.4, size = 3, color = "black") +
  labs(x = "Grau de liberdade no modelo", y = "Estimativa do erro",
       title = "Erro da Validação Cruzada vs. Grau de liberdade") +
  scale_x_continuous(breaks = 1:10, labels = 1:10)+
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5))

Aqui foi usado o método k-Fold Cross-Validation, que produziu o resultado apresentado pelo gráfico. Podemos perceber que o menor erro foi no modelo com grau de liberdade igual à 10.

O método usado acima funciona de forma a separar k grupos de observações de tamanho igual. O primeiro grupo é usado como um conjunto validador e o método se ajusta para os demais grupos, processo que é repetido k vezes, isolando cada grupo uma vez. Comparando com os outros métodos, ele possui uma certa vantagem pois reduz a variância, mas é inferior por ser mais enviesado, já que separa as observações em grupos maiores do que o método LOOCV.