##
## Exemplo: mostrando regressões lineares com interceptos aleatórios e efeito visual 
##          ignorando o efeito de grupo.
rm(list=ls())

df2 <- read.csv("http://www.leg.ufpr.br/ce092/dados/ModMix2.csv", header = TRUE)
df2$grupo <- factor(df2$grupo)  

n.grupos <- with(df2, length(table(grupo)))
table(df2$grupo)

with(df2, plot(x, y))
with(df2, cor(x, y))

## regressão linear simples
fit2.lm <- lm(y ~ x, data = df2)
abline(fit2.lm, col = "red", lty = 2)
summary(fit2.lm)

## um outro ajuste:
abline(c(52.72, 1.287), col = "blue", lty = 2)

## Qual o melhor ajuste? ... será?

df2
for (i in 1:n.grupos) {
    with(subset(df2, grupo == i), lines(x, y))
}

##
## ajustando modelo misto
##
## intercepto aleatório por grupo
##   y_ij = beta_0 + b_i + beta_1 x_ij + e_ij,
##   b_i ~ N(0, s2g)  indep. de  e_ij ~ N(0, s2)
##

## pacote nlme
(fit2.lme <- nlme::lme(y ~ x, random = ~ 1 | grupo, data = df2))
summary(fit2.lme)
rm(fit2.lme)

## pacote lme4
require(lme4)
fit2.lmer1 <- lmer(y ~ x + (1 | grupo), data = df2)
summary(fit2.lmer1)

coef(fit2.lmer1)
coef(fit2.lmer1)$grupo
fixef(fit2.lmer1)
ranef(fit2.lmer1)
fixef(fit2.lmer1)[1] + ranef(fit2.lmer1)$grupo

df2$predM0 <- predict(fit2.lmer1, re.form = NA)  ## predição marginal
df2$predC0 <- predict(fit2.lmer1)                ## predição condicional
df2

## reta ajustada para cada grupo
with(df2, plot(x, y, pch = 20))
abline(fixef(fit2.lmer1), col = "blue", lwd = 2)
###with(df2, lines(x, predM0, col = "blue", lwd = 2)) ## x não ordenado

for (i in 1:n.grupos) {
    with(subset(df2, grupo == i), {
        lines(x, y);
        lines(x, predC0, col = "blue", lty = 2)
    })
}

## agora ajustando interceptos e inclinações aleatórios (ind.) por grupo
##   (x || grupo)  <=>  (1 | grupo) + (0 + x | grupo)
fit2.lmer2 <- lmer(y ~ x + (x || grupo), data = df2)

summary(fit2.lmer2)
fixef(fit2.lmer2)
ranef(fit2.lmer2)
ranef(fit2.lmer2)$grupo[,1]
coef(fit2.lmer2)
coef(fit2.lmer2)$grupo

cbind(
    fixef(fit2.lmer2)[1] + ranef(fit2.lmer2)$grupo[,1],
    fixef(fit2.lmer2)[2] + ranef(fit2.lmer2)$grupo[,2]
)

df2$predM1 <- predict(fit2.lmer2, re.form = NA)  ## predição marginal
df2$predC1 <- predict(fit2.lmer2)                ## predição condicional

## reta ajustada para cada grupo
with(df2, plot(x, y, pch = 20))
abline(fixef(fit2.lmer1), col = "blue", lwd = 2)
abline(fixef(fit2.lmer2), col = "darkorange", lwd = 2)

for (i in levels(df2$grupo)) {
    with(subset(df2, grupo == i), {
        lines(x, y)
        lines(x, predC0, col = "blue", lty = 2)
        lines(x, predC1, col = "darkorange", lty = 2)
    })
}

##
## Interceptos e inclinações aleatórios correlacionados
##
fit2.lmer3 <- lmer(y ~ x + (x | grupo), data = df2)
summary(fit2.lmer3)

## CUIDADO:
## verificar convergência / ajuste singular antes de interpretar
sapply(list(fit2.lmer1 = fit2.lmer1,
            fit2.lmer2 = fit2.lmer2,
            fit2.lmer3 = fit2.lmer3),
       isSingular)

##
## Comparando modelos
##
## (a) Ajuste com REML:
##     válido aqui porque os quatro modelos têm a MESMA parte fix
##     e diferem só na estrutura aleatória.
##     lm() aceita o argumento REML = TRUE em logLik().
## (b) Ajuste com ML:
##     usado quando a parte FIXA diferir entre modelos.
##     refitML() reajusta um modelo lme4 por máxima verossimilhança.
##     Opção neste cso: ajustes com ML para escolher modelo,
##            ajuste final com REML do modelo escolhido
##
c(logLik(fit2.lm, REML = TRUE),
  logLik(fit2.lmer1),
  logLik(fit2.lmer2),
  logLik(fit2.lmer3)
)

## [CORRIGIDO] 'fit2.lmer' não existe -- o objeto foi renomeado para
## 'fit2.lmer1' acima; as duas linhas abaixo ainda usavam o nome antigo.
c(attr(logLik(fit2.lm), "df"),
  attr(logLik(fit2.lmer1), "df"),
  attr(logLik(fit2.lmer2), "df"),
  attr(logLik(fit2.lmer3), "df"))

c(AIC(fit2.lm),
  AIC(fit2.lmer1),
  AIC(fit2.lmer2),
  AIC(fit2.lmer3))

## anova() para classe merMod reajusta como ML
## [CORRIGIDO] a ordem importa: com fit2.lm (classe lm) em primeiro,
## o R despacha para anova.lm, que não sabe lidar com objetos merMod
## (S4) e quebra. Com um merMod em primeiro lugar, anova.merMod assume
## e compara todos corretamente (reajustando por ML).
anova(fit2.lmer1, fit2.lmer2, fit2.lmer3, fit2.lm)


## Comentários:

## - verificar convergência
## - Interpretação final dos resultados
## - REML vs ML
## - Testes para valores proximos da "beirada" do espaço paramétrico
## - Resíduos e diagnósticos: verificação de pressupostos do modelo
## - Interpretações: efeitos fixos, efeitos aleatórios, predições marginais e condicionais
##