IDA_SurvLongData.Rmd

---
title: "Untitled"
author: "Yibin Feng"
date: '2022-05-26'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Dependencies

```{r}
rm(list = ls())

library(tidyverse)
library(ggpubr)
library(ggfortify)

library(moments)

library(foreach)

library(survival)
library(survminer)

library(Hmisc)
#library(ADNIMERGE)

source("function_utility.R") # Utility functions

conflicted::conflict_prefer("summarize", "dplyr")
conflicted::conflict_prefer("summarise", "dplyr")
```

# Import data

Set the surv and long data for IDA in the chunk below:

```{r}
load("./data_cleaned/adni_cleaned.RData")

```

list of covariates
```{r}

baseline.covs <- c("AGE", "PTGENDER", "PTEDUCAT", "status.bl", "APOE4") 


# Manual exclusion
vars_manual_remove <- c("TAU", "PTAU", "ABETA")
# Note: Type of TAU, PTAU and ABETA are character
# need to handle non-numerical values first. currently excluded

# Exclude irrelevant variables
vars_irrelevant <- c(
  names(df.long_censored)[grepl(".bl", names(df.long_censored))], # Exclude variables with `.bl` suffix including Years.bl and Months.bl
  "id", "RID",
  "time", "event", "status", "DX", # Survival information
  "VISCODE", "EXAMDATE", "Y", "M", "Month", # Time variables
  "AGE", "age.fup", 
  "COLPROT", "ORIGPROT", "PTID", "SITE", # Visit information
  "PTGENDER", "PTEDUCAT", "PTETHCAT", "PTRACCAT", "PTMARRY", "APOE4", # Baseline variables
  "FSVERSION", "IMAGEUID", "FLDSTRENG" # Metadata for image
)

vars_ignore <- c(vars_manual_remove, vars_irrelevant) # Variables that will not be considered as long covariates

# -----------------------------------------------------------------------------------
# Update values of time-varying covariates in surv data when landmark time > 0
# for methods pCox-bl and pCox-lm

# The original values observed at baseline i.e. VISCODE=="bl"
# are replaced by last observed value on or before landmark

# Step 1: set the covariates to update, should cover the candidate long covariates
vars_long <- names(df.surv_preds)[!(names(df.surv_preds) %in% vars_ignore)]

baseline.covs
vars_long
```


# Table: data format

```{r}
example.ids <- c(5, 41)


example.raw <- ADNIMERGE::adnimerge %>%
  filter(RID %in% example.ids) %>%
  group_by(RID) %>%
  select(RID, VISCODE, EXAMDATE, Years.bl, DX.bl, DX, AGE, PTGENDER, ADAS13, MMSE)
  

example.surv <- df.surv_preds %>%
  filter(RID %in% example.ids) %>%
  select(RID, time, event, status.bl, AGE, PTGENDER, ADAS13, MMSE)

example.long <- df.long_censored %>%
  filter(RID %in% example.ids) %>%
  group_by(RID) %>%
  select(RID, VISCODE, Years.bl, status.bl, AGE, PTGENDER, ADAS13, MMSE)


rm(example.ids)

```

```{r}
title.table <- "Example of raw data"
Hmisc::latex(example.raw, digits = 3,
      title = title.table,
      caption = title.table,
      booktabs = TRUE,
      what = "%",
      size = "small",
      na.blank = FALSE,
      legend.bottom = TRUE)

title.table <- "Example of surv data"
Hmisc::latex(example.surv, digits = 3,
      title = title.table,
      caption = title.table,
      booktabs = TRUE,
      what = "%", 
      size = "small",
      na.blank = FALSE,
      legend.bottom = TRUE)

title.table <- "Example of long data"
Hmisc::latex(example.long, digits = 3,
      title = title.table,
      caption = title.table,
      booktabs = TRUE,
      what = "%", 
      size = "small",
      na.blank = FALSE,
      legend.bottom = TRUE)
```

# Table: Baseline characteristics

```{r}
library(ADNIMERGE) # use adnimerge because it retain all labels

df.tmp_baselinechar <- adnimerge %>% 
  filter(RID %in% df.surv_preds$RID) %>%
  filter(VISCODE == "bl") %>%
  select(AGE,  DX.bl, DX,
         ORIGPROT,
         PTGENDER, PTEDUCAT, PTETHCAT, PTRACCAT, PTMARRY, APOE4)
  #modify_if(is.factor, as.character)
  
s <- summaryM(
  DX + 
    AGE +
    PTGENDER + PTEDUCAT + PTETHCAT + PTRACCAT + PTMARRY + APOE4 ~ ORIGPROT,
  data = df.tmp_baselinechar, 
  overall = TRUE, test = FALSE)

print(s, digits = 3)
```


```{r}
title.table <- "Baseline characteristics of ADNI cohorts"

latex(s, digits = 3,
      title = title.table, 
      what = "%", legend.bottom = TRUE)
```


# Table: cohort follow up time and censoring rate

```{r}
# all phases
summary(df.surv_preds$event)
sd(df.surv_preds$event)
summary(df.surv_preds$time)
sd(df.surv_preds$time)
```


```{r}
# phase specific
lapply(unique(df.surv_preds$ORIGPROT), function(x) {
  df.tmp <- df.surv_preds %>%
    filter(ORIGPROT == x)
  summary <- df.tmp %>%
    select(event, time) %>%
    summary()
  sd.event <- sd(df.tmp$event)
  sd.fup <- sd(df.tmp$time)
  
  list(x, summary, event.sd = sd.event, fup.sd = sd.fup)
})


```

## number of visits 
```{r}
# by original protocol 
df.tmp <- df.long_censored %>%
  group_by(ORIGPROT, id) %>%
  dplyr::summarize(
    visits.subject = n(), .groups = "drop"
  )

df.tmp %>%
  group_by(ORIGPROT) %>%
  dplyr::summarise(
    visits.avg = mean(visits.subject),
    visits.sd = sd(visits.subject),
    visits.median = median(visits.subject),
    visits.max = max(visits.subject),
    visits.min = min(visits.subject)
  )
```


```{r}
# combined
df.tmp <- df.long_censored %>%
  group_by(id) %>%
  dplyr::summarize(
    visits.subject = n(), .groups = "drop"
  )

df.tmp %>%
  dplyr::summarise(
    visits.avg = mean(visits.subject),
    visits.sd = sd(visits.subject),
    visits.median = median(visits.subject),
    visits.max = max(visits.subject),
    visits.min = min(visits.subject)
  )
```

# Table: summary of observed events

```{r}
table.eventrate_lm <- data.frame(do.call(
  "rbind",
  lapply(
    2:6, function(t_l) {
      df.surv_preds %>%
        filter(time > t_l) %>%
        select(status) %>%
        table() %>%
        t()
    }
  )
))


table.eventrate_lm2 <- data.frame(do.call(
  "rbind",
  lapply(
    2:6, function(t_l) {
      df.surv_preds %>%
        filter(time > t_l) %>%
        summarize(
          "Avg time to event after landmark" = round(mean(time), 1) - t_l
        )
    }
  )
))
  
table.eventrate_lm$`Landmark time` <- 2:6

table.eventrate_lm <- table.eventrate_lm %>%
  mutate(
    `Event rate` = round(dementia / (censored + dementia), 3)
  ) %>%
  select("Landmark time", "censored", "dementia", "Event rate")

table.eventrate_lm <- table.eventrate_lm %>%
  cbind(table.eventrate_lm2)


table.eventrate_lm

Hmisc::latex(
  object = table.eventrate_lm,
  title = "Event rate at different landmark times",
  label = "tab:events_lms",
  caption = "Event rate at different landmark times",
  rowlabel = NULL,
  booktabs = TRUE
)

```


# Table: Number of repeated measurements at different landmark times

```{r}
table.rep_measure_lm <- do.call("rbind", lapply(
  2:6, function(t_l) {
  df.long_censored %>%
    filter(time > t_l) %>%
    filter(Years.bl < t_l) %>%
    group_by(id) %>%
    summarize(
      across(all_of(vars_long), function(x) sum(!is.na(x)))
    ) %>%
    summarize(
      across(-id, mean)
    )
  }
))

table.rep_measure_lm <- table.rep_measure_lm %>%
  round(digits = 1)


table.rep_measure_lm <- data.frame(t(table.rep_measure_lm))

table.rep_measure_lm <- data.frame(
  row.names(table.rep_measure_lm), table.rep_measure_lm, 
  row.names = NULL
)

colnames(table.rep_measure_lm) <- c("Longitudinal covariate", 2:6)

table.rep_measure_lm <- table.rep_measure_lm %>%
  arrange(`Longitudinal covariate`)

table.rep_measure_lm

Hmisc::latex(
  object = table.rep_measure_lm,
  title = "Number of repeated measurements before the landmark time",
  label = "tab:rep_measure_lms",
  caption = "Number of repeated measurements before the landmark time",
  rowlabel = NULL,
  booktabs = TRUE
)

```

## Figure: Number of repeated measurements at different landmark times

```{r, fig.height=5, fig.width=4}
vars_long_21 <- c(
  "ADAS11",
  "ADAS13",
  "ADASQ4",
  "CDRSB",  
  "Entorhinal",
  "FAQ",
  'Fusiform',
  "Hippocampus",
  "ICV",
  "LDELTOTAL",
  "MidTemp",
  "MMSE",
  "mPACCtrailsB",
  "mPACCdigit",
  "RAVLT.forgetting",
  "RAVLT.immediate",
  "RAVLT.learning",
  "RAVLT.perc.forgetting",
  "TRABSCOR",
  "WholeBrain",
  "Ventricles"
)


g.rep_measure_lm <-  table.rep_measure_lm %>%
  filter(`Longitudinal covariate` %in% vars_long_21) %>%
  pivot_longer(
    cols = -`Longitudinal covariate`, 
    names_to = "Landmark time", 
    values_to = "Avg repeated measurements"
    ) %>%
  mutate(
    `Landmark time` = factor(as.numeric(`Landmark time`))
  ) %>%
  ggplot() + 
    geom_point(
      aes(
        x = `Avg repeated measurements`,
        y = `Longitudinal covariate`,
        group = interaction(`Landmark time`, `Longitudinal covariate`),
        color = `Landmark time`
      ),
#      position = position_dodge2(0.7), 
      size = 2, alpha = 0.6
    ) +
    scale_x_continuous(breaks = 0:7, limits = c(0, 7)) +
    ggsci::scale_color_lancet() +
    theme_pubr() +
    theme(
      legend.position = "top",
      legend.box = "vertical",
      panel.grid.major.x = element_line(color = "grey80"),
#      panel.grid.minor.x = element_line(color = "grey60"),
      panel.grid.major.y = element_line(color = "grey80")
    )
  

ggsave(
    filename = paste0("repeated_measure_lms.png"),
    plot = g.rep_measure_lm,
    device = "png",
    path = "./figure/",
    width = 7, height = 7,
    bg = "white"
)

```


# Checking prepared datasets

```{r, eval=FALSE}
# double check values of var X.bl and var X at VISCODE == bl

#df.surv_to_check <- df.clean
df.surv_to_check <- df.surv_preds_prc

# select variables containing suffix ".bl" for checking
vars_with_bl <- names(df.surv_to_check)[grepl(".bl", names(df.surv_to_check), fixed=TRUE)]
vars_with_bl <- vars_with_bl[!vars_with_bl %in% c("Years.bl", "Month.bl")] # exclude

# logic check
res.compare <- data.frame(vars_with_bl)
for (x in vars_with_bl){
  x.no_bl <- gsub(".bl", "", x)
  tmp <- df.surv_to_check %>% filter(VISCODE == "bl") %>% select(RID, x, x.no_bl)
  tmp <- tmp %>% mutate(isSame = identical(get(x), get(x.no_bl)))
  res.compare$allSame[which(x == vars_with_bl)] <- all(tmp$isSame)
}

View(res.compare)

# visual check
for (x in vars_with_bl){
  x.no_bl <- gsub(".bl", "", x)
  tmp <- df.surv_to_check %>% filter(VISCODE == "bl") %>% select(RID, x, x.no_bl) 
  g <- tmp %>% ggplot +
    geom_point(aes(get(x), get(x.no_bl))) +
    labs(title = "compare value with var X.bl and var X at VISCODE == bl") +
    xlab(x) + ylab(x.no_bl)
  print(g)
}
```

# Check frequency of different protocol in prepared data

```{r}
# Cross-table of dementia / censoring
print("[Summary] cross table for df.surv_preds:")
table(df.surv_preds %>% select(ORIGPROT, event))
```

# Check missing values


## Figure: missing in long
```{r, eval=FALSE}
# Check missing data

# Difficult to read
missingness.surv <- df.surv_preds %>%
  select(all_of(baseline.covs)) %>%
  naniar::vis_miss(sort_miss = TRUE)

missingess.long <- df.long_censored %>%
  select(all_of(vars_long)) %>%
  naniar::vis_miss(sort_miss = TRUE)

# adjust margin
missingess.long <- missingess.long +
  theme(plot.margin = margin(t = 0, r = 70, b = 0, l = 0, unit = "pt"))
```

```{r}
# save missing plot for baseline cov
#ggsave(
#  filename = "missingness_surv.png", plot = missingness.surv, device = "png", path = "./figure/", width = 8)

# save missing plot for long covs
ggsave(
  filename = "missingness_long.png", 
  plot = missingess.long, 
  device = "png", 
  path  = "./figure/", 
  width = 8, height = 6,
  bg = "white")
```


### NAs in repeated measurement

#### Table: missing in long
```{r}
load(file = "./output/temp/output_folds_template_scenario2_seed721_lm2_transformed.RData")

missing_proportions <- folds[[1]]$missing_proportions
colnames(missing_proportions) <- c("Covariate", "Proportion of subjects without any observations", "Average number of observation per subject")
missing_proportions[,2:3] <- missing_proportions[,2:3] %>%
  round(digits = 2)


title.table <- "Missingness in longitudinal data"
latex(missing_proportions, 
      digits = 3,
      title = title.table,
      caption = title.table,
      booktabs = TRUE,
      what = "%",
      size = "small",
      na.blank = FALSE,
      legend.bottom = TRUE)
```


#### Number of candidate long covariates
```{r}
# Obtain long covariates
# Manual exclusion ----------------------------------------------------------------------------------------
# Type of TAU, PTAU and ABETA are character, need to handle non-numerical values first. currently excluded
vars_manual_remove <- c("TAU", "PTAU", "ABETA")

# Exclude irrelevant variables
vars_irrelevant <- c(
  names(data.long)[grepl(".bl", names(data.long))], # Exclude baselines, Years.bl, Month.bl
  "RID", "time", "event", "status", "DX",
  "VISCODE", "EXAMDATE", "Y", "M", "Month",
  "AGE", "age.fup", "COLPROT", "ORIGPROT", "PTID", "SITE",
  "PTGENDER", "PTEDUCAT", "PTETHCAT", "PTRACCAT", "PTMARRY", "APOE4",
  "FSVERSION", "IMAGEUID", "FLDSTRENG"
)

y.names <- names(data.long)[!(names(data.long) %in% c("id", vars_manual_remove, vars_irrelevant))]

# ----------------------------------------------------------------------------------------



missing_proportions_all <- Compute_missing_proportions(
  data.long = df.long, 
  vars_filter = c(vars_irrelevant, vars_manual_remove),
  is_include = FALSE # Set FALSE to exclude vars_filter
)

# Recalculate the missing proportions when landmark at year 2
missing_proportions_yr2 <- Compute_missing_proportions(
  data.long = df.long %>% filter(Years.bl <= 2), 
  vars_filter = c(vars_irrelevant, vars_manual_remove),
  is_include = FALSE # Set FALSE to exclude vars_filter
)

missing_proportions_compare <- missing_proportions_all %>%
  left_join(missing_proportions_yr2, by = "var_name", suffix = c("_all", "_yr2")) %>%
  mutate(diff = proportion_y_missing_yr2 - proportion_y_missing_all)


missing_proportions_compare %>%
  mutate(round(across(-var_name), digits = 3))
# This comparison shows that selecting landmark at 2 does not affect the pool of candidate long covariates
# When we set the cut-off to 0.1, we select 21 long covariates.
# But for some covariates, the proportion did increase, meaning more subjects are without any observation
```





#### Cognitive assessment

```{r, fig.width=12}
check_covariates <- c("ADAS13", "MMSE", "RAVLT.immediate", "RAVLT.learning", "FAQ")

plots <- lapply(check_covariates, function(feature){
  df.long %>%
  group_by(RID) %>%
  dplyr::summarize(
    n.visits = n(),
    n.missing = sum(is.na(get(feature))),
    p.missing = n.missing / n.visits,
    time = unique(time)
    ) %>%
  ggplot() +
    geom_histogram(aes(p.missing), binwidth = 0.05) +
    coord_cartesian(xlim = c(0, 1)) +
    labs(
      title = paste("Percentage of missing value(s) in", feature),
      x = "Ratio of NAs to # visits per subject"
    ) + 
    theme(
      title = element_text(size = 8)
    )
})

ggarrange(plotlist = plots)
```



```{r, fig.width=12}
check_covariates <- c("ADAS13", "MMSE", "RAVLT.immediate", "RAVLT.learning", "FAQ")

plots <- lapply(check_covariates, function(feature){
  df.long %>%
  group_by(RID) %>%
  dplyr::summarize(
    n.visits = n(),
    n.missing = sum(is.na(get(feature))),
    p.missing = n.missing / n.visits,
    time = unique(time),
    protocol = unique(COLPROT)
    ) %>%
  ggplot() + 
    geom_jitter(aes(time, p.missing, color = protocol),
                alpha = 0.2) +
    coord_cartesian(ylim = c(0, 1)) +
    labs(
      title = paste("Covariate =", feature),
      x = "Time to event (year)",
      y = "Ratio of NAs to # visits per subject"
    ) + 
    theme(
      title = element_text(size = 8)
    )
})

p <- ggarrange(plotlist = plots, common.legend = TRUE)
annotate_figure(
  p, 
  top = text_grob("Per subject NAs ratio against follow-up time by data collection protocol", 
                  face = "bold")
  )
```

#### Imaging

```{r, fig.width=12}
check_covariates <- c("Ventricles", "Hippocampus", "WholeBrain", "Entorhinal", "Fusiform", "MidTemp")

plots <- lapply(check_covariates, function(feature){
  df.long %>%
  group_by(RID) %>%
  dplyr::summarize(
    n.visits = n(),
    n.missing = sum(is.na(get(feature))),
    p.missing = n.missing / n.visits,
    time = unique(time)
    ) %>%
  ggplot() +
    geom_histogram(aes(p.missing), binwidth = 0.05) +
    coord_cartesian(xlim = c(0, 1)) +
    labs(
      title = paste("Percentage of missing value(s) in", feature),
      x = "Ratio of NAs to # visits per subject"
    ) + 
    theme(
      title = element_text(size = 8)
    )
})

ggarrange(plotlist = plots)
```


```{r, fig.width=12}
check_covariates <- c("Ventricles", "Hippocampus", "WholeBrain", "Entorhinal", "Fusiform", "MidTemp")

plots <- lapply(check_covariates, function(feature){
  df.long %>%
  group_by(RID) %>%
  dplyr::summarize(
    n.visits = n(),
    n.missing = sum(is.na(get(feature))),
    p.missing = n.missing / n.visits,
    time = unique(time),
    protocol = unique(COLPROT)
    ) %>%
  ggplot() + 
    geom_jitter(aes(time, p.missing, color = protocol),
                alpha = 0.2) +
    coord_cartesian(ylim = c(0, 1)) +
    labs(
      title = paste("Covariate =", feature),
      x = "Time to event (year)",
      y = "Ratio of NAs to # visits per subject"
    ) + 
    theme(
      title = element_text(size = 8)
    )
})

p <- ggarrange(plotlist = plots, common.legend = TRUE)
annotate_figure(
  p, 
  top = text_grob("Per subject NAs ratio against follow-up time by data collection protocol", 
                  face = "bold")
  )
```

#### ABETA TAU PTAU

```{r, fig.width=12}
check_covariates <- c("ABETA", "TAU", "PTAU")

plots <- lapply(check_covariates, function(feature){
  df.long %>%
  group_by(RID) %>%
  dplyr::summarize(
    n.visits = n(),
    n.missing = sum(is.na(get(feature))),
    p.missing = n.missing / n.visits,
    time = unique(time)
    ) %>%
  ggplot() +
    geom_histogram(aes(p.missing), binwidth = 0.05) +
    coord_cartesian(xlim = c(0, 1)) +
    labs(
      title = paste("Percentage of missing value(s) in", feature),
      x = "Ratio of NAs to # visits per subject"
    ) + 
    theme(
      title = element_text(size = 8)
    )
})

ggarrange(plotlist = plots)
```

```{r, fig.width=12}
check_covariates <- c("ABETA", "TAU", "PTAU")

plots <- lapply(check_covariates, function(feature){
  df.long %>%
  group_by(RID) %>%
  dplyr::summarize(
    n.visits = n(),
    n.missing = sum(is.na(get(feature))),
    p.missing = n.missing / n.visits,
    time = unique(time),
    protocol = unique(ORIGPROT)
    ) %>%
  ggplot() + 
    geom_jitter(aes(time, p.missing, color = protocol),
                alpha = 0.2) +
    coord_cartesian(ylim = c(0, 1)) +
    labs(
      title = paste("Covariate =", feature),
      x = "Time to event (year)",
      y = "Ratio of NAs to # visits per subject"
    ) + 
    theme(
      title = element_text(size = 8)
    )
})

p <- ggarrange(plotlist = plots, common.legend = TRUE)
annotate_figure(
  p, 
  top = text_grob("Per subject NAs ratio against follow-up time by data collection protocol", 
                  face = "bold")
  )
```


```{r, fig.width=12}
check_covariates <- c("ABETA", "TAU", "PTAU")

plots <- lapply(check_covariates, function(feature){
  df.long %>%
  group_by(RID) %>%
  dplyr::summarize(
    n.visits = n(),
    n.missing = sum(is.na(get(feature))),
    p.missing = n.missing / n.visits,
    time = unique(time),
    protocol = unique(COLPROT)
    ) %>%
  ggplot() + 
    geom_jitter(aes(time, p.missing, color = protocol),
                alpha = 0.2) +
    coord_cartesian(ylim = c(0, 1)) +
    labs(
      title = paste("Covariate =", feature),
      x = "Time to event (year)",
      y = "Ratio of NAs to # visits per subject"
    ) + 
    theme(
      title = element_text(size = 8)
    )
})

p <- ggarrange(plotlist = plots, common.legend = TRUE)
annotate_figure(
  p, 
  top = text_grob("Per subject NAs ratio against follow-up time by data collection protocol", 
                  face = "bold")
  )
```

#### Summary

```{r}
df.long %>%
  summarytools::descr() %>%
  t() %>%
  as.data.frame() %>%
  mutate(Pct.Valid = round(Pct.Valid / 100, 3)) %>% # Convert from [0,100] to [0,1]
  select(N.Valid, Pct.Valid) %>% 
  arrange(Pct.Valid)
```


# Transformation

## Table: Transformation

```{r}
load(file = "./output/temp/output_folds_template_scenario2_seed721_lm2_transformed.RData")

missing_proportions <- folds[[1]]$missing_proportions
colnames(missing_proportions) <- c("Covariate", "Proportion of subjects without any observations", "Average number of observation per subject")
missing_proportions[,2:3] <- missing_proportions[,2:3] %>%
  round(digits = 2)

table.transformation <- transformation_summary %>%
  filter(var_name %in% missing_proportions$Covariate[1:21]) %>%
  select(var_name, transformation, skewness_t, skewness)
table.transformation[3:4] <- table.transformation[3:4] %>%
  round(digits = 3)
colnames(table.transformation) <- c("Covariate", "Transformation", "Skewness (After)", "Skewness (Before)")

table.transformation <- table.transformation %>%
  arrange(`Skewness (Before)`)

table.transformation <- table.transformation %>%
  select(Covariate, `Skewness (Before)`, Transformation, `Skewness (After)`)

table.transformation$Transformation[is.na(table.transformation$Transformation)] <- "-"

table.transformation$`Skewness (After)`[table.transformation$Transformation == "-"] <- "-"

table.transformation

title.table <- "Transformation of longitudinal covariates"
Hmisc::latex(table.transformation, 
      digits = 3,
      title = title.table,
      caption = title.table,
      booktabs = TRUE,
      what = "%",
      size = "small",
      na.blank = FALSE,
      legend.bottom = TRUE)
```


```{r}
# Obtain long covariates
# Manual exclusion ----------------------------------------------------------------------------------------
# Type of TAU, PTAU and ABETA are character, need to handle non-numerical values first. currently excluded
vars_manual_remove <- c("TAU", "PTAU", "ABETA")

# Exclude irrelevant variables
vars_irrelevant <- c(
  names(data.long)[grepl(".bl", names(data.long))], # Exclude baselines, Years.bl, Month.bl
  "RID", "time", "event", "status", "DX",
  "VISCODE", "EXAMDATE", "Y", "M", "Month",
  "AGE", "age.fup", "COLPROT", "ORIGPROT", "PTID", "SITE",
  "PTGENDER", "PTEDUCAT", "PTETHCAT", "PTRACCAT", "PTMARRY", "APOE4",
  "FSVERSION", "IMAGEUID", "FLDSTRENG"
)

y.names <- names(data.long)[!(names(data.long) %in% c("id", vars_manual_remove, vars_irrelevant))]

```

```{r}
# Testing function for transformation


y.stats <- foreach(y = y.names, .combine = rbind) %do% {
  c(skewness(data.long[, y], na.rm = TRUE), fivenum(data.long[, y], na.rm = TRUE))
} %>% 
  data.frame(row.names = NULL)


colnames(y.stats) <- c("skewness", "minimum", "lower_hinge", "median", "upper_hinge", "maximum")
y.stats$var_name <- y.names

y.stats <- y.stats %>% 
  mutate(symmetry = (upper_hinge - median) / (median - lower_hinge)) %>%
  select(var_name, skewness, symmetry, everything())


threshold.sym <- 0.2
threshold.skew <- 0.5

# Check for values < 0 => if true, add positive constant before power transformation
# the power transformations are not monotone
# that is, not order preserving—if there are both positive and negative values 
y.transf <- y.stats %>% 
  mutate(is_neg = minimum < 0,
         start = -floor(minimum),
         is_sym_over = symmetry > 1 + threshold.sym,
         is_sym_under = symmetry < 1 - threshold.sym,
         is_skew_pos = skewness > 0 + threshold.skew,
         is_skew_neg = skewness < 0 - threshold.skew
         ) %>%
  select(var_name, is_neg, start, is_sym_over, is_sym_under, is_skew_pos, is_skew_neg, everything())


View(y.transf)
```

```{r}
# for thinking about the choice of threshold
y.transf %>%
  filter(!is.infinite(symmetry)) %>%
  ggplot(aes(x = skewness, y = symmetry)) +
    geom_point(alpha = 0.5) +
    geom_hline(yintercept = 1, color = "red", alpha = 0.5) + 
    geom_vline(xintercept = 0, color = "red", alpha = 0.5)
```


```{r, fig.width=4, fig.height=2}
y.hist_original <- lapply(y.names, function(y) {
  g <- data.long %>%
    ggplot() + 
      geom_histogram(aes(get(y)), na.rm = TRUE) +
      geom_vline(xintercept = 0, color = "red") + # Zero
      geom_vline(xintercept = y.stats$minimum[y.stats$var_name == y], color = "blue") + # Minimum
      xlab(y) +
      ggtitle(paste(y, "before transform"))
  
  return(g)
})


y.hist_original
```



```{r}

data.long %>% 
  filter(RAVLT.forgetting < -10)

data.long %>% 
  filter(RAVLT.perc.forgetting < -500)

data.long %>% filter(id == 6073)

# new issue: why is there negative RAVLT forgetting and with very large magnitude...
# already found in adnimerge
# Q1: is negative valid?
# Q2: is very negative valid?

```






# Figure: Subject at risk

```{r}
landmark <- 0:15

res.status <- do.call(
  rbind, 
  lapply(landmark, function(x) {
  df <- df.surv_preds %>%
    filter(time <= x) %>%
    group_by(status) %>%
    summarise(n = n())
  df$t <- x
  return(df)
}))
res.status$status <- res.status$status %>% as.character()

res.at_risk <- do.call(
  rbind, 
  lapply(landmark, function(x) {
  df <- df.surv_preds %>%
    filter(time > x) %>%
    summarise(n = n())
  df$t <- x
  df$status <- "at risk"
  return(df)
}))

res.status <- res.status %>% add_row(res.at_risk)

g.status_by_time <- res.status %>%
  ggplot(aes(x = t, y = n, group = status, color = status)) +
    geom_line(aes(linetype = status), size = 0.8) + geom_point() +
    scale_x_continuous(breaks = landmark) +
    scale_y_continuous(breaks = 0:16 * 100) +
    scale_color_manual(values = c("#00BA38", "#F8766D", "#619CFF")) +
    xlab("Years since baseline (t)") + ylab("Number of subjects") + 
    labs(
      color = "Status",
      linetype = "Status"
    ) +
    ggtitle("Cumulative number of subjects") +
  theme_pubclean()

g.status_by_time

ggsave(filename = "./figure/status_by_time.png", plot = g.status_by_time, 
       width = 8, height = 5)
```



```{r}
g.hist_survtime <- df.surv_preds %>%
  select(time, status) %>%
  ggplot() + 
  geom_histogram(aes(time, group = status, fill = status),
                 binwidth = 0.5, ) +
  scale_x_continuous(breaks = 0:17) + 
  scale_y_continuous(breaks = 0:8 * 50) +
  ggtitle("Distribution of survival time") +
  labs(fill = "Type") + 
  xlab("Years since baseline (t)") + ylab("Number of subjects") +
  theme_pubclean()

g.hist_survtime

ggsave(filename = "./figure/hist_survtime.png", plot = g.hist_survtime, 
       width = 8, height = 5)
```

```{r}
# Output in two panel format
g.combine <- ggarrange(plotlist = list(g.status_by_time, g.hist_survtime), 
          ncol = 2, labels = "AUTO", 
          widths = c(6, 6), heights = c(4, 4))
g.combine
ggsave(filename = "summary_outcome_time.png", plot = g.combine, device = "png",
       path = "./figure/", width = 8, height = 4)
```


```{r}
landmark <- 0:15

res.status <- do.call(
  rbind, 
  lapply(landmark, function(x) {
  df <- df.surv_preds %>%
    filter(time > x) %>%
    group_by(status) %>%
    summarise(n = n())
  df$t <- x
  return(df)
}))
res.status$status <- res.status$status %>% as.character()

res.at_risk <- do.call(
  rbind, 
  lapply(landmark, function(x) {
  df <- df.surv_preds %>%
    filter(time > x) %>%
    summarise(n = n())
  df$t <- x
  df$status <- "at risk in total"
  return(df)
}))

res.status <- res.status %>% add_row(res.at_risk)


g <- res.status %>%
  ggplot(aes(x = t, y = n, group = status, color = status)) +
    geom_line(aes(linetype = status), size = 0.8) + geom_point() +
    scale_x_continuous(breaks = landmark) +
    scale_y_continuous(breaks = 0:16 * 100) +
    xlab("Years since baseline (t)") + ylab("Number of subjects") +
    labs(
      color = "Type",
      linetype = "Type"
    ) +
    ggtitle("Number of subjects at risk at time t, by survival outcome") +
  theme_bw()

g

ggsave(filename = "./figure/status_beyond_time.png", plot = g, width = 8, height = 5)

```
# Figure: Proportions of survival outcome and subjects at risk

```{r}
# proportion of subject at risk by survival outcome

landmark <- 0:15

res.status <- do.call(
  rbind, 
  lapply(landmark, function(x) {
  df <- df.surv_preds %>%
    filter(time > x) %>%
    group_by(status) %>%
    summarise(n = n())
  df$t <- x
  return(df)
}))
res.status$status <- res.status$status %>% as.character()

res.at_risk <- do.call(
  rbind, 
  lapply(landmark, function(x) {
  df <- df.surv_preds %>%
    filter(time > x) %>%
    summarise(n = n())
  df$t <- x
  df$status <- "at risk in total"
  return(df)
}))

#res.status <- res.status %>% add_row(res.at_risk)

res.status <- res.status %>%
  left_join(res.at_risk[, c("n", "t")], by = "t", suffix = c("", ".total"))

res.status$r <- res.status$n / res.status$n.total

max.y.axis <- 1250 # scale the second axis
interval.y.axis <- 5

g.landmark_proportion <- res.status %>%
  ggplot(aes(x = t)) +
    geom_bar(aes(y = n, group = status, fill = status), 
             width = 0.8, alpha = 0.3,
             stat = "identity", position = position_dodge()) +
    geom_line(aes(y = r*max.y.axis, group = status, 
                  color = status, linetype = status), 
              size = 0.8) + 
    geom_point(aes(y = r*max.y.axis, group = status, 
                  color = status)) +
    scale_x_continuous(name = "Landmark time (year)", breaks = landmark) +
    scale_y_continuous(
      name = "Number of subjects",
      breaks = 0:interval.y.axis * max.y.axis / interval.y.axis ,
      sec.axis = sec_axis(
        trans = ~./max.y.axis,
        name = "Proportion",
        breaks = 0:interval.y.axis / interval.y.axis # range 0,1
      )
    ) +
    labs(
      color = "Proportion",
      linetype = "Proportion",
      fill = "Count"
    ) +
    ggtitle("Subjects at risk") +
  theme_pubclean()

g.landmark_proportion
```


```{r}
# proportion of events observed

landmark <- 0:15

res.status <- do.call(
  rbind, 
  lapply(landmark, function(x) {
  df <- df.surv_preds %>%
    filter(time <= x) %>%
    group_by(status) %>%
    summarise(n = n())
  df$t <- x
  return(df)
}))
res.status$status <- res.status$status %>% as.character()

res.at_risk <- do.call(
  rbind, 
  lapply(landmark, function(x) {
  df <- df.surv_preds %>%
    filter(time <= x) %>%
    summarise(n = n())
  df$t <- x
  df$status <- "at risk in total"
  return(df)
}))

#res.status <- res.status %>% add_row(res.at_risk)

res.status <- res.status %>%
  left_join(res.at_risk[, c("n", "t")], by = "t", suffix = c("", ".total"))

res.status$r <- res.status$n / res.status$n.total

max.y.axis <- 1250 # scale the second axis
interval.y.axis <- 5

g.observed_proportion <- res.status %>%
  ggplot(aes(x = t)) +
    geom_bar(aes(y = n, group = status, fill = status), 
             width = 0.8, alpha = 0.3,
             stat = "identity", position = position_dodge()) +
    geom_line(aes(y = r*max.y.axis, group = status, 
                  color = status, linetype = status), 
              size = 0.8) + 
    geom_point(aes(y = r*max.y.axis, group = status, 
                  color = status)) +
    scale_x_continuous(name = "Years since baseline", breaks = landmark) +
    scale_y_continuous(
      name = "Number of subjects",
      breaks = 0:interval.y.axis * max.y.axis / interval.y.axis ,
      sec.axis = sec_axis(
        trans = ~./max.y.axis,
        name = "Proportion",
        breaks = 0:interval.y.axis / interval.y.axis # range 0,1
      )
    ) +
    labs(
      color = "Proportion",
      linetype = "Proportion",
      fill = "Count"
    ) +
    ggtitle("Cumulative observed outcomes") +
  theme_pubclean()

g.observed_proportion
```



```{r}
# Output in two panel format
g.combine <- ggarrange(plotlist = list(g.observed_proportion, g.landmark_proportion), 
          ncol = 2, labels = "AUTO", 
          widths = c(6, 6), heights = c(4, 4),
          common.legend = TRUE, legend = "bottom"
          )
g.combine
ggsave(filename = "summary_outcome_proportion.png", plot = g.combine, device = "png",
       path = "./figure/", width = 10, height = 4)
```

# TEXT: baseline covariate

```{r}
df.surv_preds %>%
  group_by(PTGENDER) %>%
  dplyr::summarise(
    event = sum(event),
    total = n(),
    rate = event / total
  )

df.surv_preds %>%
  group_by(APOE4) %>%
  dplyr::summarise(
    event = sum(event),
    total = n(),
    rate = event / total
  )


df.surv_preds %>%
  group_by(status.bl) %>%
  dplyr::summarise(
    event = sum(event),
    total = n(),
    rate = event / total
  )
```


# Survival curves

should use adjusted for age instead. fit cox model.

- Subjects must stay in the same group over their entire observation time i.e. non-time-dependent
- msm routine estimates time until progress occured. survfit estimates time until progress was observed.
- Aalen-Johansen method used by survfit does not account for interval censoring.

## [fig:thesis] KM cumhaz

https://rpkgs.datanovia.com/survminer/
```{r}
list.cumhaz <- vector(mode = "list", length = 3)

list.cumhaz[[1]] <- ggsurvplot(
    survfit(Surv(time, status=="dementia") ~ PTGENDER,
                     data = df.surv_preds),
    fun = "cumhaz",
    xlab = "Follow up time in years", title = "Gender",
    break.time.by = 3,
    ncensor.plot = FALSE, ncensor.plot.height = 0.25,
    conf.int = TRUE,
    risk.table = FALSE, cumevents = TRUE,
    palette = "aaas")

list.cumhaz[[2]] <- ggsurvplot(
    survfit(Surv(time, status=="dementia") ~ APOE4,
                     data = df.surv_preds),
    fun = "cumhaz",
    xlab = "Follow up time in years", title = "No. of APOE4 allele",
    break.time.by = 3,
    ncensor.plot = FALSE, ncensor.plot.height = 0.25,
    conf.int = TRUE,
    risk.table = FALSE, cumevents = TRUE,
    palette = "lancet")

list.cumhaz[[3]] <- ggsurvplot(
    survfit(Surv(time, status=="dementia") ~ status.bl,
                     data = df.surv_preds),
    fun = "cumhaz",
    xlab = "Follow up time in years", title = "Baseline diagnosis",
    break.time.by = 3,
    ncensor.plot = FALSE, ncensor.plot.height = 0.25,
    conf.int = TRUE,
    risk.table = FALSE, cumevents = TRUE,
    palette = "lancet")

g <- survminer::arrange_ggsurvplots(
  x = list.cumhaz,
  ncol = 2, nrow = 2
)

ggsave(
  filename = "cumhaz.png",
  plot = g, 
  device = "png",
  path = "./figure/",
  width = 12, height = 14
)

# Stratified version
# list.cumhaz[[1]] <- ggsurvplot(
#     survfit(Surv(time, status=="dementia") ~ PTGENDER,
#                      data = df.surv_preds),
#     fun = "cumhaz", facet.by = "status.bl",
#     conf.int = TRUE,
#     risk.table = FALSE, cumevents = TRUE,
#     legend.labs = c("Female", "Male"), 
#     palette = "lancet")
# 
# list.cumhaz[[2]] <- ggsurvplot(
#     survfit(Surv(time, status=="dementia") ~ APOE4,
#                      data = df.surv_preds),
#     fun = "cumhaz", facet.by = "status.bl",
#     conf.int = TRUE,
#     risk.table = FALSE, cumevents = TRUE,
#     palette = "lancet")
# 
# list.cumhaz[[3]] <- ggsurvplot(
#     survfit(Surv(time, status=="dementia") ~ status.bl,
#                      data = df.surv_preds),
#     fun = "cumhaz", facet.by = "APOE4",
#     conf.int = TRUE,
#     risk.table = FALSE, cumevents = TRUE,
#     palette = "lancet") + 
#   ggtitle("Cumulative hazard ")

```

## types of curve

```{r}
# Check visually
ggsurvplot(survfit(Surv(time, status == "dementia") ~ 1,
                   data = df.surv_preds))

ggsurvplot(survfit(Surv(time, status == "dementia") ~ 1,
                   data = df.surv_preds),
           fun = "pct") # Surv prob in percentage

ggsurvplot(survfit(Surv(time, status == "dementia") ~ 1,
                   data = df.surv_preds),
           fun = "cumhaz") # Cumulative hazard
```

## by gender

```{r, fig.height=6}
x <- "PTGENDER"

ggsurvplot(
  survfit(Surv(time, status=="dementia") ~ get(x),
                   data = df.surv_preds),
  fun = "cumhaz",
  conf.int = TRUE,
  risk.table = TRUE,
  legend.labs = c("Female", "Male"),
  palette = "lancet")

rm(x)
```

## by APOE4

```{r, fig.height=6}
# Known genetic risk factor
ggsurvplot(survfit(Surv(time, status=="dementia") ~ APOE4,
                   data = df.surv_preds),
           conf.int = TRUE,
           risk.table = TRUE)
```





```{r, fig.height=8}
ggsurvplot(survfit(Surv(time, status=="dementia") ~ status.bl,
                   data = df.surv_preds),
           conf.int = TRUE,
           risk.table = TRUE,
           title = "KM curve by baseline diagnosis")
```

## [fig:slide] KM surv

```{r, fig.width=8, fig.height=4.5}
ggsurvplot(
    survfit(Surv(time, status=="dementia") ~ status.bl,
                     data = df.surv_preds),
    fun = "pct",
    xlab = "Follow up time in years",
    break.time.by = 2,
#    ncensor.plot = FALSE, ncensor.plot.height = 0.25,
    conf.int = TRUE,
    risk.table = TRUE, cumevents = FALSE, 
    palette = "lancet", cumcensor = FALSE,
    legend = "right", 
    ggtheme = theme_pubclean(),
    tables.theme = theme_survminer()
)
```

```{r, fig.width=10, fig.height=5}
ggsurvplot(
    survfit(Surv(time, status=="dementia") ~ status.bl,
                     data = df.surv_preds),
    fun = "pct",
    xlab = "Follow up time in years",
    break.time.by = 2,
#    ncensor.plot = FALSE, ncensor.plot.height = 0.25,
    conf.int = TRUE,
    risk.table = FALSE, cumevents = TRUE, 
    palette = "lancet", cumcensor = TRUE)

```

# Progression of ADAS13 compared to diagnosis

```{r, eval=FALSE}
# to do: add spaghetti
Plot.progression <- function(df=df.clean, rid, var.name){
  
  subplots <- list(
    df %>% filter(RID == rid) %>% 
      ggplot() + geom_point(aes(M, DX)) + 
      labs(title = "DX over visits"),
    df %>% filter(RID == rid) %>% 
      ggplot() + geom_point(aes(M, get(var.name))) + 
      labs(title = paste(var.name, "over visits"),
           y = var.name)
    )
  
  ggarrange(plotlist = subplots)
  
}
```


```{r}
# Plot.progression(df.clean, 30, "ADAS13")
# Plot.progression(df.clean, 2, "ADAS13")
```

# Distribution of ADAS13

```{r}
df.clean %>% ggplot() +
  geom_violin(aes(DX, ADAS13,
                  fill = DX))
```

```{r}
df.clean %>% ggplot() +
  geom_violin(aes(DX, EcogPtTotal,
                  fill = DX))
```

```{r, fig.width=12}
df.long %>% 
  ggplot() +
    geom_line(aes(Years.bl, ADAS13, group = RID, color = status),
              alpha=0.4) +
    ggtitle("Spaghetti chart of ADAS13") +
    facet_wrap(vars(DX.bl), nrow = 1) + 
    theme(
      legend.position = "top"
    )
```


# Visualize longitudinal data

### Literature covariates

```{r, fig.width=12}
df.long <- df.long_censored
# Spaghetti plots of selected covariates
# Broken lines due to NAs
plot_covariates <- c("ADAS13", "MMSE", "RAVLT.immediate", "RAVLT.learning", "FAQ")

plots <- lapply(plot_covariates, function(feature){
  df.long %>%
    filter(time > Years.bl) %>% # Display measurement before event only i.e. time to event > follow up time
    ggplot() +
      geom_line(aes(M, get(feature), group = RID, color = status),
                alpha = 0.4) +
      facet_wrap(vars(status.bl), nrow = 1) + 
      labs(
        title = paste("Spaghetti chart of", feature, "over time by baseline diagnosis"),
        x = "Years after baseline",
        y = feature,
        color = "Outcome"
      ) +
      theme(
        legend.position = "bottom"
      )
})

plots

```

```{r, fig.width=12}
# Spaghetti plots of selected covariates
# Broken lines due to NAs
plot_covariates <- c("ADAS13", "MMSE", "RAVLT.immediate", "RAVLT.learning", "FAQ")

plots <- lapply(plot_covariates, function(feature){
  df.long_censored %>%
    filter(time > Years.bl) %>% # Display measurement before event only i.e. time to event > follow up time
    ggplot() +
      geom_line(aes(M, get(feature), group = RID, color = status),
                alpha = 0.4) +
      facet_wrap(vars(status.bl), nrow = 1) + 
      labs(
        title = paste("Spaghetti chart of", feature, "over time by baseline diagnosis"),
        x = "Years after baseline",
        y = feature,
        color = "Outcome"
      ) +
      theme(
        legend.position = "bottom"
      )
})

plots

```

#### Random subset of subjects

```{r, fig.width=12}

df.surv <- df.surv_preds
df.long <- df.long_censored
# Spaghetti plots of selected covariates
# Broken lines due to NAs
set.seed(721)

size <- 100
random_RID <- df.surv$RID %>% sample(size)

plot_covariates <- c("ADAS13", "MMSE", "RAVLT.immediate", "RAVLT.learning", "FAQ")

plots <- lapply(plot_covariates, function(feature){
  df.long %>%
    filter(RID %in% random_RID) %>%
    filter(time > Years.bl) %>% # Display measurement before event only i.e. time to event > follow up time
    ggplot() +
      geom_line(aes(Years.bl, get(feature), group = RID, color = status),
                alpha = 0.8) +
      facet_wrap(vars(status.bl), nrow = 1) + 
      labs(
        title = paste("Spaghetti chart of", feature, "over time by baseline diagnosis for", size, "random samples"),
        x = "Years after baseline",
        y = feature,
        color = "Outcome"
      ) +
      theme(
        legend.position = "bottom"
      )
})

plots

```
### [fig:thesis] example of long covs

```{r, fig.width=8, fig.height=4}
# Spaghetti plots of selected covariates
# Broken lines due to NAs
set.seed(721)

size <- 50
random_RID <- df.surv_preds$RID %>% sample(size)

plot_covariates <- c("ADAS13", "MMSE", "RAVLT.immediate", "RAVLT.learning", "FAQ", "MidTemp")

list.g.trajectories <- lapply(plot_covariates, function(feature){
  df.long_censored %>%
    filter(RID %in% random_RID) %>%
    filter(time > Years.bl) %>% # Display measurement before event only i.e. time to event > follow up time
    ggplot() +
      # Base plot all subjects
      geom_line(aes(Years.bl, get(feature), group = RID), 
                alpha = 0.04,
                data = df.long_censored) + 
      geom_line(aes(Years.bl, get(feature), group = RID, linetype = status, color = status),
                alpha = 0.7) +
      geom_point(aes(Years.bl, get(feature), group = RID, shape = status, color = status),
                alpha = 1) +
      #scale_y_continuous(limits = c(0, NA)) + 
      facet_wrap(vars(status.bl), nrow = 1) + 
      labs(
        title = feature,
#        title = paste("Spaghetti chart of", feature, "over time by baseline diagnosis for", size, "random subjects"),
#        caption = "Baseline diagnosis - cognitive normal (left panel), mild cognitive impairment (right panel)",
        x = "Years after baseline",
        y = feature,
        shape = "Outcome",
        color = "Outcome",
        linetype = "Outcome"
      ) +
      theme(
        legend.position = "bottom"
      ) +
      theme_pubr()
})
```

save outputs, show only 1 covariate
```{r}
ggsave(
  filename = "spaghetti_separate.png",
  plot = list.g.trajectories[[1]], device = "png", path = "./figure/", width = 10, height = 4, bg = "white"
)
```


save outputs, show only 2 covariates, combined plots
```{r}
g <- ggpubr::ggarrange(
  plotlist = list(list.g.trajectories[[1]], list.g.trajectories[[6]]), 
  ncol = 1, widths = 10,
  labels = "AUTO"
)

ggsave(
  filename = "spaghetti_combined.png",
  path = "./figure/", 
  plot = g, device = "png", 
  width = 10, height = 12, 
  bg = "white"
)
```

### [fig:slide] ADAS13
```{r, fig.width=8, fig.height=4}
# Spaghetti plots of selected covariates
# Broken lines due to NAs
set.seed(721)

size <- 50
random_RID <- df.surv_preds$RID %>% sample(size)

plot_covariates <- c("ADAS13", "MMSE", "RAVLT.immediate", "RAVLT.learning", "FAQ", "MidTemp")

list.g.trajectories <- lapply(plot_covariates, function(feature){
  df.long_censored %>%
    filter(RID %in% random_RID) %>%
    filter(time > Years.bl) %>% # Display measurement before event only i.e. time to event > follow up time
    ggplot() +
      # Base plot all subjects
      geom_line(aes(Years.bl, get(feature), group = RID), 
                alpha = 0.08,
                data = df.long_censored) + 
      geom_line(aes(Years.bl, get(feature), group = RID, linetype = status, color = status),
                alpha = 0.5,
                size = 0.7) +
      geom_point(aes(Years.bl, get(feature), group = RID, shape = status, color = status),
                alpha = 1) +
      #scale_y_continuous(limits = c(0, NA)) +
      scale_x_continuous(breaks = 0:8 * 2) +
      ggsci::scale_color_lancet() +
      facet_wrap(vars(status.bl), nrow = 1) + 
      labs(
        title = feature,
#        title = paste("Spaghetti chart of", feature, "over time by baseline diagnosis for", size, "random subjects"),
#        caption = "Baseline diagnosis - cognitive normal (left panel), mild cognitive impairment (right panel)",
        x = "Years after baseline",
        y = feature,
        shape = "Outcome",
        color = "Outcome",
        linetype = "Outcome"
      ) +
      theme(
        legend.position = "bottom"
      ) +
      theme_pubr()
})
```

```{r, fig.width=7.5, fig.height=3}
list.g.trajectories[[1]] + 
  theme(title = element_blank())
```



### Imaging covariates


```{r, fig.width=12}
# Spaghetti plots of selected covariates
# Broken lines due to NAs

plot_covariates <- c("Ventricles", "Hippocampus", "WholeBrain", "Entorhinal", "Fusiform", "MidTemp")

plots <- lapply(plot_covariates, function(feature){
  df.long %>%
    filter(time > Years.bl) %>% # Display measurement before event only i.e. time to event > follow up time
    ggplot() +
      geom_line(aes(Years.bl, get(feature), group = RID, color = status),
                alpha = 0.8) +
      facet_wrap(vars(status.bl), nrow = 1) + 
      labs(
        title = paste("Spaghetti chart of", feature, "over time by baseline diagnosis"),
        x = "Years after baseline",
        y = feature,
        color = "Outcome"
      ) +
      theme(
        legend.position = "bottom"
      )
})

plots

```


```{r, fig.width=12}
# Spaghetti plots of selected covariates
# Broken lines due to NAs
set.seed(721)
random_RID <- df.surv$RID %>% sample(100)

plot_covariates <- c("Ventricles", "Hippocampus", "WholeBrain", "Entorhinal", "Fusiform", "MidTemp")

plots <- lapply(plot_covariates, function(feature){
  df.long %>%
    filter(RID %in% random_RID) %>%
    filter(time > Years.bl) %>% # Display measurement before event only i.e. time to event > follow up time
    ggplot() +
      geom_line(aes(Years.bl, get(feature), group = RID, color = status),
                alpha = 0.8) +
      facet_wrap(vars(status.bl), nrow = 1) + 
      labs(
        title = paste("Spaghetti chart of", feature, "over time by baseline diagnosis"),
        x = "Years after baseline",
        y = feature,
        color = "Outcome"
      ) +
      theme(
        legend.position = "bottom"
      )
})

plots

```

### ABETA TAU PTAU

```{r, fig.width=12}
# Spaghetti plots of selected covariates
# Broken lines due to NAs

plot_covariates <- c("ABETA", "TAU", "PTAU")

plots <- lapply(plot_covariates, function(feature){
  df.long %>%
    filter(time > Years.bl) %>% # Display measurement before event only i.e. time to event > follow up time
    ggplot() +
      geom_line(aes(Years.bl, get(feature), group = RID, color = status),
                alpha = 0.8) +
      facet_wrap(vars(status.bl), nrow = 1) + 
      labs(
        title = paste("Spaghetti chart of", feature, "over time by baseline diagnosis"),
        x = "Years after baseline",
        y = feature,
        color = "Outcome"
      ) +
      theme(
        legend.position = "bottom"
      )
})

plots

```

```{r, fig.width=12}
# Spaghetti plots of selected covariates
# Broken lines due to NAs
set.seed(721)
random_RID <- df.surv$RID %>% sample(100)

plot_covariates <- c("ABETA", "TAU", "PTAU")

plots <- lapply(plot_covariates, function(feature){
  df.long %>%
    filter(RID %in% random_RID) %>%
    filter(time > Years.bl) %>% # Display measurement before event only i.e. time to event > follow up time
    ggplot() +
      geom_line(aes(Years.bl, get(feature), group = RID, color = status),
                alpha = 0.8) +
      facet_wrap(vars(status.bl), nrow = 1) + 
      labs(
        title = paste("Spaghetti chart of", feature, "over time by baseline diagnosis"),
        x = "Years after baseline",
        y = feature,
        color = "Outcome"
      ) +
      theme(
        legend.position = "bottom"
      )
})

plots

```



```{r}
transformation_summary
```


# Text: Age at event

```{r}
summary(df.surv_preds$AGE)
```


```{r}
df.temp <- df.surv_preds %>%
  mutate(event.age = time + AGE)

median(df.temp$event.age[df.temp$event == 0])
median(df.temp$event.age[df.temp$event == 1])

mean(df.temp$event.age[df.temp$event == 0])
mean(df.temp$event.age[df.temp$event == 1])


# max min of age with event
max(df.temp$event.age[df.temp$event == 1])
min(df.temp$event.age[df.temp$event == 1])

# follow up time for the max age with event
df.temp$time[df.temp$event.age == max(df.temp$event.age[df.temp$event == 1])]

df.temp%>%
  group_by(status) %>%
  dplyr::summarise(
    event.age.mean = mean(event.age),
    base.age.mean = mean(AGE)
  )
```

time to event
```{r}
mean(df.temp$time[df.temp$event == 0])
mean(df.temp$time[df.temp$event == 1])
```

# Text: drop out due to death

```{r}
ADNIMERGE::treatdis %>% 
#  filter(RID %in% df.surv_preds$RID) %>% 
  nrow() # Check n size

# Count death frequency
sapply(
  ADNIMERGE::treatdis %>% 
    filter(RID %in% df.surv_preds$RID) %>%  # include the final data only
    select(WDREASON), 
  function(reason){
   grepl("Death", reason)
  }) %>% sum()
```