# Efecto del stroke (ACV) en el deterioro cognitivo
# Datos del English Longitudinal Study of Ageing (ELSA)
load('datos_stroke.RData')

# edad del stroke (para centrar en esa edad)
ed_str<- NULL
ids <- unique(d2$idauniq)
d2$age_stroke <- NA
for (i in seq_along(ids)){
  strk_i <- d2$str[which(d2$idauniq == ids[i])]
  edad_i <- d2$age[which(d2$idauniq == ids[i])]
  k <- min(which(strk_i==1))
  if(k==1) edad<-edad_i[1] else edad<-mean(edad_i[c(k-1,k)],na.rm=TRUE)
  ed_str <- c(ed_str,edad)
  d2$age_stroke[which(d2$idauniq == ids[i])] <- edad
}

# metrica tiempo desde/hasta el stroke
d2$time <- d2$age - d2$age_stroke

# centro age_stroke en 70 para mejorar la interpretacion en el modelo
d2$age_stroke70 <- d2$age_stroke - 70

# decripcion de los datos
library(ggplot2)
ggplot(d2, aes(x = age, y = verbf, group = idauniq)) +
  geom_point(col = gray(0.7)) +
  geom_line(col = gray(0.7)) +
  theme_bw() +
  theme(text = element_text(size = 15)) +
  xlab('Age (years)') +
  ylab('Verbal Fluency')

ggplot(d2, aes(x = time, y = verbf, group = idauniq)) +
  geom_point(col = gray(0.7)) +
  geom_line(col = gray(0.7)) +
  theme_bw() +
  geom_vline(xintercept = 0, linewidth = 2, linetype = 'dashed', color = 'red') +
  theme(text = element_text(size = 15)) +
  xlab('Tiempo desde/hasta el stroke') +
  ylab('Verbal Fluency')

# podriamos filtrar aquellos con menos de 50
d2 |>
  dplyr::filter(edad>50) -> d2

# Numero de visitas de c/participante
library(dplyr)
d2 |>
  select(idauniq, verbf) |>
  filter(!is.na(verbf)) |>
  group_by(idauniq) |> 
  summarise(n_i = n()) |> 
  count(n_i)

library(nlme)
# modelo con metrica temporal = edad
mod0 <- lme(verbf ~ age, 
            random = ~age|idauniq,
            method = 'ML',
            data = na.omit(d2[,c('verbf','age','idauniq')]))

# estimaciones
summary(mod0)

# efectos aleatorios
ranef(mod0) |> 
  as.data.frame() |>
  ggplot(aes(x = `(Intercept)`, y = age)) +
    geom_point(size = 1) +
    theme_bw() +
    geom_vline(xintercept = 0, linetype = 'dashed', size = 0.5) +
    geom_hline(yintercept = 0, linetype = 'dashed', size = 0.5)

# ¿Tienen sentido las estimaciones?

# test sobre los componentes de varianza
# ¿Es necesario incorporar efecto aleatorio para el tiempo?
mod0_H0 <- update(mod0, random = ~1|idauniq)

# una alternativa es utiilzando el AIC
AIC(mod0, mod0_H0)

# otra es haciendo una prueba de hipotesis
library(varTestnlme)
varCompTest(mod0, mod0_H0, pval.comp = 'bounds')



# modelo con metrica temporal = tiempo desde/hasta el stroke
mod1 <- lme(verbf ~ time, 
            random = ~time|idauniq,
            method = 'ML',
            data = na.omit(d2[,c('verbf','time','idauniq')]))

# estimaciones
summary(mod1)

# efectos aleatorios
ranef(mod1) |> 
  as.data.frame() |>
  ggplot(aes(x = `(Intercept)`, y = time)) +
    geom_point(size = 1) +
    theme_bw() +
    geom_vline(xintercept = 0, linetype = 'dashed', size = 0.5) +
    geom_hline(yintercept = 0, linetype = 'dashed', size = 0.5)

# test sobre el efecto aleatorio de la pendiente
mod1_H0 <- update(mod1, random = ~1|idauniq)
varCompTest(mod1, mod1_H0, pval.comp = 'bounds')



# El problema con este modelo es que no incorpora la edad al momento del ACV, 
# no es lo mismo sufrir un ACV a los 60 que a los 80, volvamos al pdf
d2$edad_acv70 <- d2$age_stroke - 70

# modelo con metrica tiempo d/h stroke y con efecto de la edad al stroke
mod2 <- lme(verbf ~ time*edad_acv70, 
            random = ~time|idauniq,
            method = 'ML',
            data = na.omit(d2[,c('verbf','time','idauniq','edad_acv70')]))

# estimaciones
summary(mod2)

# test sobre el efecto aleatorio de la pendiente
mod2_H0 <- update(mod1, random = ~1|idauniq)
varCompTest(mod2, mod2_H0, pval.comp = 'bounds')


# efecto de la edad del ACV
fixef(mod2) -> parametros
times <- seq(-10,10,0.1)
edad_acv70 <- seq(-5,5)

# efectos fijos en funcion de la edad al ACV
pi_0 <- cbind(1, edad_acv70)%*%parametros[c(1, 3)]
pi_1 <- cbind(1, edad_acv70)%*%parametros[c(2, 4)]

# efectos fijos sobre E(Y)
cbind(1,times) %*% t(cbind(pi_0, pi_1)) -> mu_Y

# juntamos todo para graficar
todo <- data.frame(tiempo     = rep(times, length(edad_acv70)),
                   edad_ACV70 = rep(edad_acv70, each = length(times)) + 70,
                   VerbFlu    = as.numeric(mu_Y))
ggplot(todo, aes(x = tiempo, y = VerbFlu, group = edad_ACV70)) +
  geom_line(aes(col = edad_ACV70)) +
  theme_bw() +
  ylab('Verbal Fluency') +
  xlab('Tiempo desde/hasta el stroke') +
  theme(text = element_text(size = 15)) +
  labs(col = 'Edad al ACV')



# modelo con efecto de la edad del acv y del acv
mod3 <- lme(verbf ~ time*(edad_acv70 + str), 
            random = ~time|idauniq,
            method = 'ML',
            data = na.omit(d2[,c('verbf','time','idauniq','edad_acv70','str')]))

# estimaciones
summary(mod3)

# test sobre el efecto aleatorio de la pendiente
mod3_H0 <- update(mod1, random = ~1|idauniq)
varCompTest(mod3, mod3_H0, pval.comp = 'bounds')


# 'str' parece no tener un efecto sobre la tasa de deterioro, quitemos ese termino del modelo
mod4 <- lme(verbf ~ time*edad_acv70 + str, 
            random = ~time|idauniq,
            method = 'ML',
            data = na.omit(d2[,c('verbf','time','idauniq','edad_acv70','str')]))

# estimaciones
summary(mod4)

# si quisieramos comparar los modelos:

# agregamos variables
anova(mod1, mod2, mod3)

# quitamos variables
anova(mod4, mod3)

# o directamente comparamos todo
AIC(mod1, mod2, mod3, mod4)

# A partir del modelo 4:
# Veamos la trayectoria de la fluidez verbal entre 10 años antres y despues del ACV
# para personas que lo sufren a los 60, 65, 70, 75 y 80

# mas facil usando 'predict'
p60 <- predict(mod4, newdata = data.frame(time = times, edad_acv70 = -10, str = ifelse(times<0,0,1)), level = 0)
p65 <- predict(mod4, newdata = data.frame(time = times, edad_acv70 = -5, str = ifelse(times<0,0,1)), level = 0)
p70 <- predict(mod4, newdata = data.frame(time = times, edad_acv70 = 0, str = ifelse(times<0,0,1)), level = 0)
p75 <- predict(mod4, newdata = data.frame(time = times, edad_acv70 = 5, str = ifelse(times<0,0,1)), level = 0)
p80 <- predict(mod4, newdata = data.frame(time = times, edad_acv70 = 10, str = ifelse(times<0,0,1)), level = 0)

# juntamos todo para graficar
todo <- data.frame(tiempo     = rep(times, 5),
                   edad_ACV70 = rep(c(-10,-5,0,5,10), each = length(times)) + 70,
                   VerbFlu    = c(p60, p65, p70, p75, p80))

ggplot(todo, aes(x = tiempo, y = VerbFlu, group = edad_ACV70)) +
  geom_line(aes(col = edad_ACV70)) +
  theme_bw() +
  ylab('Verbal Fluency') +
  xlab('Tiempo desde/hasta el stroke') +
  theme(text = element_text(size = 15)) +
  labs(col = 'Edad al ACV')

# alternativamente
ggplot(todo, aes(x = tiempo, y = VerbFlu)) +
  geom_line() +
  theme_bw() +
  ylab('Verbal Fluency') +
  xlab('Tiempo desde/hasta el stroke') +
  theme(text = element_text(size = 15)) +
  facet_grid(~edad_ACV70)