##rm(list=ls())
##set.seed(2026)

##
## Exemplo: medidas repetidas com interceptos aleatórios
##
##Tmedia <- 50
##Tsigma2 <- 3
##Tsigma2g <- 3
##Tsigma2g/(Tsigma2g + Tsigma2)
##n.grupos <- 10
##(ni <- rpois(n.grupos, 2.5) + 1)
##medias <- round(rnorm(n.grupos, Tmedia, sqrt(Tsigma2g)), dig = 1)
##(nT <- sum(ni))
##df1 <- data.frame(grupo = factor(rep(1:n.grupos, ni)))
##df1$y <- round(medias[rep(1:n.grupos, ni)] + rnorm(nT, 0, sqrt(Tsigma2)), dig = 1)
##df1
##write.csv(df1, "ModMix1.csv", row.names = FALSE)

rm(list=ls())

df1 <- read.csv("http://www.leg.ufpr.br/ce092/dados/ModMix1.csv", header = TRUE)
head(df1)
Tmedia <- 50          ## média verdadeira
(n.grupos <- with(df1, length(table(grupo))))
(ni <- table(df1$grupo))
##
## Visualização dos dados
##
with(df1, plot.default(grupo, y, pch = 20,
                       xlab = "grupo", ylab = "y"))
abline(h = Tmedia, lty = 3)

(mgeral <- with(df1, mean(y)))
abline(h = mgeral, col = "orange")

(fit1.lm <- lm(y ~ 1, data = df1))
coef(fit1.lm)
abline(h = coef(fit1.lm), lty = 2, col = "red")

## médias por grupo
(mediasgr <- with(df1, tapply(y, grupo, mean)))
(fit1.lmgr <- lm(y ~ -1 + factor(grupo), data = df1))
coef(fit1.lmgr)  ## comparar com mediasgr
points(1:n.grupos, mediasgr, pch = "x", col = "red", cex = 1.5)


##
## Dois pacotes diferentes para modelos lineares mistos: nlme e lme4
## Retornam as mesmas estimativas.
##
(fit1.lme <- nlme::lme(y ~ 1, random = ~ 1 | grupo, data = df1))

library(lme4)
(fit1.lmer <- lmer(y ~ 1 + (1 | grupo), data = df1))
summary(fit1.lmer)
fixef(fit1.lmer)
ranef(fit1.lmer)
coef(fit1.lmer)
coef(fit1.lmer)$grupo
## comparar com coef(fit1.lmer)$grupo:
fixef(fit1.lmer) + ranef(fit1.lmer)$grupo

abline(h = fixef(fit1.lmer), lty = 2, col = "blue")
points(1:n.grupos, coef(fit1.lmer)$grupo[,1], pch = "x", col = "blue", cex = 1.5)

## repetindo comandos para gráfico
with(df1, plot.default(grupo, y, pch = 20))
abline(h = Tmedia, lty = 3)

abline(h = coef(fit1.lm), lty = 2, col = "red")
abline(h = fixef(fit1.lmer), lty = 2, col = "blue")

points(1:n.grupos, mediasgr, pch = "x", col = "red", cex = 1.5)
points(1:n.grupos, coef(fit1.lmer)$grupo[,1], pch = "x", col = "blue", cex = 1.5)

segments(x0 = 1:n.grupos, y0 = mediasgr,
         y1 = coef(fit1.lmer)$grupo[, 1], col = "gray40", lty = 3)
legend("bottomright", bty = "n", cex = 0.8,
       legend = c("média verdadeira", "média geral (lm)",
                  "médias por grupo (lm)", "efeito fixo (lmer)",
                  "preditos por grupo (lmer)"),
       lty = c(3, 2, NA, 2, NA), pch = c(NA, NA, "x", NA, "x"),
       col = c("black", "red", "red", "blue", "blue"))
##
## Extraindo outras informações do ajuste do modelo com lme4
##
VarCorr(fit1.lmer)
VarCorr(fit1.lmer)$grupo

## coeficiente de correlação intracalsse (ICC)
(s2g <- lme4::VarCorr(fit1.lmer)$grupo[1])
(s2 <- sigma(fit1.lmer)^2)
(ICC <- s2g / (s2g + s2))

## [CORRIGIDO] este trecho usava s2g/s2 antes de serem definidos
## (ficava logo após "ranef", mas s2g e s2 só existem daqui pra baixo);
## movido para cá, onde as duas variáveis já estão disponíveis.
##
## encolhimento ("shrinkage"): o predito de cada grupo é a média
## ponderada entre a média do próprio grupo e o efeito fixo (média
## geral), com peso w_i = s2g / (s2g + s2/n_i) -- grupos com poucas
## observações (n_i pequeno) são puxados mais para a média geral.
(w <- s2g / (s2g + s2 / as.vector(ni)))
w * mediasgr + (1 - w) * fixef(fit1.lmer)   ## comparar com coef(fit1.lmer)$grupo

class(fit1.lmer)
isREML(fit1.lmer)        ## default de lmer() é REML
nobs(fit1.lmer)
ngrps(fit1.lmer)
formula(fit1.lmer)
methods(class = "lmerMod") 

summary(fit1.lmer)$coefficients    
## lme4 não reporta valor-p para os efeitos fixos: em dados
## desbalanceados os graus de liberdade do denominador não têm forma
## fechada única. Alternativas:
##   - lmerTest::lmer()  (aproximações de Satterthwaite / Kenward-Roger)
##   - intervalos de confiança (abaixo)
##   - nlme::lme(), que reporta p-valores com g.l. "containment"
summary(fit1.lme)$tTable

##
## Estimativas e intervalos para média geral (populacional)
## Comparando diferentes formas de estimar e obter IC para a média
##

## (1) intervalo "ingênuo" (ignora grupos)
ic1 <- t.test(df1$y)$conf.int ## idêntico ao IC de confint(fit1.lm)

## (2) efeitos fixos, "no pooling": médias de grupo tratadas como
##     parâmetros fixos e distintos; a média geral é a MÉDIA NÃO
##     PONDERADA das n.grupos médias (cada grupo com peso 1/n.grupos,
##     não peso proporcional a n_i).
K <- matrix(rep(1 / n.grupos, n.grupos), nrow = 1)
est2 <- as.numeric(K %*% coef(fit1.lmgr))
se2  <- as.numeric(sqrt(K %*% vcov(fit1.lmgr) %*% t(K)))
ic2  <- est2 + c(-1, 1) * qt(0.975, df.residual(fit1.lmgr)) * se2

## (3) modelo misto (lme4), Wald: rápido, mas assume normalidade
##     assintótica também para os componentes de variância
ic3 <- confint(fit1.lmer, method = "Wald")["(Intercept)", ]

## (4) modelo misto (lme4), perfil da verossimilhança: mais confiável
##     que Wald, sobretudo com poucos grupos
ic4 <- confint(fit1.lmer, method = "profile")["(Intercept)", ]

##
## Organizando tudo em uma única tabela
##
(icTab <- data.frame(
     metodo = c("Ingênuo",
                "Efeitos fixos",
                "Misto, Wald",
                "Misto, perfil"),
     estimativa = c(mean(df1$y), est2, fixef(fit1.lmer), fixef(fit1.lmer)),
     li = c(ic1[1], ic2[1], ic3[1], ic4[1]),
     ls = c(ic1[2], ic2[2], ic3[2], ic4[2])
 ))

ypos <- nrow(icTab):1
opar <- par(mar = c(4, 5, 4, 2))
with(icTab, {
    plot(estimativa, ypos,
         xlim = range(icTab[, c("li", "ls")]) + c(-0.5, 0.5),
         ylim = c(0.5, nrow(icTab) + 0.5),
         yaxt = "n", pch = 19,
         xlab = expression(mu), ylab = "",
         main = "Estimativas e IC(95%) da média geral")
    axis(2, at = ypos, labels = metodo, las = 1, cex.axis = 0.7)
    segments(x0 = li, x1 = ls, y0 = ypos, y1 = ypos, lwd = 2)
    abline(v = Tmedia, lty = 3, col = "black")
    legend("bottomright", bty = "n", cex = 0.8, lty = 3,
           legend = "média teórica")
})
par(opar)

## Algum Diagnóstico
par(mfrow = c(1, 3))
plot(fitted(fit1.lmer), resid(fit1.lmer), pch = 20,
     xlab = "valores ajustados", ylab = "resíduos"); abline(h = 0, lty = 2)
qqnorm(resid(fit1.lmer), main = "resíduos"); qqline(resid(fit1.lmer))
qqnorm(ranef(fit1.lmer)$grupo[, 1], main = "efeitos aleatórios")
qqline(ranef(fit1.lmer)$grupo[, 1])
par(mfrow = c(1, 1))