################################################################################
### Modelos generalizados aditivos para locação, escala e forma
### Exemplo: valores de aluguel de imóveis em Munique
### CE092 - Extensões do modelo de regressão
################################################################################


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

library(gamlss)
library(ggplot2)

options(device = "x11")
options(scipen = 999)


################################################################################
### Leitura e preparação dos dados
################################################################################

data(rent)

head(rent)
summary(rent)
help(rent)

## Variáveis consideradas:
## R   - valor do aluguel
## Fl  - área do imóvel
## A   - ano de construção
## H   - presença de aquecimento central
## loc - classificação da localização do imóvel

rent$H <- factor(rent$H, labels = c("Sim", "Não"))
rent$loc <- factor(rent$loc, labels = c("Abaixo", "Média", "Acima"))

levels(rent$H)
levels(rent$loc)


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

## Distribuição dos valores de aluguel
ggplot(data = rent, aes(x = R)) +
    geom_histogram(aes(y = after_stat(density)), bins = 30, fill = 'grey80', color = 'black') +
    geom_density(linewidth = 1) +
    theme_bw(base_size = 14) +
    labs(x = "Valor do aluguel", y = "Densidade")

## Valor do aluguel segundo a área do imóvel
ggplot(data = rent, aes(x = Fl, y = R)) +
    geom_point(alpha = 0.4) +
    geom_smooth(se = FALSE) +
    scale_x_continuous(breaks = seq(0, 200, 20)) +
    scale_y_continuous(breaks = seq(0, 3000, 500)) +
    theme_bw(base_size = 14) +
    labs(x = "Área do imóvel", y = "Valor do aluguel")

## Valor do aluguel segundo o ano de construção
ggplot(data = rent, aes(x = A, y = R)) +
    geom_point(alpha = 0.4) +
    geom_smooth(se = FALSE) +
    scale_x_continuous(breaks = seq(1800, 2000, 20)) +
    scale_y_continuous(breaks = seq(0, 3000, 500)) +
    theme_bw(base_size = 14) +
    labs(x = "Ano de construção", y = "Valor do aluguel")

## Valor do aluguel segundo a presença de aquecimento central
ggplot(data = rent, aes(x = H, y = R)) +
    geom_boxplot() +
    scale_y_continuous(breaks = seq(0, 3000, 500)) +
    theme_bw(base_size = 14) +
    labs(x = "Aquecimento central", y = "Valor do aluguel")

## Distribuição dos valores de aluguel segundo o aquecimento central
ggplot(data = rent, aes(x = R)) +
    geom_histogram(aes(y = after_stat(density)), bins = 30) +
    geom_density(linewidth = 1) +
    facet_wrap(~H) +
    theme_bw(base_size = 14) +
    labs(x = "Valor do aluguel", y = "Densidade")

## Valor do aluguel segundo a localização
ggplot(data = rent, aes(x = loc, y = R)) +
    geom_boxplot() +
    scale_y_continuous(breaks = seq(0, 3000, 500)) +
    theme_bw(base_size = 14) +
    labs(x = "Classificação da localização", y = "Valor do aluguel")

## Distribuição dos valores de aluguel segundo a localização
ggplot(data = rent, aes(x = R)) +
    geom_histogram(aes(y = after_stat(density)), bins = 30) +
    geom_density(linewidth = 1) +
    facet_wrap(~loc) +
    theme_bw(base_size = 14) +
    labs(x = "Valor do aluguel", y = "Densidade")

## Ano de construção segundo a presença de aquecimento central
ggplot(data = rent, aes(x = H, y = A)) +
    geom_boxplot() +
    theme_bw(base_size = 14) +
    labs(x = "Aquecimento central", y = "Ano de construção")


################################################################################
### Modelo de regressão linear
################################################################################

## Distribuição normal, com média linear nas covariáveis e variância constante
mod1 <- gamlss(R ~ Fl + A + H + loc, family = NO, data = rent)

summary(mod1)
coef(mod1)


################################################################################
### Diagnóstico do modelo de regressão linear
################################################################################

plot(mod1)

## O diagnóstico indica limitações da distribuição normal com variância
## constante para representar adequadamente os valores de aluguel.


################################################################################
### Gráficos de efeitos
################################################################################

term.plot(
    mod1,
    pages = 1,
    ask = FALSE,
    what = "mu",
    las = 1,
    xlabs = c(
        "Área do imóvel",
        "Ano de construção",
        "Aquecimento central",
        "Classificação da localização"
    )
)


################################################################################
### Modelo linear generalizado com resposta gama
################################################################################

## A distribuição gama permite acomodar respostas positivas e assimétricas.
mod2 <- gamlss(R ~ Fl + A + H + loc, family = GA, data = rent)

summary(mod2)
coef(mod2)


################################################################################
### Diagnóstico do modelo linear generalizado
################################################################################

plot(mod2)


################################################################################
### Gráficos de efeitos
################################################################################

term.plot(
    mod2,
    pages = 1,
    ask = FALSE,
    what = "mu",
    las = 1,
    xlabs = c(
        "Área do imóvel",
        "Ano de construção",
        "Aquecimento central",
        "Classificação da localização"
    )
)


################################################################################
### Modelo generalizado aditivo
################################################################################

## Os efeitos de área e ano de construção são representados por P-splines.
mod3 <- gamlss(R ~ pb(Fl) + pb(A) + H + loc, family = GA, data = rent)

summary(mod3)


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

plot(mod3)


################################################################################
### Gráficos de efeitos
################################################################################

## Efeito da área
term.plot(
    mod3,
    pages = 1,
    ask = FALSE,
    what = "mu",
    terms = 1,
    las = 1,
    ylab = "Efeito parcial",
    xlab = "Área do imóvel"
)

## Efeito do ano de construção
term.plot(
    mod3,
    pages = 1,
    ask = FALSE,
    what = "mu",
    terms = 2,
    las = 1,
    ylab = "Efeito parcial",
    xlab = "Ano de construção"
)

## Efeitos de todas as covariáveis
term.plot(
    mod3,
    pages = 1,
    ask = FALSE,
    what = "mu",
    las = 1,
    xlabs = c(
        "Área do imóvel",
        "Ano de construção",
        "Aquecimento central",
        "Classificação da localização"
    )
)


################################################################################
### Modelo generalizado aditivo duplo
################################################################################

## Além do parâmetro de locação, permitimos que o parâmetro de escala varie
## de acordo com as covariáveis.

## Modelo com distribuição gama
mod4 <- gamlss(
    R ~ pb(Fl) + pb(A) + H + loc,
    sigma.fo = ~pb(Fl) + pb(A) + H + loc,
    family = GA,
    data = rent
)

summary(mod4)

## Modelo com distribuição normal inversa
mod4_ig <- gamlss(
    R ~ pb(Fl) + pb(A) + H + loc,
    sigma.fo = ~pb(Fl) + pb(A) + H + loc,
    family = IG,
    data = rent
)

summary(mod4_ig)


################################################################################
### Comparação das distribuições
################################################################################

GAIC(mod4, mod4_ig, k = 2)

## O AIC permite comparar as duas distribuições mantendo a mesma estrutura
## de preditores para os parâmetros de locação e escala.


################################################################################
### Gráficos de efeitos do modelo gama duplo
################################################################################

## Parâmetro de locação
term.plot(
    mod4,
    pages = 1,
    ask = FALSE,
    what = "mu",
    las = 1,
    xlabs = c(
        "Área do imóvel",
        "Ano de construção",
        "Aquecimento central",
        "Classificação da localização"
    )
)

## Parâmetro de escala
term.plot(
    mod4,
    pages = 1,
    ask = FALSE,
    what = "sigma",
    las = 1,
    xlabs = c(
        "Área do imóvel",
        "Ano de construção",
        "Aquecimento central",
        "Classificação da localização"
    )
)


################################################################################
### Diagnóstico do modelo generalizado aditivo duplo
################################################################################

plot(mod4)


################################################################################
### Modelo generalizado aditivo para locação, escala e forma
################################################################################

## A distribuição Box-Cox Cole e Green permite modelar três parâmetros:
## mu, sigma e nu. Nesta parametrização, nu é um parâmetro de forma.

help("BCCGo")

## Algumas funções associadas à distribuição BCCGo
dBCCGo(x = 1.2, mu = 1, sigma = 0.1, nu = 2.5)
pBCCGo(q = 1.2, mu = 1, sigma = 0.1, nu = 2.5)
qBCCGo(p = 0.7, mu = 1, sigma = 0.1, nu = 2.5)


################################################################################
### Simulação da distribuição BCCGo
################################################################################

set.seed(87)

x <- rBCCGo(n = 500, mu = 1, sigma = 0.1, nu = 2.5)

hist(
    x,
    probability = TRUE,
    breaks = 25,
    xlab = "x",
    ylab = "Densidade",
    main = ""
)

curve(
    dBCCGo(x, mu = 1, sigma = 0.1, nu = 2.5),
    add = TRUE,
    lwd = 2
)


################################################################################
### Ajuste do modelo GAMLSS
################################################################################

## Inicialmente, nu é mantido constante. Os parâmetros mu e sigma variam
## de acordo com as covariáveis.
mod6 <- gamlss(
    R ~ pb(Fl) + pb(A) + H + loc,
    sigma.fo = ~pb(Fl) + pb(A) + H + loc,
    family = BCCGo,
    data = rent
)

summary(mod6)

## Em seguida, permitimos que mu, sigma e nu variem de acordo com as
## covariáveis. Essa especificação é deliberadamente flexível e tem
## finalidade didática.
mod7 <- gamlss(
    R ~ pb(Fl) + pb(A) + H + loc,
    sigma.fo = ~pb(Fl) + pb(A) + H + loc,
    nu.fo = ~pb(Fl) + pb(A) + H + loc,
    family = BCCGo,
    data = rent
)

summary(mod7)


################################################################################
### Comparação dos modelos GAMLSS
################################################################################

GAIC(mod4, mod6, mod7, k = 2)

## A comparação deve considerar não apenas o AIC, mas também o diagnóstico,
## a complexidade e a interpretabilidade dos modelos.


################################################################################
### Seleção do preditor para o parâmetro de forma
################################################################################

## Avaliação da contribuição dos termos associados ao parâmetro nu
drop1(mod7, what = "nu")

## Modelo mais parcimonioso para nu
mod7_alt <- update(mod7, ~H, what = "nu")

GAIC(mod7, mod7_alt, k = 2)

summary(mod7_alt)


################################################################################
### Gráficos de efeitos do modelo selecionado
################################################################################

## Parâmetro de locação
term.plot(mod7_alt, pages = 1, ask = FALSE, what = "mu", las = 1)

## Parâmetro de escala
term.plot(mod7_alt, pages = 1, ask = FALSE, what = "sigma", las = 1)

## Parâmetro de forma
term.plot(mod7_alt, pages = 1, ask = FALSE, what = "nu", las = 1)


################################################################################
### Diagnóstico do modelo selecionado
################################################################################

plot(mod7_alt)

## Worm plot
wp(mod7_alt, ylim.all = 0.6)

## O worm plot permite avaliar desvios sistemáticos entre a distribuição
## ajustada e a distribuição dos dados.


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

## Dois perfis de imóveis para ilustração
novos_dados <- data.frame(
    Fl = c(52, 82),
    A = c(1940, 1975),
    H = factor(c("Não", "Sim"), levels = levels(rent$H)),
    loc = factor(c("Abaixo", "Acima"), levels = levels(rent$loc))
)

rownames(novos_dados) <- c("Imóvel 1", "Imóvel 2")

novos_dados


################################################################################
### Predições para o parâmetro de locação
################################################################################

## Escala do preditor
predict(mod7_alt, newdata = novos_dados, what = "mu")

## Escala do parâmetro
pmu <- predict(
    mod7_alt,
    newdata = novos_dados,
    what = "mu",
    type = "response"
)

pmu


################################################################################
### Predições para o parâmetro de escala
################################################################################

## Escala do preditor
predict(mod7_alt, newdata = novos_dados, what = "sigma")

## Escala do parâmetro
psigma <- predict(
    mod7_alt,
    newdata = novos_dados,
    what = "sigma",
    type = "response"
)

psigma


################################################################################
### Predições para o parâmetro de forma
################################################################################

pnu <- predict(
    mod7_alt,
    newdata = novos_dados,
    what = "nu",
    type = "response"
)

pnu


################################################################################
### Distribuições estimadas para os dois perfis
################################################################################

curve(
    dBCCGo(x, mu = pmu[1], sigma = psigma[1], nu = pnu[1]),
    from = 0,
    to = 3500,
    xlab = "Valor do aluguel",
    ylab = "Densidade",
    lwd = 2
)

curve(
    dBCCGo(x, mu = pmu[2], sigma = psigma[2], nu = pnu[2]),
    from = 0,
    to = 3500,
    lwd = 2,
    lty = 2,
    add = TRUE
)

legend(
    "topright",
    legend = c("Imóvel 1", "Imóvel 2"),
    lty = c(1, 2),
    lwd = 2,
    bty = "n"
)


################################################################################
### Probabilidades estimadas
################################################################################

## Probabilidade de aluguel inferior a 500
pBCCGo(500, mu = pmu[1], sigma = psigma[1], nu = pnu[1])
pBCCGo(500, mu = pmu[2], sigma = psigma[2], nu = pnu[2])

## Probabilidade de aluguel superior a 1000
pBCCGo(
    1000,
    mu = pmu[1],
    sigma = psigma[1],
    nu = pnu[1],
    lower.tail = FALSE
)

pBCCGo(
    1000,
    mu = pmu[2],
    sigma = psigma[2],
    nu = pnu[2],
    lower.tail = FALSE
)


################################################################################
### Quantis estimados
################################################################################

## Medianas
qBCCGo(0.50, mu = pmu[1], sigma = psigma[1], nu = pnu[1])
qBCCGo(0.50, mu = pmu[2], sigma = psigma[2], nu = pnu[2])

## Quantis de ordens 0,75 e 0,90
qBCCGo(c(0.75, 0.90), mu = pmu[1], sigma = psigma[1], nu = pnu[1])
qBCCGo(c(0.75, 0.90), mu = pmu[2], sigma = psigma[2], nu = pnu[2])


################################################################################
### Comparação dos modelos ajustados
################################################################################

GAIC(mod1, mod2, mod3, mod4, mod7_alt, k = 2)