# packages ----

library(tidyverse)

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

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

# 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 <- read.table("https://stats.idre.ucla.edu/wp-content/uploads/2016/02/tolerance1.txt", sep=",", header=T) # WIDE

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

# Corrrelaciones 2 a 2 entre puntuaciones de tol en 5 ocasiones.

round(cor(tolerance[,2:6]),3)


tolerance.pp$idf <- factor(tolerance.pp$id) #factor!

## Gráficos ----
# Gráfico de puntos
ggplot(data = tolerance.pp) +
  geom_point(aes(x = time, y = tolerance, color = idf)) +
  labs(y = "Tolerance") + 
  theme(legend.position = "none")

# Gráfico de líneas - spaghetti plot
ggplot(data = tolerance.pp) +
  geom_line(aes(x = time, y = tolerance, color = idf)) +
  labs(y = "Tolerance") + 
  theme(legend.position = "none")

# Gráfico de puntos para cada individuo

ggplot(data = tolerance.pp) +
  geom_point(aes(x = time, y = tolerance, color = idf)) +
  facet_wrap(~idf)  +
  labs(y = "Tolerance") + 
  theme(legend.position = "none")

# gráfico de líneas con smooth
# change method smooth
method = 'lm'
#method = NULL
ggplot(data = tolerance.pp) +
  geom_point(aes(x = time, y = tolerance, color = idf)) +
  geom_smooth(aes(x = time, y = tolerance), method="lm", se=FALSE, color="black") +
  facet_wrap(~idf) +
  labs(y = "Tolerance") + 
  theme(legend.position = "none")



## Diferencias cambio entre individuos (covariables) ----


summary(lm(tolerance ~ time, data = tolerance.pp))
coef(lm(tolerance ~ time, data = tolerance.pp))

# Modelo lineal estimado para cada individuo
# estima el modelo para c/individuo y generar vector de betas

indiv <- unique(tolerance.pp$id)
beta0 <- NULL
beta1 <- NULL
for (j in 1:16) {
  beta0[j] <- coef(lm(tolerance ~ time,
                      data = tolerance.pp %>%
                        filter(id==indiv[j])))[1]
  beta1[j] <- coef(lm(tolerance ~ time,
                      data = tolerance.pp %>%
                        filter(id==indiv[j])))[2]
}




# Estadísticas descriptivas de parámetros de trayectorias individuales  

mean(beta0)
var(beta0)
sd(beta0)
cor(beta0,beta1)

# en una tabla

data.frame(beta0=beta0,beta1=beta1) %>%
  summarise(
    int_m =mean(beta0),
    int_d =sd(beta0),
    pen_m =mean(beta1),
    pen_d =sd(beta1),
    cor   =cor(data.frame(beta0=beta0,beta1=beta1))[1,2]) %>%
  mutate_if(is.numeric,~ round(., 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")