Análise Descritiva

Decidimos trabalhar com as variáveis mpg (milhas por galão de combustível) e disp (volume do motor), do conjunto de dados mtcars, por conta da estrutura não linear entre elas, e da possibilidade de ajustar uma regressão segmentada aos dados (aparente mudança de comportamento em disp=200).

# Scatterplot
plot(mtcars$mpg ~ mtcars$disp, xlab = 'Volume do Motor', ylab = 'Milhas por Galão')

# Medidas Descritivas
kable(summary(mtcars[,c(3,5)]))

Regressão Linear

O primeiro modelo a ser testado é o mais comum, a regressão linear.

Equação do modelo:

\[y = \beta_0 + \beta_1x + \epsilon\] Onde:

\[E(y/x) = \beta_0 +\beta_1x\]

\[Var(y/x) = Var(\epsilon) = \sigma^2\]

Assim:

\[y/x \sim N(\beta_0+\beta_1x, \sigma^2)\]

# Ajuste 1 - Regressão Linear
ajuste1 <- lm(mpg ~ disp, data= mtcars)

preds <- data.frame(disp = seq(0,500, length=101))
preds$pred1 <- predict(ajuste1, newdata = preds)

plot(mpg~disp, data=mtcars, xlab = 'Volume do Motor', ylab = 'Milhas por Galão')
lines(pred1~disp, data=preds, col = 'steelblue', lwd = 2)

# Descrição do Ajuste
summary(ajuste1)

## 
## Call:
## lm(formula = mpg ~ disp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8922 -2.2022 -0.9631  1.6272  7.2305 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 29.599855   1.229720  24.070  < 2e-16 ***
## disp        -0.041215   0.004712  -8.747 9.38e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.251 on 30 degrees of freedom
## Multiple R-squared:  0.7183, Adjusted R-squared:  0.709 
## F-statistic: 76.51 on 1 and 30 DF,  p-value: 9.38e-10

# Análise de Resíduos
par(mfrow = c(2,2))
plot(ajuste1, 1:4)

Apesar da reta se ajustar aos dados de maneira razoável, os resíduos variam entre -6 e 8 e apresentam um padrão não aleatório. Além disso, o gráfico quantil-quantil indica uma certa fuga da normalidade na calda superior. Por fim, a distância de Cook (identificação de observações influentes) não é alarmante para nenhum ponto em específico (valores muito baixos).

Transformação Box-Cox

Para contornar os pontos indicados no tópico anterior, decidimos aplicar uma transformação na variável resposta que induzisse a normalidade (pressuposto do modelo de regressão linear). A transformação de Box-Cox baseia-se em encontrar o valor ótimo de \(\lambda\) (equação abaixo), de maneira que \(y^*\) tenha distribuição normal.

\[y^* = \frac{y^\lambda - 1}{\lambda}\] Principais transformações:

Como o intervalo indicado no gráfico inclui o valor zero, vamos ajustar um novo modelo aplicando uma transformação logarítmica à variável resposta.

# Ajuste 2 - Log(y)
ajuste2 <- lm(log(mpg) ~ disp, data= mtcars)

preds$pred2 <- predict(ajuste2, newdata = preds)

plot(log(mpg)~disp, data=mtcars, xlab = 'Volume do Motor', ylab = 'Log(Milhas por Galão)')
lines(pred2~disp, data=preds, col = 'steelblue', lwd = 2)

# Descrição do Ajuste
summary(ajuste2)

## 
## Call:
## lm(formula = log(mpg) ~ disp, data = mtcars)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.21183 -0.10837 -0.04732  0.08251  0.35546 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.445548   0.054290   63.47  < 2e-16 ***
## disp        -0.002115   0.000208  -10.17  3.1e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1435 on 30 degrees of freedom
## Multiple R-squared:  0.7751, Adjusted R-squared:  0.7676 
## F-statistic: 103.4 on 1 and 30 DF,  p-value: 3.095e-11

# Análise de Resíduos
par(mfrow = c(2,2))
plot(ajuste2, 1:4)

Apesar da escala dos resíduos diminuir de forma considerável, alguns pontos não estão sendo bem representados pelo ajuste. Por conta da estrutura dos dados, acreditamos que é necessário trabalhar com alternativas não lineares.

Polinômio de 2° Grau

Equação do modelo:

\[y = \beta_0 + \beta_1x + \beta_2x^2 + \epsilon\]

# Ajuste 3 - Polinômio de 2° Grau
ajuste3 <- lm(mpg ~ disp + I(disp^2), data= mtcars)

preds$pred3 <- predict(ajuste3, newdata = preds)

plot(mpg~disp, data=mtcars, xlab = 'Volume do Motor', ylab = 'Milhas por Galão')
lines(pred3~disp, data=preds, col = 'steelblue', lwd = 2)

# Descrição do Ajuste
summary(ajuste3)

## 
## Call:
## lm(formula = mpg ~ disp + I(disp^2), data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9112 -1.5269 -0.3124  1.3489  5.3946 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.583e+01  2.209e+00  16.221 4.39e-16 ***
## disp        -1.053e-01  2.028e-02  -5.192 1.49e-05 ***
## I(disp^2)    1.255e-04  3.891e-05   3.226   0.0031 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.837 on 29 degrees of freedom
## Multiple R-squared:  0.7927, Adjusted R-squared:  0.7784 
## F-statistic: 55.46 on 2 and 29 DF,  p-value: 1.229e-10

# Análise de Resíduos
par(mfrow = c(2,2))
plot(ajuste3, 1:4)

Visualmente, o ajuste do modelo com a adição de um polinômio de segundo grau parece representar os dados com maior precisão. Como esperado, como estamos utilizando a escala original da variável resposta, com a intenção de facilitar a interpretação do modelo, os problemas com os resíduos voltam a aparecer.

Spline Cúbico

Ainda no intuito de capturar a estrutura não linear entre as variáveis, decidimos ajustar um spline cúbico aos dados. Na função smooth.spline, definimos os graus de liberdade para controlar o número de nós. Em resumo, graus de liberdade mais baixos resultam em um ajuste mais próximo de uma regressão linear, enquanto um número maior de graus de liberdade proporciona uma suavização maior.

# Ajuste 4 - Spline Cúbico
ajuste4 <- smooth.spline(mtcars$disp, mtcars$mpg, df = 5)

plot(mtcars$disp, mtcars$mpg,
     xlab = 'Volume do Motor', ylab = 'Milhas por Galão')

lines(ajuste4, col = 'steelblue', lwd = 2)

Visualmente, o modelo se ajusta muito bem aos dados.

Regressão Segmentada

Como os dados possuem uma certa mudança de comportamento em disp=200 (comportamento linear antes do ponto, não linear depois do ponto), vamos ajustar uma regressão segmentada da seguinte forma:

\[ y = \begin{cases} \beta_0 + \beta_1 x & \text{se } x < 200 \\ \beta_0 + \beta_2 x + \beta_3( x -200)^2 + \beta_4 (x-200)^3 & \text{se } x \geq 200 \end{cases} \]

# Ajuste 5 - Regressão Segmentada

x <- mtcars$disp
y <- mtcars$mpg

dados <- data.frame(x = x, y = y)

ajuste_seg <- lm(y ~ I(x * (x < 200)) + I((x) * (x >= 200)) + I((x - 200)^2 * (x >= 200)) + I((x - 200)^3 * (x >= 200)), data = dados)

plot(mtcars$disp, mtcars$mpg,
     xlab = 'Volume do Motor', ylab = 'Milhas por Galão')

pred <- data.frame(x = seq(0, 500, length = 1000))

pred$y_pred <- predict(ajuste_seg, newdata = pred)

pred$x[pred$x < 200 & pred$x >= 185] <- NA
pred$y_pred[pred$x < 200 & pred$x >= 185] <- NA

lines(pred$x, pred$y_pred, col = 'steelblue', lwd = 2)

KNN

O KNN (K-Nearest Neighbors) é um modelo de regressão não paramétrico, que estima o valor de cada uma das observações com base na média das k observações mais próximas delas. Utilizando o pacote FNN, com k=2 já encontramos um ajuste satisfatório, e decidimos seguir assim (para evitar o overfitting do modelo).

# Ajuste 6 - KNN
ajuste_knn <- FNN::knn.reg(train = mtcars$disp, y = mtcars$mpg, k = 2)

plot(mtcars$disp, mtcars$mpg,
     xlab = 'Volume do Motor', ylab = 'Milhas por Galão')


sorted_disp <- sort(mtcars$disp)
sorted_preds <- ajuste_knn$pred[order(mtcars$disp)]

lines(sorted_disp, sorted_preds, col = 'steelblue', lwd = 2)

Comparação dos Ajustes

# Regressão Linear x Polinômio de 2º Grau
logLik(ajuste1)

## 'log Lik.' -82.10469 (df=3)

logLik(ajuste3)

## 'log Lik.' -77.19785 (df=4)

# Regressão Linear x Polinômio de 2º Grau
AIC(ajuste1)

## [1] 170.2094

AIC(ajuste3)

## [1] 162.3957

# Regressão Linear x Transformação Box-Cox x Polinômio de 2º Grau
summary(ajuste1)$r.squared

## [1] 0.7183433

summary(ajuste2)$r.squared

## [1] 0.7751119

summary(ajuste3)$r.squared

## [1] 0.7927323

func_rmse <- function(modelo, log = FALSE, spline = FALSE, knn = FALSE) {
    if (knn == TRUE) {
        # Para modelo KNN ajustado 
        pred <- modelo$pred
        mse <- mean((mtcars$mpg - pred)^2)
        rmse <- sqrt(mse)
        
    } else if (spline == TRUE) {
        # Caso o modelo seja um spline
        pred <- predict(modelo, mtcars$disp)$y
        mse <- mean((mtcars$mpg - pred)^2)
        rmse <- sqrt(mse)
        
    } else if (log == TRUE) { 
        # Caso seja na escala logarítmica
        pred_log <- predict(modelo, newdata = mtcars)
        pred <- exp(pred_log)
        mse <- mean((mtcars$mpg - pred)^2)
        rmse <- sqrt(mse)
        
    } else {
        # Predições na escala original
        pred <- predict(modelo, newdata = mtcars)
        mse <- mean((mtcars$mpg - pred)^2)
        rmse <- sqrt(mse)
    }
    
    return(rmse)
}

# RMSE dos Modelos
rmses <- c(
    func_rmse(ajuste1), 
    func_rmse(ajuste2, log = TRUE), 
    func_rmse(ajuste3),
    func_rmse(ajuste4, spline = TRUE), 
    func_rmse(ajuste_seg), 
    func_rmse(ajuste_knn, knn =TRUE)
)

# Nomes dos modelos
modelos <- c(
    'Regressão Linear', 
    'Transformação Box-Cox', 
    'Polinômio 2º Grau', 
    'Spline Cúbico', 
    'Regressão Segmentada', 
    'KNN'
)

# Tabela com RMSEs
knitr::kable(
data.frame('Modelos' = modelos, 'RMSE' = rmses))