################################################################################
### Modelos aditivos generalizados
### Sessão prática: ajuste e interpretação de um GAM
### Exemplo: emissões de NOx na combustão de etanol
### CE092 - Extensões do modelo de regressão
################################################################################


################################################################################
### Pacotes
################################################################################

library(ggplot2)
library(lattice)
library(mgcv)
library(patchwork)

setwd(this.path::here())
options(device = "x11")


################################################################################
### Leitura dos dados
################################################################################

data(ethanol, package = "lattice")

dados <- ethanol[, c("E", "NOx")]

str(dados)
summary(dados)


################################################################################
### Análise descritiva
################################################################################

ggplot(dados, aes(x = E, y = NOx)) +
    geom_point(size = 2) +
    labs(x = "Razão de equivalência", y = "Emissão de NOx") +
    theme_bw(base_size = 14)+
    scale_x_continuous(breaks = seq(0.5, 1.2, 0.1))


################################################################################
### Ajuste do GAM
################################################################################

## A função s(E) representa um efeito suave da razão de equivalência.
## A base "tp" corresponde a thin plate regression splines.

## Thin plate regression splines: representação de funções suaves derivada dos thin 
## plate splines, em que a complexidade da função é controlada pela penalização de 
## sua curvatura. Tem a vantagem de não exigir a escolha explícita da posição de 
## nós (knots).

## O argumento k define a dimensão da base disponível para o termo suave.
## O método REML é utilizado para estimar o parâmetro de suavização.

ajuste <- gam(
    NOx ~ s(E, bs = "tp", k = 10),
    data = dados,
    method = "REML"
)

ajuste

## k = 10 significa que a dimensão da base é 10, limitando a complexidade
## máxima disponível para representar s(E). Neste caso, o EDF máximo é
## aproximadamente k - 1 = 9.

################################################################################
### Resumo do ajuste
################################################################################

summary(ajuste)

## Elementos importantes da saída:
## - edf: graus de liberdade efetivos do termo suave;
## - Ref.df: graus de liberdade de referência usados no teste aproximado;
## - F e p-value: teste aproximado para o termo suave;
## - R-sq.(adj): coeficiente de determinação ajustado;
## - Deviance explained: proporção do desvio explicada pelo modelo;
## - REML: critério utilizado na estimação do parâmetro de suavização.

tab_suave <- summary(ajuste)$s.table
tab_suave

EDF <- tab_suave[1, "edf"]
EDF


################################################################################
### Efeito suave estimado
################################################################################

## O gráfico padrão de plot.gam() apresenta o efeito parcial s(E).
## Por construção, o termo suave é centrado aproximadamente em zero.

plot(
    ajuste,
    shade = TRUE,
    residuals = TRUE,
    pch = 16,
    cex = 0.7,
    xlab = "Razão de equivalência",
    ylab = "Efeito parcial estimado"
)


################################################################################
### Curva ajustada na escala da resposta
################################################################################

## Agora são obtidas as médias estimadas de NOx na escala original da resposta.

novo <- data.frame(E = seq(min(dados$E), max(dados$E), length.out = 300))

pred <- predict(ajuste, newdata = novo, se.fit = TRUE)

plot(NOx ~ E, data = dados, pch = 16,
     xlab = "Razão de equivalência", ylab = "Emissão de NOx", las = 1)

lines(novo$E, pred$fit, lwd = 2)
lines(novo$E, pred$fit - 1.96 * pred$se.fit, lty = 2)
lines(novo$E, pred$fit + 1.96 * pred$se.fit, lty = 2)


################################################################################
### Diagnóstico do ajuste
################################################################################

## gam.check() apresenta diagnósticos gerais do modelo e avalia se a dimensão
## escolhida para a base (k) parece suficiente para representar o termo suave.

gam.check(ajuste)

## p=0.005. Há evidência de que k=10 pode estar restringindo demais o smooth.

################################################################################
### Resíduos
################################################################################

dados$residuo <- residuals(ajuste, type = "deviance")
dados$ajustado <- fitted(ajuste)

g_resid <- ggplot(dados, aes(x = ajustado, y = residuo)) +
    geom_point(size = 2) +
    geom_hline(yintercept = 0, linetype = 2) +
    geom_smooth(method = "loess", se = FALSE) +
    labs(x = "Valores ajustados", y = "Resíduos de desvio") +
    theme_bw(base_size = 14)

g_qq <- ggplot(dados, aes(sample = residuo)) +
    stat_qq() +
    stat_qq_line() +
    labs(x = "Quantis teóricos", y = "Quantis observados") +
    theme_bw(base_size = 14)

g_resid + g_qq


################################################################################
### Avaliação da dimensão da base
################################################################################

## Como verificação, podemos ajustar um modelo com uma base mais ampla.

ajuste_k20 <- gam(
    NOx ~ s(E, bs = "tp", k = 20),
    data = dados,
    method = "REML"
)

summary(ajuste_k20)
gam.check(ajuste_k20)

## p=0.07. Não há evidência de que k=20 pode estar restringindo demais o smooth.

################################################################################
### Comparação entre k = 10 e k = 20
################################################################################

novo$pred_k10 <- predict(ajuste, newdata = novo, type = "response")
novo$pred_k20 <- predict(ajuste_k20, newdata = novo, type = "response")

ggplot(dados, aes(x = E, y = NOx)) +
    geom_point(size = 1.8) +
    geom_line(
        data = novo,
        aes(x = E, y = pred_k10, linetype = "k = 10"),
        inherit.aes = FALSE,
        linewidth = 1
    ) +
    geom_line(
        data = novo,
        aes(x = E, y = pred_k20, linetype = "k = 20"),
        inherit.aes = FALSE,
        linewidth = 0.9
    ) +
    scale_linetype_manual(
        values = c("k = 10" = 1, "k = 20" = 2),
        name = NULL
    ) +
    labs(x = "Razão de equivalência", y = "Emissão de NOx") +
    scale_x_continuous(breaks = seq(0.5, 1.2, 0.1)) +
    theme_bw(base_size = 14) +
    theme(legend.position = "bottom")


################################################################################
### Síntese
################################################################################

## - s(E) permite representar uma relação não linear sem especificar previamente
##   uma forma paramétrica para a relação entre E e NOx.
##
## - O EDF resume a complexidade efetivamente utilizada pelo termo suave.
##
## - O gráfico de s(E) representa um efeito parcial, enquanto a predição com
##   type = "response" fornece a média estimada na escala da resposta.
##
## - O diagnóstico deve avaliar tanto a adequação usual do modelo quanto se a
##   dimensão da base k foi suficiente para representar o padrão dos dados.