################################################################################
### Modelos aditivos generalizados
### Exemplo: cifose após cirurgia da coluna
### CE092 - Extensões do modelo de regressão
################################################################################


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

library(ggplot2)
library(mgcv)
library(patchwork)
library(visreg)
library(DHARMa)

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


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

data("kyphosis", package = "rpart")

str(kyphosis)
summary(kyphosis)
head(kyphosis)


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

## Distribuição das covariáveis segundo a ocorrência de cifose.

g1 <- ggplot(kyphosis, aes(x = Age, fill = Kyphosis)) +
    geom_histogram(
        colour = "grey50",
        alpha = 0.5,
        position = "identity",
        bins = 10
    ) +
    labs(x = "Idade", y = "Frequência", fill = "Cifose") +
    theme_bw(base_size = 14) +
    theme(legend.position = "bottom")

g2 <- ggplot(kyphosis, aes(x = Number, fill = Kyphosis)) +
    geom_bar(
        colour = "grey50",
        alpha = 0.5,
        position = "identity"
    ) +
    scale_x_continuous(breaks = 0:10) +
    labs(x = "Número de vértebras", y = "Frequência", fill = "Cifose") +
    theme_bw(base_size = 14) +
    theme(legend.position = "bottom")

g3 <- ggplot(kyphosis, aes(x = Start, fill = Kyphosis)) +
    geom_histogram(
        colour = "grey50",
        alpha = 0.5,
        position = "identity",
        bins = 10
    ) +
    scale_x_continuous(breaks = seq(0, 20, 2)) +
    labs(x = "Vértebra inicial", y = "Frequência", fill = "Cifose") +
    theme_bw(base_size = 14) +
    theme(legend.position = "bottom")

g1 + g2 + g3


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

## Modelo logístico aditivo:
##
## logit[P(Kyphosis = "present")] =
##     beta0 + s(Age) + beta1*Number + s(Start)
##
## Os termos s(Age) e s(Start) permitem relações não lineares com a
## probabilidade de ocorrência de cifose.

ajuste <- gam(
    Kyphosis ~ s(Age) + Number + s(Start),
    family = binomial(link = "logit"),
    data = kyphosis,
    method = "REML"
)

summary(ajuste)


################################################################################
### Graus de liberdade efetivos
################################################################################

## O EDF resume a complexidade efetivamente utilizada pelos termos suaves.
## Valores próximos de 1 sugerem relações aproximadamente lineares.

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

tab_suave[, "edf"]


################################################################################
### Efeitos parciais estimados
################################################################################

## Efeitos na escala do preditor linear.

g1 <- visreg(ajuste, "Age",
    scale = "linear",
    points = list(size = 1.5)
)

g2 <- isreg(ajuste, "Start",
    scale = "linear",
    points = list(size = 1.5)
)

g3 <- visreg(ajuste, "Number",
    scale = "linear",
    points = list(size = 1.5)
)

g1 + g2 + g3


################################################################################
### Efeitos na escala da resposta
################################################################################


g1 <- visreg(ajuste, "Age",
             scale = "response",
             points = list(size = 1.5)
)

g2 <- visreg(ajuste, "Start",
            scale = "response",
            points = list(size = 1.5)
)

g3 <- visreg(ajuste, "Number",
             scale = "response",
             points = list(size = 1.5)
)

g1 + g2 + g3


################################################################################
### Interpretação dos termos
################################################################################

## Nos termos suaves, a interpretação principal é feita pela forma estimada da
## função, seu intervalo de confiança e seu EDF.
##
## Para Number, o efeito permanece paramétrico e pode ser interpretado pelo
## coeficiente estimado na escala do logito.

coef(ajuste)
summary(ajuste)$p.table


################################################################################
### Modelo com Age linear
################################################################################

ajuste_age_linear <- gam(
    Kyphosis ~ Age + Number + s(Start),
    family = binomial(link = "logit"),
    data = kyphosis,
    method = "REML"
)

summary(ajuste_age_linear)


################################################################################
### Modelo com Start linear
################################################################################

ajuste_start_linear <- gam(
    Kyphosis ~ s(Age) + Number + Start,
    family = binomial(link = "logit"),
    data = kyphosis,
    method = "REML"
)

summary(ajuste_start_linear)


################################################################################
### Modelo logístico linear
################################################################################

ajuste_linear <- gam(
    Kyphosis ~ Age + Number + Start,
    family = binomial(link = "logit"),
    data = kyphosis,
    method = "REML"
)

summary(ajuste_linear)


################################################################################
### Comparação de estruturas via ML
################################################################################

## Para comparar modelos com diferentes estruturas de termos, os candidatos
## são ajustados por máxima verossimilhança.

ajuste_ml <- gam(
    Kyphosis ~ s(Age) + Number + s(Start),
    family = binomial(link = "logit"),
    data = kyphosis,
    method = "ML"
)

ajuste_age_linear_ml <- gam(
    Kyphosis ~ Age + Number + s(Start),
    family = binomial(link = "logit"),
    data = kyphosis,
    method = "ML"
)

ajuste_start_linear_ml <- gam(
    Kyphosis ~ s(Age) + Number + Start,
    family = binomial(link = "logit"),
    data = kyphosis,
    method = "ML"
)

ajuste_linear_ml <- gam(
    Kyphosis ~ Age + Number + Start,
    family = binomial(link = "logit"),
    data = kyphosis,
    method = "ML"
)

AIC(
    ajuste_ml,
    ajuste_age_linear_ml,
    ajuste_start_linear_ml,
    ajuste_linear_ml
)


################################################################################
### Predições
################################################################################

novos_dados <- data.frame(
    Age = c(60, 30, 90),
    Number = c(3, 2, 5),
    Start = c(5, 10, 13)
)

## Predições na escala do preditor linear.

pred_link <- predict(
    ajuste,
    newdata = novos_dados,
    type = "link",
    se.fit = TRUE
)

pred_link


################################################################################
### Predições na escala da resposta
################################################################################

## Predições da probabilidade de ocorrência de cifose.

pred_resp <- predict(
    ajuste,
    newdata = novos_dados,
    type = "response",
    se.fit = TRUE
)

pred_resp


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

## gam.check() apresenta diagnósticos usuais e verifica a suficiência da
## dimensão das bases utilizadas nos termos suaves.

gam.check(ajuste)


################################################################################
### Ajuste final
################################################################################

## Após a definição da estrutura, o modelo escolhido deve ser reajustado por
## REML. A fórmula abaixo mantém os dois termos suaves.

ajuste_final <- gam(
    Kyphosis ~ s(Age) + Number + s(Start),
    family = binomial(link = "logit"),
    data = kyphosis,
    method = "REML"
)

summary(ajuste_final)


################################################################################
### Diagnóstico do ajuste usando resíduos quantílicos normalizados aleatorizados
################################################################################

gam.check(ajuste_final)


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

## - A distribuição da resposta e a função de ligação seguem a estrutura de
##   um modelo logístico binomial.
##
## - O GAM flexibiliza a forma como Age e Start entram no preditor linear.
##
## - Os termos suaves são interpretados pela forma estimada, pelo intervalo
##   de confiança e pelo EDF.
##
## - O efeito de Number permanece paramétrico e pode ser interpretado pelo
##   coeficiente correspondente.
##
## - As curvas podem ser visualizadas tanto na escala do logito quanto na
##   escala da probabilidade.
##
## - A comparação de estruturas pode ser feita com modelos ajustados por ML,
##   enquanto o modelo final pode ser reajustado por REML.
##
## - O diagnóstico deve avaliar tanto o comportamento dos resíduos quanto a
##   suficiência da dimensão das bases.