# packages ----

library(tidyverse)

library(lme4)
library(lmerTest)
library(nlme)

library(viridis)
library(qwraps2) #for mean_ci
library(here)

library(gridExtra)

# themes ----
theme_set(theme_bw())

# data tolerance ----

# Encuesta Nacional sobre la Juventud (NYS). En cada año con 11, 12, 13, 14 y 15 años, se releva sobre la tolerancia hacia comportamientos "desviados". El grado en que se consideran aceptables los comportamientos diferentes o inusuales dentro de un contexto particular.


tolerance.pp <- read.table("https://stats.idre.ucla.edu/wp-content/uploads/2016/02/tolerance1_pp.txt", sep=",", header=T) # LONG

tolerance.pp <- 
  tolerance.pp %>% 
  mutate(
    idf = factor(id),
    group = factor(ifelse(exposure < median(exposure), 0, 1)),
    malef = factor(male)
  )

levels(tolerance.pp$group) <- c("Low exposure","High exposure")
levels(tolerance.pp$malef) <- c("Female","Male")






## Gráficos ----
ggplot(tolerance.pp, aes(time, tolerance, colour=idf)) + 
  geom_point() + 
  facet_wrap(~idf) + 
  geom_smooth(method=lm,se=FALSE, color="black") +
  labs(y = "Tolerance") 

ggplot(tolerance.pp, aes(time, tolerance, group=id, colour=idf)) +
  geom_line() + 
  labs(y = "Tolerance", title = "Spaghetti plot") +
  theme(legend.position = "none")



## Modelo multinivel ----

### Modelo RE intercepto y pendiente incondicional ----
model1_lm<- lm(tolerance ~ time, data=tolerance.pp)
summary(model1_lm)

model1<- lme(tolerance~time, data=tolerance.pp, random= ~time | id, method="ML")
summary(model1)

anova(model1,model1_lm) # RE int&slope vs LM incond

# Extra
model1_lme4 <- lme4::lmer(tolerance ~ time + (1 + time | id), data = tolerance.pp, REML = FALSE)
summary(model1_lme4)


### Modelo RE intercepto y pendiente condicional (sexo) ----
model2_lm<- lm(tolerance ~ time*malef,
               data=tolerance.pp)
summary(model2_lm)

model2<- lme(tolerance~time*malef, data=tolerance.pp, random= ~time | id, method="ML")
summary(model2)

anova(model2,model2_lm) # RE int&slope vs LM cond male

# Extra
model2_lme4 <- lme4::lmer(tolerance ~ time * malef + (1 + time | id), data = tolerance.pp, REML = FALSE)
summary(model2_lme4)


### Modelo RE intercepto y pendiente condicional (sexo y exposure) ----
model3_lm<- lm(tolerance ~ time*malef + time*group,
               data=tolerance.pp)
summary(model3_lm)

model3<- lme(tolerance~time*malef + time*group, data=tolerance.pp, random= ~time | id, method="ML")
summary(model3)

anova(model3,model3_lm) #RE int&slope vs LM cond male&expo

# comparar la inclusión de covariables invariantes
anova(model3,model2) # expo&male vs male 
anova(model3,model1) # expo&male vs incond








### Modelo RE intercepto y pendiente condicional (exposure) ----
model4_lm<- lm(tolerance ~ time*group,
               data=tolerance.pp)
summary(model4_lm)

model4<- lme(tolerance~time*group, data=tolerance.pp, random= ~time | id, method="ML")
summary(model4)

anova(model4,model4_lm) # RE int&slope vs LM cond expo




### Modelo RE solo intercepto condicional (exposure) ----
model4_lm<- lm(tolerance ~ time*group,
               data=tolerance.pp)
summary(model4_lm)

model5<- lme(tolerance~time*group, data=tolerance.pp, random= ~1 | id, method="ML")
summary(model5)

anova(model5,model4_lm) # RE int vs LM cond expo
anova(model4,model5) # RE int&slope vs LM cond expo







### Modelo RE solo intercepto incondicional ----
model1_lm<- lm(tolerance ~ time,data=tolerance.pp)
summary(model1_lm)

model6<- lme(tolerance~time, data=tolerance.pp, random= ~1 | id, method="ML")
summary(model6)

### Modelo RE solo pendiente incondicional ----
model1_lm<- lm(tolerance ~ time, data=tolerance.pp)
summary(model1_lm)

model7<- lme(tolerance~time, data=tolerance.pp, random= ~0 + time | id, method="ML")
summary(model7)


anova(model6,model1) #RE int vs RE int&slope
anova(model7,model1) #RE slope vs int&slope




# Resultados postestimation del modelo incondicional ----
# Veremos el model1 (RE int&slope), model6 (RE solo int) y model7 (RE solo pendiente)
model1<- lme(tolerance~time, data=tolerance.pp, random= ~time | id, method="ML")
summary(model1)

model6<- lme(tolerance~time, data=tolerance.pp, random= ~1 | id, method="ML")
summary(model6)

model7<- lme(tolerance~time, data=tolerance.pp, random= ~0 + time | id, method="ML")
summary(model7)

# Predicción de RE individuales
ranef(model1)
ranef(model6)
ranef(model7)
# Resultados de FE 
fixef(model1)
fixef(model6)
fixef(model7)
# Coeficiones individuales del modelo
coef(model1)
coef(model6)
coef(model7)

ranef(model1) + matrix(fixef(model1), byrow = TRUE, nrow = 16, ncol=2)

cbind(
  ranef(model6) + matrix(fixef(model6)[1], byrow = TRUE, nrow = 16, ncol=1),
  fixef(model6)[2])

cbind(
  coef(model7)[1],
  ranef(model7) + matrix(fixef(model7)[2], byrow = TRUE, nrow = 16, ncol=1))


# Valores ajustados según el modelo (Ygorro)
tolerance.pp$ygorro1 <- predict(model1)
tolerance.pp$ygorro6 <- predict(model6)
tolerance.pp$ygorro7 <- predict(model7)


# El mejor es ???
# Hacemos visualización del incondicional LM, RE int, RE slope, RE int&slope para comparar
anova(model1,model6) 
anova(model1,model7)

# Gráfico de dispersión para c/individuo con sus valores ajustados por LME con RE en int&slope y la recta de regresión

ggplot(data = tolerance.pp) +
  geom_line(aes(time, ygorro1, color = idf)) +
  geom_point(aes(time, tolerance, color = idf)) +
  geom_abline(intercept = fixef(model1)[1], slope = fixef(model1)[2], color = "black") +
  facet_wrap(~ idf) + 
  theme(legend.position = "none") + 
  labs(y = "Tolerance")


# Gráfico de dispersión para c/individuo con sus valores ajustados por LME con RE en int y la recta de regresión

ggplot(data = tolerance.pp) +
  geom_line(aes(time, ygorro6, color = idf)) +
  geom_point(aes(time, tolerance, color = idf)) +
  geom_abline(intercept = fixef(model1)[1], slope = fixef(model1)[2], color = "black") +
  facet_wrap(~ idf) + 
  theme(legend.position = "none") + 
  labs(y = "Tolerance")

# Gráfico de dispersión para c/individuo con sus valores ajustados por LME con RE en int&slope y la recta de regresión

ggplot(data = tolerance.pp) +
  geom_line(aes(time, ygorro7, color = idf)) +
  geom_point(aes(time, tolerance, color = idf)) +
  geom_abline(intercept = fixef(model1)[1], slope = fixef(model1)[2], color = "black") +
  facet_wrap(~ idf) + 
  theme(legend.position = "none") + 
  labs(y = "Tolerance")

# sin wrap los 4 modelos

p <- ggplot(data = tolerance.pp) +
  geom_point(aes(time, tolerance, color = idf)) +
  geom_abline(intercept = coef(model1_lm)[1], slope = coef(model1_lm)[2], color = "black") +
  theme(legend.position = "none") + 
  labs(y = "Tolerance")

p1 <- ggplot(data = tolerance.pp) +
  geom_point(aes(time, tolerance, color = idf)) +
  geom_line(aes(time, ygorro1, color = idf)) +
  theme(legend.position = "none") + 
  labs(y = "Tolerance")

p6 <- ggplot(data = tolerance.pp) +
  geom_point(aes(time, tolerance, color = idf)) +
  geom_line(aes(time, ygorro6, color = idf)) +
  theme(legend.position = "none") + 
  labs(y = "Tolerance")

p7 <- ggplot(data = tolerance.pp) +
  geom_point(aes(time, tolerance, color = idf)) +
  geom_line(aes(time, ygorro7, color = idf)) +
  theme(legend.position = "none") + 
  labs(y = "Tolerance")

gridExtra::grid.arrange(p,p1,p6,p7, nrow=2, ncol=2)







# data sleepstudy ----
# Tiempos de reacción en un estudio de privación del sueño

sleepstudymod <- read_csv(here('data/sleepstudymod.csv'), col_select = 2:5)
datos <- sleepstudymod %>% 
  mutate(SubjectF=factor(Subject),
         FemaleF=factor(Female))
levels(datos$FemaleF) <- c("Male","Female")
  
ggplot(data = datos, aes(x = Days, y = Reaction, color = SubjectF)) +
  geom_point() +
  theme_bw() +
  facet_wrap(~ Subject) + 
  labs(y = "Reaction time") + 
  theme(legend.position = "none")