## Overview

This document contains plots for safety data as well as the R code that generates these graphs. These plots are not strictly exploratory plots as data from a PopPK model are used to generate some of the plots.

## Setup

```{r, echo = FALSE, warning=FALSE, message=FALSE}
# remove reference to home directory in libPaths
.libPaths(grep("home", .libPaths(), value=TRUE, invert=TRUE))
.libPaths(grep("usr", .libPaths(), value=TRUE, invert=TRUE))

# add localLib to libPaths for locally installed packages
.libPaths(c("localLib", .libPaths()))

# will load from first filepath first, then look in .libPaths for more packages not in first path
# version matches package in first filepath, in the case of multiple instances of a package

# library(rmarkdown)
library(gridExtra)
library(grid)
library(ggplot2)
library(dplyr)
library(zoo)
library(caTools)
source("../R/xgx_packages_functions.R")

theme_set(theme_bw(base_size=12))

```

## Load Dataset
The plots presented here are based on merged safety data with a popPK model generated parameters data ([download dataset](Data/AE_xgx.csv)), as well as a dataset containing model-generated AUCs ([download dataset](Data/AUC_Safety.csv)). Data specifications can be accessed on [Datasets](Datasets.html) and Rmarkdown template to generate this page can be found on [Rmarkdown-Template](Rmarkdown/Adverse_Events.Rmd).


```{r, warning=FALSE, message=FALSE}
auc <- read.csv("../Data/AUC_Safety.csv")
my.data <- read.csv("../Data/AE_xgx.csv")
my.data$DOSE_label <- paste(my.data$Dose,"mg")
my.data <- arrange(my.data, -Dose)
my.data$DOSE_label <- factor(my.data$DOSE_label,levels = c(paste(unique(my.data$Dose),"mg")))
# add binary 0 and 1 colun to the dataset
my.data$AE_binary <- as.numeric(as.character( plyr::mapvalues(my.data$AE,
                                                               c('Yes','No'), 
                                                               c(1, 0)) ))
auc_my_data <- auc %>%
  left_join(my.data  ,by="SUBJID")

# Function to add number of subjects in each group
n_fun <- function(x){
  return(data.frame(y = median(x)*0.01, label = paste0("n = ",length(x))))
}
```


## Exposure-safety plots

These plots are looking at the AUC (or Cmax and Ctrough) for patients with the AE and for those who don’t have the AE at each dose group. The number of the patients at each dose level is also included in the plots. With these type of plots we are looking to see if there is a correlation (or absence a correlation) between exposure and AE.

### AUC vs Dose group by AE 

```{r, cache = TRUE, fig.width = 8, fig.height = 4, warning = FALSE, message= FALSE} 
gg <- ggplot(data= my.data,aes(x=factor(Dose), y=AUCAVE, fill=factor(AE))) +theme_bw()
gg <- gg + geom_boxplot(data=my.data,aes(x=factor(Dose),y=AUCAVE, fill=factor(AE)),outlier.shape = NA)
gg <- gg + geom_point(color="black",binaxis='y', stackdir='center',dotsize=0.5,
                   position=position_jitterdodge(0.27))
gg <- gg + guides(color=guide_legend(""))
gg <- gg + stat_summary(data=my.data, fun.data = n_fun, geom = "text",
                     fun.y=mean, position = position_dodge(width = 0.75)) 
gg <- gg + scale_fill_discrete(name="AE occurrence:")
gg <- gg + ylab("AUCtau (ng.h/ml)")
gg <- gg + xlab("Dose (mg)")
gg <- gg + theme(text = element_text(size=15))
gg
```



### Explore AE vs AUC

Plot binary AE occurrence against AUC or Cmax.

```{r, cache = TRUE, fig.width = 8, fig.height = 4, warning = FALSE, message= FALSE}
gg <- ggplot(data = my.data, aes(y=AUCAVE,x=AE))+theme_bw()
gg <- gg + geom_jitter(data = my.data, 
                       aes(color = DOSE_label), shape=19,  width = 0.1, height = 0, alpha = 0.5)
gg <- gg + geom_boxplot(width = 0.5, fill = NA, outlier.shape=NA) 
gg <- gg + guides(color=guide_legend(""),fill=guide_legend(""))
gg <- gg + coord_flip() 
gg <- gg + xlab("AE") + ylab("AUCtau (ng.h/ml)")
gg
```



## Probability of AE for each dose group

Plot binary AE occurrence against dose. Using summary statistics can be helpful, e.g. Mean +/- SE, or median, 5th & 95th percentiles. For binary response data, plot the percent responders along with binomial confidence intervals.

Here are some questions to ask yourself when looking at Dose-safety plots: Do you see any relationship? Does AE increase (decrease) with increasing dose? 


```{r, cache = TRUE, fig.width = 8, fig.height = 4, warning = FALSE, message= FALSE} 
gg <- ggplot(data = my.data, aes(x=Dose,y=AE_binary))+theme_bw()
gg <- gg + stat_summary(geom="errorbar", 
                          fun.data = function(y){
                              data.frame(ymin = binom::binom.exact(sum(y), length(y), conf.level = 0.95)$lower,
                                         ymax = binom::binom.exact(sum(y), length(y), conf.level = 0.95)$upper)
                            }, alpha = 0.5, size = 1, width= 0.2)
gg <- gg + stat_summary(geom="point", fun.y=mean, shape = 0)
gg <- gg + guides(color=guide_legend(""),fill=guide_legend(""))
gg <- gg + scale_y_continuous(labels=scales::percent) + ylab("Probability of AE (%)")
gg <- gg + xlab("Dose (mg)") 
gg <- gg + geom_smooth( method = "glm",method.args=list(family=binomial(link = logit)), color = "black")

pp1 <- gg

## Same plot but on a log scale
pp2 <- gg + scale_x_log10(breaks=unique(my.data$Dose)) 

grid.arrange(pp1,pp2,nrow=1)
```



### Probability of AE by AUC 

Plot AE against exposure. Include a logistic regression for binary data to help determine the shape of the exposuresafety relationship. Summary information such as mean and 95% confidence intervals by quartiles of exposure can also be plotted. Exposure and Cmax metric that you use in these pltots could be either be raw concentrations, or NCA or model-derived exposure metrics (e.g. Cmin, Cmax, AUC), and may depend on the level of data that you have available.


```{r, cache = TRUE, fig.width = 8, fig.height = 4, warning = FALSE, message= FALSE}
data_to_plot <- my.data

data_to_plot$AUC_bins <- cut(data_to_plot$AUCAVE, quantile(data_to_plot$AUCAVE, na.rm=TRUE), na.rm=TRUE)
data_to_plot <- data_to_plot[!is.na(data_to_plot$AUC_bins), ]
data_to_plot <- data_to_plot %>%
  group_by(AUC_bins) %>%
  mutate(AUC_midpoints = median(AUCAVE))

gg <- ggplot(data = data_to_plot, aes(x=AUCAVE,y=AE_binary))+theme_bw()
gg <- gg + geom_jitter(aes( color = DOSE_label), width = 0, height = 0.05, alpha = 0.5) 
gg <- gg + geom_smooth( method = "glm",method.args=list(family=binomial(link = logit)), color = "black")
gg <- gg + stat_summary(aes(x=AUC_midpoints, y=AE_binary),geom = "errorbar",
                            fun.data = function(y){
                              data.frame(ymin = binom::binom.exact(sum(y), length(y), conf.level = 0.95)$lower,
                                         ymax = binom::binom.exact(sum(y), length(y), conf.level = 0.95)$upper)
                            })
gg <- gg + stat_summary(aes(x=AUC_midpoints, y=AE_binary),shape=0,size=5, geom="point", fun.y = mean)
gg <- gg + guides(color=guide_legend(""),fill=guide_legend(""))
gg <- gg + ylab("AE")  + scale_y_continuous(breaks=c(0,1)) + coord_cartesian(ylim=c(-0.2,1.2))
gg <- gg + xlab("AUCtau (ng.h/mL)")
gg
```


## Boxplots of Exposure over time, with AE highlighted

The plot below shows time vs AUC from a pop PK model at each day. The colored dots corresponding to adverse events. It is not strictly an exploratory plot. But once you have a PopPK model, it is simple to generate this plot. Just plot the time vs the predicted AUC and color the time points red on days that an adverse event occurred.


```{r, cache = TRUE, fig.width = 8, fig.height = 4, warning = FALSE, message= FALSE}
data_to_plot <- my.data[!is.na(my.data$AETOXGRDN), ]
gg <- ggplot()+theme_bw()
gg <- gg + geom_jitter(data=auc[auc$AUC_day<30,], 
                       aes(x=AUC_day, y=AUC_popPK/AUC_day, group=AUC_day),
                       color="grey", alpha = 0.5, 
                       position=position_jitter(width=.1, height=0))
gg <- gg + geom_boxplot(data=auc[auc$AUC_day<30,], 
                        aes(x=AUC_day, y=AUC_popPK/AUC_day, group=AUC_day),
                        outlier.shape = NA, fill = NA)
gg <- gg + geom_point(data=data_to_plot[data_to_plot$DAY<30,], 
                      aes(x=DAY, y=AUC/DAY, color=factor(AETOXGRS)),size=2)
gg <- gg + scale_color_manual(breaks = c("GR1", "GR2", "GR3"),
                        values=c(rgb(1,0.5,0.5), rgb(0.75,0.25,0.25), rgb(0.5,0,0)))
gg <- gg + guides(color=guide_legend(""),fill=guide_legend(""))
gg <- gg + scale_x_continuous(breaks = seq(0, 30, by = 3))
gg <- gg + ylab("AUCtau (ng.day/ml)")
gg <- gg + xlab("Time (day)")
gg
```

## Oncology individual plots of percent change from baseline, including dosing history, labeled by "Overall Response" and AE grade

These plots allow one to look for subtle trends in the individual trajectories with respect to the dosing history, safety events and efficacy as percent change of tunor size from baseline. The plots below make use of the following oncology datasets: RECIST and nmpk dataset [(download here)](Oncology_Efficacy_Data.csv) and dose record dataset [(download here)](Oncology_Efficacy_Dose.csv)

```{r, cache = TRUE, fig.width = 10, fig.height = 14, warning = FALSE, message= FALSE}
# Read oncology efficacy data from the oncology efficacy 
# page and combine them with safety data in this page

safety_data <- read.csv("../Data/AE_xgx.csv")
efficacy_data <- read.csv("../Data/Oncology_Efficacy_Data.csv")
dose_record <- read.csv("../Data/Oncology_Efficacy_Dose.csv")

efficacy_data$DOSE_label <- paste(efficacy_data$DOSE_ABC,"mg")
efficacy_data$DOSE_label <- factor(efficacy_data$DOSE_label,levels = c(paste(unique(efficacy_data$DOSE_ABC),"mg")))

efficacy_data.monotherapy = efficacy_data %>% filter(COMB=="Single")
efficacy_data.combo = efficacy_data %>% filter(COMB=="Combo")


# Dose record data preparation for making step function plot
# in order to show any dose reduction during the study
# the idea is that at each time point, you have both the current dose and the previous dose
# the dplyr::lag command implements this
data_areaStep <- bind_rows(old = dose_record,
                           new = dose_record %>% 
                 group_by(IDSHORT) %>% 
                 mutate(DOSE = lag(DOSE)),
                        .id  = "source") %>%
                 arrange(IDSHORT, TIME, source) %>%
  ungroup() %>%
  mutate(DOSE = ifelse(lag(IDSHORT)!=IDSHORT, NA, DOSE), 
          TIME = TIME/24) #convert to days

data_areaStep.monotherapy = filter(data_areaStep,COMB=="Single")

# calculate average dose intensity up to the first assessment:
# "TIME==57"" is the first assessment time in this dataset
first.assess.time = 57
dose_record <- dose_record %>%
  group_by(IDSHORT) %>%
  mutate(ave_dose_intensity = mean(DOSE[TIME/24 < first.assess.time]))

dose_intensity <- dose_record[c("IDSHORT","COMB","ave_dose_intensity")]
dose_intensity <- subset(dose_intensity, !duplicated(IDSHORT))


# This part is optional to label "OR" in the plot
# "OR" can be substituted with other information, such as non-target, new target lesions
#  make the OR label for the plot


safety_label <- safety_data %>%
  select(SUBJID, DAY, AETOXGRS, Dose)

colnames(safety_label)[2] <- "TIME"
colnames(safety_label)[4] <- "DOSE_ABC"
safety_label$AETOXGRS <- as.character(safety_label$AETOXGRS)
safety_label <- safety_label[!safety_label$AETOXGRS =="",]

efficacy_AE_label <- efficacy_data %>%
  select(SUBJID, TIME, psld, DOSE_ABC)
efficacy_AE_label <- merge(safety_label,efficacy_AE_label, by = c("SUBJID", "TIME","DOSE_ABC"),
                            all.x=T, all.y=T)

subj <- efficacy_AE_label  %>% 
  subset(!is.na(psld)) %>%
  group_by(SUBJID) %>%
  mutate(CountNonNa = length(psld))

subj <- c(unique(subset(subj, CountNonNa>1, "SUBJID")))

efficacy_AE_label <- efficacy_AE_label  %>% 
  subset(SUBJID%in%subj$SUBJID)%>%
  group_by(SUBJID) %>%
  mutate(ValueInterp = na.approx(psld,TIME, na.rm=FALSE))

efficacy_AE_label <- efficacy_AE_label[!is.na(efficacy_AE_label$AETOXGRS),]
efficacy_AE_label <- efficacy_AE_label[!is.na(efficacy_AE_label$ValueInterp),]
efficacy_AE_label <- subset( efficacy_AE_label, select = -psld )
colnames(efficacy_AE_label)[5] <- "psld"
colnames(efficacy_AE_label)[1] <- "IDSHORT"



efficacy_data.label <- efficacy_data %>%
  group_by(SUBJID) %>%
  mutate(label_psld = as.numeric(ifelse(TIME==TIME_OR , psld,""))) %>%
  filter(!(is.na(label_psld) | label_psld==""))

dose.shift = 50
dose.scale = 1.2
data_areaStep.monotherapy = data_areaStep.monotherapy %>%
  mutate(DOSE.shift = DOSE/dose.scale+dose.shift)

dose.unique = c(0,unique(efficacy_data.monotherapy$DOSE_ABC))

gg <- ggplot(data = efficacy_data.monotherapy)
gg <- gg + geom_point(mapping = aes(y= psld, x= TIME))
gg <- gg + geom_text(data= efficacy_data.label,aes(y= label_psld, x= TIME_OR, label=OR), vjust=-.5)
gg <- gg + geom_hline(aes(yintercept = 0),size=0.1, colour="black")
gg <- gg + geom_line(mapping = aes(y= psld, x= TIME))
gg <- gg + geom_ribbon(data= data_areaStep.monotherapy,
                       aes( ymin = 50, ymax = DOSE.shift , x= TIME),
                       fill="palegreen2", color = "black", alpha=0.5 )
gg <- gg + geom_text(data= efficacy_AE_label,
                     aes(y= psld, x= TIME, label=AETOXGRS), colour="red",fontface=2,
                     size=5, show.legend = F, hjust=-0.05, vjust=2)
gg <- gg + geom_vline(data= efficacy_AE_label,
                      aes(xintercept= TIME), 
                      size=1, linetype="dashed", colour="red")
gg <- gg + facet_wrap(~IDSHORT, ncol=6)
gg <- gg + scale_y_continuous(
  sec.axis = sec_axis(~(.-dose.shift)*dose.scale, name = "Dose(mg)", breaks = dose.unique))
gg <- gg + labs(y = "Percent Change in\nSum of Longest Diameters",
                x = "Time (months)",
                colour = "Parameter")
gg <- gg + scale_x_units(units.input = "day",units.output="month",t.start = 0,t.end = 15, increment = 3)
gg <- gg + theme(text = element_text(size=15))
gg
```

## R Session Info

```{r}
sessionInfo()
```
