Introdução

A regressão linear é um dos métodos mais utilizados em modelagem estatística. Ela parte do pressuposto de que existe uma relação linear entre as variáveis explicativas e a variável resposta, o que implica que a relação entre elas pode ser explicada por uma reta. Em termos práticos, espera-se uma relação aditiva e constante entre as variáveis.

No entanto, essa premissa frequentemente é quebrada, principalmente em situações onde os dados apresentam comportamentos não lineares. Um exemplo claro disso pode ser observado ao estudar a relação entre as variáveis medv (valor mediano das propriedades) e indus (proporção de acres destinados a negócios não varejistas) no famoso Boston Housing Dataset. Sabe-se de antemão que essa relação não é linear, o que pode levar a falhas significativas na modelagem.

Problema da Linearidade

Quando a linearidade dos dados não é alcançada, ocorrem problemas na distribuição dos resíduos, o que indica que o modelo não está capturando a verdadeira relação entre as variáveis. Isso pode gerar previsões enviesadas e imprecisas, onde o modelo tende a subestimar ou superestimar a variável resposta em diferentes regiões dos dados. Além disso, a falha em capturar relações não lineares pode resultar em uma modelagem inadequada para diferentes subgrupos dos dados.

Soluções para o Problema

Neste projeto, o objetivo é identificar essa falta de linearidade entre as variáveis medv e indus, demonstrando como é possível detectar esse problema e aplicar soluções. Entre as estratégias, podemos usar técnicas como transformações de variáveis, modelos polinomiais e splines para ajustar a relação de forma mais precisa e minimizar o erro nas previsões.

“Não existe modelo de regressão perfeito, existe o péssimo, o ruim e o menos pior.”

Essa frase nos lembra que, em análise de dados, a busca por um modelo ideal pode ser ilusória. O objetivo é encontrar uma solução que ofereça o melhor ajuste possível, minimizando os erros e as distorções.

Análise Exploratória Inicial

O primeiro passo na análise é verificar algumas informações fundamentais sobre as variáveis envolvidas no modelo. Começamos investigando as correlações (Pearson e Spearman) entre as variáveis indus e medv, com o objetivo de identificar se há uma relação linear ou não linear entre elas. A correlação de Pearson foi de 0.48 e a correlação de Spearman foi de 0.59. Esses valores indicam que existe uma correlação moderada entre as variáveis, sugerindo que, embora a relação não seja perfeitamente linear, há espaço para ajustes que poderiam melhorar a linearidade.

Distribuição da Variável medv

Além das correlações, também investigamos a distribuição da variável resposta medv. Embora a distribuição não impacte diretamente a linearidade da relação entre indus e medv, pode indicar problemas como heterocedasticidade (variação não constante dos resíduos). Nesse caso, observamos uma alta assimetria à direita na distribuição de medv.

Essa assimetria não afeta necessariamente a modelagem linear diretamente, mas pode sugerir que o erro do modelo não será distribuído uniformemente ao longo dos valores preditos, o que pode levar a previsões enviesadas ou inconsistências nos resíduos. Isso é algo que devemos ter em mente ao avançar com a modelagem e buscar soluções para melhorar o ajuste.

#Casas mais caras estão sitiadas em regiões de alto comportamento industrial?

#Primeiro vamos entender a correlação e a distribuição de Y

y = Boston$medv #variável dependente Medv - preço médio das casas
x = Boston$indus #variável independente Indus - área industrial da cidade


#Correlações
pearson_corr <- round(cor(y, x, method = 'pearson'), 2)
spearman_corr <- round(cor(y, x, method = 'spearman'), 2)

cat("Correlação Pearson entre X e Y:", pearson_corr, "\n")
## Correlação Pearson entre X e Y: -0.48
cat("Correlação Spearman entre X e Y:", spearman_corr, "\n")
## Correlação Spearman entre X e Y: -0.58
#Dimensão de Y
cat("Comprimento de y:", length(y), "\n")
## Comprimento de y: 506
# Visualizando a distribuição de y
library(e1071)
hist(y, breaks = 20, freq = FALSE, main = "Histograma da Distribuição de Y", xlab = "Y")
lines(density(y), col = "blue", lwd = 2)

# Teste de normalidade e estatísticas de assimetria e curtose
shapiro_result <- shapiro.test(y)
skewness_value <- skewness(y)
kurtosis_value <- kurtosis(y)

# Exibindo os resultados
cat("Resultado do Teste de Shapiro-Wilk:", shapiro_result$p.value, "\n")
## Resultado do Teste de Shapiro-Wilk: 4.941386e-16
cat("Assimetria (Skewness):", skewness_value, "\n")
## Assimetria (Skewness): 1.101537
cat("Curtose (Kurtosis):", kurtosis_value, "\n")
## Curtose (Kurtosis): 1.450984
# Scatter plot com formatação
plot(Boston$indus, Boston$medv, pch = 19, col = "blue",
     xlab = 'Indus', ylab = 'MedV',
     main = 'Scatter Plot da relação entre \nMedian Property Value e Non-Retail Business Acres por Town')

Avaliação da Linearidade do Modelo

Ao avaliar o gráfico de resíduos versus valores ajustados (Residuals vs Fitted), podemos observar que a forma dos resíduos parece se aproximar de uma curva em U. Esse padrão indica que a relação entre as variáveis medv e indus não pode ser adequadamente explicada por uma reta linear, reforçando a ideia de que há uma não linearidade significativa na interação entre essas variáveis.

Além disso, identificamos a presença de outliers e uma distribuição assimétrica à direita nos resíduos, sugerindo valores elevados de assimetria (skewness) e curtose (kurtosis). Essa observação é consistente com a análise preliminar da distribuição da variável medv. Outro ponto importante é o gráfico Scale-Location, que demonstra um padrão nos resíduos. Isso indica que pode haver problemas relacionados à heterocedasticidade, o que reforça a ideia de que estamos interpretando incorretamente a relação entre medv e indus.

Embora o estudo tenha identificado esses problemas de linearidade, é importante notar que nossa análise está focada especificamente no primeiro aspecto — a falta de linearidade entre as variáveis. Questões como a presença de outliers e o comportamento dos resíduos em relação à homocedasticidade são importantes, mas estão fora do escopo desta investigação.

Outras Formas de Identificar a Não Linearidade

Existem várias maneiras de identificar a falta de linearidade em um modelo de regressão. Uma das minhas formas favoritas é através do gráfico de resíduos parciais, especialmente quando trabalhamos com modelos de regressão linear múltipla.

A fórmula para o resíduo parcial é dada por:

\[ \text{Resíduo Parcial} = \text{Resíduo} + \hat{\beta}_i X_i \]

Esse gráfico permite isolar o relacionamento entre uma variável independente (X_i) e a variável dependente (Y), levando em consideração todas as outras variáveis do modelo. No nosso caso, como estamos lidando apenas com uma variável explicativa (indus), ele nos ajuda a entender a contribuição individual dessa variável para o modelo.

Em termos práticos, um resíduo parcial pode ser interpretado como o valor resultante da combinação da previsão de uma variável isolada (X_i) com o resíduo da equação completa. Esse gráfico é particularmente útil para verificar se a relação entre uma variável específica e a resposta é linear, enquanto os efeitos das demais variáveis são controlados.

Outra abordagem útil é o uso do gráfico crPlot(modelo, variável = '') do pacote car no R. Esse gráfico visualiza a relação entre a variável de interesse e o resíduo ajustado, permitindo uma inspeção mais detalhada da linearidade.

Alternativas para Modelagem: Splines e Regressão Polinomial

Com base na nossa análise, comprovamos que não existe uma relação linear adequada entre medv e indus. Portanto, não podemos utilizar uma regressão linear simples para explicar essa relação de forma precisa.

Diante desse cenário, existem várias alternativas. Uma das abordagens mais comuns seria aplicar transformações matemáticas em X ou Y. No entanto, essas transformações raramente resolvem completamente o problema da não linearidade. Outra abordagem seria o uso de regressões polinomiais.

O Problema das Regressões Polinomiais

Embora uma regressão polinomial possa capturar a não linearidade, ela apresenta um problema significativo: esses modelos são construídos de forma global, ou seja, aplicam a mesma função polinomial a todos os dados. Isso pode ser problemático se houver diferentes regiões nos dados onde a variação da relação muda. Nesse caso, o polinômio pode ainda falhar em explicar a complexidade dos dados.

Splines: Uma Solução para Não Linearidade

Para lidar com essa limitação, podemos recorrer a uma abordagem que permite segmentar o modelo em regiões distintas ou nós (knots). Essa técnica é conhecida como Splines, e uma das mais comuns é a Piecewise Splines.

Piecewise Splines

A ideia por trás das Piecewise Splines é aplicar diferentes regressões (lineares, quadráticas ou cúbicas) em cada região dos dados, separadas por nós. A fórmula para um Piecewise Cubic Spline é dada por:

\[ y_i = \begin{cases} \beta_{01} + \beta_{11} x_i + \beta_{21} x_i^2 + \beta_{31} x_i^3 + \epsilon_i & \text{se} \ x_i < c \\ \beta_{02} + \beta_{12} x_i + \beta_{22} x_i^2 + \beta_{32} x_i^3 + \epsilon_i & \text{se} \ x_i \geq c \end{cases} \]

Onde c define o nó (ou a região de quebra no plano cartesiano).

Por convenção, utilizamos o Piecewise Cubic porque ele geralmente oferece melhor desempenho em relação ao polinômio quadrático, capturando de forma mais suave as transições entre as regiões.

Em termos práticos, aplicamos duas funções polinomiais diferentes aos dados. Uma função é usada para o subconjunto onde x_i < c e a outra para x_i ≥ c. Cada uma dessas funções é ajustada separadamente usando o método de mínimos quadrados, o que permite maior flexibilidade para capturar as variações em diferentes partes dos dados.

Restrições de Suavidade nos Nós

Para evitar que o modelo apresente descontinuidades nas transições entre as regiões, podemos impor restrições que garantam a suavidade nos nós. Isso significa que as funções polinomiais nas diferentes regiões devem se encontrar nos pontos de quebra. Podemos impor até três restrições nas Piecewise Splines:

  1. Continuidade: Garantir que os valores das duas funções polinomiais sejam os mesmos no nó. \[ f_1(x = c) = f_2(x = c) \]
  2. Primeira Derivada: Garantir que as inclinações das duas funções sejam as mesmas no nó. \[ f'_1(x = c) = f'_2(x = c) \]
  3. Segunda Derivada: Garantir que a curvatura das duas funções também seja a mesma no nó. \[ f''_1(x = c) = f''_2(x = c) \]

Essas restrições fazem com que o modelo seja mais suave e contínuo, evitando quebras abruptas entre as regiões.

Outras Abordagens: Natural Splines e Smoothing Splines

Além das Piecewise Splines, temos outras variações como as Natural Splines e Smoothing Splines.

Natural Splines

As Natural Splines são Piecewise Splines cúbicas, mas com restrições nas pontas (extremos). Isso é importante porque, nas Piecewise Splines, geralmente observamos variações muito grandes nas extremidades, o que pode gerar intervalos de confiança extensos e instabilidade nos modelos. As Natural Splines adicionam restrições para que a função spline seja quase linear nos extremos, melhorando a estabilidade.

Smoothing Splines

As Smoothing Splines são mais complexas e poderosas, mas não as abordaremos em detalhes aqui. Elas funcionam criando quebras nos dados por meio de uma penalização na função de custo, o que pode dificultar a interpretação dos intervalos de confiança e dos erros padrão. Embora sejam eficientes, prefiro evitá-las em favor de métodos mais simples e fáceis de interpretar, como as Piecewise Splines.

sort.x = sort(x)

cuts <- summary(x)[c(2,3,5)] #Em primeiro momento vamos usar a famosa convenção teórica de definir os nós nos quartis 25, 50 e 75 da distribuição de X


#Para construir Piecewise Splines usamos o bs() para alterar a forma de X -> uma funçao em X
#Knots definimos as quebras e degree os graus polinomiais a serem usados
spline.bs = lm(y ~ bs(x, knots = cuts, degree = 3))
pred.bs = predict(spline.bs, newdata = list(x = sort.x), se = TRUE)
se.bands.bs = cbind(pred.bs$fit + 2 * pred.bs$se.fit, 
                    pred.bs$fit - 2 * pred.bs$se.fit)

y.lab = 'Median Property Value'
x.lab = 'Non-retail business acres per town'


plot(x, y, cex.lab = 1.1, col="darkgrey", xlab = x.lab, ylab = y.lab, 
     main = "Cubic Spline", bty = 'l')
lines(sort.x, pred.bs$fit, lwd = 2, col = "red")
matlines(sort.x, se.bands.bs, lwd = 2, col = "red", lty = 3)

#Veja como a modelagem (o encaixe do modelo) parece um pouco melhor (visualmente)

Escolha dos Nós (Knots) e Graus de Liberdade

Você deve estar se perguntando: Onde devemos posicionar os nós (knots) e quantos devemos usar, certo?

Posicionando os Nós

A escolha dos nós é fundamental para capturar as variações na relação entre X e Y. Devemos posicioná-los onde acreditamos que a variação entre as variáveis começa a mudar significativamente, ou seja, em pontos onde há uma curva acentuada ou uma expansão na relação. Esses pontos representam mudanças estruturais no comportamento dos dados e, portanto, são locais ideais para aplicar os nós.

Definindo Graus de Liberdade em vez de Nós

Embora possamos definir manualmente onde posicionar os nós, existe uma abordagem mais moderna e sofisticada que permite ao algoritmo ajustar automaticamente a posição dos nós com base na complexidade do modelo que desejamos.

Em vez de escolher as posições exatas dos nós, podemos definir a quantidade de graus de liberdade (degrees of freedom, ou df) que queremos que o modelo utilize. A quantidade de nós será então determinada automaticamente pela fórmula:

\[ \text{Knots} = df - \text{graus} + 1 \]

O modelo aplicará os nós em pontos de mudança da variação, de acordo com o número de graus de liberdade especificado. Isso permite que o modelo se ajuste de maneira mais dinâmica às mudanças nos dados.

Trade-off entre Flexibilidade e Overfitting

Quanto maior o número de graus de liberdade utilizado, mais flexível será o modelo. No entanto, essa flexibilidade adicional também pode aumentar o risco de overfitting, onde o modelo se ajusta demais aos dados de treino, perdendo capacidade de generalização.

Uma maneira de encontrar o número ideal de graus de liberdade, especialmente em cenários de previsão, é avaliar o trade-off entre viés e variância. Podemos testar diferentes valores de df e observar como o modelo se comporta em termos de erro de previsão e estabilidade.

Aplicação em Piecewise e Natural Splines

Essa abordagem de definir graus de liberdade em vez de nós pode ser aplicada tanto em Piecewise Splines quanto em Natural Splines. Em ambos os casos, a definição dos df permite que o modelo ajuste automaticamente a posição dos nós de forma otimizada, capturando as variações nos dados sem a necessidade de intervenções manuais.

# Definir uma semente para reprodutibilidade
set.seed(123)

# Definir o número de folds para cross-validation
k <- 20

# Definir uma faixa de graus de liberdade a serem testados
df_values <- seq(1, 24, by = 2)  # Graus de liberdade de 3 a 30
mean_squared_errors_test <- numeric(length(df_values))  # Para armazenar os erros de teste
mean_squared_errors_train <- numeric(length(df_values))  # Para armazenar os erros de treinamento

# Realizar k-fold cross-validation para cada número de graus de liberdade
for (i in seq_along(df_values)) {
  df <- df_values[i]
  
  # Criar uma estrutura para os folds
  folds <- createFolds(Boston$medv, k = k, list = TRUE)
  
  fold_errors_test <- numeric(k)  # Armazenar os erros de cada fold para teste
  fold_errors_train <- numeric(k)  # Armazenar os erros de cada fold para treinamento
  
  # Loop sobre cada fold
  for (j in 1:k) {
    # Criar conjunto de treino e teste
    train_set <- Boston[-folds[[j]], ]
    test_set <- Boston[folds[[j]], ]
    
    # Ajustar o modelo spline com o número atual de graus de liberdade
    modelo <- lm(medv ~ bs(indus, df = df), data = train_set)
    
    # Fazer previsões no conjunto de teste
    predictions_test <- predict(modelo, newdata = test_set)
    predictions_train <- predict(modelo, newdata = train_set)
    
    # Calcular o erro quadrático médio (MSE) para este fold
    fold_errors_test[j] <- mean((test_set$medv - predictions_test)^2)
    fold_errors_train[j] <- mean((train_set$medv - predictions_train)^2)
  }
  
  # Calcular o erro médio para este número de graus de liberdade
  mean_squared_errors_test[i] <- mean(fold_errors_test)
  mean_squared_errors_train[i] <- mean(fold_errors_train)
}

# Criar um data frame para armazenar os resultados
results <- data.frame(
  DF = df_values,
  MSE_Test = mean_squared_errors_test,
  MSE_Train = mean_squared_errors_train
)


print(results)
##    DF MSE_Test MSE_Train
## 1   1 61.51387  61.02982
## 2   3 61.88965  61.03319
## 3   5 62.10340  60.70588
## 4   7 62.00885  59.61618
## 5   9 61.06686  58.88041
## 6  11 58.47850  54.63621
## 7  13 55.80849  53.15543
## 8  15 55.80756  52.72706
## 9  17 55.64548  51.97902
## 10 19 55.09152  51.80888
## 11 21 55.73629  52.04513
## 12 23 80.15632  50.78438
# Identificar os índices dos dois menores erros de teste
min_indices <- order(results$MSE_Test)[1:2] 

# Criar o gráfico
ggplot(results) +
  geom_line(aes(x = DF, y = MSE_Test, color = "Teste")) +
  geom_line(aes(x = DF, y = MSE_Train, color = "Treinamento")) +
  geom_point(aes(x = DF, y = MSE_Test, color = "Teste")) +
  geom_point(aes(x = DF, y = MSE_Train, color = "Treinamento")) +
  geom_point(data = results[min_indices, ], aes(x = DF, y = MSE_Test), color = "tomato3", size = 5) +  
  geom_text(data = results[min_indices, ], aes(x = DF, y = MSE_Test, label = paste("DF:", DF, "\nMSE:", round(MSE_Test, 2))), vjust = -1, color = "red", size=3) +  
  labs(title = "Erro Quadrático Médio por Graus de Liberdade",
       x = "Graus de Liberdade",
       y = "Erro Quadrático Médio (MSE)") +
  scale_color_manual(name = "Conjunto", values = c("Teste" = "red", "Treinamento" = "blue")) +
  theme_minimal()

Construindo o Modelo com Graus de Liberdade Ideais

Agora que entendemos a importância dos graus de liberdade e como eles impactam a flexibilidade do modelo, vamos construir o nosso modelo de Splines usando um número “ideal” de df. Vamos também explorar os pontos onde o modelo decidiu aplicar os nós.

Definindo o Modelo

Ao definir o modelo, podemos controlar diretamente os graus de liberdade, e o algoritmo ajustará os nós automaticamente com base nos pontos de maior variação em X.

Depois de construir o modelo, podemos verificar onde os nós foram posicionados usando a função attr(). Isso nos permite acessar os valores de X onde as quebras ocorreram.

# Ajustar o modelo de spline
#Optei com DF = 17 (focamos em deixar o menos flexível e com menor erro)
spline.bs <- lm(y ~ bs(x, df=17, degree=3))
pred.bs <- predict(spline.bs, newdata = list(x = sort.x), se = TRUE)
se.bands.bs <- cbind(pred.bs$fit + 2 * pred.bs$se.fit, 
                     pred.bs$fit - 2 * pred.bs$se.fit)

# Criar data frame para ggplot
data_plot <- data.frame(
  x = sort.x,
  y_fit = pred.bs$fit,
  lower_ci = se.bands.bs[, 2],
  upper_ci = se.bands.bs[, 1]
)

# Obter os nós do spline
knots <- attr(bs(x, df = 17, degree = 3), "knots")

# Calcular RMSE
rmse <- sqrt(mean((y - predict(spline.bs))^2))

# Criar o gráfico
ggplot(data_plot, aes(x = x)) +
  geom_point(aes(y = y), color = "darkgrey") +  
  geom_line(aes(y = y_fit), color = "red", size = 1) + 
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), fill = "red", alpha = 0.2) + 
  geom_vline(xintercept = knots, linetype = "dashed", color = "blue", size = 0.8) + 
  labs(
    title = "Cubic Spline with Confidence Interval",
    subtitle = paste("Modelo de Spline cúbico ajustado com DF = 17 RMSE =", round(rmse, 2)),
    x = x.lab,
    y = y.lab
  ) + theme_minimal()

Avaliação do Ajuste e Limitações

Após o ajuste com splines, conseguimos resolver o problema de linearidade que observamos anteriormente. O modelo agora reflete melhor a relação entre medv e indus.

No entanto, ainda enfrentamos desafios relacionados à heterocedasticidade e à presença de outliers. Esses problemas indicam que:

  • Heterocedasticidade: A variabilidade dos resíduos pode sugerir que nosso modelo ainda não está capturando todas as nuances da relação entre as variáveis. É possível que existam preditores adicionais que explicam essa variação.

  • Outliers: A presença de valores extremos sugere que os dados de medv podem estar mal interpretados, ou que precisamos incluir mais variáveis explicativas no modelo.

Contudo, não vamos abordar esses aspectos neste estudo, pois nosso foco principal foi na correção da linearidade entre medv e indus.

par(mfrow = c(2,2))
plot(spline.bs)

#podemos confirmar, com os gráficos que nos resta um problemas de heterocedastividade
#provável que por conta de alguma variável ou relação que precisamos capturar
#para explicar a ponta desse modelo (onde a variância aumenta)

Vamos, por fim, apenas verificar a distribuição dos resíduos e podemos notar uma distribuição provável quasipoisson ou binomial negativa, talvez usar um glm() da familia dessas com o splines fosse o ideal no nosso cenário

# Calcular os resíduos do modelo
residuals <- residuals(spline.bs)

# Calcular média, skewness e kurtosis
mean_residuals <- mean(residuals)
skewness_residuals <- skewness(residuals)
kurtosis_residuals <- kurtosis(residuals)

# Criar data frame dos resíduos
residuals_df <- data.frame(residuals)


ggplot(residuals_df, aes(x = residuals)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "lightblue", color = "black", alpha = 0.7) + 
  geom_density(color = "red", size = 1) +  # Curva de densidade
  geom_vline(aes(xintercept = mean_residuals), linetype = "dashed", color = "blue", size = 1) + 
  labs(
    title = "Histograma dos Resíduos do Modelo",
    x = "Resíduos",
    y = "Densidade"
  ) +
  annotate("text", x = mean_residuals + 2.2, y = 0.08, label = paste("Média:", round(mean_residuals, 4)), color = "blue", size = 4, vjust = -1) + 
  annotate("text", x = mean_residuals + 20.8, y = 0.1, label = paste("Skewness:", round(skewness_residuals, 4)), color = "red", size = 4) +  
  annotate("text", x = mean_residuals + 20.25, y = 0.09, label = paste("Kurtosis:", round(kurtosis_residuals, 4)), color = "red", size = 4) +  
  theme_minimal(base_family = "sans") +
  theme(
    plot.title = element_text(size = 16, face = "bold"), 
    axis.title = element_text(size = 14), 
    axis.text = element_text(size = 12) 
  ) + 
  theme_classic()

Por fim, vamos verificar os resultados desse modelo, vamos? Podemo notar bons resultados mesmo em uma modelo bem complicada

# Ajustar o modelo de spline
spline.bs <- lm(y ~ bs(x, df=17, degree=3))

# Previsões
pred.bs <- predict(spline.bs, newdata = list(x = sort.x), se = TRUE)

# Cálculo das métricas
mse <- mean((y - predict(spline.bs))^2)
rmse <- sqrt(mse)
rss <- sum((y - predict(spline.bs))^2)  
mape <- mean(abs((y - predict(spline.bs)) / y)) * 100 
r_squared_adj <- summary(spline.bs)$adj.r.squared

# Criar tabela de métricas
metrics_table <- tibble(
  Metric = c("MSE", "RMSE", "RSS", "MAPE", "R² Ajustado"),
  Value = round(c(mse, rmse, rss, mape, r_squared_adj), 3) 
)

# Criar tabela para IC das previsões
se.bands.bs <- cbind(pred.bs$fit + 2 * pred.bs$se.fit, 
                     pred.bs$fit - 2 * pred.bs$se.fit)
ic_table <- data.frame(
  DF = sort.x,
  Fit = pred.bs$fit,
  Lower_CI = se.bands.bs[,2],
  Upper_CI = se.bands.bs[,1]
)

# Visualizar a tabela de métricas
metrics_table %>%
  kable(caption = "Tabela de Métricas do Modelo") %>%
  kable_styling(full_width = F)
Tabela de Métricas do Modelo
Metric Value
MSE 52.097
RMSE 7.218
RSS 26361.226
MAPE 25.133
R² Ajustado 0.361 
# Exibir apenas o cabeçalho da tabela de IC das previsões
head(ic_table) %>%
  kable(caption = "Cabeçalho da Tabela de Intervalo de Confiança das Previsões") %>%
  kable_styling(full_width = F)
Cabeçalho da Tabela de Intervalo de Confiança das Previsões
DF Fit Lower_CI Upper_CI
0.46 47.13581 34.14091 60.13072
0.74 37.75682 30.64060 44.87304
1.21 30.01824 25.60682 34.42967
1.22 29.94003 25.55568 34.32439
1.25 29.72324 25.42106 34.02542
1.25 29.72324 25.42106 34.02542

Exemplos de Natural Splines e Smoothing Splines

O Natural Splines segue a mesma lógica que os Piecewise Splines, sendo uma alternativa viável para o ajuste de modelos. Assim como os Piecewise Splines, os Natural Splines também minimizam o Residual Sum of Squares (RSS), o que permite que ambos sejam utilizados com a função lm() em R para ajustes de modelos lineares.

Natural Splines

Os Natural Splines introduzem restrições adicionais nas extremidades, garantindo que as funções spline se aproximem de uma linha reta fora dos nós. Isso ajuda a evitar flutuações excessivas nas bordas, melhorando a estabilidade do modelo, especialmente em casos onde os dados podem apresentar variabilidade elevada.

Smoothing Splines

Os Smoothing Splines, por outro lado, aplicam uma abordagem diferente. Eles introduzem uma penalização na função de custo, permitindo que o modelo ajuste a suavidade da curva. Essa suavização pode ser poderosa, mas torna a interpretação dos intervalos de confiança, erros padrão e testes estatísticos mais complexos.

Ambos os métodos oferecem flexibilidade na modelagem de relações não lineares e são opções a serem consideradas quando se trabalha com dados que não se encaixam bem em um modelo de regressão linear simples.

# Natural cubic spline.
spline.ns = lm(y ~ ns(x, df=17))
pred.ns = predict(spline.ns, newdata = list(x = sort.x), se = TRUE)
se.bands.ns = cbind(pred.ns$fit + 2 * pred.ns$se.fit, 
                    pred.ns$fit - 2 * pred.ns$se.fit)

# Smoothing spline
spline.smooth = smooth.spline(x, y, df = 17)

par(mfrow = c(1,3))

# Natural cubic spline.
plot(x, y, cex.lab = 1.1, col="darkgrey", xlab = x.lab, ylab = y.lab, 
     main = "Natural Cubic Spline", bty = 'l')
lines(sort.x, pred.ns$fit, lwd = 2, col = "orange")
matlines(sort.x, se.bands.ns, lwd = 2, col = "orange2", lty = 3)

# Smoothing spline.
plot(x, y, cex.lab = 1.1, col="darkgrey", xlab = x.lab, ylab = y.lab, 
     main = "Smoothing Spline (13 df)", bty = 'l')
lines(spline.smooth, lwd = 2, col = "brown")


# Piecewise Cubic Spline
plot(x, y, cex.lab = 1.1, col="darkgrey", xlab = x.lab, ylab = y.lab, 
     main = "Piecewise Cubic Spline \nNosso modelo", bty = 'l')
lines(sort.x, pred.bs$fit, lwd = 2, col = "steelblue")
matlines(sort.x, se.bands.bs, lwd = 2, col = "steelblue", lty = 3)

# Conclusões

Neste documento, exploramos a relação entre as variáveis medv e indus utilizando abordagens de regressão não linear, especificamente através do uso de splines. A análise inicial dos resíduos revelou a presença de não linearidade, sugerindo que um modelo linear simples não seria suficiente para capturar a complexidade da relação entre as variáveis.

A utilização de gráficos de resíduos parciais permitiu isolar o efeito de cada variável, proporcionando uma visão mais clara de suas interações. Além disso, a consideração de outliers e a heterocedasticidade demonstraram a necessidade de ajustes adicionais no modelo, como a introdução de variáveis preditoras ou transformações.

A introdução de splines, especialmente as Piecewise Splines, mostrou-se uma solução eficaz para modelar a não linearidade, permitindo uma flexibilidade maior na representação das relações complexas entre medv e indus. A discussão sobre a escolha do número de nós com base nos graus de liberdade enfatiza a importância do trade-off entre viés e variância, um aspecto crucial na construção de modelos preditivos robustos.

Para aprofundar a compreensão sobre esses tópicos, as seguintes referências são recomendadas:

  1. Gareth James, Daniela Witten, Trevor Hastie, e Robert Tibshirani. An Introduction to Statistical Learning with Applications in R. Springer, 2013.
  2. Peter Bruce e Andrew Bruce. Practical Statistics for Data Scientists: 50 Essential Concepts. O’Reilly Media, 2016.
  3. Écio Souza. Modelos de Regressão em R. Editora Atlas, 2018.

Esses materiais oferecem uma base sólida para entender não apenas as técnicas discutidas, mas também outras abordagens estatísticas e práticas que podem ser aplicadas em análises de dados complexas.