Além da linearidade

Regressões polinomial e com splines

Author

Rodrigo M. R. de Medeiros

Pacotes necessários

library(gamlss.data) ## Pacote do conjunto de dados
library(tidyverse)   ## Manipulação e visualização de dados
library(splines)     ## Splines
theme_set(theme_classic())  ## Definição do tema dos gráficos ggplot

1 Conjunto de dados

## Conjunto de dados dbbmi
glimpse(dbbmi)

# Rows: 7,294
# Columns: 2
# $ age <dbl> 0.03, 0.04, 0.04, 0.04, 0.04, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, …
# $ bmi <dbl> 13.23529, 12.43877, 14.54177, 11.77395, 15.32561, 13.21439, 16.834…

1.1 Divisão entre conjunto de treinamento e validação

## Divisão aleatória (70% treinamento e 30% validação)
set.seed(123); id <- sample(1:nrow(dbbmi), 0.7 * nrow(dbbmi))
train <- dbbmi[id, ]
val <- dbbmi[-id, ]

1.2 Estatística descritiva no conjunto de treinamento

## Estatísticas descritivas
summary(dbbmi)

#       age              bmi       
#  Min.   : 0.030   Min.   :11.17  
#  1st Qu.: 1.863   1st Qu.:15.96  
#  Median :10.450   Median :17.45  
#  Mean   : 9.291   Mean   :18.03  
#  3rd Qu.:15.130   3rd Qu.:19.60  
#  Max.   :21.700   Max.   :35.42

## Histograma da idade
ggplot(train, aes(x = age)) +
  geom_histogram(fill = "white", col = 1, bins = 20) +
  labs(x = "Idade", y = "Frequência")

## Histograma do IMC
ggplot(train, aes(x = bmi)) +
  geom_histogram(fill = "white", col = 1, bins = 20) +
  labs(x = "IMC", y = "Frequência")

## Gráfico de dispersão
ggplot(train, aes(x = age, y = bmi)) +
  geom_point(col = "darkgray") +
  labs(x = "Idade", y = "IMC")

2 Regressão linear simples

O modelo de regressão linear simples explica o IMC do \(i\)-ésimo indivíduo (\(Y_i\)) em função da sua idade (\(x_i\)) como

\[ Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \]

em que \(\varepsilon_i\) é uma variável aleatória com média zero e variância constante.

## Ajuste do modelo de regressão linear simples
fit0 <- lm(bmi ~ age, train)

## Reta ajustada
data.frame(age = train$age, bmi = train$bmi, pred = predict(fit0)) %>% 
ggplot() +
  geom_point(aes(x = age, y = bmi), col = "darkgray") +
  geom_line(aes(x = age, y = pred), col = "blue", lwd = 1) +
  labs(x = "Idade", y = "IMC")

3 Regressão polinomial

O modelo de regressão polinomial para \(Y_i\) é definido por

\[ Y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \cdots + \beta_d x_i^d \varepsilon_i, \]

em que \(\varepsilon_i\) é uma variável aleatória com média zero e variância constante.

## Ajustes da regressão polinomial com grau 9 e 10
fitp9 <- lm(bmi ~ poly(age, 9, raw = TRUE), train)
fitp10 <- lm(bmi ~ poly(age, 10, raw = TRUE), train)

## Curva ajustada com d = 9
data.frame(age = train$age, bmi = train$bmi, pred = predict(fitp9)) %>% 
  ggplot() +
  geom_point(aes(x = age, y = bmi), col = "darkgray") +
  geom_line(aes(x = age, y = pred), col = "blue", lwd = 1) +
  labs(x = "Idade", y = "IMC")

## Curva ajustada com d = 10
data.frame(age = train$age, bmi = train$bmi, pred = predict(fitp10)) %>% 
  ggplot() +
  geom_point(aes(x = age, y = bmi), col = "darkgray") +
  geom_line(aes(x = age, y = pred), col = "blue", lwd = 1) +
  labs(x = "Idade", y = "IMC")

4 Regressão com splines

Sejam \(c_1 \leqslant c_2 \leqslant \cdots \leqslant c_K\) os nós. O modelo de regressão com splines para \(Y_i\) é definido por

\[ \begin{array}{ll}\begin{split}Y_i = & \; \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \cdots + \beta_d x_i^d + \\& \hspace{0.85cm} \; \beta_{d+1} g(x_i, c_1) + \beta_{d+2} g(x_i, c_2) \cdots + \beta_{d+K} g(x_i, c_K) + \varepsilon_i,\end{split}\end{array} \]

em que \(\varepsilon_i\) é uma variável aleatória com média zero e variância constante e

\[ g(x_i, c_j) = \begin{cases}(x_i - c_j)^d, & \text{se } x_i > c_j,\\0, & \text{caso contrário}. \end{cases} \]

## Nós
nos <- c(2, 5, 10, 18)

## Ajustes da regressão com cpline linear, quadrática e cúbica
fitsp1 <- lm(bmi ~ bs(age, knots = nos, degree = 1), train)
fitsp2 <- lm(bmi ~ bs(age, knots = nos, degree = 2), train)
fitsp3 <- lm(bmi ~ bs(age, knots = nos, degree = 3), train)

## Curva ajustada considerando o spline linear
data.frame(age = train$age, bmi = train$bmi, pred = predict(fitsp1)) %>% 
  ggplot() +
  geom_point(aes(x = age, y = bmi), col = "darkgray") +
  geom_vline(xintercept = nos, lty = 2) +
  geom_line(aes(x = age, y = pred), col = "blue", lwd = 1) +
  labs(x = "Idade", y = "IMC")

## Curva ajustada considerando o spline quadrático
data.frame(age = train$age, bmi = train$bmi, pred = predict(fitsp2)) %>% 
  ggplot() +
  geom_point(aes(x = age, y = bmi), col = "darkgray") +
  geom_vline(xintercept = nos, lty = 2) +
  geom_line(aes(x = age, y = pred), col = "blue", lwd = 1) +
  labs(x = "Idade", y = "IMC")

## Curva ajustada considerando o spline cúbico
data.frame(age = train$age, bmi = train$bmi, pred = predict(fitsp3)) %>% 
  ggplot() +
  geom_point(aes(x = age, y = bmi), col = "darkgray") +
  geom_vline(xintercept = nos, lty = 2) +
  geom_line(aes(x = age, y = pred), col = "blue", lwd = 1) +
  labs(x = "Idade", y = "IMC")

5 Comparação

## EQM da regressão linear simples no conjunto de validação
rmse0 <- mean((val$bmi - predict(fit0, newdata = val))^2)

## EQM da regressão polinomial no conjunto de validação
rmsep9 <- mean((val$bmi - predict(fitp9, newdata = val))^2)
rmsep10 <- mean((val$bmi - predict(fitp10, newdata = val))^2)

## EQM da regressão com splines no conjunto de validação
rmses1 <- mean((val$bmi - predict(fitsp1, newdata = val))^2)
rmses2 <- mean((val$bmi - predict(fitsp2, newdata = val))^2)
rmses3 <- mean((val$bmi - predict(fitsp3, newdata = val))^2)

## Tabela
tab <- data.frame(Modelo = c("RLS", paste0("Polinomial - ", 9:10),
                      "Spline linear", "Spline quadrático", "Spline cúbico"),
           EQM = round(c(rmse0, rmsep9, rmsep10,
                         rmses1, rmses2, rmses3), 3))

tab

#              Modelo   EQM
# 1               RLS 5.308
# 2    Polinomial - 9 4.350
# 3   Polinomial - 10 4.336
# 4     Spline linear 4.517
# 5 Spline quadrático 4.412
# 6     Spline cúbico 4.334

## Ajuste no conjunto de validação
data.frame(age = val$age, bmi = val$bmi, 
           predsp = predict(fitsp3, val),
           predpl = predict(fitp10, val)) %>% 
  pivot_longer(cols = 3:4, names_to = "model", values_to = "prediction") %>% 
  mutate(model = factor(model, labels = c("Polinomial (d = 10)",
                                          "Spline cúbico"))) %>% 
  ggplot() +
  geom_point(aes(x = age, y = bmi), col = "darkgray", cex = 1) +
  geom_vline(xintercept = nos, lty = 2) +
  geom_line(aes(x = age, y = prediction, col = model), lwd = 1) +
  scale_color_manual(values = c("red2", "blue")) +
  labs(x = "Idade", y = "IMC", col = "Modelo") +
  theme(legend.position = "top")