################################################################################
### Modelos aditivos generalizados
### Exemplo: desempenho educacional e indicadores socioeconômicos
### CE092 - Extensões do modelo de regressão
################################################################################


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

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

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


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

load("pisa.Rdata")

str(pisa)
summary(pisa)
head(pisa)


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

## Relações marginais entre Overall e cada uma das covariáveis.

g1 <- ggplot(pisa, aes(x = Income, y = Overall)) +
    geom_point(size = 1.5) +
    geom_smooth(se = FALSE) +
    labs(x = "Renda", y = "Desempenho geral") +
    theme_bw(base_size = 14)

g2 <- ggplot(pisa, aes(x = Edu, y = Overall)) +
    geom_point(size = 1.5) +
    geom_smooth(se = FALSE) +
    labs(x = "Educação", y = "Desempenho geral") +
    theme_bw(base_size = 14)

g1 + g2


################################################################################
### Modelo linear aditivo
################################################################################

## Modelo com efeitos lineares de Income e Edu.

ajuste_linear <- gam(
    Overall ~ Income + Edu,
    data = pisa,
    method = "REML"
)

summary(ajuste_linear)


################################################################################
### GAM aditivo
################################################################################

## O modelo aditivo permite relações não lineares separadas para Income e Edu.

ajuste_aditivo <- gam(
    Overall ~ s(Income) + s(Edu),
    data = pisa,
    method = "REML"
)

summary(ajuste_aditivo)


################################################################################
### Efeitos parciais do modelo aditivo
################################################################################

## Cada curva representa o efeito parcial da respectiva covariável,
## controlando pelo efeito da outra covariável.

g1 <- visreg(ajuste_aditivo, "Income",
             points = list(size = 1.5))

g2 <- visreg(ajuste_aditivo, "Edu",
             points = list(size = 1.5))

g1 + g2


################################################################################
### Superfície aditiva estimada
################################################################################

## Mesmo sem interação, a soma dos dois efeitos suaves produz uma superfície
## estimada sobre o plano definido por Income e Edu.

vis.gam(
    ajuste_aditivo,
    view = c("Income", "Edu"),
    plot.type = "contour",
    color = "heat",
    too.far = 0.10
)


################################################################################
### GAM com efeito bidimensional
################################################################################

## O termo te(Income, Edu) representa um efeito suave conjunto das duas
## covariáveis, permitindo que o efeito de Income dependa de Edu e vice-versa.
##
## Tensor product smooths são especialmente úteis quando as covariáveis possuem
## escalas ou unidades distintas.

ajuste_interacao <- gam(
    Overall ~ te(Income, Edu),
    data = pisa,
    method = "REML"
)

summary(ajuste_interacao)


################################################################################
### Visualização do efeito bidimensional
################################################################################

plot(
    ajuste_interacao,
    select = 1,
    scheme = 2
)

vis.gam(
    ajuste_interacao,
    view = c("Income", "Edu"),
    plot.type = "persp",
    theta = -45,
    phi = 30,
    color = "heat",
    too.far = 0.10
)


################################################################################
### Comparação conceitual dos modelos
################################################################################

## Modelo linear:
## E(Overall | Income, Edu) = beta0 + beta1*Income + beta2*Edu
##
## Modelo aditivo:
## E(Overall | Income, Edu) = beta0 + s1(Income) + s2(Edu)
##
## Modelo bidimensional:
## E(Overall | Income, Edu) = beta0 + s(Income, Edu)


################################################################################
### Comparação dos modelos via ML
################################################################################

## Para comparar estruturas diferentes, ajustamos os modelos por máxima
## verossimilhança.

ajuste_linear_ml <- gam(
    Overall ~ Income + Edu,
    data = pisa,
    method = "ML"
)

ajuste_aditivo_ml <- gam(
    Overall ~ s(Income) + s(Edu),
    data = pisa,
    method = "ML"
)

ajuste_interacao_ml <- gam(
    Overall ~ te(Income, Edu),
    data = pisa,
    method = "ML"
)

AIC(
    ajuste_linear_ml,
    ajuste_aditivo_ml,
    ajuste_interacao_ml
)


################################################################################
### Comparação entre modelo aditivo e efeito conjunto
################################################################################

## A comparação deve considerar conjuntamente:
## - AIC;
## - complexidade dos termos suaves;
## - forma dos efeitos estimados;
## - plausibilidade substantiva;
## - comportamento dos resíduos.

summary(ajuste_aditivo)
summary(ajuste_interacao)


################################################################################
### Diagnóstico do modelo aditivo
################################################################################

gam.check(ajuste_aditivo)


################################################################################
### Diagnóstico do modelo bidimensional
################################################################################

gam.check(ajuste_interacao)


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

## Após a comparação das estruturas, o modelo escolhido deve ser reajustado
## por REML.

## Exemplo: se o modelo bidimensional for escolhido.

ajuste_final <- gam(
    Overall ~ te(Income, Edu),
    data = pisa,
    method = "REML"
)

summary(ajuste_final)

vis.gam(
    ajuste_final,
    view = c("Income", "Edu"),
    plot.type = "contour",
    color = "heat",
    too.far = 0.10
)

gam.check(ajuste_final)


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

## - O modelo linear supõe efeitos lineares e aditivos de Income e Edu.
##
## - O GAM aditivo permite efeitos não lineares separados para cada covariável.
##
## - O termo te(Income, Edu) permite um efeito suave conjunto, de modo que o
##   efeito de uma covariável possa variar conforme o valor da outra.
##
## - A superfície do modelo aditivo resulta da soma de dois efeitos univariados.
##
## - A superfície bidimensional é mais flexível e permite representar padrões
##   que não podem ser decompostos em efeitos aditivos separados.
##
## - A comparação entre os modelos deve considerar ajuste, complexidade,
##   interpretação, diagnóstico e plausibilidade substantiva.