pacman::p_load(tidyverse, dplyr, ggplot2, quantreg, qgam, splines)

Introdução

Este documento apresenta a aplicação e comparação de três modelos de regressão quantílica para predizer o volume expiratório forçado no primeiro segundo (FEV1) em mulheres, segundo os quantis 0.05, 0.5 e 0.95 da distribuição.

Os modelos utilizados foram:

rq linear: regressão quantílica tradicional, log-linear.

rq + splines: regressão quantílica com splines cúbicos para capturar não linearidades.

qgam: modelo aditivo quantílico generalizado com suavização penalizada.

O desempenho preditivo foi avaliado pelo Erro Absoluto Mediano (MAE) fora da amostra, métrica robusta que indica a mediana das diferenças absolutas entre valores preditos e observados. Menor MAE indica melhor precisão.

Carregamento e transformação dos dados

setwd("C:/Users/Usuario/Downloads")
a <- read.csv2("jones_espiro.csv", sep = ",")

# Transformações iniciais
a <- a %>%
  mutate(
    id = as.factor(id),
    cidade = as.factor(cidade),
    sex = factor(sex, levels = c(1, 2), labels = c("Masculino", "Feminino")),
    ethnic = as.factor(ethnic),
    age = as.numeric(age),
    height = as.numeric(height) / 100,  # cm para metros
    weight = as.numeric(weight),
    fvc = as.numeric(fvc),
    fev1 = as.numeric(fev1),
    fev1fvc = as.numeric(fev1fvc),
    fef2575 = as.numeric(fef2575),
    branco = factor(ifelse(ethnic == 1, 1, 0))
  )

# Filtro para mulheres
dfm <- filter(a, sex == "Feminino")

Ajuste dos Modelos Quantílicos e Predições

Regressão quantílica

# Definindo função para calcular MAE fora da amostra
# Ajuste do modelo quantílico linear para quantis 5%, 50% e 95%
mod_rq <- rq(log(fev1) ~ log(age) + log(height) + branco, data = dfm, tau = c(0.05, 0.5, 0.95))

# Criar grade para predição, fixando idade na mediana, branco = 1
newdata <- data.frame(
  age = median(dfm$age, na.rm = TRUE),
  height = seq(min(dfm$height, na.rm = TRUE), max(dfm$height, na.rm = TRUE), length.out = 100),
  branco = factor(1, levels = levels(dfm$branco))
)

# Predições nos quantis para a grade de alturas
preds <- predict(mod_rq, newdata = newdata)

# Organizar dados para plotagem
preds_long <- as.data.frame(preds) %>%
  mutate(height = newdata$height) %>%
  pivot_longer(cols = starts_with("tau"), names_to = "quantil", values_to = "log_predito") %>%
  mutate(predito = exp(log_predito))  # Retorna à escala original

# Visualização
p1 = ggplot(preds_long, aes(x = height, y = predito, color = quantil)) +
  geom_line(size = 1.2) +
  geom_point(data = dfm, aes(x = height, y = fev1), inherit.aes = FALSE, alpha = 0.3, color = "grey30") +
  labs(title = "Predição do FEV1 pelos quantis 0.05, 0.5 e 0.95 (Regressão Quantílica Linear)",
       x = "Altura (m)",
       y = "FEV1 (L)",
       color = "Quantil") +
  scale_color_manual(values = c("tau= 0.05" = "blue", "tau= 0.50" = "black", "tau= 0.95" = "red"),
                     labels = c("Q05", "Q50", "Q95")) +
  theme_classic()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
p1
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

Especificação: Modelo linear clássico, assume relação log-linear entre preditores (log(age), log(height), branco) e log(fev1).

Flexibilidade: Baixa, captura apenas efeitos lineares na escala log.

Limitação: Incapaz de modelar efeitos não lineares, o que pode levar a viés se a relação verdadeira for curva.

Regressão quantílica aditiva generalizada

# Ajustes para os quantis 5%, 50% e 95% com termos suavizados para idade e altura


mod_qgam_05 <- qgam(log(fev1) ~ s(log(age)) + s(log(height)) + branco, data = dfm, qu = 0.05)
## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.05.........done
mod_qgam_50 <- qgam(log(fev1) ~ s(log(age)) + s(log(height)) + branco, data = dfm, qu = 0.5)
## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.5.........done
mod_qgam_95 <- qgam(log(fev1) ~ s(log(age)) + s(log(height)) + branco, data = dfm, qu = 0.95)
## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.95..............done

# Predição para a mesma grade
newdata$pred_log_05 <- predict(mod_qgam_05, newdata)
newdata$pred_log_50 <- predict(mod_qgam_50, newdata)
newdata$pred_log_95 <- predict(mod_qgam_95, newdata)

newdata <- newdata %>%
  mutate(
    pred_05 = exp(pred_log_05),
    pred_50 = exp(pred_log_50),
    pred_95 = exp(pred_log_95)
  )

# Filtrar para idade próxima da mediana para visualização
grid_filtered <- newdata %>% filter(abs(age - median(dfm$age, na.rm = TRUE)) < 0.1)

df_plot <- grid_filtered %>%
  select(height, pred_05, pred_50, pred_95) %>%
  pivot_longer(cols = starts_with("pred_"), names_to = "quantil", values_to = "FEV1") %>%
  mutate(quantil = factor(quantil, levels = c("pred_05", "pred_50", "pred_95"), labels = c("Q05", "Q50", "Q95")))

# Plotagem dos modelos GAM quantílicos
ggplot(df_plot, aes(x = height, y = FEV1, color = quantil)) +
  geom_line(size = 1) +
  geom_point(data = dfm, aes(x = height, y = fev1), alpha = 0.3, color = "grey30") +
  scale_color_manual(values = c("Q05" = "blue", "Q50" = "black", "Q95" = "red")) +
  labs(title = "Curvas Quantílicas GAM para FEV1 vs Altura (Idade ~ Mediana)",
       x = "Altura (m)",
       y = "FEV1 (L)",
       color = "Quantil") +
  theme_minimal()

Especificação: Estima os quantis da resposta incluindo termos suavizados (splines) para log(age) e log(height); modelo semi-paramétrico.

Flexibilidade: Alta, permite capturar relações não lineares e mais complexas entre preditores e resposta.

Vantagem: Ajusta curvas suaves, evita suposições rígidas de linearidade.

Limitação: Mais complexo, pode superajustar em amostras pequenas ou com ruído.

Regressão quantílica com splines cúbicos

# Ajustes para os quantis 5%, 50% e 95% com termos suavizados para idade e altura
mod_qgam_05 <- qgam(log(fev1) ~ s(log(age)) + s(log(height)) + branco, data = dfm, qu = 0.05)
## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.05.........done
mod_qgam_50 <- qgam(log(fev1) ~ s(log(age)) + s(log(height)) + branco, data = dfm, qu = 0.5)
## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.5.........done
mod_qgam_95 <- qgam(log(fev1) ~ s(log(age)) + s(log(height)) + branco, data = dfm, qu = 0.95)
## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.95..............done

# Predição para a mesma grade
newdata$pred_log_05 <- predict(mod_qgam_05, newdata)
newdata$pred_log_50 <- predict(mod_qgam_50, newdata)
newdata$pred_log_95 <- predict(mod_qgam_95, newdata)

newdata <- newdata %>%
  mutate(
    pred_05 = exp(pred_log_05),
    pred_50 = exp(pred_log_50),
    pred_95 = exp(pred_log_95)
  )

# Filtrar para idade próxima da mediana para visualização
grid_filtered <- newdata %>% filter(abs(age - median(dfm$age, na.rm = TRUE)) < 0.1)

df_plot <- grid_filtered %>%
  select(height, pred_05, pred_50, pred_95) %>%
  pivot_longer(cols = starts_with("pred_"), names_to = "quantil", values_to = "FEV1") %>%
  mutate(quantil = factor(quantil, levels = c("pred_05", "pred_50", "pred_95"), labels = c("Q05", "Q50", "Q95")))

# Plotagem dos modelos GAM quantílicos
ggplot(df_plot, aes(x = height, y = FEV1, color = quantil)) +
  geom_line(size = 1) +
  geom_point(data = dfm, aes(x = height, y = fev1), alpha = 0.3, color = "grey30") +
  scale_color_manual(values = c("Q05" = "blue", "Q50" = "black", "Q95" = "red")) +
  labs(title = "Curvas Quantílicas GAM para FEV1 vs Altura (Idade ~ Mediana)",
       x = "Altura (m)",
       y = "FEV1 (L)",
       color = "Quantil") +
  theme_minimal()
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

Especificação: Modelo linear quantílico, porém com splines cúbicos de base para log(age) e log(height), controlando a flexibilidade via graus de liberdade.

Flexibilidade: Moderada a alta, similar ao qgam, mas modelagem explícita via splines cúbicos.

Vantagem: Permite não linearidade com estrutura mais controlada (df fixos), facilita interpretação dos componentes spline.

Limitação: Menos flexível que qgam, depende da escolha dos graus de liberdade para evitar under/overfitting.

Avaliação Comparativa dos Modelos Quantílicos

A capacidade preditiva dos modelos quantílicos foi avaliada com base na divisão aleatória da base em treino (80%) e teste (20%). Três modelos foram comparados:

rq linear: regressão quantílica linear,
qgam: regressão quantílica aditiva generalizada,
rq + splines: regressão quantílica com splines cúbicos.

A métrica utilizada foi o erro absoluto mediano (MAE), calculado sobre a escala original (litros). Os resultados para os quantis 5%, 50% e 95% são apresentados a seguir:

library(knitr)
library(purrr)

# Dividir base
n <- nrow(dfm)
set.seed(123)  # para reprodutibilidade
idx_treino <- sample(1:n, size = 0.8 * n)
  treino <- dfm[idx_treino, ]
teste <- dfm[-idx_treino, ]

# Função de erro
calc_mae <- function(mod, teste, log_transform = TRUE) {
  preds <- predict(mod, newdata = teste)
  if (log_transform) preds <- exp(preds)
  median(abs(preds - teste$fev1), na.rm = TRUE)
}

# Quantis e modelos
quantis <- c(0.05, 0.5, 0.95)
resultados <- map_dfr(quantis, function(tau) {
  mod_rq <- rq(log(fev1) ~ log(age) + log(height) + branco, data = treino, tau = tau)
  mae_rq <- calc_mae(mod_rq, teste)
  
  mod_qgam <- qgam(log(fev1) ~ s(log(age)) + s(log(height)) + branco, data = treino, qu = tau)
  mae_qgam <- calc_mae(mod_qgam, teste)
  
  mod_rq_bs <- rq(log(fev1) ~ bs(log(age), df = 4) + bs(log(height), df = 4) + branco, data = treino, tau = tau)
  mae_rq_bs <- calc_mae(mod_rq_bs, teste)
  
  tibble(
    Quantil = paste0("Q", sprintf("%02d", tau * 100)),
    Modelo = c("rq linear", "qgam", "rq + splines"),
    MAE = c(mae_rq, mae_qgam, mae_rq_bs)
  )
})

Estimating learning rate. Each dot corresponds to a loss evaluation. qu = 0.05………..done Estimating learning rate. Each dot corresponds to a loss evaluation. qu = 0.5……..done Estimating learning rate. Each dot corresponds to a loss evaluation. qu = 0.95………..done



resultados %>%
  arrange(Quantil, MAE) %>%
  mutate(MAE = round(MAE, 3)) %>%
  kable(
    format = "markdown",
    caption = "Erro absoluto mediano (MAE) por modelo e quantil.",
    align = c('l', 'l', 'c'),
    col.names = c("Quantil", "Modelo", "MAE (L)")
  )

Erro absoluto mediano (MAE) por modelo e quantil.
Quantil	Modelo	MAE (L)
Q05	rq + splines	0.333
Q05	rq linear	0.338
Q05	qgam	0.355
Q50	rq linear	0.151
Q50	qgam	0.153
Q50	rq + splines	0.163
Q95	rq + splines	0.379
Q95	qgam	0.401
Q95	rq linear	0.412

Desempenho preditivo dos modelos quantílicos A Tabela apresenta os valores do erro absoluto mediano (MAE) para três abordagens de regressão quantílica (linear, com splines cúbicos e aditiva generalizada via qgam), avaliadas nos quantis 0,05, 0,50 e 0,95. O MAE foi calculado na escala original do desfecho (FEV1, em litros), com base em um conjunto de teste independente.

No quantil inferior (Q05), que representa indivíduos com valores mais baixos de FEV1, o modelo com splines cúbicos apresentou o melhor desempenho (MAE = 0,333 L), seguido pela regressão linear (0,338 L), enquanto o modelo qgam apresentou o maior erro (0,355 L). Isso sugere que, para a cauda inferior da distribuição, ganhos marginais de precisão são alcançados com o uso de splines, embora as diferenças sejam pequenas.

No quantil mediano (Q50), os modelos apresentaram desempenho bastante semelhante. A regressão linear obteve o menor erro (0,151 L), ligeiramente melhor que o qgam (0,153 L), e ambos superaram o modelo com splines (0,163 L). Isso indica que, no centro da distribuição do FEV1, a relação entre as covariáveis e o desfecho pode ser adequadamente capturada por uma estrutura linear simples.

Já no quantil superior (Q95), os resultados se invertem: o modelo com splines novamente foi o mais preciso (0,379 L), seguido pelo qgam (0,401 L), e a regressão linear obteve o pior desempenho (0,412 L). Essa tendência reforça a hipótese de que nas extremidades da distribuição — onde padrões mais complexos ou não lineares podem emergir — modelos com maior flexibilidade (como splines) oferecem vantagens preditivas.

Em síntese, não há um modelo universalmente superior: o desempenho depende do quantil de interesse. Modelos mais flexíveis (splines, qgam) tendem a melhorar o ajuste nas extremidades da distribuição, ao custo de desempenho marginalmente inferior no centro. Para aplicações que envolvam identificação de casos extremos de função pulmonar, os modelos com splines cúbicos mostram-se mais apropriados.

Algumas renderizações não permitiram omitir o texto de saída de alguns tipos de regressão quantílica.

Equações para predição - Espirometria

Invíduos do sexo feminino

Marcos Brum

2025-06-01