---
title: "Supplementary code for: Topological network analysis of patient similarity for precision manAgement of acute blood pressure in spinal cord injury"
author: "A Torres-Espin, J Haefeli, R Ehsanian, D Torres, CA Almeida, JR Huie, A Chou, D Morozov, N Sanderson, B Dirlikov, CG Suen, JL Nielson, N Kyritsis, DD Hemmerle, JF Talbott, GT Manley, SS Dhall, WD Whetstone, JC Bresnahan, MS Beattie, SL McKenna, JZ Pan, AR Ferguson and the TRACK-SCI Investigators"
date: "09/06/2021"
output: 
  html_document:
    toc: yes
    toc_float: yes
    collapsed: false
    smooth_scroll: no
    number_sections: true
    toc_depth: 3
editor_options: 
  chunk_output_type: console
---

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

# Setting up

## Libraries

Load the libreries used through the entire analysis

```{r}
library(igraph) #to build and manAge graph data  
library(tidyverse) #data cleaning capabilities
library(reshape2) #data formatting
library(visNetwork) #for interactive netwrok visualization
library(vegan) #it implements a nice isomap interface
library(RColorBrewer) #for managing colors
library(TDAstats) #it implements persistent homology and other TDA methods
library(cccd) #to compute the k-NN graph
library(table1)#to create the demographic table
library(glmnet) #to fit the LASSO
library(glmnetUtils) #to fit the LASSO
library(caret) #for cross-validate the logistic regression
library(splines) #to fit the spline models
library(boot) ## for estimating LOOCV in glm models
library(viridis) ## color schema
library(ROCR) ## ROC processing
library(patchwork) ## for figures
library(png) ## for imAge processing for figures
library(grid) ## for figures
library(webshot) ## to take a shot at the graphs for figures
library(glmulti) ## for model selection
```

## Custom helper functions

This script contains the custom functions used as helpers for the paper. Please note these functions are specific to the paper and have been coded to run with the main script of the publication. Running them in different settings might fail or produce unwanted results.

**plot_legend:** Custom function for plotting the legends of the networks

```{r}
plot_legend<-function(values, colors, legend_title, breaks, labels){
  #Custom function for plotting the legends of the networks
  #
  #@param values: values to be plotted
  #@param colors: the vector of colors to plot
  #@param legend_title: the string value for the legend title
  #@param breaks: the numeric vector for the x axis breaks
  #@param labels: the string vector for the x axis breaks
  #
  #@returns ggplot object with network gradient legend
  
  if (!is.numeric(values)){
    value<-as.numeric(as.factor(values))
  }else{
    value<-values
  }
  # plot(value, y=rep(0.5, length(value)),yaxt="n", ylab = '',
  #      type='n',frame.plot = F, xlim=c(min(value),max(value)))
  # legend_imAge <- as.raster(matrix(colors, nrow=1))
  # rasterImAge(legend_imAge, min(value), 0.3, max(value),0.5)
  
  
  # values<-mapping$mean.map
  # color<-net.color.Pal(250)
  # legend_title<-"MAP (mmHg)"
  df<-data.frame(values=seq(min(values), max(values), length.out = length(colors)))
  legend_plot<-ggplot(df, aes(values, y=1, fill=values))+
    geom_tile()+
    scale_fill_gradientn(colours = colors)+
    scale_x_continuous(breaks = breaks, labels = labels)+
    ggtitle(legend_title)+
    theme_minimal()+
    theme(plot.title = element_text(hjust = 0.5), axis.text.y = element_blank(),
          axis.title = element_blank(),
          axis.ticks.x = element_line(color="black"),
          axis.ticks.length.x = unit(3,"mm"),
          legend.position = "none",
          panel.grid = element_blank())
  
  return(legend_plot)
}
```

<a href="#header">Back to top</a>

**extract_cumfreq:** Function to extract the cummulative frequency distribution for each patient by binning. The bin width is specified by the per_step argument determining the percentAge respect the range of the specific measure the bin will cover. 1% by default

```{r }
extract_cumfreq<-function(df,per_step=1){
  
  # Function to extract the cummulative frequency distribution for each patient by binning. The bin width is specified by the per_step argument determining the percentAge respect the range of the specific measure the bin will cover. 1% by default
  
  #@param df: the dataframe to extract the cummulative frequency dist in percentAge step
  #@param per_step: the percentAge width of binning
  
  fmin<-floor(min(df, na.rm = T))
  fmax<-ceiling(max(df, na.rm = T))
  step<-per_step*(fmax-fmin)/100
  range.t<-seq(fmin,fmax,step)
  
  cumfreq_int<-function(x, range.int=range.t){
    temp<-table(cut(as.numeric(x), breaks = range.t))
    temp<-temp/sum(temp)
    return(cumsum(temp))
  }
  
  result<-as.matrix(t(apply(df, 1, cumfreq_int)))
  
  return(result)
}
```

<a href="#header">Back to top</a>

**phom.dist2:** Wrap function of the phom.dist function from the TDAstats packAge to perform the wasserstein distance with power of 2 (see TDAstats packAge for details)

```{r}
phom.dist2<-function (phom1, phom2){
  
#   Wrap function of the phom.dist function from the TDAstats packAge
# to perform the wasserstein distance with power of 2 (see TDAstats packAge for details)
  
  max.dim1 <- max(phom1[, 1])
  max.dim2 <- max(phom2[, 1])
  TDAstats:::wass_mat_calc(phom1, phom2, pow.val = 2, dim = max(max.dim1, 
                                                                max.dim2))
}
```

<a href="#header">Back to top</a>

**Net.viz:** This function plots the network. VisNetwork is an interface to vis.js packAge in JavaScript. The randomSeed of the network is established for reproducibility proposes so that we always get the same network shape.

```{r }
Net.viz<-function(colors=NULL, network_name="", shot=T){
  
  # This function plots the network. VisNetwork is an interface to vis.js packAge in JavaScript. The randomSeed of the network is established for reproducibility proposes so that we always get the same network shape.
  
  #@param colors: colors to use in the network
  #@param network_name: name of the network to save file snapshot
  
  g.data$nodes$color<-colors
  net_vis<-visNetwork(nodes = g.data$nodes, edges = g.data$edges, height = "800px", width = "800px")%>%
    visIgraphLayout(randomSeed = 78,smooth = F, physics = T,
                    layout='layout_with_fr')%>%
    visEdges(color = list('inherit'='both', 'opacity'=0.5), width = 10)%>%
    visNodes(size=25, mass=3)%>%
    visPhysics(stabilization = T)
  
  if (shot){
    ## Save a png version of the network
    html_name <- tempfile(fileext = ".html") #set temporary file
    visSave(net_vis, html_name) #save the html network to temporary file
    webshot(html_name,file = paste0(network_name,".png"), delay = 1,
          vwidth = 800, vheight = 800,cliprect = c(80,100,700,700)) #webshot the html network
    # dev.off()
  }
  
  return(net_vis)
}
```

<a href="#header">Back to top</a>

**cNet.viz:** Visualize the contracted network. The CNet.viz function will plot the network

```{r }
cNet.viz<-function(colors=NULL, network_name=""){
  
  #Visualize the contracted network. The CNet.viz function will plot the network
  
  #@param colors: colors to use in the network
  #@param network_name: name of the network to save file snapshot
  
  gcon.data$nodes$color<-colors
  net_vis_con<-visNetwork(nodes = gcon.data$nodes, edges = gcon.data$edges, height = "800px", width = "800px")%>%
    visIgraphLayout(randomSeed = 74,smooth = F, physics = T,
                    layout='layout_with_fr')%>%
    visEdges(color = list('inherit'='both', 'opacity'=0.5),length = 50, width = 3)%>%
    visNodes(scaling = list('max'=35, 'min'=10), mass=1.5)%>%
    visPhysics(stabilization = T)
  
  ## Save a png version of the network
  html_name <- tempfile(fileext = ".html") #set temporary file
  visSave(net_vis_con, html_name) #save the html network to temporary file
  webshot(html_name,file = paste0(network_name,"_con.png"), delay = 1,
        vwidth = 800, vheight = 800,cliprect = c(260,240,410,360)) #webshot the html network
  # dev.off()
  return(net_vis_con)
}
```

<a href="#header">Back to top</a>

**plot_net:** Following function is used to plot the network plots in the paper and cluster plots with the specified variable and aggregating function (mean, median or proportion)

```{r }
plot_net<-function(var, aggre.fun, y_axe_text, value.proportion,
                   prop.colors, prop.labels=c("A","B","C","D","E","NA"),
                   ranges=c(1,100), network_name){
  
  #Following function is used to plot the network plots in the paper and cluster plots with the specified variable and aggregating function (mean, median or proportion)
  
  #@param var: the variable to map onto the network
  #@param aggre.fun: the aggregated function to use: "mean", "median", "proportion"
  #@param y_axe_text: the text label in y axis to plot
  #@param value.proportion: value to compute the proportion for.
  #@param prop.colors: sequence of colors for the proportions
  #@param prop.labels: labels for the ggplot legend when aggre.fun = "proportion"
  #@param ranges:
  #@param network_name: name of the network to save file snapshot

  proportions<-function(x, ...){sum(x==value.proportion, na.rm=T)/length(na.omit(x))}
  
  if (ranges[1]!=1|ranges[2]!=100){
    values<-as.numeric(mapping[,var])
    
    minimal<-min(values, na.rm = T)+min(values, na.rm = T)*(ranges[1]/100)
    maximal<-max(values,na.rm = T)*(ranges[2]/100)
    values2<-ifelse(values<minimal, minimal, values)
    values2<-ifelse(values2>maximal, maximal, values2)
    values<-values2
  }else{
    values<-mapping[,var]
  }
  
  agg_cluster<-lapply(split(values, clusters$membership), get(aggre.fun), na.rm=T)
  agg_cluster2<-lapply(split(values, mapping$cluster), get(aggre.fun),na.rm=T)
  
  agg_cluster<-unlist(agg_cluster)
  agg_cluster2<-unlist(agg_cluster2)
  node.col.agg <- net.color.Pal(250)[as.numeric(cut(agg_cluster,breaks = 250))]
  node.col2.agg <- net.color.Pal(250)[as.numeric(cut(agg_cluster2,breaks = 250))]
  
  
  if(aggre.fun=='proportions'){
    map_colors<-as.numeric(mapping[,var])
    map_colors[is.na(map_colors)]<-which(prop.labels=="NA")
    node.col<-prop.colors[map_colors]
  }else{
    node.col<-net.color.Pal(250)[as.numeric(cut(values,breaks = 250))] 
  }
  
  net<-Net.viz(node.col, network_name) #replot the similarity network with the mapping var
  cnet<-cNet.viz(node.col.agg, network_name) #replot the contracted network with the mapping var
  
  if (aggre.fun%in%c('mean','median')){
    plot_c<-ggplot(mapping, aes_string('cluster', var,fill='cluster'))+
      geom_boxplot(size=0.6, outlier.size = 1.5)+
      scale_fill_manual(values=node.col2.agg)+
      theme_minimal(base_size = 11)+
      xlab('Patient cluster')+ylab(y_axe_text)+
      theme(legend.position = 'none', axis.text.x = element_text(size = 8))+
      scale_x_discrete(labels=paste('C', 1:11, sep=''))
  }else if(aggre.fun=='proportions'){
    prop.df<-mapping
    prop.df[,var]<-as.factor(prop.df[,var])
    plot_c<-ggplot(prop.df,aes_string('cluster', fill=var))+
      geom_bar(position='fill')+
      scale_fill_manual(values=prop.colors, na.value='grey50',
                        labels=prop.labels)+
      theme_minimal(base_size = 11)+
      xlab('Patient cluster')+ylab(y_axe_text)+
      theme(legend.position = 'none', 
            panel.grid.minor = element_blank(), 
            panel.grid.major.x = element_blank())+
      theme(legend.position = 'bottom', 
            legend.title = element_blank(),
            legend.key.size = unit(4, "mm"),
            legend.box = "horizontal",
            axis.text.x = element_text(size = 8))+
      guides(fill=guide_legend(nrow=1,byrow=TRUE))+
      scale_x_discrete(labels=paste('C', 1:11, sep=''))
  }
  
  return(list(net=net,con.net=cnet, plot=plot_c, 
              agg_values=agg_cluster,network_name=network_name))
}

```

<a href="#header">Back to top</a>

**prediction_fun**: function to extract prediction metrics and prediction threshold

```{r}
prediction_fun<-function(prob_prediction, real_class, ref){
  ## Function to extract LOOCV prediction metrics and prediction threshold
  #
  #@param prob_prediction: probability of prediction
  #@param real_class: the class to predict
  
  predictions<-prediction(prob_prediction, real_class)
  
  #Extract sensitivity and specificity into a dataframe
  roc_df<- data.frame(x=unlist(performance(predictions, "sens")@x.values),
                      sens_y=unlist(performance(predictions, "sens")@y.values),
                      spec_y=unlist(performance(predictions, "spec")@y.values))
  
  #calculate threshold as the point that balance sensitivity and specificity
  roc_df<-roc_df%>%
    mutate(diff_y = abs(sens_y-spec_y))
  
  #best prediction threshold
  pred_threshold<-roc_df$x[which.min(roc_df$diff_y)]
  
  #Predicting with the best threshold
  prediction_best<-ifelse(prob_prediction>pred_threshold, 1,0)
  conf_mat<-caret::confusionMatrix(reference=factor(real_class),
                                   data=factor(prediction_best), positive=ref)
  
  #Sensitivity-Specificity plot
  sense_spec_plot<-ggplot(roc_df, aes(x, sens_y))+
    geom_line(color="orange")+
    geom_line(data=roc_df, aes(x, spec_y), color = "steelblue")+
    geom_vline(xintercept = pred_threshold, linetype="dashed", 
               color="red", alpha=0.5)+
    theme_minimal()+
    ylab("Value (spec or sense)")+
    xlab("Prediction threshold")+
    annotate("text", label=round(pred_threshold, 3), x=0.55, y=0.1)+
    annotate("text", label="Sensitivity", x=0.1, y=0.85, color="orange")+
    annotate("text", label="Specificity", x=0.8, y=0.85, color="steelblue")+
    xlim(c(0,1))
  
  #AUC
  auc<-unlist(performance(predictions, "auc")@y.values)
  
  return(list("conf_mat"=conf_mat, "sense_spec_plot"=sense_spec_plot, 
              "auc"=auc, "roc_df"=roc_df))
}
```

<a href="#header">Back to top</a>

## global settings

```{r}
set.seed(666) #set a random seed for reproducibility of the random iterations in the code

# set color palette color-blind 'proof' for some graphs #from http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/
cbbPalette <- c("#000000", "#E69F00", "#56B4E9","#009E73", "#F0E442", 
                "#0072B2", "#D55E00", "#CC79A7")

# Colors and pallet for the network plots
# net.color <- c('blue2','blue', 'green2','yellow','#ffa500','red2')
net.color <- plasma(6)
net.color.Pal <- colorRampPalette(net.color)
```

<a href="#header">Back to top</a>

# Reading data

Reading files for both centers data and joining them. See paper.
To reproduce the entire analysis and figures, make sure you download the dataset for ODC-SCI and specify the path.

```{r}
## path to downloaded datasets
data1<-read.csv('../data/odc_sci_245.csv')
data2<-read.csv('../data/odc_sci_246.csv')

raw.df<-rbind(data1, data2)
colnames(raw.df)[1] <- "Subject_ID"

# raw.df<-read.csv('data/OR data publication.csv')
mapping<-raw.df
```

<a href="#header">Back to top</a>

# Demographic table

Construct demographics table 1 in the paper

```{r}
demographics.df<-raw.df[,c('Age', 'AIS_ad','AIS_dis', 'length_surg',
                           'AIS_is_improved','days_post_surgery_dish',
                           'Injury_level')]
demographics.df$AIS_is_improved<-as.character(demographics.df$AIS_is_improved)
demographics.df$AIS_is_improved[is.na(demographics.df$AIS_is_improved)]<-'N/A'
demographics.df$days_post_surgery_dish<-as.numeric(as.character(demographics.df$days_post_surgery_dish))
label(demographics.df$Age)<-'Age'
label(demographics.df$AIS_ad)<-'AIS admission'
label(demographics.df$AIS_dis)<-'AIS discharge'
label(demographics.df$length_surg)<-'Surgery duration'
label(demographics.df$days_post_surgery_dish)<-'Surgery to discharge'
label(demographics.df$AIS_is_improved)<-'AIS improvement'
label(demographics.df$Injury_level)<-'NLI admission'

units(demographics.df$Age)<-'years'
units(demographics.df$length_surg)<-'minutes'
units(demographics.df$days_post_surgery_dish)<-'days'

#Plot table 1
table1(~Age+AIS_ad+AIS_dis+length_surg+
         Injury_level+
         days_post_surgery_dish|AIS_is_improved, data=demographics.df, overall = F)
```

## Univariate test

Univariate test for demographics between improvers and non-improvers. T-test for numerical variables, Chi square test for categorical variables.

```{r }
p.value.df<-demographics.df[demographics.df$AIS_is_improved!='N/A',]
p.value.df<-droplevels(p.value.df)
x<-p.value.df$AIS_is_improved

p.value.results<-list()
for (var in colnames(demographics.df)){
  y<-p.value.df[,var]
  if (isTRUE(all.equal(x,y))){next()}
  
  if (is.numeric(y)){
    p.value<-t.test(y~x)$p.value
  }else{
    p.value<-fisher.test(x,y)$p.value
  }
  p.value.results[[var]]<-p.value
}

p.value.results<-as.data.frame(do.call(rbind, p.value.results))

p.value.results
```

<a href="#header">Back to top</a>

# Data preprocessing

This section performs the data preprocessing steps necessary to conduct the rest of the analysis.

-   Extraction of the OR q5 values
-   Extract the cummulative frequency distribution for each patient and measure

```{r}
#Extraction of the OR q5 values
map.columns<-colnames(raw.df)[str_detect(colnames(raw.df), 'OR_MAP_Q')]
map.wide<-raw.df[,c('Subject_ID',map.columns)]

#Format in long for some applications
map_df_long<-melt(map.wide)

HR.columns<-colnames(raw.df)[str_detect(colnames(raw.df), 'OR_HeartRate_Q')]
HR.wide<-raw.df[,c('Subject_ID',HR.columns)]

#Format in long for some applications
hr_df_long<-melt(HR.wide)

#Extract the cumulative frequency distribution for each patient and measure
MAP.freq<-extract_cumfreq(map.wide[,-1],1) #the bin width is 1% of the MAP range
HR.freq<-extract_cumfreq(HR.wide[,-1],1) #the bin width is 1% of the HR range
```

**Table Example**

```{r}
knitr::kable(map.wide[1:10,1:8])
```

## Figure 1

Code for the construction of figure 1

```{r fig.width = 10, fig.height = 12}
## Load consort diagram
consort<-readPNG("or consort.png") ## This is an external imAge included as supplementary file. Make sure to specify the path location

consort_grob<-rasterGrob(consort)

or_map_labels<-c("OR_MAP_Q120","OR_MAP_Q480","OR_MAP_Q955")
or_HR_labels<-c("OR_HeartRate_Q120","OR_HeartRate_Q480","OR_HeartRate_Q955")

#MAP time plot
map_patient_plot<-ggplot(map_df_long, aes(x=variable, y=Subject_ID, fill=value))+
  geom_tile()+
  xlab('Time (Q5 minutes)')+ylab('SCI patients')+
  scale_fill_gradientn(colours = c('blue3','steelblue','lightgreen','orange','firebrick3'),
                       limits=c(40,140), na.value = 'grey')+
  scale_x_discrete(breaks=or_map_labels, labels=c('2h', '8h','16h'))+
  geom_vline(xintercept = c(24,96, 191), linetype=2)+
  labs(fill = "mmHg")+
  theme(legend.position = c(0.89,0.2), axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.background = element_blank(),
        legend.box.background = element_blank())

#plot mean MAP
map_patient_mean_plot<-ggplot(map_df_long, aes(x=Subject_ID, y=value, group=Subject_ID))+
  stat_summary(fun.y = 'mean', geom='bar', fill='orange', width=0.5)+
  stat_summary(fun.y = 'mean', geom='point')+
  coord_flip(ylim = c(55, 110))+
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        panel.grid = element_blank(), panel.background = element_blank())+
  ylab('aMAP (mmHg)')+xlab(NULL)+
  scale_y_continuous(breaks = c(60,80,100))

#HR time plot
hr_patient_plot<-ggplot(hr_df_long, aes(x=variable, y=Subject_ID, fill=value))+
  geom_tile()+
  xlab('Time (Q5 minutes)')+ylab('SCI patients')+
  scale_fill_gradientn(colors=c('yellow4','yellow','purple','purple4'),
                       values=c(0,0.2,0.5,0.9,1),na.value = 'grey',
                       limits=c(0,160))+
  scale_x_discrete(breaks=or_HR_labels, labels=c('2h', '8h','16h'))+
  geom_vline(xintercept = c(24,96, 191), linetype=2, alpha=0.6)+
    labs(fill = "bpm")+
  theme(legend.position = c(0.89,0.2), axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.background = element_blank(),
        legend.box.background = element_blank())

#Plot mean HR
hr_patient_mean_plot<-ggplot(hr_df_long, aes(x=Subject_ID, y=value))+
  stat_summary(fun.y = 'mean', geom='bar', fill='orange', width=0.5)+
  stat_summary(fun.y = 'mean', geom='point', fill='red')+
  coord_flip(ylim = c(50, 110))+
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        panel.grid = element_blank(), panel.background = element_blank())+
  ylab('aHR (beats per min.)')+xlab(NULL)+
  scale_y_continuous(breaks = c(60,80,100))

#Layout of the figure
layout <- c(
  #Line 1
  area(t=1, l=1, r=10, b=4),
  # Line 2
  area(t=5, l=1, r=3, b=7),area(t=5, l=4, r=5, b=7),area(t=5, l=6, r=8, b=7),
  area(t=5, l=9, r=10, b=7)
)

#Building the figure
figure1<-wrap_elements(consort_grob)+
  map_patient_plot+ggtitle("Mean Arterial Pressure (MAP)")+
  map_patient_mean_plot+
  hr_patient_plot+ggtitle("Heart Rate (HR)")+
  hr_patient_mean_plot

#Specify figure options
figure1<-figure1+plot_layout(design = layout)+
  plot_annotation(tag_levels = "a")

#Save the figure
ggsave(figure1,file="figure 1.png", width = 10, height = 12,dpi = "retina")

figure1
```

<a href="#header">Back to top</a>

## Example plots

Example plots for visualizing the data and the preprocessing.

Generate the dataframe for ploting 4 raw examples and 4 processed examples (see paper for details)

```{r fig.width=10, fig.height=3}
#Raw map data for 4 patient examples
df.plot.raw<-data.frame(patient=c(rep('Patient 100', length(map.wide[100,-1])),
                                  rep('Patient 1', length(map.wide[1,-1])),
                                  rep('Patient 10', length(map.wide[10,-1])),
                                  rep('Patient 86', length(map.wide[86,-1]))),
                        map=c(unlist(map.wide[100,-1]),unlist(map.wide[1,-1]),
                              unlist(map.wide[10,-1]),unlist(map.wide[86,-1])),
                        q5=c(seq(1,length(map.wide[100,-1])),
                             seq(1, length(map.wide[1,-1])),
                             seq(1, length(map.wide[10,-1])),
                             seq(1, length(map.wide[86,-1]))))

# cummulative frequency distribution of map data for 4 patient examples
df.plot.cumfreq<-data.frame(patient=c(rep('Patient 100', length(MAP.freq[100,])),
                                      rep('Patient 1',length(MAP.freq[1,])),
                                      rep('Patient 10',length(MAP.freq[10,])),
                                      rep('Patient 86', length(MAP.freq[86,]))),
                            freq=c(unlist(MAP.freq[100,]), unlist(MAP.freq[1,]),
                                   unlist(MAP.freq[10,]),unlist(MAP.freq[86,])),
                            bin=c(1:length(MAP.freq[100,]), 1:length(MAP.freq[1,]),
                                  1:length(MAP.freq[10,]), 1:length(MAP.freq[86,])))

#Panel A
fig2_sup1_a<-ggplot(df.plot.raw,aes(q5, map, color=patient))+
  geom_line()+
  theme_minimal(base_size = 10)+
  ylab('MAP(mmHg)')+xlab('Time (Q5)')+
  theme(legend.title = element_blank(), legend.position = c(0.7, 0.9))+
  guides(color=guide_legend(nrow=2,byrow=TRUE))

#Panel B
fig2_sup1_b<-ggplot(df.plot.cumfreq,aes(bin, freq, color=patient))+
  geom_line()+
  theme_minimal(base_size = 10)+
  ylab('MAP cummulative\ndist. function')+xlab('Bin of MAP (1% binwidth)')+
  theme(legend.title = element_blank(), legend.position = c(0.8, 0.3))

fig2_sup1_a+fig2_sup1_b
```

<a href="#header">Back to top</a>

## Optimization of isomap solutions

The following section iterates through the procedure to optimize the parameter selection for isomap according to the residual variance and the persistence homology criteria (see methods in paper for details)

```{r}
#Combine MAP and HR cumulative freq. distributions. Each row is a patient
processed_df<-dist(cbind(MAP.freq,HR.freq), method = 'euclidean')

#'Iterates through the calculations of RV and wasserstein distances for
#'k 3 to 7 and dimensions kept from 1 to 12. (see methods in the paper for details)
opt.df<-list()
counter<-1
for (k in 3:7){
  gdist.k<-isomapdist(processed_df, k=k)
  ph.gdist.k<-calculate_homology(as.matrix(gdist.k),format = 'distmat')
  
  for (i in 1:12){
    isomap.k<-isomap(processed_df, k=k,ndim = i)
    ph.emd.k<-calculate_homology(as.matrix(dist(isomap.k$points)), 
                                 format = 'distmat')
    RV<-1-(cor(c(gdist.k), c(dist(isomap.k$points)))^2)
    WD<-phom.dist2(ph.gdist.k,ph.emd.k)
    opt.df[[counter]]<-data.frame(k=k, d=i, RV=RV,WD0=WD[1],WD1=WD[2])
    counter<-counter+1
  }
}
isomap_solution<-do.call(rbind, opt.df)

knitr::kable(isomap_solution)
```

### Optimization plots

```{r fig.width=10, fig.height=3}
# Plot RV, panel C
fig2_sup1_c<-ggplot(isomap_solution, aes(d, RV, color=factor(k), shape=factor(k)))+
  geom_point()+
  geom_line()+
  scale_x_continuous(breaks=1:12)+
  theme_minimal(base_size = 10)+
  theme(panel.grid.minor.x = element_blank(),
        legend.key.size = unit(0.4,'cm'),
        legend.position = c(0.8,0.7))+
  scale_color_manual(values=cbbPalette)+
  xlab('# dimensions kept')+ylab('Residual variance')+
  labs(color='Isomap k', shape='Isomap k')

# Plot WD 0 dimension, panel D
fig2_sup1_d<-ggplot(isomap_solution, aes(d, WD0, color=factor(k), shape=factor(k)))+
  geom_line()+
  geom_point()+
  scale_x_continuous(breaks = 1:12)+
  xlab('# dimensions kept')+ylab('Wasserstein 0-d distance')+
  theme_minimal(base_size = 10)+
  scale_color_manual(values=cbbPalette)+
  theme(panel.grid.minor = element_blank(), legend.position = "none")+
  labs(color='Isomap k', shape='Isomap k')

# Plot WD 1 dimension, panel E
fig2_sup1_e<-ggplot(isomap_solution, aes(d, WD1, color=factor(k), shape=factor(k)))+
  geom_line()+
  geom_point()+
  scale_x_continuous(breaks = 1:12)+
  xlab('# dimensions kept')+ylab('Wasserstein 1-d distance')+
  theme_minimal(base_size = 10)+
  scale_color_manual(values=cbbPalette)+
  theme(panel.grid.minor = element_blank(),legend.position = "none")+
  labs(color='Isomap k', shape='Isomap k')

fig2_sup1_c+fig2_sup1_d+fig2_sup1_e
```

<a href="#header">Back to top</a>

# Exploratory network analysis suggests an upper and lower limit of intraoperative MAP for recovery

The following section of the code extracts the similarity network as specified in the paper.

## Extract Network

```{r}
iso.all<-isomap(processed_df,k=6, ndim = 4)#Solution from the isomap optimization
D<-iso.all$points

#' Iteration until to find the minimal k for K-nn mutual graph with 
#' one single component (each patient connected to at least one other)
k<-1
g<-nng(x=D, k = k, mutual = T)
g<-igraph::simplify(g)
g<-as.undirected(g)

while (max(components(g)$membership)>1) {
  k<-k+1
  g<-nng(x=D, k = k, mutual = T)
  g<-igraph::simplify(g)
  g<-as.undirected(g)
}

g<-igraph::simplify(g)#simplify the graph to eliminate self referencing edges
g<-add_layout_(g, nicely())#add layout to the graph (only for visualization using plot())
# plot(g)# raw plot of the graph (see igraph packAge for more details)
```

<a href="#header">Back to top</a>

## Plot the 'naked network'

```{r fig.height=5, fig.width=5}
g.data<-toVisNetworkData(g) #transform g into a VisNetwork object for visualization
Net.viz(shot=T)

if (!is.null(dev.list())) {
  dev.off()
}
```

<a href="#header">Back to top</a>

## Network clustering

The following code finds the community structure (clusters) using the walktrap algorithm.

```{r }
clusters<-c()
member.temp<-list()
modularity.df<-c()

#iterates through an arbitrary range of randoms steps (from 1 to 100)
for (i in 1:100){
  g.clusters<-cluster_walktrap(g, steps = i)
  clusters[i]<-max(g.clusters$membership)
  member.temp[[i]]<-g.clusters$membership
  modularity.df[i]<-modularity(g.clusters,g.clusters$membership)
}

#Choose the solution that maximise modularity
steps.t<-which.max(modularity.df) 

#Construct data frame from the results for visualization
cluster.df<-data.frame(steps=c(1:length(clusters)), 
                       n_clusters=clusters, modu=modularity.df)
```

<a href="#header">Back to top</a>

### Random walk plots

```{r fig.width=10, fig.height=3}
#Panel F
fig2_sup1_f<-ggplot(cluster.df, aes(x=steps, y=modu))+
  geom_line(color='blue')+
  xlab('Number of random steps')+
  ylab('Modularity (Q)')+
  theme_minimal(base_size = 10)+
  geom_vline(xintercept = steps.t, linetype='dashed')

# ggsave('modu-steps.png', dpi='print', height = 2, width = 2.4)

#Panel G
fig2_sup1_g<-ggplot(cluster.df, aes(x=steps, y=n_clusters))+
  geom_line(color='blue')+
  xlab('Number of random steps')+
  ylab('Number of clusters')+
  theme_minimal(base_size = 10)+
  geom_vline(xintercept = steps.t, linetype='dashed')

fig2_sup1_f+fig2_sup1_g

# ggsave('clustes-steps.png', dpi='print', height = 2, width = 2.4)
```

<a href="#header">Back to top</a>

## Collapsed network

The following code extract the collapsed network based on the clustering

```{r}
clusters<-cluster_walktrap(g, steps =steps.t)
gcon<-contract.vertices(g, clusters$membership,vertex.attr.comb=toString) #contract the network
gcon<-igraph::simplify(gcon)#Eliminate the self-referencing edges

gcon.data <- toVisNetworkData(gcon)
n_patients_cluster<-as.vector(table(clusters$membership)) #calculate number of patients per cluster

#Extraction of colors for the dendogram
cluster_colors<-brewer.pal(11, 'Paired')
```

### Plot the dendogram

```{r fig.width=10, fig.height=3}
dendPlot(clusters, mode = 'p',direction='downwards',show.tip.label=F,
         palette = cluster_colors)

# save the plot
png("dendogram.png")
  dendPlot(clusters, mode = 'p',direction='downwards',show.tip.label=F,
         palette = cluster_colors)
dev.off()
```

<a href="#header">Back to top</a>

## Network plots

The following section include the network plots.

```{r}
#Here the labels for the clusters are override to those provided in the paper
mapping$cluster<-factor(clusters$membership)
levels(mapping$cluster)<-c(9,5,7,2,3,4,6,11,10,8,1)
mapping$cluster<-factor(mapping$cluster, levels=1:11)

gcon.data$nodes$label<-paste('C',c(9,5,7,2,3,4,6,11,10,8,1), sep='')
gcon.data$nodes$value<-n_patients_cluster*2 #Set the size of the node proportional to the number of patients on it
```

### Extract mapped networks

This section creats the networks mapped to specified variables using the custom helper functions at the begining of this document.

```{r}
#MAP net
MAP.gcon<-plot_net('mean.map','mean','AverAge MAP (mmHg)', network_name = "map_net")

#HR net
HR.gcon<-plot_net('mean.HR','mean','AverAge HR (bpm)', network_name = "HR_net")

# AIS improved
mapping$AIS_is_improved<-as.factor(mapping$AIS_is_improved)
mapping$AIS_dis<-as.factor(mapping$AIS_dis)
mapping$AIS_is_improved_num<-as.numeric(mapping$AIS_is_improved)
mapping$AIS_is_improved_num[mapping$AIS_is_improved_num==3]<-NA
mapping$AIS_dis_num<-as.numeric(mapping$AIS_dis)
#2 is the code for improvers, this tells the function to call proportion of improvers

AIS_is_improved.gcon<-plot_net(var='AIS_is_improved_num','proportions','# of patients',2,prop.colors=c('#0000EE', '#FF3030', '#7F7F7F'),network_name = "AIS_is_improved_net", prop.labels = c("No","Yes", "NA"))

#AIS at discharge
proportion_colors<-c('blue','#329ce5', 'green','#ffa500','red','grey50')

# Extract proportions of AIS discharge by cluster (A=1, E=5)
AIS_dis1<-plot_net(var='AIS_dis_num', 'proportions','# of patients',1,
                             prop.colors=proportion_colors,
                   prop.labels = c("A","B","C","D","E", "NA"),
                   network_name = "AIS_disA")
AIS_dis2<-plot_net('AIS_dis_num', 'proportions','# of patients',2,
                              prop.colors=proportion_colors,
                   prop.labels = c("A","B","C","D","E", "NA"),
                   network_name = "AIS_disB")
AIS_dis3<-plot_net('AIS_dis_num', 'proportions','# of patients',3,
                              prop.colors=proportion_colors,
                   prop.labels = c("A","B","C","D","E", "NA"),
                   network_name = "AIS_disC")
AIS_dis4<-plot_net('AIS_dis_num', 'proportions','# of patients',4,
                              prop.colors=proportion_colors,
                   prop.labels = c("A","B","C","D","E", "NA"),
                   network_name = "AIS_disD")
AIS_dis5<-plot_net('AIS_dis_num', 'proportions','# of patients',5,
                              prop.colors=proportion_colors,
                   prop.labels = c("A","B","C","D","E", "NA"),
                   network_name = "AIS_disE")
```

<a href="#header">Back to top</a>

**Data frame for aggregated cluster scatter plots**

```{r fig.width=5, fig.height=5}
# data frame for Scatter plots of the aggregated values for each cluster
agg_plot_df<-data.frame(map=MAP.gcon$agg_values, 
                        AIS=AIS_is_improved.gcon$agg_values,
                        HR=HR.gcon$agg_values,
                        AIS_dis_A=AIS_dis1$agg_values,
                        AIS_dis_B=AIS_dis2$agg_values,
                        AIS_dis_C=AIS_dis3$agg_values,
                        AIS_dis_D=AIS_dis4$agg_values,
                        AIS_dis_E=AIS_dis5$agg_values)

agg_MAP_net_plot<-ggplot(agg_plot_df, aes(map, AIS, fill=AIS))+
  geom_point(shape=21, size=3)+
  scale_fill_gradientn(colors=net.color)+
  xlab('Mean cluster AverAge MAP (mmHg)')+
  ylab('Proportion of\n AIS improve.')+labs(fill='Proportion\nAIS improve.')+
  theme_minimal(base_size = 11)+
  geom_smooth(span=1, se=F)

agg_HR_net_plot<-ggplot(agg_plot_df, aes(HR, AIS, fill=AIS))+
  geom_point(shape=21, size=3)+
  scale_fill_gradientn(colors=net.color)+
  xlab('Mean cluster AverAge HR (bpm)')+
  ylab('Proportion of\n AIS improve.')+labs(fill='Proportion\nAIS improve.')+
  theme_minimal(base_size = 11)+
  geom_smooth(span=1, se=F)

AIS_disA_plot<-ggplot(agg_plot_df, aes(map, AIS_dis_A, fill=AIS_dis_A))+
  geom_point(shape=21, size=2)+
  scale_fill_gradientn(colors=net.color)+
  xlab('Mean cluster AverAge MAP (mmHg)')+
  ylab('Dis. proportion of\n AIS A')+labs(fill='% AIS A')+
  theme_minimal(base_size = 11)+
  geom_smooth(span=1, se=F)

AIS_disB_plot<-ggplot(agg_plot_df, aes(map, AIS_dis_B, fill=AIS_dis_B))+
  geom_point(shape=21, size=2)+
  scale_fill_gradientn(colors=net.color)+
  xlab('Mean cluster AverAge MAP (mmHg)')+
  ylab('Dis. proportion of\n AIS B')+labs(fill='% AIS B')+
  theme_minimal(base_size = 11)+
  geom_smooth(span=1, se=F)

AIS_disC_plot<-ggplot(agg_plot_df, aes(map, AIS_dis_C, fill=AIS_dis_C))+
  geom_point(shape=21, size=2)+
  scale_fill_gradientn(colors=net.color)+
  xlab('Mean cluster AverAge MAP (mmHg)')+
  ylab('Dis. proportion of\n AIS C')+labs(fill='% AIS C')+
  theme_minimal(base_size = 11)+
  geom_smooth(span=1, se=F)

AIS_disD_plot<-ggplot(agg_plot_df, aes(map, AIS_dis_D, fill=AIS_dis_D))+
  geom_point(shape=21, size=2)+
  scale_fill_gradientn(colors=net.color)+
  xlab('Mean cluster AverAge MAP (mmHg)')+
  ylab('Dis. proportion of\n AIS D')+labs(fill='% AIS D')+
  theme_minimal(base_size = 11)+
  geom_smooth(span=1, se=F)

AIS_disE_plot<-ggplot(agg_plot_df, aes(map, AIS_dis_E, fill=AIS_dis_E))+
  geom_point(shape=21, size=2)+
  scale_fill_gradientn(colors=net.color)+
  xlab('Mean cluster AverAge MAP (mmHg)')+
  ylab('Dis. proportion of\n AIS E')+labs(fill='% AIS E')+
  theme_minimal(base_size = 11)+
  geom_smooth(span=1, se=F)
```

<a href="#header">Back to top</a>

### Figure 2

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

## Flow

flow<-readPNG("Figure 2 flow.png") #external file added as supplement in the paper
flow<-wrap_elements(rasterGrob(flow, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")))

#Map net
map_net <- readPNG('map_net.png')
map_net_grob<-wrap_elements(rasterGrob(map_net, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1.2, "native")),clip = F)
#Map contracted net
map_net_con <- readPNG("map_net_con.png")
map_net_con_grob<-wrap_elements(rasterGrob(map_net_con, interpolate = T,
                                           width = unit(1, "native"), 
                       height = unit(1.2, "native")),clip = F)
#HR net
HR_net <- readPNG('HR_net.png')
HR_net_grob<-wrap_elements(rasterGrob(HR_net, interpolate = T,
                                      width = unit(1, "native"), 
                       height = unit(1.2, "native")), clip = F)
#HR contracted net
HR_net_con <- readPNG("HR_net_con.png")
HR_net_con_grob<-wrap_elements(rasterGrob(HR_net_con, interpolate = T,
                                          width = unit(1, "native"), 
                       height = unit(1.2, "native")),clip = F)

#AIS_is_improved net
AIS_is_improved_net <- readPNG('AIS_is_improved_net.png')
AIS_is_improved_net_grob<-wrap_elements(rasterGrob(AIS_is_improved_net, interpolate = T,
                                      width = unit(1, "native"), 
                       height = unit(1.2, "native")), clip = F)
#AIS_is_improved contracted net
AIS_is_improved_net_con <- readPNG("AIS_is_improved_net_con.png")
AIS_is_improved_net_con_grob<-wrap_elements(rasterGrob(AIS_is_improved_net_con, interpolate = T,width = unit(1, "native"), 
                       height = unit(1.2, "native")),clip = F)

#Legends
legend.map<-plot_legend(mapping$mean.map, net.color.Pal(250), "AverAge MAP\n(mmHg)", c(70, 85, 100), c("70", "85","100"))

legend.HR<-plot_legend(mapping$mean.HR, net.color.Pal(250), "AverAge HR\n(bpm)", c(60, 80, 100), c("60", "80","100"))

legend.AIS<-plot_legend(25:75, net.color.Pal(250), "Proportion\n AIS improve. (%)", c(25, 50, 75), c("<25", "50",">75"))

#temporal AIS ggplot to extract the legend
temp_AIS<-AIS_is_improved.gcon$plot+
  labs(fill="AIS improve.")+
  theme(legend.position = "right", 
        legend.title = element_text(), 
        legend.direction = "vertical")+
  guides(fill=guide_legend(nrow=3,byrow=TRUE))

legend_AIS_cat<-wrap_elements(ggplotGrob(temp_AIS)$grobs[[15]])

#Build figure
figure<-flow+
  map_net_grob+
  map_net_con_grob+inset_element(legend.map, left=-0.2, bottom = -0.1, 
                                   right = 0.3, top=0.4, ignore_tag = T)+
  MAP.gcon$plot+
  HR_net_grob+
  HR_net_con_grob+inset_element(legend.HR, left=-0.2, bottom = -0.1, 
                                   right = 0.3, top=0.4, ignore_tag = T)+
  HR.gcon$plot+
  AIS_is_improved_net_grob+
  inset_element(legend_AIS_cat, left=0.1, bottom = -0.1, 
                                   right = 0.3, top=0.4, ignore_tag = T)+
  AIS_is_improved_net_con_grob+inset_element(legend.AIS, left=-0.2, bottom = -0.1, 
                                   right = 0.3, top=0.4, ignore_tag = T)+
     AIS_is_improved.gcon$plot+theme(legend.position = "none")+
  wrap_elements(agg_MAP_net_plot,clip = F)+
  wrap_elements(agg_HR_net_plot, clip=F)

#Figure layout
layout <- c(
  # line 1
  area(t=1, l=1, r=6, b=2),
  # line 2
  area(t=3, l=1, r=2, b=4),area(t=3, l=3, r=4, b=4),area(t=3, l=5, r=6, b=4),
  # line 3
  area(t=5, l=1, r=2, b=6),area(t=5, l=3, r=4, b=6),area(t=5, l=5, r=6, b=6),
  # line 4
  area(t=7, l=1, r=2, b=8),area(t=7, l=3, r=4, b=8),area(t=7, l=5, r=6, b=8),
  # line 5
  area(t=9, l=1, r = 3, b=10),area(t=9, l=4, r = 6, b=10)
)

#Add figure options
figure2<-figure+plot_layout(design = layout)+
  plot_annotation(tag_levels = 'a',theme=theme(text = element_text(size=13)))

## Save figure
ggsave(figure2,file="figure 2.png", width = 12, height = 15.2, dpi = "retina")

figure2
```

<a href="#header">Back to top</a>

### Figure 2 supplement 1
```{r fig.height=8, fig.width=10}
## read dendogram for panel h
dendo <- readPNG('dendogram.png')
dendo_grob<-rasterGrob(dendo, interpolate = T,width = unit(1.2, "native"), 
                       height = unit(1, "native"))

## compose figure
figure2sup1<-(fig2_sup1_a+ggtitle("Data processing")|fig2_sup1_b)/
  (fig2_sup1_c+ggtitle("Optimization Isomap solution")|fig2_sup1_d|fig2_sup1_e)/
  (fig2_sup1_f+ggtitle("Walktrap solution")|fig2_sup1_g|dendo_grob)+
  plot_annotation(tag_levels = 'a')

## Save figure
ggsave(figure2sup1,file="figure2sup1.png", width = 10, height = 8)

figure2sup1
```

<a href="#header">Back to top</a>

### Figure 2 supplement 2

```{r fig.height=8,fig.width=12}
#AIS dis A
AIS_disA_con_net <- readPNG('AIS_disA_con.png')
AIS_disA_con_net_grob<-wrap_elements(rasterGrob(AIS_disA_con_net, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")),clip = F)
#AIS dis B
AIS_disB_con_net <- readPNG('AIS_disB_con.png')
AIS_disB_con_net_grob<-wrap_elements(rasterGrob(AIS_disB_con_net, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")),clip = F)
#AIS dis C
AIS_disC_con_net <- readPNG('AIS_disC_con.png')
AIS_disC_con_net_grob<-wrap_elements(rasterGrob(AIS_disC_con_net, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")),clip = F)
#AIS dis D
AIS_disD_con_net <- readPNG('AIS_disD_con.png')
AIS_disD_con_net_grob<-wrap_elements(rasterGrob(AIS_disD_con_net, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")),clip = F)
#AIS dis E
AIS_disE_con_net <- readPNG('AIS_disE_con.png')
AIS_disE_con_net_grob<-wrap_elements(rasterGrob(AIS_disE_con_net, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")),clip = F)

## Compose figure
figure<-AIS_disA_con_net_grob+wrap_elements(AIS_disA_plot, ignore_tag = T)+
  AIS_disB_con_net_grob+wrap_elements(AIS_disB_plot, ignore_tag = T)+
  AIS_disC_con_net_grob+wrap_elements(AIS_disC_plot, ignore_tag = T)+
  AIS_disD_con_net_grob+wrap_elements(AIS_disD_plot, ignore_tag = T)+
  AIS_disE_con_net_grob+wrap_elements(AIS_disE_plot, ignore_tag = T)+
  wrap_elements(AIS_dis1$plot+labs(fill="AIS dis.")+
    guides(fill=guide_legend(ncol = 1,byrow=F))+
    theme(legend.position = "right", legend.title = element_text()))

## Figure layout
layout <- c(
  ## line 1
  area(t=1, l=1, r=2, b=2),area(t=1, l=3, r=4, b=2),area(t=1, l=5, r=6, b=2),
  area(t=1, l=7, r=8, b=2),
  ## line 2
  area(t=3, l=1, r=2, b=4),area(t=3, l=3, r=4, b=4),area(t=3, l=5, r=6, b=4),
  area(t=3, l=7, r=8, b=4),
  ## line 3
  area(t=5, l=1, r=2, b=6),area(t=5, l=3, r = 4, b=6),area(t=5, l=5, r = 8, b=6)
)

## Add figure options
figure2sup2<-figure+plot_layout(design = layout)+
  plot_annotation(tag_levels = 'a')

## Save figure
ggsave(figure2sup2,file="figure2sup2.png", width = 12, height = 8, dpi = "retina")

figure2sup2
```

<a href="#header">Back to top</a>

# MAP has a non-linear association with probability of recovery
## Logistic regression

We build several models and compare them by Likelihood ratio test and cross validation errors. See methods in the paper for considered models.

```{r}
## build dataframe for regression analysis
regression.df<-raw.df[!is.na(raw.df$AIS_is_improved),]
regression.df$AIS_is_improved<-as.character(regression.df$AIS_is_improved)

## specify outcome metric as 0: NO and 1: YES
regression.df$AIS_is_improved[regression.df$AIS_is_improved=='NO']<-0
regression.df$AIS_is_improved[regression.df$AIS_is_improved=='YES']<-1

regression.df$AIS_is_improved<-as.numeric(regression.df$AIS_is_improved)
regression.df$days_post_surgery_dish<-as.numeric(as.character(regression.df$days_post_surgery_dish))

## Fit the models
fit0<-glm(AIS_is_improved~1, data=regression.df, family = 'binomial', na.action = 'na.exclude')
fit1<-glm(AIS_is_improved~poly(mean.map,1), data=regression.df, family = 'binomial', na.action = 'na.exclude')
fit2<-glm(AIS_is_improved~poly(mean.map,2), data=regression.df, family = 'binomial', na.action = 'na.exclude')
fit3<-glm(AIS_is_improved~poly(mean.map,3), data=regression.df, family = 'binomial', na.action = 'na.exclude')
fit.ns2<-glm(AIS_is_improved~splines::ns(mean.map,2), data=regression.df, family = 'binomial', na.action = 'na.exclude')
fit.ns3<-glm(AIS_is_improved~splines::ns(mean.map,3), data=regression.df, family = 'binomial', na.action = 'na.exclude')
```

<a href="#header">Back to top</a>

### Likelihood ratio test

```{r}
anova(fit0, fit1, fit2, fit3, test ='LRT') # naive, linear and polynomials
anova(fit1, fit.ns2, test='LRT') # linear vs. n-spline 2
anova(fit2,fit.ns3, test='LRT') # quadratic vs. n-spline 3
```

### LOOCV cross validation

Leave-one-out cross validation

```{r }
set.seed(666)
cv.fit0<-cv.glm(regression.df, fit0) #the default K-CV is the number of observations, which means leave-one-out CV
cv.fit1<-cv.glm(regression.df, fit1)
cv.fit2<-cv.glm(regression.df, fit2)
cv.fit3<-cv.glm(regression.df, fit3)
cv.fit.ns2<-cv.glm(regression.df, fit.ns2)
cv.fit.ns3<-cv.glm(regression.df, fit.ns3)

CVerror0<-cv.fit0$delta[1] #Cross validation error
CVerror1<-cv.fit1$delta[1]#Cross validation error
CVerror2<-cv.fit2$delta[1]#Cross validation error
CVerror3<-cv.fit3$delta[1]#Cross validation error
CVerrorns2<-cv.fit.ns2$delta[1]#Cross validation error
CVerrorns3<-cv.fit.ns3$delta[1]#Cross validation error

model_names<-c('naive','linear','quadratic','cubic','n splin 2','n splin 3')
poly_cv_df<-data.frame(model=factor(model_names, levels=model_names),
                       CVerror=c(CVerror0,CVerror1,CVerror2,CVerror3,CVerrorns2,CVerrorns3))

poly_cv_df
```

<a href="#header">Back to top</a>

# Factors influencing MAP association with recovery

## FULL model LOOCV

Fitting the full model with LOOCV
```{r }
## process the data for regression
regression.df$AIS_ad<-factor(regression.df$AIS_ad,
                             levels=c("A","B","C","D","E"))
regression.df$AIS_ad<-relevel(regression.df$AIS_ad, "E")

regression.df.f<-regression.df%>%select(mean.map, mean.HR, length_surg, days_post_surgery_dish, Age, AIS_ad,AIS_is_improved, Injury_level)

## format outcome variable as factor
regression.df.f$AIS_is_improved<-as.factor(regression.df.f$AIS_is_improved)

regression.df.f<-regression.df.f%>%
  mutate(AIS_is_improved_c = ifelse(AIS_is_improved==0,"NO","YES"))

## fit full model with LOOCV
ctrl <- trainControl(method = "LOOCV",savePredictions = "all",classProbs = T,
                     number = 1)

mod_fit <- train(AIS_is_improved_c ~ poly(mean.map,2)+mean.HR+length_surg+days_post_surgery_dish+Age+AIS_ad,
                 data=regression.df.f, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

#Prediction probability
prob_predict_improve<-mod_fit$pred[,4]

summary(mod_fit)
```

<a href="#header">Back to top</a>

## Only poly MAP LOOCV

Fitting the model that only contains MAP (no covariates)
```{r}
mod_fit_map <- train(AIS_is_improved_c ~ poly(mean.map,2),
                 data=regression.df.f, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

#Prediction probability
prob_predict_improve_map<-mod_fit_map$pred[,4]

summary(mod_fit_map)
```

<a href="#header">Back to top</a>

## Covariates only

Fitting the covariate only model
```{r}
mod_fit_cov <- train(AIS_is_improved_c ~mean.HR+length_surg+days_post_surgery_dish+Age+AIS_ad,
                 data=regression.df.f, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

#Prediction probability
prob_predict_improve_cov<-mod_fit_cov$pred[,4]

summary(mod_fit_cov)
```

<a href="#header">Back to top</a>

### Figure 3
**Plot LOOCV estimated error all models and polynomial**

```{r fig.width=4, fig.height=6}
## LOOCV model plot
CV_models_plot<-ggplot(poly_cv_df, aes(model, CVerror))+
  geom_point()+
  geom_line(aes(as.numeric(model)))+
  theme_minimal()+
  ylab('Estimated LOOCV error')

## Plot polynomial

poly_plot<-regression.df%>%
  ggplot(aes(mean.map, AIS_is_improved))+
  geom_point(shape=16, color='black', alpha=0.2)+
  geom_smooth(method = "glm",method.args = list(family = "binomial"),
              se = T, formula = y~poly(x,2))+
  xlab('Mean MAP during surgery (mmHg)')+ylab(expression('Pr('*Delta*'AIS grade > 0)'))+
  theme_minimal(base_size = 11)

## add options to the figure
figure3<-CV_models_plot/poly_plot+plot_annotation(tag_levels = "a")

## save figure
ggsave(figure3, file="figure3.png", width = 4, height = 6)

figure3

```

<a href="#header">Back to top</a>

## Only AIS-A

Models fitted to the AIS-A subcohort

### Full LOOCV model A

```{r }
## process the data for regression
regression.df.fA<-regression.df.f%>%
  filter(AIS_ad=="A")

mod_fit_A_poly <- train(AIS_is_improved_c ~ poly(mean.map,2)+mean.HR+length_surg+days_post_surgery_dish+Age,
                 data=regression.df.fA, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')
summary(mod_fit_A_poly)
```

<a href="#header">Back to top</a>

## NLI

Models fitted with the NLI split cohorts

### Model LOOCV Cervical

```{r }
## process the data for regression
regression.df.fCer<-regression.df.f%>%
  filter(Injury_level=="Cervical")

regression.df.fCer<-droplevels(regression.df.fCer)

mod_fit_Cer <- train(AIS_is_improved_c ~ poly(mean.map,2)+mean.HR+length_surg+
                     days_post_surgery_dish+Age+AIS_ad,
                 data=regression.df.fCer, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')
summary(mod_fit_Cer)
```

### Model LOOCV Other

```{r }
## process the data for regression
regression.df.fOther<-regression.df.f%>%
  filter(Injury_level=="Other")

regression.df.fOther<-droplevels(regression.df.fOther)

mod_fit_Other <- train(AIS_is_improved_c ~ poly(mean.map,2)+mean.HR+length_surg+
                     days_post_surgery_dish+Age+AIS_ad,
                 data=regression.df.fOther, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')
summary(mod_fit_Other)
```

<a href="#header">Back to top</a>

## Days to discharge

### AIS improvement by days to discharge

```{r}
summary(glm(AIS_is_improved ~ days_post_surgery_dish,
                 data=regression.df.f,family="binomial"))
```

#### Interaction with MAP

```{r}
summary(glm(AIS_is_improved ~ days_post_surgery_dish*poly(mean.map, 2),
                 data=regression.df.f,family="binomial"))
```


### Days to discharge by quadratic aMAP

```{r}
summary(lm(days_post_surgery_dish ~ poly(mean.map, 2), data = regression.df.f))
```

### Full model - days to discharge

```{r}
regression.df.fnoDays<-regression.df.f%>%
  select(AIS_is_improved_c, mean.map, mean.HR, length_surg, Age, AIS_ad, days_post_surgery_dish)%>%
  .[complete.cases(.),]

full_model<-train(AIS_is_improved_c ~ poly(mean.map,2)+mean.HR+days_post_surgery_dish+
                    length_surg+Age+AIS_ad,
                 data=regression.df.fnoDays, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

mod_fit_noDays <- train(AIS_is_improved_c ~ poly(mean.map,2)+mean.HR+length_surg+Age+AIS_ad,
                 data=regression.df.fnoDays, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

## Likelihood ratio test
anova(full_model$finalModel, mod_fit_noDays$finalModel, test = "LRT")

summary(full_model)
summary(mod_fit_noDays)
```


# Intraoperative MAP range from 76-[104-117] mmHg predicts recovery

The following section of the code works the MAP threshold determination as in the paper

## Symetric threshold
### Generate moving window vars

This section generates variables of the time inside/outside a increasing window centered on the mean MAP with higher prob. of improvement (90). It increments the window by 1mmHg each side at the time, so the first window is 89-91, second would be 88-92 and so on.

```{r fig.width=5, fig.height=3}
temp.list<-list()

## Calculate time out of MAP range for each range
for (t in 1:30){
  temp.df<-map_df_long%>%
    filter(!is.na(value))%>%
    group_by(Subject_ID)%>%
    summarise(num=sum(between(value, -Inf,90-t),between(value, 90+t,Inf)), time=num*5)
  colnames(temp.df)<-c('Subject_ID', 
                       paste('q5_MAP_out_',t, sep=''),
                       paste('Range ',90-t,'-',90+t, sep=''))
  temp.list[[t]]<-temp.df
}

## build a single dataframe for analysis
OR_MAP_threshold_out<-Reduce(function(x, y) merge(x, y, by='Subject_ID'), temp.list)
OR_MAP_threshold_out$AIS_is_improved<-raw.df$AIS_is_improved
OR_MAP_threshold_out_long<-OR_MAP_threshold_out%>%
  pivot_longer(-c(Subject_ID,AIS_is_improved))

# plot for window range schema
schema.plot.df<-map_df_long%>%filter(Subject_ID=='TSretro110')
schema.plot.df$n_q5<-1:dim(schema.plot.df)[1]

middle<-90
mw_symetric_plot<-ggplot(schema.plot.df, aes(n_q5, value))+
  geom_line()+
  coord_cartesian(xlim = c(0,120), ylim=c(69,115))+
  theme_minimal()+
  xlab('Time bin (Q5)')+ylab('MAP (mmHg)')+
  annotate("rect", xmin=85, xmax=130, ymin=65, ymax=115, fill="white")+
  annotate("rect", xmin=90, xmax=120, ymin=75, ymax=105, fill="grey85")+
  annotate("polygon", x=c(90, 120, 120, 90),y=c(89, 75, 105, 91), fill="white")+ 
  geom_hline(yintercept = middle, linetype='dashed')+
  annotate("segment", x=0, xend=90, y=middle+1, yend=middle+1, color="orange")+
  annotate("segment", x=0, xend=90, y=middle-1, yend=middle-1, color="orange")+
  annotate("segment", x=0, xend=100, y=middle+5, yend=middle+5, color="red")+
  annotate("segment", x=0, xend=100, y=middle-5, yend=middle-5, color="red")+
  annotate("segment", x=0, xend=110, y=middle+10, yend=middle+10, color="purple")+
  annotate("segment", x=0, xend=110, y=middle-10, yend=middle-10, color="purple")+
  annotate("segment", x=90, xend=90, y=89, yend=91, color="orange")+
  annotate("segment", x=100, xend=100, y=85, yend=95, color="red")+
  annotate("segment", x=110, xend=110, y=80, yend=100, color="purple")+
  annotate("segment", x=120, xend=120, y=75, yend=105, 
           color="black", linetype="dashed", 
           arrow=arrow(ends = "both",type = "closed", length = unit(2,"mm")))+
  annotate("text", x=105, y=110, label="Increase\n1 mmHg", size=4)+
  annotate("text", x=90, y=72, label="Time\noutside range", size=3, hjust=0)+
  annotate("rect",xmin=105, xmax=110, ymin=72, ymax=74, fill="grey85")+
  scale_x_continuous(breaks = seq(0,70,20))

mw_symetric_plot
```

<a href="#header">Back to top</a>

### Assortativity of the range

Here we measure the assortativity of the network for the time outside each range

```{r }
A_time_MAP_range<-data.frame(variable='', A=NA)
vars<-colnames(OR_MAP_threshold_out)[
  grepl(x = colnames(OR_MAP_threshold_out),pattern = 'Range')]

for(i in 1:length(vars)){
  temp<-assortativity(g,OR_MAP_threshold_out[,vars[i]],directed=F)
  temp.df<-data.frame(variable=vars[i],A=temp)
  A_time_MAP_range<-rbind(A_time_MAP_range, temp.df)
}

A_time_MAP_range<-A_time_MAP_range[-1,]
A_time_MAP_range$variable_num<-as.numeric(as.factor(A_time_MAP_range$variable))
A_time_MAP_range$variable<-str_replace(A_time_MAP_range$variable, 
                                       'Range ','')

labels_range<-A_time_MAP_range$variable
labels_range[seq(2,length(labels_range),2)]<-''
```

<a href="#header">Back to top</a>

### Plot assortativity coeff

```{r fig.width=4, fig.height=3}
ac_range_plot<-ggplot(A_time_MAP_range, aes(variable_num, rev(A)))+
  geom_point()+
  scale_x_continuous(breaks =rev(A_time_MAP_range$variable_num),
                     labels = labels_range)+
  geom_vline(xintercept = 17, linetype='dashed')+
  theme_minimal()+
  geom_smooth(se=F, color='red', span=0.4, size=0.5)+
  ylab('Assortativity coeff. (Ar)')+xlab('MAP range (mmHg)')+
  theme(axis.text.x = element_text(angle = 90, 
                                   hjust = 0.5, vjust = 0.5, size=9),
        axis.ticks.x = element_line(color='black'))

# ggsave('Ar range.png', width = 3, height = 2.5, dpi='print')
ac_range_plot
```

Network for range 73-107 found by assortativity analysis

```{r fig.width=4, fig.height=3}
mapping$range_73_107<-OR_MAP_threshold_out$`Range 73-107`
range_73_107.net<-plot_net('range_73_107','mean','xx', ranges=c(1,30),
                           network_name = "range_73_107_net")

agg_plot_df$range_73_107<-range_73_107.net$agg_values

agg_plot_73_107<-ggplot(agg_plot_df, aes(range_73_107, AIS, fill=range_73_107))+
  geom_point(shape=21, size=2)+
  scale_fill_gradientn(colors=net.color)+
  xlab('Mean cluster time out of MAP\nrange 73-107 mmHg')+
  ylab('Proportion of\n AIS improve.')+labs(fill='Mean\ntime out')+
  theme_minimal(base_size = 11)+
  geom_smooth(span=1, se=F)

agg_plot_73_107
```

### LASSO

logistic regression with LASSO regularization and cross-validation to find the range of MAP that is most predictive of AIS improvement.

#### Data processing for LASSO

```{r }
lasso.df<-OR_MAP_threshold_out%>%
  filter(!is.na(AIS_is_improved))

vars<-colnames(lasso.df)[grepl(x = colnames(lasso.df), 
                               pattern = 'Range')]
lasso.x<-as.matrix(lasso.df%>%dplyr::select(vars))
```

#### Fit logistic LASSO regression

Lambda is determined by cross validation. N folds controls for the number of folds CV. Setting nfolds for 103 (number of observations) implies doing LOOCV for lambda (see glmnet and glmnetUtils packAges for details)

```{r fig.width=4, fig.height=3}
logLasso<-cv.glmnet(x=lasso.x, y=lasso.df$AIS_is_improved,family='binomial',
                    alpha=1,nfolds=103)
#Coefficients
coef(logLasso$glmnet.fit,s=logLasso$lambda.1se)
```

#### LASSO range cross validation
Logistic cross LASSO range solution validation

```{r}
log.CV.df<-OR_MAP_threshold_out%>%
  select(`Range 76-104`,AIS_is_improved)%>%
  filter(!is.na(AIS_is_improved))

ctrl <- trainControl(method = "LOOCV",savePredictions = TRUE, classProbs = T)

mod_fit_sym <- train(AIS_is_improved ~ `Range 76-104`,log.CV.df, method="glm", family="binomial",trControl = ctrl)

summary(mod_fit_sym)
mod_fit_sym
```

## Upper limit
### Generate asymetric moving window

This section generates variables of the time inside/outside a increasing window fixed on the mean MAP 76 on the lower bound. It increments the window by 1mmHg each time, so the first window is 76-77, second would be 76-78 and so on.

```{r fig.width=5, fig.height=3}
temp.list<-list()
# map_df_long<-map.wide%>%pivot_longer(-Subject_ID,names_to = 'variable', values_to = 'value')

## calculate time outside MAP range for each range
for (t in 1:44){
  temp.df<-map_df_long%>%
    filter(!is.na(value))%>%
    group_by(Subject_ID)%>%
    summarise(num=sum(between(value, -Inf,76),between(value, 76+t,Inf)), time=num*5)
  colnames(temp.df)<-c('Subject_ID', 
                       paste('q5_MAP_out_',t, sep=''),
                       paste('Range ',76,'-',76+t, sep=''))
  temp.list[[t]]<-temp.df
}

## build single dataframe for analysis
OR_MAP_threshold_out_upper<-Reduce(function(x, y) merge(x, y, by='Subject_ID'), temp.list)
OR_MAP_threshold_out_upper$AIS_is_improved<-raw.df$AIS_is_improved
OR_MAP_threshold_out_upper_long<-OR_MAP_threshold_out_upper%>%
  pivot_longer(-c(Subject_ID,AIS_is_improved))

# plot for window range schema
schema.plot.df_upper<-map_df_long%>%filter(Subject_ID=='TSretro110')
schema.plot.df_upper$n_q5<-1:dim(schema.plot.df_upper)[1]

middle<-76
mw_asymetric_plot<-ggplot(schema.plot.df_upper, aes(n_q5, value))+
  geom_line()+
  coord_cartesian(xlim = c(0,120), ylim=c(69,115))+
  theme_minimal()+
  xlab('Time bin (Q5)')+ylab('MAP (mmHg)')+
  annotate("rect", xmin=85, xmax=130, ymin=76, ymax=115, fill="white")+
  annotate("rect", xmin=90, xmax=120, ymin=65, ymax=91, fill="grey85")+
  annotate("polygon", x=c(90, 120, 120, 90),y=c(76, 76, 91, 77), fill="white")+ 
  geom_hline(yintercept = middle, linetype='dashed')+
  annotate("segment", x=0, xend=90, y=middle+1, yend=middle+1, color="orange")+
  annotate("segment", x=0, xend=100, y=middle+5, yend=middle+5, color="red")+
  annotate("segment", x=0, xend=110, y=middle+10, yend=middle+10, color="purple")+
  annotate("segment", x=90, xend=90, y=76, yend=77, color="orange")+
  annotate("segment", x=100, xend=100, y=76, yend=81, color="red")+
  annotate("segment", x=110, xend=110, y=76, yend=86, color="purple")+
  annotate("segment", x=120, xend=120, y=76, yend=91, 
           color="black", linetype="dashed", 
           arrow=arrow(ends = "both",type = "closed", length = unit(2,"mm")))+
  annotate("text", x=105, y=96, label="Increase\n1 mmHg", size=4)+
  annotate("text", x=60, y=72, label="Time\noutside range", size=3, hjust=0)+
  annotate("rect",xmin=75, xmax=80, ymin=72, ymax=74, fill="grey85")+
  scale_x_continuous(breaks = seq(0,70,20))

mw_asymetric_plot
```

### Assortativity of the range

```{r }
A_time_MAP_range_upper<-data.frame(variable='', A=NA)
vars_upper<-colnames(OR_MAP_threshold_out_upper)[
  grepl(x = colnames(OR_MAP_threshold_out_upper),pattern = 'Range')]

for(i in 1:length(vars_upper)){
  temp<-assortativity(g,OR_MAP_threshold_out_upper[,vars_upper[i]],directed=F)
  temp.df<-data.frame(variable=vars_upper[i],A=temp)
  A_time_MAP_range_upper<-rbind(A_time_MAP_range_upper, temp.df)
}

A_time_MAP_range_upper<-A_time_MAP_range_upper[-1,]
A_time_MAP_range_upper$variable_num<-1:44
A_time_MAP_range_upper$variable<-str_replace(A_time_MAP_range_upper$variable, 
                                       'Range ','')

labels_range_upper<-A_time_MAP_range_upper$variable
labels_range_upper[seq(2,length(labels_range_upper),2)]<-''
```

### Plot assortativity coeff

```{r fig.width=4, fig.height=3}
ac_range_plot_upper<-ggplot(A_time_MAP_range_upper, aes(variable_num, A))+
  geom_point()+
  scale_x_continuous(breaks =A_time_MAP_range_upper$variable_num,
                     labels = labels_range_upper)+
  geom_vline(xintercept = 40, linetype='dashed')+
  theme_minimal()+
  geom_smooth(se=F, color='red', span=0.4, size=0.5)+
  ylab('Assortativity coeff. (Ar)')+xlab('MAP range (mmHg)')+
  theme(axis.text.x = element_text(angle = 90, 
                                   hjust = 0.5, vjust = 0.5, size=9),
        axis.ticks.x = element_line(color='black'))

# ggsave('Ar range.png', width = 3, height = 2.5, dpi='print')
ac_range_plot_upper
```

Network for range 76-116 found by assortativity analysis

```{r fig.width=4, fig.height=3}
mapping$range_76_116<-OR_MAP_threshold_out_upper$`Range 76-116`
range_76_116.net<-plot_net('range_76_116','mean','xx', ranges=c(1,30),
                           network_name = "range_76_116_net")

agg_plot_df$range_76_116<-range_76_116.net$agg_values

agg_plot_76_116<-ggplot(agg_plot_df, aes(range_76_116, AIS, fill=range_76_116))+
  geom_point(shape=21, size=2)+
  scale_fill_gradientn(colors=net.color)+
  xlab('Mean cluster time out of MAP\nrange 76-116 mmHg')+
  ylab('Proportion of\n AIS improve.')+labs(fill='Mean\ntime out')+
  theme_minimal(base_size = 11)+
  geom_smooth(span=1, se=F)

agg_plot_76_116
```

### LASSO upper limit

#### Data processing for LASSO

```{r }
lasso.df_upper<-OR_MAP_threshold_out_upper%>%
  filter(!is.na(AIS_is_improved))

vars_upper<-colnames(lasso.df_upper)[grepl(x = colnames(lasso.df_upper), 
                               pattern = 'Range')]
lasso.x_upper<-as.matrix(lasso.df_upper%>%dplyr::select(vars_upper))
```

#### Fit logistic LASSO regression

Lambda is determined by cross validation. N folds controls for the number of folds CV. Setting nfolds for 103 (number of observations) implies doing LOOCV for lambda (see glmnet and glmnetUtils packAges for details)

```{r }
logLasso_upper<-cv.glmnet(x=lasso.x_upper, y=lasso.df_upper$AIS_is_improved,family='binomial',
                    alpha=1,nfolds=103)

#Coefficients
coef(logLasso_upper$glmnet.fit,s=logLasso_upper$lambda.1se)
```

#### LASSO range cross validation

```{r}
log.CV.df_upper<-OR_MAP_threshold_out_upper%>%
  select(`Range 76-117`,AIS_is_improved)%>%
  filter(!is.na(AIS_is_improved))

ctrl <- trainControl(method = "LOOCV",savePredictions = TRUE, classProbs = T)

mod_fit_asym <- train(AIS_is_improved ~ `Range 76-117`,log.CV.df_upper, method="glm", family="binomial",trControl = ctrl)

summary(mod_fit_asym)
mod_fit_asym
```

## Figure 4

```{r fig.width=12, fig.height=7}
OR_MAP_threshold_out$AIS_is_improved<-as.character(OR_MAP_threshold_out$AIS_is_improved)
OR_MAP_threshold_out$AIS_is_improved[OR_MAP_threshold_out$AIS_is_improved=='NO']<-0
OR_MAP_threshold_out$AIS_is_improved[OR_MAP_threshold_out$AIS_is_improved=='YES']<-1
OR_MAP_threshold_out$AIS_is_improved<-as.numeric(OR_MAP_threshold_out$AIS_is_improved)

range_log_plot<-OR_MAP_threshold_out%>%filter(!is.na(AIS_is_improved))%>%
  select(`Range 76-104`,AIS_is_improved)%>%
  ggplot(aes(`Range 76-104`, AIS_is_improved))+
  geom_point(shape=16, color='black', alpha=0.2)+
  geom_smooth(method = "glm",method.args = list(family = "binomial"),
              se = T)+
  xlab('Estimated time outside\n 76-104mmHg MAP range (minutes)')+ylab(expression('Pr('*Delta*'AIS grade > 0)'))+
  annotate("text", label = "Time MAP outside\n76-104mmHg", x=300, y=0.8)+
  ggtitle("Confirmatory Logistic LASSO")+
  theme_minimal()

OR_MAP_threshold_out_upper$AIS_is_improved<-as.character(OR_MAP_threshold_out_upper$AIS_is_improved)
OR_MAP_threshold_out_upper$AIS_is_improved[OR_MAP_threshold_out_upper$AIS_is_improved=='NO']<-0
OR_MAP_threshold_out_upper$AIS_is_improved[OR_MAP_threshold_out_upper$AIS_is_improved=='YES']<-1
OR_MAP_threshold_out_upper$AIS_is_improved<-as.numeric(OR_MAP_threshold_out_upper$AIS_is_improved)

range_log_plot_upper<-OR_MAP_threshold_out_upper%>%filter(!is.na(AIS_is_improved))%>%
  select(`Range 76-117`,AIS_is_improved)%>%
  ggplot(aes(`Range 76-117`, AIS_is_improved))+
  geom_point(shape=16, color='black', alpha=0.2)+
  geom_smooth(method = "glm",method.args = list(family = "binomial"),
              se = T)+
  xlab('Estimated time outside\n76-117mmHg MAP range (minutes)')+ylab(expression('Pr('*Delta*'AIS grade > 0)'))+
  annotate("text", label = "Time MAP outside\n76-117mmHg", x=300, y=0.8)+
  ggtitle("Confirmatory Logistic LASSO")+
  theme_minimal()

figure4<-(mw_symetric_plot+ggtitle("Symmetric MAP threshold")|ac_range_plot+ggtitle("Exploratory network MAP range")|range_log_plot)/
  (mw_asymetric_plot+ggtitle("Upper limit threshold")|ac_range_plot_upper+ggtitle("Exploratory network MAP range")|range_log_plot_upper)

figure4<-figure4+plot_annotation(tag_levels = 'a')

ggsave(figure4,file="figure 4.png", width = 12, height = 7)

figure4
```

### Fig4 supplement 1

```{r fig.width=12, fig.height=7}
fig4s1_a<-readPNG("range_73_107_net.png")
fig4s1_a<-wrap_elements(rasterGrob(fig4s1_a, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")))

fig4s1_b<-readPNG("range_73_107_net_con.png")
fig4s1_b<-wrap_elements(rasterGrob(fig4s1_b, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")))

fig4s1_d<-readPNG("range_76_116_net.png")
fig4s1_d<-wrap_elements(rasterGrob(fig4s1_d, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")))

fig4s1_e<-readPNG("range_76_116_net_con.png")
fig4s1_e<-wrap_elements(rasterGrob(fig4s1_e, interpolate = T,
                                       width = unit(1, "native"), 
                       height = unit(1, "native")))

fig4s1<-fig4s1_a+fig4s1_b+wrap_elements(agg_plot_73_107)+
           fig4s1_d+fig4s1_e+wrap_elements(agg_plot_76_116)

layout <- c(
  # line 1
  area(t=1, l=1, r=2, b=2),area(t=1, l=3, r=4, b=2),area(t=1, l=5, r=7, b=2),
  # line 2
  area(t=3, l=1, r=2, b=4),area(t=3, l=3, r=4, b=4),area(t=3, l=5, r=7, b=4)
)

fig4s1<-fig4s1+plot_layout(design = layout)+plot_annotation(tag_levels = 'a',theme = theme(text = element_text(size=13)))

ggsave(fig4s1,file="fig4 supplement 1.png", width = 12, height = 6)

fig4s1
```

### Fig4 supplement 2

```{r fig.width=8, fig.height=10}
#Plot LASSO solutions
fig4s2_a<-rasterGrob({
  png(filename = "temp.png",res = 150)
    plot(logLasso)
  dev.off()
  readPNG("temp.png")
},interpolate = T,width = unit(1, "native"), height = unit(1, "native"))

fig4s2_b<-rasterGrob({
  png(filename = "temp.png",res = 150)
    plot(logLasso$glmnet.fit, xvar='lam',label = T)
  dev.off()
  readPNG("temp.png")
},interpolate = T,width = unit(1, "native"), height = unit(1, "native"))


# Plotting the individual panels
fig4s2_c<-OR_MAP_threshold_out%>%
  select(colnames(OR_MAP_threshold_out)[
    str_detect(colnames(OR_MAP_threshold_out),"Range")],
    AIS_is_improved)%>%
  pivot_longer(-AIS_is_improved)%>%
  mutate(name = factor(name, levels = unique(name)))%>%
  ggplot(aes(value, AIS_is_improved))+
    geom_point(shape=16, color='black', alpha=0.2)+
    geom_smooth(method = "glm",method.args = list(family = "binomial"),
              se = T)+
    xlab('Estimated time outside MAP')+
    ylab(expression('Pr('*Delta*'AIS grade > 0)'))+
  theme_minimal(base_size = 11)+
  facet_wrap(~name)

fig4s2<-wrap_elements(fig4s2_a)+wrap_elements(fig4s2_b)+fig4s2_c

layout <- c(
  # line 1
  area(t=1, l=1, r=2, b=2),area(t=1, l=3, r=4, b=2),
  # line 2
  area(t=3, l=1, r=4, b=6)
)

fig4s2<-fig4s2+plot_layout(design = layout)+plot_annotation(tag_levels = 'a',theme = theme(text = element_text(size=13)))

ggsave(fig4s2,file="fig4 supplement 2.png", width = 8, height = 10)

fig4s2
```

### Fig4 supplement 3
```{r fig.width=8, fig.height=10}
#Plot LASSO solutions
fig4s3_a<-rasterGrob({
  png(filename = "temp.png",res = 150)
    plot(logLasso_upper)
  dev.off()
  readPNG("temp.png")
},interpolate = T,width = unit(1, "native"), height = unit(1, "native"))

fig4s3_b<-rasterGrob({
  png(filename = "temp.png",res = 150)
    plot(logLasso_upper$glmnet.fit, xvar='lam',label = T)
  dev.off()
  readPNG("temp.png")
},interpolate = T,width = unit(1, "native"), height = unit(1, "native"))


# Plotting the individual panels
fig4s3_c<-OR_MAP_threshold_out_upper%>%
  select(colnames(OR_MAP_threshold_out_upper)[
    str_detect(colnames(OR_MAP_threshold_out_upper),"Range")],
    AIS_is_improved)%>%
  pivot_longer(-AIS_is_improved)%>%
  mutate(name = factor(name, levels = unique(name)))%>%
  ggplot(aes(value, AIS_is_improved))+
    geom_point(shape=16, color='black', alpha=0.2)+
    geom_smooth(method = "glm",method.args = list(family = "binomial"),
              se = T)+
    xlab('Estimated time outside MAP')+
    ylab(expression('Pr('*Delta*'AIS grade > 0)'))+
  theme_minimal(base_size = 11)+
  facet_wrap(~name)

fig4s3<-wrap_elements(fig4s3_a)+wrap_elements(fig4s3_b)+fig4s3_c

layout <- c(
  # line 1
  area(t=1, l=1, r=2, b=2),area(t=1, l=3, r=4, b=2),
  # line 2
  area(t=3, l=1, r=4, b=6)
)

fig4s3<-fig4s3+plot_layout(design = layout)+plot_annotation(tag_levels = 'a',theme = theme(text = element_text(size=13)))

ggsave(fig4s3,file="fig4 supplement 3.png", width = 8, height = 10)

fig4s3
```


# Prediction models

## Predicting AIS improvement
### Data preprocessing

Process the data for prediction model
```{r }
## build dataframe for analysis
prediction.df<-regression.df%>%
  left_join(OR_MAP_threshold_out_upper[,c("Subject_ID","Range 76-117")])%>%
  left_join(OR_MAP_threshold_out[,c("Subject_ID","Range 76-104")])

prediction.df.f<-prediction.df%>%
  select(mean.map, mean.HR, length_surg, days_post_surgery_dish, Age,
         AIS_ad,AIS_is_improved, Injury_level, 
        `Range 76-117`,`Range 76-104`)%>%
  mutate(AIS_is_improved_c = ifelse(AIS_is_improved==0,"NO","YES"))%>%
  rename(range_76_104 = `Range 76-104`, range_76_117 = `Range 76-117`)%>%
  filter(complete.cases(.))
```

### Search for best model

Search for the best model as for AIC. Single terms  and pairwise interaction are allowed.

```{r}
model_search<-glmulti(AIS_is_improved~poly(mean.map,2)+mean.HR+length_surg+ days_post_surgery_dish+Age+AIS_ad+Injury_level+range_76_104+range_76_117, 
                 data = prediction.df.f,family="binomial", level = 1, crit = "aicc", plotty = F, report = F)

## Extract the model that minimizes AIC
best_model<-model_search@objects[[1]]
```

### Best model

```{r}
best_model
```

### Prediction models

Next we build a prediction model with the terms from the best model in the previous search.

```{r}
## Best model
mod_fit_best <- train(AIS_is_improved_c~poly(mean.map,2)+AIS_ad, 
                 data = prediction.df.f, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

#Prediction probability
prob_predict_improve_best<-mod_fit_best$pred[,4]

model_prediction_best<-prediction_fun(prob_predict_improve_best,
                        prediction.df.f$AIS_is_improved,ref = "1")



##No MAP
mod_fit_best_noMap <- train(AIS_is_improved_c~AIS_ad, 
                 data = prediction.df.f, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

#Prediction probability
prob_predict_improve_best_noMAP<-mod_fit_best_noMap$pred[,4]

model_prediction_best_noMAP<-prediction_fun(prob_predict_improve_best_noMAP,
                        prediction.df.f$AIS_is_improved,ref = "1")

model_prediction_best$conf_mat
model_prediction_best_noMAP$conf_mat

coef(mod_fit_best$finalModel)
```

## Predicting AIS discharge

Process the data for prediction model
```{r }
prediction.df.f_dis<-prediction.df%>%
  select(mean.map, mean.HR, length_surg, days_post_surgery_dish, Age,
         AIS_ad,AIS_is_improved, Injury_level, AIS_dis,
        `Range 76-117`,`Range 76-104`)%>%
  mutate(A_dis = ifelse(AIS_dis == "A", 1, 0),
         D_dis = ifelse(AIS_dis == "D", 1, 0),
         A_dis_f = ifelse(AIS_dis == "A", "A", "nonA"),
         D_dis_f = ifelse(AIS_dis == "D", "D", "nonD"))%>%
  rename(range_76_104 = `Range 76-104`, range_76_117 = `Range 76-117`)%>%
  filter(complete.cases(.))
```

### Predicting A discharge

#### Search for best model

Search for the best model as for AIC.

```{r}
model_search_disA<-glmulti(A_dis~poly(mean.map,2)+mean.HR+length_surg+ days_post_surgery_dish+Age+AIS_ad+Injury_level+range_76_104+range_76_117, 
                 data = prediction.df.f_dis,family="binomial", level = 1,
                 crit = "aicc", plotty = F, report = F)

## Extract the model that minimizes AIC
best_model_disA<-model_search_disA@objects[[1]]
```

#### Best model

```{r}
best_model_disA
```

#### Prediction models

Next we build a prediction model with the terms from the best model in the previous search. In addition, we augment that model with adding the time ranges.

```{r}
## Best model
mod_fit_disA_best <- train(A_dis_f~poly(mean.map,2)+AIS_ad+Injury_level+range_76_117, 
                 data = prediction.df.f_dis, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

#Prediction probability
prob_predict_disA_best<-mod_fit_disA_best$pred[,3]

model_prediction_disA_best<-prediction_fun(prob_predict_disA_best,
                        prediction.df.f_dis$A_dis,ref = "1")

## noMap
mod_fit_disA_best_noMap <- train(A_dis_f~AIS_ad+Injury_level, 
                 data = prediction.df.f_dis, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

#Prediction probability
prob_predict_disA_best_noMap<-mod_fit_disA_best_noMap$pred[,3]

model_prediction_disA_best_noMap<-prediction_fun(prob_predict_disA_best_noMap,
                        prediction.df.f_dis$A_dis,ref = "1")

model_prediction_disA_best$conf_mat
model_prediction_disA_best_noMap$conf_mat

coef(mod_fit_disA_best$finalModel)
```

### Predicting D discharge

#### Search for best model

Search for the best model as for AIC. Single terms  and pairwise interaction are allowed.

```{r}
model_search_disD<-glmulti(D_dis~poly(mean.map,2)+mean.HR+length_surg+ days_post_surgery_dish+Age+AIS_ad+Injury_level+range_76_104+range_76_117, 
                 data = prediction.df.f_dis,family="binomial", level = 1,
                 crit = "aicc", plotty = F, report = F)

## Extract the model that minimizes AIC
best_model_disD<-model_search_disD@objects[[1]]
```

#### Best model

```{r}
best_model_disD
```

#### Prediction models

Next we build a prediction model with the terms from the best model in the previous search.

```{r}
## Best model
mod_fit_disD_best <- train(D_dis_f~AIS_ad+length_surg+Age, 
                 data = prediction.df.f_dis, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

#Prediction probability
prob_predict_disD_best<-mod_fit_disD_best$pred[,3]

model_prediction_disD_best<-prediction_fun(prob_predict_disD_best,
                        prediction.df.f_dis$D_dis,ref = "1")

## add MAP
mod_fit_disD_best_addMap <- train(D_dis_f~poly(mean.map, 2)+AIS_ad+length_surg+Age, 
                 data = prediction.df.f_dis, method="glm", family="binomial",
                 trControl = ctrl,na.action = 'na.omit')

#Prediction probability
prob_predict_disD_best_addMap<-mod_fit_disD_best_addMap$pred[,3]

model_prediction_disD_best_addMap<-prediction_fun(prob_predict_disD_best_addMap,
                        prediction.df.f_dis$D_dis,ref = "1")

model_prediction_disD_best$conf_mat
model_prediction_disD_best_addMap$conf_mat

coef(mod_fit_disD_best$finalModel)
```

## Figure 5

```{r fig.width=8, fig.height=8}
sense_spec_plot_improve<-model_prediction_best$sense_spec_plot+
  ggtitle("AIS impro.")

sense_spec_plot_disA<-model_prediction_disA_best$sense_spec_plot+
  ggtitle("A at disch.")

sense_spec_plot_disD<-model_prediction_disD_best$sense_spec_plot+
  ggtitle("D at disch.")

# ROC plot
auc_improve<-model_prediction_best$auc
auc_disA<-model_prediction_disA_best$auc
auc_disD<-model_prediction_disD_best$auc

roc_df_improve<-model_prediction_best$roc_df
roc_df_disA<-model_prediction_disA_best$roc_df
roc_df_disD<-model_prediction_disD_best$roc_df

roc_plot<-ggplot(roc_df_improve, aes(1-spec_y, sens_y))+
  geom_line(data=roc_df_disA, aes(1-spec_y, sens_y), color="green2")+
  geom_line(data=roc_df_disD, aes(1-spec_y, sens_y), color="steelblue1")+
  geom_line(color="firebrick")+
  theme_minimal()+
  scale_x_continuous(breaks=seq(0,1,0.25), 
                     labels = c("1", "0.75","0.5","0.25","0"))+
  xlab("Specificity")+
  ylab("Sensitivity")+
  geom_segment(aes(x=0, xend=1, y=0, yend=1), linetype="dashed")+
  
  annotate("text", label=paste("AIS impro. =",round(auc_improve, 3)),x=0.66, y=0.16)+
  geom_segment(aes(x=0.3, xend=0.36, y=0.15, yend=0.15), color="firebrick")+
  
  annotate("text", label=paste("A at disch. =",round(auc_disA, 3)),
           x=0.66, y=0.1)+
  geom_segment(aes(x=0.3, xend=0.36, y=0.1, yend=0.1), color="green2")+
  
  annotate("text", label=paste("D at disch. =",round(auc_disD, 3)),
           x=0.66, y=0.04)+
  geom_segment(aes(x=0.3, xend=0.36, y=0.03, yend=0.03), color="steelblue1")+
  #Annotate points of best model
  annotate("point", x=1-roc_df_improve$spec_y[which.min(roc_df_improve$diff_y)],
           y=roc_df_improve$sens_y[which.min(roc_df_improve$diff_y)], color="firebrick",
           size=4, alpha=0.3)+
  annotate("point", x=1-roc_df_disA$spec_y[which.min(roc_df_disA$diff_y)],
           y=roc_df_disA$sens_y[which.min(roc_df_disA$diff_y)],color="green2",
           size=4, alpha=0.3)+
  annotate("point", x=1-roc_df_disD$spec_y[which.min(roc_df_disD$diff_y)],
           y=roc_df_disD$sens_y[which.min(roc_df_disD$diff_y)],color="steelblue1",
           size=4, alpha=0.3)+
  annotate("point", x=0.05, y=0.94,color="grey", size=4, alpha=0.5)+
  annotate("text", label = "Threshold", x=0.22, y=0.95)

## Compose the model prediction figure
figure5<-(sense_spec_plot_improve+sense_spec_plot_disA+sense_spec_plot_disD+
            roc_plot+ggtitle("ROC"))+
  plot_annotation(tag_levels = 'a')+plot_layout(ncol=2)

ggsave(figure5, file="figure 5.png", width = 8, height = 8)

figure5
```

# Session info

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