---
title: "Approaches to analyze experimental *R. solanacearum* infections:\n*Supplementary Material 4:* Datasets from the Prior Lab\n Survival Data"
author: "Niklas Schandry"
output: 
  html_document:
    toc: true
    toc_float: true
    theme: yeti
    highlight: kate
---
# General information on this file
This document is part of the supplementary material of "A practical guide to descriptive and statistical analysis of R. solanacearum infection data using R". This file explores more complex survival analysis, that contain multiple levels of Plants (genotype) and Strains. This unpublished work by Santatra Ravelomanantsoa and Philippe Prior, both with CIRAD / INRA UMR PVBMT (Réunion Island), has been pseudonymized.

# Data Import

## Specifying the data

Initially, data needs to be read into R. 

```{r eval = T}
###Name of the file to be read
table <- c("S4_Ravelomanantsoa_Prior.csv") 
````


## Reading data and formatting
````{r, eval=T}
disease_index <- as.data.frame(read.table(table, header=T,
                                          sep=";" , ###Sets the seperator of the csv file
                                          dec=","), ###Sets the decimal operatore of the csv file
                                          stringsAsFactors=T) ###
````

Using str() one can see if the table was properly imported.

````{r}
str(disease_index)
````

This table is a survival table already, go directly into surv_table.

````{r, message=F}
library("tidyr")
```
````{r}
surv_table<- gather(disease_index, key=DPI, 
                     value=Diseased,
                     c(X1,X8,X11,X15,X18,X22,X25,X32, X43),
                     na.rm=T)
str(surv_table)
###Do stuff to surv_table that makes it more useable
surv_table$DPI <- as.numeric(
  na.omit(
    as.numeric(
      unlist(
        strsplit(
          as.character( 
            surv_table$DPI
          ),
          "X" )
      )
    )
  )
)
surv_table <- surv_table[(surv_table$Diseased %in% c(0,1)),] ##Drop Disease recordings that are not binary
surv_table$Strain <- as.factor(surv_table$Strain)
surv_table$Plant <- as.factor(surv_table$Plant)
surv_table$ROTOPLAN <- as.factor(surv_table$ROTOPLAN)
surv_table$StrainPlant <- interaction(surv_table$Strain,surv_table$Plant)
surv_table$StrainRoto <- interaction(surv_table$Strain, surv_table$ROTOPLAN)
surv_table$PlantRoto <- interaction(surv_table$Plant, surv_table$ROTOPLAN)
contrasts(surv_table$Strain) <- "contr.treatment"
contrasts(surv_table$Plant) <- "contr.treatment"
contrasts(surv_table$ROTOPLAN) <- "contr.sum"
contrasts(surv_table$StrainPlant) <- "contr.treatment"

```

In this dataset, Strain is the bacterial strain, Plant is the plant genotype that was infected, and ROTOPLAN, is the rotoplan (a kind of greenhouse) that this infection was carried out in. In other words, ROTOPLAN is an indicator that descriminates replicates. 

Above, i am generating interaction columns. In the R notation interactions are denoted using \* in a formula. Strain\*Plant in a formula, is equivalent to interaction(Strain,Plant). For me, predifining the interaction terms, makes it easier to tangle them apart for visualization.I also find it easier to interpret.

# Survival Analysis

Initially, one may enjoy generating Kaplan-Meier estimates of the survival data. This is done using the survfit function. Here, the rotoplan is not included as a predictor.

````{r}
library("survival")
surv_fit <- survfit(Surv(DPI, as.numeric(Diseased), type="right") ~Strain + Plant, data=surv_table)

###Make Dataframe for plotting
strata_dummy <-NULL
for(i in 1:length(surv_fit$strata)){
      # add vector for one strata according to number of rows of strata
      strata_dummy <- c(strata_dummy, rep(names(surv_fit$strata)[i], surv_fit$strata[i]))
}
###Data frame generation inspired by a post by Hadley Wickham to the ggplot2 Googlegroup
surv_fit.df <- data.frame( 
  time = surv_fit$time, 
  n.risk = surv_fit$n.risk, 
  n.event = surv_fit$n.event, 
  surv = surv_fit$surv, 
  strata = strata_dummy, 
  upper = surv_fit$upper, 
  lower = surv_fit$lower 
) 
zeros <- data.frame(time = 0, surv = 1, strata = names((surv_fit$strata)), 
  upper = 1, lower = 1)

surv_fit.df <- plyr::rbind.fill(zeros, surv_fit.df) ###I dont want to load plyr because i guess it will interfere with dplyr...
rm(strata_dummy)
rm(zeros)
stratas <- matrix( nrow=length(surv_fit.df$strata),ncol=2, unlist(strsplit(as.character(surv_fit.df$strata),", ")), byrow=T )
surv_fit.df$Plant <- as.factor(matrix(nrow=length(surv_fit.df$strata), ncol=2,unlist(strsplit(stratas[,2],"=")), byrow=T )[,2])
surv_fit.df$Strain <- as.factor(matrix(nrow=length(surv_fit.df$strata), ncol=2,unlist(strsplit(stratas[,1],"=")), byrow=T )[,2])

###End of data frame generation
library("ggplot2")
###Start plotting
ggplot(surv_fit.df,aes(time, surv, colour = Strain)) + 
  facet_grid(Strain~Plant) +
  stat_summary(aes(y = surv), fun.y=mean, geom="step") +
  ggtitle("All KM survival estimates")
```

## Pairwise logrank testing

If one is interested between all pairwise differences of ALL StrainPlant combinations (6216 pairs), run the below. This will perform a logrank test for every possible combination of strain and plant and test against all others. This is stratified by ROTOPLAN.
This is not evaluated in the rmarkdown file because the output is huge and not particularily helpful. 
The below does not correct for multiple comparisons. 

```{r eval=F}
###Pairwise--- http://stackoverflow.com/questions/17338774/r-formula-how-to-constrain-calculations-to-two-groups-using-formula
library("gtools")

levels<-(unique(surv_table$StrainPlant))

groups<-combinations(length(levels),2,levels)

pairWise_StrainPlant <- plyr::alply(groups, 1, 
    function(pair) {
        survdiff(Surv(DPI, as.numeric(Diseased), type="right")~StrainPlant  + frailty(ROTOPLAN), surv_table , subset=surv_table$StrainPlant %in% c(pair[[1]], pair[[2]]) )
    })
names(pairWise_StrainPlant) <- plyr::alply(groups, 1, 
    function(pair) {
        c(pair[[1]], pair[[2]]) 
    })
```

## Hazards

```{r}
library("survcomp")
###Cox-Proportional hazards####
#Build model
srv_coxph <- coxph(Surv(DPI, as.numeric(Diseased),type="right") ~Strain+Plant+ROTOPLAN, data=surv_table) 
###Check porportionality of hazards
cox.zph(srv_coxph, transform = "log")
###Hazard ration analysis, not appropriate if cox.zph above returns a low pvalue
hazrat <- hazard.ratio(surv_table$Strain, surv_table$DPI, surv.event = as.numeric(surv_table$Diseased))
````

The proportional hazards assumption is violated.

## Survival Regression

As the hazards are non-proportional, and generating a ton of pairwise tests as done above without correcting for multiple comparisons is a really bad idea, survival regression could be handy.

```{r message = F}
library("rms")
```
```{r}

psm_gaus <- psm(Surv(DPI, as.numeric(Diseased), type="right") ~ Strain * Plant + ROTOPLAN, surv_table, dist = "gaussian")
psm_logistic <- psm(Surv(DPI, as.numeric(Diseased), type="right") ~ Strain * Plant + ROTOPLAN, surv_table, dist = "logistic")
psm_lnorm <- psm(Surv(DPI, as.numeric(Diseased), type="right") ~ Strain * Plant + ROTOPLAN, surv_table, dist = "lognormal")
psm_wei <- psm(Surv(DPI, as.numeric(Diseased), type="right") ~ Strain * Plant + ROTOPLAN, surv_table, dist = "weibull")

aic.scores.psm <- rbind(
  extractAIC(psm_wei),
  extractAIC(psm_gaus),
  extractAIC(psm_logistic),
  extractAIC(psm_lnorm))
###Make useable AIC table
rownames(aic.scores.psm) <- c("Weibull", "Gaussian", "Logist", "Lognorm")
colnames(aic.scores.psm) <- c("df", "AIC")
aic.scores.psm
```

Gaussian is the best fit. Explore the gaussian psm

```{r}
library("multcomp")
summary(glht(psm_gaus, linfct=mcp(Strain="Tukey")))
summary(glht(psm_gaus, linfct=mcp(Plant="Tukey")))
```

## By plant genotype analysis 

Maybe, it would be interesting to only know which strain performed how within each Plant genotype used.

```{r}
options(warn = -1)
plants_array <- array(levels(surv_table$Plant))
plant_within <- sapply(plants_array, 
                    function(plant) {
                      psm(Surv(DPI, as.numeric(Diseased), type="right") ~ Strain + ROTOPLAN, data=surv_table[surv_table$Plant==plant,], dist="gaussian")
                    }, simplify=F)
names(plant_within) <- levels(surv_table$Plant)

plant_within_cld <- sapply(plants_array, 
                    function(plant) {
                      cld( glht( 
                        psm(Surv(DPI, as.numeric(Diseased), type="right") ~ Strain + ROTOPLAN, data=surv_table[surv_table$Plant==plant,], dist="gaussian"),
                        linfct=mcp(Strain="Tukey")
                      )
                      )
                    }
                    , simplify=F)
names(plant_within_cld) <- levels(surv_table$Plant)

plant_within_cld
```
## By strain analysis
```{r}

strain_array <- array(levels(surv_table$Strain))
strain_within <- sapply(strain_array, 
                    function(strain) {
                      psm(Surv(DPI, as.numeric(Diseased), type="right") ~ Plant + ROTOPLAN, data=surv_table[surv_table$Strain==strain,], dist="gaussian")
                    }, simplify=F)
names(strain_within) <- levels(surv_table$Strain)

###Strain4 / Rs14, can not be fit. This is probably because it is flatlining (not causing symptoms) for some plants. This strain is excluded.
###I would consider this a relevant finding, but it is also quite obvious from all analysis in this file that this strain is different from the others

strain_within_cld <- sapply(strain_array[c(1:5,7:length(strain_array))], 
                    function(strain) {
                      cld( glht( 
                        psm(Surv(DPI, as.numeric(Diseased), type="right") ~ Plant + ROTOPLAN, data=surv_table[surv_table$Strain==strain,], dist="gaussian"),
                        linfct=mcp(Plant="Tukey")
                      )
                      )
                    }
                    , simplify=F)
names(strain_within_cld) <- levels(surv_table$Strain)[ c(1:5,7:length(strain_array)) ]

strain_within_cld
contrasts(surv_table$Strain)
```

Finally, check for a rotoplan effect across the full dataset, using the same distribution as earlier

## Rotoplan effect?
```{r}
cld(glht(
  survreg(Surv(DPI, as.numeric(Diseased), type="right") ~ Strain + ROTOPLAN, data=surv_table, dist="gaussian"),
                        linfct=mcp(ROTOPLAN="Tukey")
))
```
 
 As can be seen, ROTOPLAN has a signficant effect. Rotoplans 2 and 3 are in the same signficance group, meaning they are not significantly different. The same is true for 1 and 4. One could explore this further by subsetting the dataset into two groups, one made of ROTOPLANs (1,4) and one of (2,3) and analyze those individually.
 
## Plotting survival regression 
 
```{r generate curves of survreg}
library("stringr")
library("dplyr")
###The StrainPlant variable was generated earlier.
###The below assumes interactions between Strain and Plant, meaning we could write
### Surv(DPI, as.numeric(Diseased), type="right") ~ Strain * Plant
###and obtain the same results. I prefer the way it is currently done, because i can split it back into the original things more conveniently.
s_reg_gaus <- survreg(Surv(DPI, as.numeric(Diseased), type="right") ~ StrainPlant, surv_table, dist="gaussian")
s_reg_logistic <- survreg(Surv(DPI, as.numeric(Diseased), type="right") ~ StrainPlant, surv_table, dist="logistic")
s_reg_lnorm <- survreg(Surv(DPI, as.numeric(Diseased), type="right") ~ StrainPlant, surv_table, dist="lognormal")
s_reg_wei <- survreg(Surv(DPI, as.numeric(Diseased), type="right") ~ StrainPlant, surv_table, dist="weibull")

###Step 1, extract the coefficients. These are relative to StrainPlant combination RS1.A because StrainPlant is treatment contrasted.
surv_fit.df$StrainPlant <- interaction(surv_fit.df$Strain,surv_fit.df$Plant)
for (i in 1:nlevels(surv_fit.df$StrainPlant)) { #For loop through strains
  if(i==1) { #Strain1 is relative to itself, so no change
  coef_wei <- list()
  coef_logistic <- list()
  coef_gaus <- list()
  coef_lnorm <- list()
  coef_wei[i] <- coef(s_reg_wei)[i]
  coef_logistic[i] <- coef(s_reg_logistic)[i]
  coef_gaus[i] <- coef(s_reg_gaus)[i]
  coef_lnorm[i] <- coef(s_reg_lnorm)[i]
  } else { ###Other strains are relative to 1
  coef_wei[i] <- coef(s_reg_wei)[1] + coef(s_reg_wei)[i]
  coef_logistic[i] <- coef(s_reg_logistic)[1] + coef(s_reg_logistic)[i]
  coef_gaus[i] <- coef(s_reg_gaus)[1] + coef(s_reg_gaus)[i]
  coef_lnorm[i] <- coef(s_reg_lnorm)[1] + coef(s_reg_lnorm)[i]
  }
}
##Step 2
####Store the coefficients and the scale in a new data frame, of parameters
### Keep in mind that survreg.distributions$weibull is different from rweibull, hence the difference in names.
sregparams <- data.frame(
  StrainPlant = rep(levels(surv_table$StrainPlant),4 ), #Fill with strains
  scale.wei = exp(unlist(coef_wei)), #weibull fit scale parameters
  scale.logistic = rep(s_reg_logistic$scale, nlevels(surv_table$StrainPlant)), #fill with logis scales
  scale.gaus = rep(s_reg_gaus$scale, nlevels(surv_table$StrainPlant)), #fill with gaus scales
  scale.lnorm = rep(s_reg_lnorm$scale, nlevels(surv_table$StrainPlant)), #fill with lnorm scale
  shape.wei =  rep(1/s_reg_wei$scale, nlevels(surv_table$StrainPlant)), #shape for weibull
  shape.logistic = unlist(coef_logistic), #shape for logistic
  shape.gaus =  unlist(coef_gaus), #shape for gaus
  shape.lnorm =   unlist(coef_lnorm) #shape for lnorm
  )
##Step 3
###Calculate the "daily" value of each curve
for (i in 1:nlevels(surv_fit.df$StrainPlant)){
  if(i==1) {
    wei <- list()
    logis <- list()
    gaus <- list()
    lnorm <- list()
  }
  x <- levels(surv_fit.df$StrainPlant)[i]
  n <- c(1:max(surv_table$DPI))
  data <- dplyr::filter(sregparams, StrainPlant==x) ###Watch out here. It is crucial to use dplyr::filter and not base::filter
  time <- n
  wei <- cbind(wei, pweibull(
    q=n,
    scale=data$scale.wei,
    shape=data$shape.wei,
    lower.tail=FALSE))
  logis <- cbind(logis,plogis(
    q=n,
    scale=data$scale.logistic,
    location=data$shape.logistic,
    lower.tail=FALSE  ))
  gaus <- cbind(gaus,pnorm(
   q=n,
   sd=data$scale.gaus,
   mean=data$shape.gaus,
   lower.tail = F))
  lnorm <- cbind(lnorm, plnorm(
    q=n,
    sd=data$scale.lnorm,
    mean=data$shape.lnorm,
    lower.tail=F))
}

##Step 4
###Put all the curves into a data.frame that contains information on "time" and also "Strain", for compatibility with other data.frames
sreg_curves <- data.frame(
  wei.sreg = cbind(unlist(wei)),
  logis.sreg = cbind(unlist(logis)),
  gaus.sreg = cbind(unlist(gaus)),
  lnorm.sreg = cbind(unlist(lnorm)),
  StrainPlant = rep(unlist(levels(surv_fit.df$StrainPlant)),each=max(surv_table$DPI)),
  time = rep(c(1:max(surv_table$DPI)), nlevels(surv_fit.df$StrainPlant))
)
sreg_curves$Strain <- str_split_fixed(sreg_curves$StrainPlant, "\\.",2)[,1]
sreg_curves$Plant <- str_split_fixed(sreg_curves$StrainPlant, "\\.",2)[,2]
##Step 5
###Turn that data.frame into a long data.frame (not used here but for other figures.)
sreg_long <- sreg_curves %>% gather(., key="Distribution",values = c(lnorm.sreg, wei.sreg,gaus.sreg,logis.sreg) )
sreg_long$Distribution <- as.factor(sreg_long$Distribution)
##Levels: gaus.sreg wei.sreg
levels(sreg_long$Distribution) <- c("Gaussian","Lognormal","Loglogistic","Weibull")

sreg_long$Strain <- str_split_fixed(sreg_long$StrainPlant, "\\.",2)[,1]
sreg_long$Plant <- str_split_fixed(sreg_long$StrainPlant, "\\.",2)[,2]
```
Now, these can be plotted and inspected visually..
```{r plot KM and curve}
###Plot of KM+Weibull
  ggplot(surv_fit.df, aes(time, surv, colour = Strain)) +
  geom_step() +
  geom_line(data=sreg_curves,aes(y=wei.sreg),color="black") +
  facet_grid(Plant~Strain)+ 
  ggtitle("Kaplan-Meier estimates and fit to\nWeibull distribution")
###Plot of KM+Gaussian  
  ggplot(surv_fit.df, aes(time, surv, colour = Strain)) +
    geom_step() +
    geom_line(data=sreg_curves,aes(y=gaus.sreg),color="black") +
    facet_grid(Plant~Strain) + 
    ggtitle("Kaplan-Meier estimates and fit to\nGaussian distribution")
##Logis
    ggplot(surv_fit.df, aes(time, surv, colour = Strain)) +
    geom_step() +
    geom_line(data=sreg_curves,aes(y=logis.sreg),color="black") +
    facet_grid(Plant~Strain) + 
    ggtitle("Kaplan-Meier estimates and fit to\nLogistic distribution")
  ##Lnorm
    ggplot(surv_fit.df, aes(time, surv, colour = Strain)) +
    geom_step() +
    geom_line(data=sreg_curves,aes(y=lnorm.sreg),color="black") +
    facet_grid(Plant~Strain) + 
    ggtitle("Kaplan-Meier estimates and fit to\nLognormal distribution")
```

#Session Info
```{r session info}
sessionInfo()
```