## Fonte:
## Rats data (Verbeke & Mohlenbergs, 2006, Linear Mixed Models for Longitudinal Data. Springer)
## https://gbiomed.kuleuven.be/english/research/50000687/50000696/geertverbeke/documents/potsdam.pdf
## 
## Como o crescimento craniofacial depende da produção de testosterona?
## Crescimento craniofacial em ratos: efeito do Decapeptyl (inibidor de testosterona) em 3 doses
## Comparação de 3 grupos: controle, dose baixa (de decapeptyl 0.4 mg/kg) e dose alta (de decapeptyl 4 mg/kg)
## Tratamento começa aos 45 dias de idade; medidas tomadas a cada 10 dias, a partir do dia 50
## A variável resposta é o crescimento craniofacial, medido como a distância entre pontos anatômicos bem definidos em radiografias do crânio de cada rato


## Modelo misto
## - dado longitudinal
## - 50 indivíduos em 3 grupos de tratamentos observados até 7 vezes no tempo
## - presença "dropouts"
##
#
## Modelo:
## Y_i = \beta_0 _ b_{0i} + (\beta_1 I(trat=1) + \beta_2 I(trat=2) + \beta_3 I(trat=3) + b_{1i})*t{ij} + \epsilon{e_{ij}}
## 

(rats <- read.table("http://www.leg.ufpr.br/ce092/dados/rats.txt",
                    header = TRUE, row.names = 1, stringsAsFactors = TRUE))
dim(rats)
head(rats)
rats <- transform(rats,
                  time = log(1 + (age - 45)/10),
                  rat = as.factor(rat)
                  )
str(rats)

## fator ordenado (con < low < hig):
rats$trat <- factor(rats$trat, levels = c("con", "low", "hig"))

require(ggplot2)

gg <- ggplot(data = rats) +
    theme_bw() +
    facet_wrap(vars(trat))

ggori <- gg +
    geom_line(aes(x = age, y = growth, group = rat), alpha = 0.3) +
    geom_point(aes(x = age, y = growth, group = rat))
    labs(x = "Idade (dias)", y = "Crescimento craniofacial (pixels)",
         title = "Escala original de tempo")
ggori

ggtransf <- gg +
    geom_line(aes(x = time, y = growth, group = rat), alpha = 0.3) +
    geom_point(aes(x = time, y = growth, group = rat))
    labs(x = "Idade (dias)", y = "Crescimento craniofacial (pixels)",
         title = "Escala transformada de tempo")
ggtransf

ggmedias <- gg +
    stat_summary(aes(x = time, y = growth, group = trat, color = trat),
                   fun = mean, geom = "line", lwd = 1.5) +
    labs(x = "Idade (dias)", y = "Crescimento craniofacial (pixels)",
         title = "Médias por grupo de tratamento")
ggmedias

##
## Ajuste de 3 modelos
##
require(lme4)
## Modelo somente com intecepto aleatório
fit1 <- lmer(growth ~ time:trat + (1|rat), data = rats)
## Modelo com 2 efeitos aleatórios (intercept e inclinação) independentes
fit2Ind <- lmer(growth ~ time:trat + (1|rat) + (0 + time|rat), data = rats)
##fit2Ind <- lmer(growth ~ time:trat + (1 + time||rat), data = rats)
## Modelo com 2 efeitos aleatórios (intercept e inclinação) correlacionados
fit2Cor <- lmer(growth ~ time:trat + (1 + time|rat), data = rats)

## "time:trat" inclui APENAS a interação, sem os efeitos principais
## de 'trat' e 'time'. Isso é equivalente a um modelo com:
##   - intercepto global único (Intercept)
##   - inclinação DIFERENTE por tratamento (não uma diferença,
##     mas uma inclinação separada para cada grupo)
## Interpretação: cada grupo tem sua própria taxa de crescimento
## ("slope" em 'time'), mas compartilham o mesmo intercepto médio.
## Isso é razoável pois os grupos eram equivalentes ao início
## (age=50, antes do efeito do tratamento se consolidar),
## mas pode ser restritivo em outros contextos.

fit1
summary(fit1)
faraway::sumary(fit1)
VarCorr(fit1)
logLik(fit1)

fit2Ind
summary(fit2Ind)
faraway::sumary(fit2Ind)
VarCorr(fit2Ind)
logLik(fit2Ind)

fit2Cor
summary(fit2Cor)
faraway::sumary(fit2Cor)
VarCorr(fit2Cor)
logLik(fit2Cor)

dim(coef(fit1)$rat)
dim(coef(fit2Ind)$rat)
dim(coef(fit2Cor)$rat)

##
## Escolha do modelo
##
fitsML <- list()
fitsML$fit1 <- lmer(growth ~ time:trat + (1 | rat),
                    data = rats, REML = FALSE)
fitsML$fit2Ind <- lmer(growth ~ time:trat + (1 | rat) + (0 + time | rat),
                       data = rats, REML = FALSE)
fitsML$fit2Cor <- lmer(growth ~ time:trat + (1 + time | rat),
                       data = rats, REML = FALSE)

with(fitsML, {
    anova(fit1, fit2Ind, fit2Cor)
})
with(fitsML, {
    AIC(fit1, fit2Ind, fit2Cor)
})
with(fitsML, {
    BIC(fit1, fit2Ind, fit2Cor)
})

## modelo escolhido: fit1
fit <- fit1
rm(fit1, fit2Ind, fit2Cor, fitsML)

summary(fit)
faraway::sumary(fit)
VarCorr(fit)
logLik(fit)

## Visualizar os BLUPs (shrinkage):
lattice::dotplot(ranef(fit, condVar = TRUE))


## Preditos marginais/populacionais (sem efeitos aleatórios)
## método Delta para intervalos de confiança

## predição nos pontos observados
rats$pred <- predict(fit, re.form = NA)
Xpred <- model.matrix(~ time:trat, data = rats)
dim(Xpred)
Vbeta <- as.matrix(vcov(fit))
rats$se <- sqrt(diag(Xpred %*% Vbeta %*% t(Xpred)))
rats <- transform(rats,
                lower = pred - 1.96*se,
                upper = pred + 1.96*se)
head(rats)

ggtransf +
    geom_line(aes(x = time, y = pred, group = trat), data = rats, col = "orange", lwd = 2) +
    geom_ribbon(aes(x = time, ymin = lower, ymax = upper, group = trat), data = rats, col = "orange", fill = "orange", alpha = 0.1)


##options(width = 250)
## predição em um grid
nd <- expand.grid(age = seq(50, 110, by = 5),
                  trat = with(rats, levels(trat)))
nd <- transform(nd, time = log(1 + (age - 45)/10))
dim(nd)
head(nd)
nd$pred <- predict(fit, newdata = nd, re.form = NA) ## marginal por grupos
head(nd, n = 20)

## intervalos de predição via método delta
X_pred <- model.matrix(~ time:trat, data = nd)
V_beta <- as.matrix(vcov(fit))
nd$se <- sqrt(diag(X_pred %*% V_beta %*% t(X_pred)))
nd <- transform(nd,
                lower = pred - 1.96*se,
                upper = pred + 1.96*se)
dim(nd)
head(nd, n = 20)

## alternativa "direta"
fit.eff <- as.data.frame(
    effects::Effect(c("trat", "time"), fit,
                    xlevels = list(time = nd$time)))
fit.eff
fit.eff <- transform(fitInterc.eff,
                     age = 10 * (exp(time) - 1) + 45)
head(fit.eff)

##
gg <- ggplot(data = rats) +
    theme_bw() +
    facet_wrap(vars(trat))

ggori <- gg +
    geom_line(aes(x = age, y = growth, group = rat), alpha = 0.3) +
    geom_point(aes(x = age, y = growth, group = rat))
    labs(x = "Idade (dias)", y = "Crescimento craniofacial (pixels)",
         title = "Escala original de tempo")
ggori

ggtransf <- gg +
    geom_line(aes(x = time, y = growth, group = rat), alpha = 0.3) +
    geom_point(aes(x = time, y = growth, group = rat))
    labs(x = "Idade (dias)", y = "Crescimento craniofacial (pixels)",
         title = "Escala transformada de tempo")
ggtransf
##

## geral
gpred1 <- ggplot(data = nd) +
    theme_bw() +
    facet_wrap(vars(trat)) +
    geom_line(aes(x = time, y = pred, group = trat, color = trat), lwd = 2) +
    geom_ribbon(aes(x = time, ymin = lower, ymax = upper, group = trat, fill = trat),
                alpha = 0.2)

gpred2 <- ggplot(data = nd) +
    theme_bw() +
    facet_wrap(vars(trat)) +
    geom_line(aes(x = age, y = pred, group = trat, color = trat), lwd = 2) +
    geom_ribbon(aes(x = age, ymin = lower, ymax = upper, group = trat, fill = trat),
                alpha = 0.2)

## por tratamento
gpred3 <- ggplot(data = nd) +
    theme_bw() +
    geom_line(aes(x = time, y = pred, group = trat, color = trat),
              lwd = 1.5) +
    geom_ribbon(aes(x = time, ymin = lower, ymax = upper,
                    group = trat, fill = trat),
                alpha = 0.2) +
    scale_color_manual(values = c("blue", "orange", "green")) +
    scale_fill_manual(values = c("blue", "orange", "green"))

gpred4 <- ggplot(data = nd) +
    theme_bw() +
    geom_line(aes(x = age, y = pred, group = trat, color = trat),
              lwd = 1.5) +
    geom_ribbon(aes(x = age, ymin = lower, ymax = upper,
                    group = trat, fill = trat),
                alpha = 0.2) +
    scale_color_manual(values = c("blue", "orange", "green")) +
    scale_fill_manual(values = c("blue", "orange", "green"))

require(patchwork)
gpred1 / gpred2
gpred3 + gpred4



fit0 <- lmer(growth ~ time + (1|rat), data = rats)
anova(fit0, fit)
AIC(fit0, fit)
BIC(fit0, fit)