---
title: "Simulations Sylvestre 2009"
author: "Alexis van STRAATEN"
date: "`r Sys.Date()`"
output:
  pdf_document:
    latex_engine: xelatex
  html_document: default
---

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

# Sylvestre 2009 - supplementary material

```{r, echo=FALSE, message=FALSE, warning=FALSE}

# Packages

library(ggplot2)
library(dplyr)
library(purrr)
library(survival)
library(WCE)
library(parallel)

```

*1. Design of simulations and data generation*

To validate our model and investigate its properties, we simulated a hypothetical prospective
cohort study of the association between exposure to a single drug and time to an adverse
event. The cohort consisted of 500 new users of the drug, and time zero was defined as the
first day of the drug use [1]. Individuals could interrupt and resume their treatment repeatedly
thereafter. The follow-up was limited to one year. In the next three sub-sections, we describe:
(i) the generation of a matrix of individual drug use patterns, that was kept fixed across the
simulations; (ii) the selection of alternative weight functions and parameters used in different
data-generation scenarios; and (iii) the generation of event times conditional on the weighted
cumulative doses.

```{r}
# Pre-processing
n = 500 # subjects
m = 365 # days

# Generate the matrix of covariate, in a 'long' format.
Xmat = matrix(ncol = 1,
              nrow = n * m)
```

*1.1. Generation of time-dependent exposures*

The drug treatment was assumed to vary in dose and duration both
between individuals and over time within an individual. Since
individuals could interrupt and then resume the treatment repeatedly, we
generated consecutive periods of drug use, each followed by a period
during which the treatment was interrupted. The duration, in weeks, of
the initial treatment for individual i, was generated from a lognormal
distribution with mean of log(0.5) and standard deviation of log(0.8),
and rounded up to the next week. At the end of that week, the subject
was assumed to start the first period of interruption, when the subject
was assumed not to be exposed to the drug. The duration of the first
interruption was generated from the same lognormal distribution. Then,
the subsequent, alternating periods of use and interruptions were
generated similarly, until the end of follow-up was reached. For each
uninterrupted period of use, the standardized daily dose [1] was assumed
to remain constant across the entire period. The daily dose $X(t)$ in
equation (2), was randomly assigned the values of 0.5, 1, 1.5, 2, 2.5,
or 3, with the recommended daily dose corresponding to 1.


```{r}

# Function to generate an individual time-dependent exposure history
# e.g. generate prescriptions of different durations and doses.

TDhist <- function(m){
  
  # Duration : lognormal distribution(0.5,0.8)
  duration <-  7 + 7*round(rlnorm(1, meanlog = 0.5, sdlog = 0.8), 0) # in weeks
  
  # Dose : random assigment of values 0.5, 1, 1.5, 2, 2.5 and 3
  dose <-  sample(seq(from = 0.5, to = 3, by = 0.5), size = 1)
  
  # Start with drug exposure
  vec <- rep(dose, duration)
  
  # Repeat until the vector larger than m
  while (length(vec) <= m){
    intermission <- 7 + 7*round(rlnorm(1, meanlog = 0.5, sdlog = 0.8), 0) # in weeks
    duration <-  7 + 7*round(rlnorm(1, meanlog = 0.5, sdlog = 0.8), 0) # in weeks
    dose <-  sample(c(0.5, 1, 1.5, 2, 2.5, 3), size = 1)
    vec <- append(vec, c(rep(0, intermission), rep(dose, duration)))}
  return(vec[1:m])
  }


# Create time-dependent variable - Dose
Xmat[,1] <- do.call("c", lapply(1:n, function(i) TDhist(m)))
dim(Xmat) <- c(365, 500) 
```


*1.2. Weight functions*

We considered 6 different scenarios, each corresponding to a different true weight function,
each defined over a $[0,1]$ interval, corresponding to one year. As in (4), each weight function
described how the relative weight of past exposure on the risk at current time u changed
with increasing time $(u − t)$ elapsed since the exposure. Alternatively, each function can be
interpreted as showing how the impact of current exposure changes with increasing time since
exposure.*
*The six weight functions are shown as white curves in the six panels of Figure 1. The first
two scenarios assumed that weights decreased monotonically as time since exposure increased.
For scenario 1, we used an exponential decay function with $w(u-t) = 7e^{-7(u-t)}$ while for
scenario 2, we used a Bi-linear function $w(u-t) = 1- \frac{u-t}{50/365}$ for $u-t \leq \frac{50}{365}$
and 0 otherwise.
Scenarios 3 and 4 corresponded to non-monotonic functions where the weights first increased
and then decreased. The Early peak function of scenario 4 corresponded to the density of a
N [0.04; 0.05] distribution, while the Inverted U function for scenario 3, to the density function
of a N [0.2; 0.06] distribution, both left-truncated at t = 0. In comparison with the Inverted
U, in the Early peak scenario, the maximum weight was assigned to more recent doses and
the weights declined more sharply afterward. Scenario 5 corresponded to a Constant weight function with $w(u-t) = \frac{1}{180}$ for $0 \leq u-t \leq \frac{180}{365}$ , so that the resulting WCE was in fact a standard un-weighted time-dependent cumulative dose variable $\sum X(t)$ calculated over the
previous 6-month period. In the scenario 6, we considered a weight function that initially
increased, reached a plateau at around u − t = 180/365 and then started to decrease to reach
0 at $u-t = \frac{240}{365}$. This function was labeled Hat and was specifically designed to investigate
the impact of imposing an incorrect a priori constraint on the weight function. Indeed, in the
analysis we constrained all weight functions to reach 0 at $u-t = \frac{180}{365}$ where the Hat function
was still at its maximum value. To enhance comparability, all the six true weight functions
were standardized so that the area under each function summed up to 1 over the interval
[u, u − 180/365].


```{r}
# Scénario 1 - Exponential
exponential_function <- function(u_t){7*exp(-7*u_t)}

# Scénario 2 - Bi-Linear
bi_linear_function <- function(u_t){ifelse((u_t) <= (50/365),
                                           1-((u_t)/(50/365)),
                                           0)}

## Scénario 3 - Early peak
early_peak_function <- do.call("c", lapply(c(1:365)/365, dnorm, 0.04, 0.05)) / 365

## Scénario 4 - Inverted U
inverted_u_function <- do.call("c", lapply(c(1:365)/365, dnorm, 0.2, 0.06)) / 365

## Scénario 5 - Constant
constant_function <- function(u_t){ifelse(u_t >= 0 & u_t <= 180/365,
                                          1/180,
                                          0)}


scenario <- c("exponential_function",
             # "bi_linear_function",
              "early_peak_function",
              "inverted_u_function",
              "constant_function")
```


*1.3. Events generation*

We generated event times conditional on the WCE, using the permutational algorithm,
specifically developed and validated for simulating event times conditional on time-dependent
covariates and/or effects [2, 3]. The permutational algorithm involves three major steps: (i)
generating individual vectors of time-dependent covariate values; (ii) generating the event
times, as well as random censoring times, from the pre-specified marginal distributions,
independent of covariates; (iii) matching individual event times with individual covariate
vectors based on pre-specified assumptions about the covariates impact on the hazard [4].
Matching at step (iii) is performed so that probability of a subject, who remained at risk until
time t, with the time-dependent covariate vector $X_i(t)$, being matched with the event at time
t is proportional to the subjects current hazard $h(t|X_i(t))$ [2, 3]. A detailed description of the
algorithm and its validation can be found in [3].
Our implementation of the permutational algorithm to generate each simulated sample
involved the following steps. First, for a given scenario with a pre-specified true weight
function, we calculated the 500 individual vectors $WCE_i(u)$, of the true values of the time-
dependent variable WCE in (2), for each day of potential follow-up $0 < u < 365$. To this
end, for each subject i = 1, .., 500, we used the individual vector of the time-dependent daily
drug exposure $X_i(t)$, 0 < t < 365 days (generated in Section 1.1), and the corresponding
weight function $w(ut)$ (as described in Section 1.2). 

```{r}
# Function to obtain WCE vector
wce_vector <- function(u, scenario, Xmat){
  
  t <- 1:u
  
  if(scenario %in% c("exponential_function", "bi_linear_function", "constant_function")){
    scenario_function <- do.call(scenario, list((u - t)/365))
  }else if(scenario == "early_peak_function"){
    scenario_function <- early_peak_function[(u - t + 1)]
  }else if(scenario == "inverted_u_function"){
    scenario_function <- inverted_u_function[(u - t + 1)]
  }
  
  wce <- scenario_function * Xmat[t,]
  
  if(u == 1){
    res <- wce
  }else{
    res <- apply(wce, 2, sum)
  }
  
  return(res)
}

```


```{r}
# Example with scenario 1 - Exponential decay function
wce_mat <- do.call("rbind", lapply(1:365, wce_vector, scenario = "exponential_function", Xmat = Xmat))
wce_mat[,1:10] %>%
  matplot() %>%
  title("WCE curves for scénario 1 - Exponential decay function")
```


We then generated N = 500 independent
event times $\tau_i$, independent of the exposure, assuming their marginal distribution is uniform
U [0, 365] days. We also generated N = 500 random right censoring times $C_i$, independent
of exposure and event times, from uniform U [0, 730] days distribution which was selected
to obtain an approximately 50% rate of right random censoring. Next, we determined the
individual observed follow-up time $T_i = min(\tau_i,C_i)$ and censoring/event status at the end
of his/her follow-up $\delta_i = I(T_i = \tau_i)$, which resulted in about 250 (un-censored) events per
simulated sample.  


```{r}
# Function to generate event times and censoring times
event_censor_generation <- function(){
  
  # Event times : Uniform[1;365] for all scenarios
  eventRandom <- round(runif(n, 1, 365), 0)
  
  # Censoring times : Uniform[1;730] for all scenarios
  censorRandom <- round(runif(n, 1,730),0)
  
  return(list(eventRandom = eventRandom,
              censorRandom = censorRandom))
}
```


The final step of the permutational algorithm involved matching each of the N individual
event or censoring times $T_i$ with one of the N time-dependent WCE vectors. To this end, we
first ranked the observed times $T_i$, i = 1, ..., 500 in increasing order: $T_i-1 < T_i < T_{i +1}$. We
then proceeded iteratively from the earliest to the latest time, by matching each consecutive
time with one of the available at risk WCE vectors, and then removing the matched WCE
vector from all subsequent risk sets [4, 3]. If the observed time $T_i$ corresponded to a censored
observation $(\delta_i = 0)$, we matched it with one among the WCE vectors in the corresponding
risk set, using simple random sampling, with all at risk vectors assigned equal probabilities
of being matched, independent of the individual exposure history. If, however, the observed
time $T_i$ corresponded to an un-censored event $(\delta_i = 1)$, we employed weighted sampling,
with probability of selecting a specific vector $WCE(u)$ we employed weighted sampling,
with probability of selecting a specific vector $T_i$ [4, 3]. Specifically, the probability of sampling a vector $WCE_s(T_i)$, from the corresponding risk set $R_i$ , to be matched with the event at time $T_i$ , was calculated as following: $P_s(T_i) = \frac{exp(ln(4))WCE_s(T_i)}{\sum exp(ln(4))WCE_p(T_i)}$


```{r}
# Function for 'the final step of the permutational algorithm'
matching_algo <- function(wce_mat){
  
  n_patient <- ncol(wce_mat)
  events_generation <- event_censor_generation()
  
  df_event <- data.frame(patient      = 1:n_patient,
                         eventRandom  = events_generation$eventRandom,
                         censorRandom = events_generation$censorRandom)
  
  df_event <- df_event %>%
    group_by(patient) %>%
    mutate(FUP_Ti = min(eventRandom, censorRandom)) %>%
    mutate(event  = ifelse(FUP_Ti == eventRandom, 1, 0)) %>%
    ungroup() %>%
    arrange(FUP_Ti)
  
  # init
  patient_order <- df_event$patient
  j = 1
  id <- 1:n_patient
  wce_mat_df <- wce_mat %>% as.data.frame()
  matching_result <- data.frame()
  
  # Iterative matching, start with the lowest FUP
  for (i in patient_order) {
    
    event      <- df_event[j, "event"]  %>% pull()
    time_event <- df_event[j, "FUP_Ti"] %>% pull()
    
    if(event == 0){
      
      # If no event, all probabilities are the same
      sample_id <- sample(id, 1)
      
    }else if(event == 1){
      
      # If event, matching with different probabilities 
      wce_matrix <- wce_mat_df %>% select(paste0("V", id)) %>% as.matrix()
      proba <- (4 * wce_matrix[time_event,]) / sum(4 * wce_matrix[time_event,])
      sample_id <- sample(id, 1, prob = proba)
    }
    
    matching_result <- rbind(matching_result, 
                             data.frame(id_patient = i,
                                        id_dose_wce = sample_id))
    
    id <- id[!id %in% sample_id]
    j = j + 1
    
    # Stop when last id of iterative algo
    if(length(id) == 1){
      
      matching_result <- rbind(matching_result,
                               data.frame(id_patient  = patient_order[n_patient],
                                          id_dose_wce = id))
      
      return(list(matching_result = matching_result,
                  df_event        = df_event,
                  patient_order   = patient_order))
    }
  }
}
  
# Function to render dataset after the matching algo
get_dataset <- function(Xmat, wce_mat){
  
  df_wce <- data.frame()
  Xmat_df <- Xmat %>%
    as.data.frame()
  
  matching_result <- matching_algo(wce_mat)
  
  for (i in matching_result$patient_order) {
    
    fu <- matching_result$df_event[matching_result$df_event$patient == i, "FUP_Ti"] %>% pull()
    event_patient <- matching_result$df_event[matching_result$df_event$patient == i, "event"] %>% pull()
    
    if(event_patient == 1){
      event_vec <- c(rep(0, (fu-1)),1)
    }else{
      event_vec <- rep(0, fu)
    }
    
    id_dose <- matching_result$matching_result[matching_result$matching_result$id_patient == i, "id_dose_wce"]
    
    df_dose <- data.frame(patient = rep(i, fu),
                          start = 0:(fu-1),
                          stop = 1:fu,
                          event = event_vec,
                          dose = Xmat_df[1:fu, paste0("V", id_dose)])
    
    df_wce <- rbind(df_wce, df_dose)
    
  }
  return(df_wce)
}

```


# Simulations


```{r}
# Function to simulate right constrained and unconstrained WCE with the same
#  dataset according to a specific scenario
Sylvestre_simulations <- function(scenario, Xmat){
  
  wce_mat <- do.call("rbind", lapply(1:365, wce_vector, scenario = scenario, Xmat = Xmat))
  df_wce <- get_dataset(Xmat = Xmat, wce_mat = wce_mat)
  
  # cutoff at 180 - right constrained and unconstrained with the same dataset
  wce_right_constrained <- WCE(df_wce, "Cox", 1:3, 180, constrained = "right",
                               id = "patient", event = "event", start = "start",
                               stop = "stop", expos = "dose")
  
  wce_unconstrained <- WCE(df_wce, "Cox", 1:3, 180, constrained = FALSE,
                           id = "patient", event = "event", start = "start",
                           stop = "stop", expos = "dose")
  
  # Best result according to BIC
  return(list(results_right_constrained = wce_right_constrained$WCEmat[which.min(wce_right_constrained$info.criterion), 1:180],
              results_unconstrained = wce_unconstrained$WCEmat[which.min(wce_unconstrained$info.criterion), 1:180]))
  
}

n_simul <- 10
#Sylvestre_simulations_results <- lapply(scenario, function(x) replicate(n_simul, Sylvestre_simulations(x, Xmat = Xmat)))

#Sylvestre_simulations_results %>% saveRDS("results.rds")
Sylvestre_simulations_results <- readRDS("results.rds")

# results <- lapply(1:2,
#                   function(i) Sylvestre_simulations("exponential_function",
#                                                     contrainte = "R"))
```


```{r}
# Function to extract results for plot correctly risk functions
extract_simulations_result <- function(result){
  
  right_constrained <- do.call("rbind", result[1,])
  unconstrained <- do.call("rbind", result[2,])
  
  return(list(right_constrained = right_constrained,
              unconstrained = unconstrained))
}

res <- lapply(Sylvestre_simulations_results, extract_simulations_result)

results_right_constrained <- map(res, "right_constrained")
results_unconstrained <- map(res, "unconstrained")

```



```{r}
# Correct shapes
expo_decay_list <- lapply((0:180)/365, exponential_function)
expo_decay <- do.call("rbind", expo_decay_list) / 365

bi_linear_list <- lapply((0:180)/365, bi_linear_function)
bi_linear <- do.call("rbind", bi_linear_list) / 25

early_peak <- early_peak_function[1:180]

inverted_u <- inverted_u_function[1:180]

constant_list <- lapply((0:180)/365, constant_function)
constant <- do.call("rbind", constant_list)

# Function to plot simulations results
wce_plot <- function(results, scenario_correct_shape){
 
  matplot(scenario_correct_shape, type = "l", col = "red",
          ylim = c(-0.02, 0.04), ylab = " ")
  matlines(t(results), type = "l", col = "black", lty = 1)
  matlines(scenario_correct_shape, type = "l", col = "red", lty = 1)
 
}

scenario_correct_shape <- list(expo_decay,
                              # bi_linear,
                               early_peak,
                               inverted_u,
                               constant)

```

\newpage
# Results

## Right constrained

```{r}
par(mfrow=c(2, 3))
plot_right_constrained <- map2(results_right_constrained, scenario_correct_shape, wce_plot)
```

\newpage
## Unconstrained

```{r}
par(mfrow=c(2, 3))
plot_unconstrained <- map2(results_unconstrained, scenario_correct_shape, wce_plot)
```

