---
title: "Approaches to analyze experimental *R. solanacearum* infections:\n*Supplementary Material 3:* Datasets from the Allen Lab\n Part III: Analysis of Lowe-Power_mBio_2016" 
author: "Niklas Schandry"
output: 
  html_document:
    toc: true
    toc_float: true
    theme: yeti
    highlight: kate
---
# General information on this document
This document is part of the supplementary material of "A practical guide to descriptive and statistical analysis of R. solanacearum infection data using R". In this document, the analysis of the Lowe_mBio_2016 dataset is performed, based on the formatted file generated in section I of material S3.

# Data Import

## Specifying the data
The below is retained so packages can be installed if necessary.
```{r workingdirectory, eval=F}
###Install all required packages:
#install.packages(c("MESS","lme4","lmerTest","multcomp","survival","rms","coxme","stargazer","survcomp","tidyverse","rcompanion"))


###Define Working Directory and set it
###Note for the Markdown version: R-Markdown cannot set the working directory
###R markdown will always use the directory the .Rmd file is located in
###In the .Rmd file this code section is not actually evaluated and only serves illustratory purpose.
wd <- c("~/My_Data/DataDirectory/")
setwd(wd)
````

Import the mBio2016 long table. This was generated in part I.

```{r input definition}
###Name of the file to be read
di_long <- read.table("S3_mBio2016_long.csv", sep=",", header=T)
di_long$DPI <- di_long$DPI
###Redifine the variables
di_long$Batch <- as.factor(di_long$Batch)
contrasts(di_long$Strain) <- "contr.Treatment"###First alphabetical strain will be the baseline!ant) <- "contr.Treatment"###First alphabetical plant will be the baseline!
contrasts(di_long$Batch) <- "contr.sum" ###Batches will be averaged to generate the baseline!
````


# Analysis of AUDPC

## Rationale and data preparation

Everything from here on out is very similar to the analysis performed on sim_data.csv.
Changes are explained.
*Change* in below code: Removed "Plant" from group_by, since Plant here has a different meaning than in sim_data.csv (here it is a batch specific subject, in the other dataset it is a plant genotype)
```{r load dplyr, message=F}
library("dplyr")
```
```{r summarize}
di_summary <- di_long %>% group_by(Strain, Batch, DPI) %>% 
  summarise(mean(DI),sd(DI),sd(DI)/sqrt(length(DI))) ##calculate within Batch mean, sd, and se for each Plant/Strain combination.
##Averages within each replicate. This summary table is mainly helpful for plotting, not used for analysis...
colnames(di_summary) <- c("Strain", "Batch", "DPI", "mean", "sd", "se") ###Assign correct columnnames
```

Using these summaries, one can take a look at the averaged disease progression, using

```{r load ggplot, message=F}
library("ggplot2")
```
```{r plot disease area from summary}
ggplot(filter(di_summary, Batch==di_summary$Batch[[1]])) + ###Use Only Batch A (too busy otherwise)
  aes(x=DPI,y=mean,color=Strain) + ###Color by Strain, specify x and y.
  geom_area(aes(fill=Strain),position="identity",alpha=0.15) + ###Area plot, colored by Strain
  ##geom_errorbar(aes(x=DPI, ymax=mean+se, ymin=mean-se), size=0.25)+ ###This line adds SE. 
  ##geom_errorbar(aes(x=DPI, ymax=mean+sd, ymin=mean-sd), size=0.25)+ ###This line adds SD. Don't use both.
  facet_wrap(~Strain) + ###One plot per strain 
  labs(x = "Days post infection", y = "Avg. Disease Index") + #Labels
  ggtitle("Figure 1\nDisease Areas,\nper strain, for batch A") #Title
```  
  
  Here, from a plot of the Areas, it seems these strains are not obviously different regarding the disease.
  
*Change*: Modified subject selection, as here the subjects are based on interaction, hence have character names. Each level of this corresponds to a number, so this can be used for iteration in the loop below.
```{r load MESS, message=F}
library("MESS")
```
```{r make auc_df}
####Build a table of AUDCPs, per subject
auc_df <- data.frame() ###Make auc_df data frame
for (i in 1:max(as.numeric(di_long$subject))) { ##Go by subject
  temp <- di_long[as.numeric(di_long$subject)==i,] ###Subset full table into the subject table
  temp <- droplevels(temp) ###Drop levels, so levels works properly below
  auc_df[i,1] <- i ###Subject number
  auc_df[i,2] <- levels(temp$Strain) ###Strain
  auc_df[i,3] <- levels(temp$Plant) ###Plant
  auc_df[i,4] <- levels(temp$Batch) ###Batch
  auc_df[i,5] <- auc(temp$DPI,temp$DI) ###AUC; i assume that trapezoid rule is fine here. 
  ### Additionally, auc calculation starts with the lowest x (DPI). I think this is sensible
  ### I assume that if one specifies "from=0", the curve is expanded by a triangle that covers the range from 
  ### 0 to whatever is the value at the first observation. I think the first observation should ideally be 0
  ### if data was recorded from the beginning..
}
colnames(auc_df) <- c("subject","Strain","Plant","Batch","AUC") ###Name columns in AUC datafarame
auc_df$Strain <- as.factor(auc_df$Strain) #refactor
auc_df$Plant <- as.factor(auc_df$Plant) #refactor
auc_df$Batch <- as.factor(auc_df$Batch) #refactor
str(auc_df)
```

The auc data.frame contains one area under the curve per subject and all other subject specific variables as stored in the original table.

##Analysis of differences in area under the disease progression curve
An initial assessment of strain or batch specific differences in AUDPC can be performed visually, for example by generating boxpots.
```{r boxplot of auc}
ggplot(auc_df) + geom_boxplot(
  aes(x=Strain, y=AUC, color=Strain),  #Plot boxplots of AUCs, by strains
  notch=F) +
  labs(y="AUDPC") +
  facet_wrap(~Batch) + ###Individual plots per batch (and plant if applicable)
  ggtitle("Figure 2\nArea Under the Disease Progression Curve per strain across batches")
```
  
  Here, it can be seen that comp_nagOP performed poorly in replicated 3.  
  
  Next, one can use the area under the disease progression curve, to build a linear model, or a linear mixed effects model.
```{r load lm4 and lmerTest, message=F}
library("lme4")
library("lmerTest")
```
```{r analysis AUC}
summary(lm(AUC~Batch,data=auc_df)) ###Can be used to identify batch effects. If there are none, including batch as a random factor below is not necessary (but also not necessarily wrong).

auc_lmer <- lmer(
  AUC ~ Strain + (1|Batch),   ### AUC modeled as a function of strain, random effects of batch 
    data=auc_df  ) #Linear mixed effects model.

auc_lm <- lm(AUC~Strain+Batch,data=auc_df) #Linear model.

AIC(auc_lm,auc_lmer) #Lower AIC, better fit, linear model is slightly better.
ggplot(data=auc_df, aes(y=AUC, x=Strain)) +geom_boxplot(aes(colour=Batch)) + ggtitle("Boxplot of AUDCP per strain by experimental batch")
```
  
  The model can be explored using various functions, such as summary.
  
```{r summary auc lmer}
library("broom")
summary(auc_lmer)  ### A model summary, containing information on the model.
tidy(auc_lmer) ### A cleaner display using tidy.

###The tidy output explained:
#Term: A discription: (Intercept) overall intercept. Intercept depends on the contrasts set initially. Here treatment contrasts are used, so Intercept = First alphabetical strain (Strain1).
#StrainStrain2: Difference in the estimate (slope), between Strain2 and the (Intercept)

##Estimate: The estimated slope

```
  
For example, we see the estimated slopes (Estimate) and standard errors, together with a t- and corresponding p-value in the output of the summary function. Note, that the above only contains information on differences between different levels and the "baseline"", which is called (Intercept).  The baseline is determined by the contrast settings that were specified earlier.

But, it may be quite relevant to know how individual strains compare to each other (and not just how each strain compares to Strain1). This can be analyzed using a generalized linear hypothesis test, while ajusting for multiple comparisons using Tukey's method.

```{r load multcomp, message=F}
library("multcomp")
library("rcompanion")
library("stringr")
```
```{r glht auc}
tidy(summary(glht(auc_lmer, linfct=mcp(Strain="Tukey"))))
```
  
 This information can, for example, be integrated into a boxplot of the individual disease areas. For example, using a compact letter display, in conbination with AUDPC, combines statistical and visual information.
First, however, the grouping letters need to be calculated. Then, these letters are added to the boxplot generated earlier.

```{r AUC model CLD plot}
auc_cld <- cld(glht(auc_lmer, linfct=mcp(Strain="Tukey"))) ###Save letters
auc_cld <- cbind(levels(auc_df$Strain),auc_cld$mcletters$Letters) ###bind letters to columns
colnames(auc_cld) <- c("Strain","Letter") ###Name columns
auc_cld <- as.data.frame(auc_cld) #Coerce to dataframe
###Integrate letters into auc_df##
###Some extra scripting to make the mean and CI plot.
auc_CI <- as.data.frame(tidy(confint(auc_lmer)))
auc_CI <- auc_CI[3:nrow(auc_CI),] ##Drop sig01, sigma
auc_CI$Strain <- levels(auc_df$Strain)
###Mean relative to Strain1 (except strain1 that one is absolute)

for (i in 1:nrow(auc_CI)){
  if (i==1){
    auc_CI$mean[i] <- mean(c(auc_CI$X2.5..[i],auc_CI$X97.5..[i]))
    auc_CI$upr[i] <- auc_CI$X97.5..[i]
    auc_CI$lwr[i] <- auc_CI$X2.5..[i]
  } else {
    auc_CI$mean[i] <- c(auc_CI$mean[1]+mean(c(auc_CI$X2.5..[i],auc_CI$X97.5..[i])))
    auc_CI$upr[i] <- c(auc_CI$mean[1]+auc_CI$X97.5..[i])
    auc_CI$lwr[i] <- c(auc_CI$mean[1]+auc_CI$X2.5..[i])
  }
}
 
 
####Generate plot of meanCI of AUDCP with significance letters and raw data as jittered points
  ggplot(aes(x=Strain, y=AUC, color=Strain),data=auc_df) + ###Plot the auc_df
  geom_crossbar(data = auc_CI, aes(x = Strain, y = mean, ymin = lwr, ymax = upr,fill=Strain), alpha=0.3) +
  geom_jitter(aes(shape=Batch)) + ###with jitter overplotted, symbol shape defined by batch
  geom_text(aes(x=Strain, y=-3, label=Letter),color="black", data=auc_cld) + ###Get the letters from auc_cld 
  #and write those to position y=-3
  labs(y="AUDPC") + #Y-Axis label
  ggtitle("AUDPC raw values and mean from the LMM\nwith 95% CI per strain and grouping letters") #Title
```
  
  One can plot the pairwise difference in AUDPC means with confidence intervals for the linear model and the linear mixed effects model.
```{r pairwise confint AUDPC} 
pairwise_confint_AUDPC_lm <- as.data.frame(confint(glht(auc_lm, mcp(Strain = "Tukey")))$confint)
pairwise_confint_AUDPC_lm$Comparison <- rownames(pairwise_confint_AUDPC_lm)

pairwise_confint_AUDPC_lmer <- as.data.frame(confint(glht(auc_lmer, mcp(Strain = "Tukey")))$confint)
pairwise_confint_AUDPC_lmer$Comparison <- rownames(pairwise_confint_AUDPC_lmer)

###Plot the comparisons, below may not be the most straight-foward way to plot this the way I want it, but it works.
ggplot(pairwise_confint_AUDPC_lm, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr, color = Comparison)) + ###Plot Comparison on x, estimate on y
  scale_x_discrete(limits = rev(levels(as.factor(pairwise_confint_AUDPC_lm$Comparison)))) + ###Rescale x, so the order is inverted 
  geom_errorbar() + geom_point() + ###Draw data
  coord_flip() +  theme(legend.position="none") + xlab("") +###Invert X and Y, hide legend
  ggtitle("Difference in means of the AUDPC \nin the linear model with 95% confidence interval") ##Add a title

####Plot of the comparisons in the LMM. Here the bars are colored by their absolute value, instead of by their name.
ggplot(pairwise_confint_AUDPC_lmer, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr, color = abs(Estimate))) + ###Plot Comparison on x, estimate on y
  scale_x_discrete(limits = rev(levels(as.factor(pairwise_confint_AUDPC_lmer$Comparison)))) + ###Rescale x, so the order is inverted 
  geom_errorbar() + geom_point() + ###Draw data
  coord_flip() +  theme(legend.position="none") + xlab("") +###Invert X and Y, hide legend
  ggtitle("Difference in means of the AUDPC \nin the LMM with 95% confidence interval") ##Add a title
```
  
  Estimates of lm and lmm are very similar. As a rule of thumb, comparisons where "0" is not part of the 95% confidence interval are likely to produce a signficant p-value (assuming significance is denoted by p<0.05).

# Analysis of Disease development

## Rationale
In this section, the data is analyzed using a linear mixed effect model. While such models have been used already in the previous section, to test for a strain specific influence on the area under the disease progression curve, different data is used to build the model in this section.
As mentioned, AUDPC summarizes disease incidence and time into a single variable, the area.
However, in certain cases the AUDCP could be very similar, while the actual disease progression is different.
Taking a look at Figure 1, we see that Strain2 and Strain3 display different disease progressions. For Strain 2, the Disease Index continually rises over time. For Strain3, the the disease appears to set on very quickly, in fact before the observations started, but from there it does not increase in severity.
However, if we compare the impact of those two strains on the AUDCP (Fig2, 3), it becomes evident, that the AUDPCs for these strains are not statistically different, both strains belong to group "a". 
Other methods may be more sensitive to such differences. Below a different approach using linear mixed effect models is taken.

## Model
This model is the same as the model in the other files.
```{r lmer DI constrasts and model}
###Define contrasts for lmer
contrasts(di_long$Strain) <- "contr.treatment"###First alphabetical strain will be the baseline!
#contrasts(di_long$Plant) <- "contr.treatment"###First alphabetical plant will be the baseline!
contrasts(di_long$Batch) <- "contr.poly" ###Batches will be averaged
###Drop things that are not "Useful"
di_long.useful <- filter(di_long, Useful=="Yes")
str(di_long.useful)
#Fixed effect from strain, and a fixed interaction between Strain and time. Random effects individual and replicate.
disease_lmer <- lmer(DI ~ Strain + Strain:DPI + (1 | subject) + (1 | Batch), di_long.useful)
```
  
 The model can be investigated using summary functions. The pairwise comparisons may be plotted with 95%CI to assess how different two strains are.
 
```{r pairwise confint plot of disease lmer}
###Check model summary
summary(disease_lmer)
###E.g. plot confints
###Make pairwise confints and plot them, with a flipped coordinate system
pairwise_confint <- as.data.frame(confint(glht(disease_lmer, mcp(Strain = "Tukey", interaction_average=T)))$confint)
pairwise_confint$Comparison <- rownames(pairwise_confint)

ggplot(pairwise_confint, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr, color = abs(Estimate))) + ###Plot Comparison on x, estimate on y
  scale_x_discrete(limits = rev(levels(as.factor(pairwise_confint$Comparison)))) + ###Rescale x, so the order is inverted 
  geom_errorbar() + geom_point() + ###Draw data
  coord_flip() +  theme(legend.position="none") + xlab("") +###Invert X and Y, hide legend
  ggtitle("Difference in means with 95% confidence interval \ncolored by absolute estimated difference") ##Add a title

confint_model <- as.data.frame(tidy(confint(disease_lmer)))

confint_slopes <- confint_model[ (1 + 3 + nlevels( di_long$Strain ) ) : ( 3 + 2*nlevels( di_long$Strain ) ) , 2:3 ]

colnames(confint_slopes) <- c("lwr","upr")
confint_slopes$Estimate <- rowMeans(confint_slopes)
confint_slopes$Strain <- levels(di_long$Strain)

###Plot the estimates, below may not be the most straight-foward way to plot this the way I want it, but it works.
ggplot(confint_slopes, aes(x = Strain, y = Estimate, ymin = lwr, ymax = upr, color = abs(Estimate))) + ###Plot Comparison on x, estimate on y
  scale_x_discrete(limits = rev(levels(as.factor(confint_slopes$Strain)))) + ###Rescale x, so the order is inverted 
  geom_errorbar() + geom_point() + ###Draw data
  coord_flip() +  theme(legend.position="none") + xlab("") +###Invert X and Y, hide legend
  ggtitle("Absolute slopes with 95%CI") ##Add a title

#Intercepts are treatment contrasted
confint_icep <- confint_model[4:( 3 + nlevels( di_long$Strain )),]

confint_icep$Strain <- levels(di_long$Strain)

for (i in 1:nrow(confint_icep)){
  if (i==1){
    confint_icep$Estimate[i] <- mean(c(confint_icep$X2.5..[i],confint_icep$X97.5..[i]))
    confint_icep$upr[i] <- confint_icep$X97.5..[i]
    confint_icep$lwr[i] <- confint_icep$X2.5..[i]
  } else {
    confint_icep$Estimate[i] <- c(confint_icep$Estimate[1]+mean(c(confint_icep$X2.5..[i],confint_icep$X97.5..[i])))
    confint_icep$upr[i] <- c(confint_icep$Estimate[1]+confint_icep$X97.5..[i])
    confint_icep$lwr[i] <- c(confint_icep$Estimate[1]+confint_icep$X2.5..[i])
  }
}
ggplot(confint_icep, aes(x = Strain, y = Estimate, ymin = lwr, ymax = upr, color = abs(Estimate))) + ###Plot Comparison on x, estimate on y
  scale_x_discrete(limits = rev(levels(as.factor(confint_icep$Strain)))) + ###Rescale x, so the order is inverted 
  geom_errorbar() + geom_point() + ###Draw data
  coord_flip() +  theme(legend.position="none") + xlab("") +###Invert X and Y, hide legend
  ggtitle("Absolute intercept with 95%CI") ##Add a title
```
  
  Using the linear model and raw data, different displays can be plotted, for example, a boxplot of the  "Useful" data-points combined with the predictions (extrapolations) of the linear model.
  
```{r predict full dataset and plot}
##Add predictions to full dataset
library("modelr")
##Add predictions to full dataset
di_long <- add_predictions(di_long,disease_lmer,var="lmer.pred")
###
ggplot(data=di_long, aes(x=DPI,y=DI))+
  geom_boxplot(aes(color=Strain,group=DPI),data=filter(di_long, Useful=="Yes"))+
  geom_smooth(aes(y=lmer.pred,color=Strain), method="lm", alpha=0.6) +
  facet_wrap(~Strain)+
  ggtitle("Boxplots and linear fit for individual strains")
```
  
  Finally, of the methods discussed here, those relying on a linear model of raw observations over time, can be used for other kinds of observations, e.g. bacterial titers. 
  
  
  
  At least briefly, plot residuals of the lmm to assess if this appears normal (random)
  
```{r residuals}
di_long$resid[di_long$Useful=="Yes"] <- resid(disease_lmer)
ggplot(di_long[di_long$Useful=="Yes",], aes(x=DPI,y=resid)) + geom_jitter(aes(color=Strain)) 
```

# Survival Analysis

## Background
Survival analysis is different in several aspects from the earlier approaches. Survival analysis, also known that time-to-event analysis, builds on a different dataset, which can  be generated from raw disease indices. Survival analysis provides a way to analyze the time-to-event recordings within one population. In the case of the data used here, what is of interest is an event that can be referred to as *death*. Correctly defining *death* in this context is crucial for the outcome of the analysis.
*Death*, here, is defined as a subject reaching a certain disease index, from which it cannot recover. In this script the disease index that defines the threashold to *death* is called "cutoff".

```{r set cutoff}
###Set the cutoff by assigning a number to the variable cutoff
cutoff <- c(4) 
```

All of the observations that are above the dotted line in the plot below are *dead*, at the day they cross that line. Those that never cross the line are alive until the end of observations, and are "right censored", meaning that their event was not observed during the time this subject was observed.
```{r plot recorded data and cutoff}
ggplot(data=di_long) +
                  geom_jitter(aes(x=DPI, y=DI,color=Strain,shape=Batch)) +
                  geom_segment(aes(x=0, xend=max(DPI)+0.5,y=cutoff, yend=cutoff), linetype="dotted") +
                  labs(x = "Days post infection", y = "Disease indices", title="Scatterplot of disease indices\n cutoff plotted as dotted line") +
                  coord_cartesian(xlim=c(-0.1,10))
```
  
##How to make survival data from disease indices
  A "survival table" can be generated using the following code. This code works on the long table, generated in the beginning, and the cutoff variable defined above. 
  
```{r disable warnings, echo=F}
oldw <- getOption("warn")
options(warn = -1)
```

```{r make survival table, message = F}
###Generate survival table
library("stringr")
surv_from_DI <- data.frame(Subject=levels(di_long$subject),
                         Strain=str_split_fixed(levels(di_long$subject), "\\.",3)[,1],
                         Plant=str_split_fixed(levels(di_long$subject), "\\.",3)[,2],
                         Batch=str_split_fixed(levels(di_long$subject), "\\.",3)[,3])
###Fill survival table based on the di_long table. This generates warnings. These can be ignored and come from the min()
for (i in 1:max(as.numeric(di_long$subject))) { #Go by subject
  dummy <- di_long[as.numeric(di_long$subject)==i,] #generate dummy for the subject
  if (is.infinite(min(dummy$DPI[which(dummy$DI >= cutoff)]))) { #If none of the DI is greater than the cutoff (this is where warnings are generated, min on an empty object returns infinite and a warning!)
    surv_from_DI$End[i] <- max(dummy$DPI) #Generate a  observation, censoring at the maximum DPI recorded
    surv_from_DI$Death[i] <- 0 #Still alive, because it did not pass the cutoff
  } else { #If more than zero DI are greater than the cutoff
    surv_from_DI$End[i] <- min(dummy$DPI[which(dummy$DI >= cutoff)]) #Use the lowest DPI where condition is met
    surv_from_DI$Death[i] <- 1 #record as dead
  }
}
rm(dummy)
```
```{r toggle warnings, echo=F}
options(warn = oldw)
```


Kaplan-Meier estimates of survival are the basic tool of survival analysis. These can be estimated using the survfit function from the "survival"" package.

```{r load survival, message=F}
library("survival")
```
```{r survfit}
surv_DI_fit <- survfit(Surv(End, Death) ~Strain +strata(Batch), data=surv_from_DI)
````

The survminer package provides the ggsurvplot() function. This works nicely on datasets with few treatments. However, for the data presented here, I think it is easier to initially generate a data frame that contains the whole fit and plot with "normal" ggplot2

```{r generate KM dataframe}
library("stringr")
###Strata dummy generation, modified from kmiddleton / rexamples 
strata_dummy <-NULL
for(i in 1:length(surv_DI_fit$strata)){
      # add vector for one strata according to number of rows of strata
      strata_dummy <- c(strata_dummy, rep(names(surv_DI_fit$strata)[i], surv_DI_fit$strata[i]))
}
###Data frame generation inspired by a post by Hadley Wickham to the ggplot2 Googlegroup
surv_DI_fit.df <- data.frame( 
  time = surv_DI_fit$time, 
  n.risk = surv_DI_fit$n.risk, 
  n.event = surv_DI_fit$n.event, 
  surv = surv_DI_fit$surv, 
  strata = strata_dummy, 
  upper = surv_DI_fit$upper, 
  lower = surv_DI_fit$lower 
) 
zeros <- data.frame(time = 0, surv = 1, strata = names((surv_DI_fit$strata)), 
  upper = 1, lower = 1)

surv_DI_fit.df <- plyr::rbind.fill(zeros, surv_DI_fit.df) ###I dont want to load plyr because i guess it will interfere with dplyr...
rm(strata_dummy)
rm(zeros)
###Some stuff to rename other stuff, this needs to be adapted if other variables are used.

surv_DI_fit.df$Batch <- as.factor( str_split_fixed(
  matrix( nrow=length(surv_DI_fit.df$strata),ncol=2, unlist(strsplit(as.character(surv_DI_fit.df$strata),", ")), byrow=T )[,2],"=",2)[,2])
surv_DI_fit.df$Strain <- as.factor( str_split_fixed(
    matrix( nrow=length(surv_DI_fit.df$strata),ncol=2, unlist(strsplit(as.character(surv_DI_fit.df$strata),", ")), byrow=T )[,1],"=",2)[,2])

###End of data frame generation
###Start plotting
ggplot(surv_DI_fit.df, aes(time, surv, colour = Strain)) + 
  geom_step(aes(y = surv*100,linetype=Batch)) +
  facet_wrap(~Strain) +
  ggtitle("Survival estimates for all Batches")

```

##Analysis of survival curves and fits
Comparing the Kaplan-Meier survival estimates can be done in different ways.
###Logrank testing
The below produces all pairwise comparisons of the Kaplan Meier estimate of survival using a logrank test.
```{r pairwise logranks}
###Make a table of pairwise chisq pvalues, for the logrank test.
#Based on a post to the R Mailing list by T. Therneau
pw_logrank_test_type <- 0 ###0 for logrank, 1 for peto and peto
pw_logrank <- matrix(0., nlevels(surv_from_DI$Strain),nlevels(surv_from_DI$Strain))
for (i in 1:nlevels(surv_from_DI$Strain)) {
  for (j in (1:nlevels(surv_from_DI$Strain))[-i]) {
    datasubset <- droplevels(subset( surv_from_DI,
      surv_from_DI$Strain %in% (unique(surv_from_DI$Strain))[c(i,j)]))
    temp <- survdiff(Surv(End, Death)~Strain+strata(Batch), data=datasubset, rho=pw_logrank_test_type)
      
    pw_logrank[i,j] <- pchisq(temp$chisq, df=1, lower=F) ##df will always be 1 because this is pairwise
    }
}
colnames(pw_logrank) <- levels(surv_from_DI$Strain)
rownames(pw_logrank) <- levels(surv_from_DI$Strain)
#Make dummy adjustment table
pw_logrank_adjBon <- pw_logrank
#Fill adjusted pvalue table.
for (i in 1:ncol(pw_logrank)) {
  pw_logrank_adjBon[,i] <- cbind(p.adjust(pw_logrank[,i], method="bonferroni"))
}
```
````{R results="asis"}
stargazer::stargazer(pw_logrank_adjBon,type="html",title="Pairwise Chisq p-values (Bonferroni adjusted)")
```


### Regressions
  Generally a survival regression does not assume proportionality of hazards. A survival regression is fit to a distribution, defined by dist="".
```{r load survcomp and rms, message=F}
library("survcomp")
library("rms")
library("modelr")
```
```{r survival regression}
####Survival Regression###
###This is done using functions from rms.
###psm is a survival::survreg wrapper. but the output is more handle-able.

ddist <- datadist(surv_from_DI)
options(datadist="ddist")
psurv_gaus <- psm(Surv(End, Death) ~Strain, data=surv_from_DI, dist="gaussian")
psurv_logistic <- psm(Surv(End, Death) ~Strain, data=surv_from_DI, dist="logistic")
psurv_lnorm <- psm(Surv(End, Death) ~Strain, data=surv_from_DI, dist="lognormal")
psurv_wei <- psm(Surv(End, Death) ~Strain, data=surv_from_DI, dist="weibull")

###Same with survreg()
s_reg_gaus <- survreg(Surv(End, Death) ~Strain, data=surv_from_DI, dist="gaussian")
s_reg_logistic <- survreg(Surv(End, Death) ~Strain, data=surv_from_DI, dist="logistic")
s_reg_lnorm <- survreg(Surv(End, Death) ~Strain, data=surv_from_DI, dist="lognormal")
s_reg_wei <- survreg(Surv(End, Death) ~Strain, data=surv_from_DI, dist="weibull")  
aic.scores.psurv <- rbind(
  extractAIC(psurv_wei),
  extractAIC(psurv_gaus),
  extractAIC(psurv_logistic),
  extractAIC(psurv_lnorm))
###Make useable AIC table
rownames(aic.scores.psurv) <- c("Weibull", "Gaussian", "Logist", "Lognorm")
colnames(aic.scores.psurv) <- c("df", "AIC")
###Call table
````
````{R results="asis"}
stargazer::stargazer(aic.scores.psurv,type="html",title="AIC Scores")
```

From the table above, the model with the lowest AIC score can be chosen. For this analysis, this is the lognormal model, but this does not have to apply to other experiments.
Then, one can inspect that model for significance
```{r glht lnorm}
summary(glht(psurv_lnorm,linfct=mcp(Strain="Tukey")))
```

Again, one can then inspect the differences in the model, for example using pairwise comparisons of means.
```{r}
pairwise_confint_sreg <- as.data.frame(confint(glht(psurv_lnorm, mcp(Strain = "Tukey")))$confint)
pairwise_confint_sreg$Comparison <- rownames(pairwise_confint_sreg)
###Plot the comparisons, below may not be the most straight-foward way to plot this the way I want it, but it works.
ggplot(pairwise_confint_sreg, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr, color = Comparison)) + ###Plot Comparison on x, estimate on y
  scale_x_discrete(limits = rev(levels(as.factor(pairwise_confint$Comparison)))) + ###Rescale x, so the order is inverted 
  geom_errorbar() + geom_point() + ###Draw data
  coord_flip() +  theme(legend.position="none") + xlab("") +###Invert X and Y, hide legend
  ggtitle("Difference in means with 95% Confidence interval") ##Add a title
```

## Plotting of parametric survival regression
It is possible, but not really easy, to plot the generated curves. These curves are the result of fitting the data to a distribution in the earlier section. 
Doing this in a manner that is compatible with ggplot2 is not straightforward. Below is code to generate plots of the KM estimates per batch and the generated regression. This is performed for the four distributions above, and can be adapted to different distributions if necessary.
```{r generate curves of survreg}
library("tidyr")
###Step 1, extract the coefficients. These are relative to Strain1 because Strain is treatment contrasted.

for (i in 1:nlevels(surv_DI_fit.df$Strain)) { #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(
  Strain = rep(levels(surv_from_DI$Strain),4 ), #Fill with strains
  scale.wei = exp(unlist(coef_wei)), #weibull fit scale parameters
  scale.logistic = rep(s_reg_logistic$scale, nlevels(surv_from_DI$Strain)), #fill with logis scales
  scale.gaus = rep(s_reg_gaus$scale, nlevels(surv_from_DI$Strain)), #fill with gaus scales
  scale.lnorm = rep(s_reg_lnorm$scale, nlevels(surv_from_DI$Strain)), #fill with lnorm scale
  shape.wei =  rep(1/s_reg_wei$scale, nlevels(surv_from_DI$Strain)), #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_DI_fit.df$Strain)){
  if(i==1) {
    wei <- list()
    logis <- list()
    gaus <- list()
    lnorm <- list()
  }
  x <- levels(surv_DI_fit.df$Strain)[i]
  n <- c(1:max(surv_from_DI$End))
  data <- filter(sregparams, Strain==x)
  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)),
  Strain = rep(unlist(levels(surv_DI_fit.df$Strain)),each=max(surv_from_DI$End)),
  time = rep(c(1:max(surv_from_DI$End)), nlevels(surv_DI_fit.df$Strain))
)
##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 lnorm.sreg logis.sreg wei.sreg
levels(sreg_long$Distribution) <- c("Gaussian","Lognormal","Loglogistic","Weibull")
```

Now, these can be plotted and inspected visually..

```{r plot KM and curve}
###Plot of KM+Weibull
  ggplot(surv_DI_fit.df, aes(time, surv, colour = Strain)) +
  geom_step(aes(linetype=Batch)) +
  geom_line(data=sreg_curves,aes(y=wei.sreg),color="black") +
  facet_wrap(~Strain) + 
  ggtitle("Kaplan-Meier estimates and fit to\nWeibull distribution")
###Plot of KM+Logistic
  ggplot(surv_DI_fit.df, aes(time, surv, colour = Strain)) +
    geom_step(aes(linetype=Batch)) +
    geom_line(data=sreg_curves,aes(y=logis.sreg),color="black") +
    facet_wrap(~Strain) + 
    ggtitle("Kaplan-Meier estimates and fit to\n Logistic distribution")
###Plot of KM+Gaussian  
  ggplot(surv_DI_fit.df, aes(time, surv, colour = Strain)) +
    geom_step(aes(linetype=Batch)) +
    geom_line(data=sreg_curves,aes(y=gaus.sreg),color="black") +
    facet_wrap(~Strain) + 
    ggtitle("Kaplan-Meier estimates and fit to\nGaussian distribution")
###Plot of KM+Lognormal  
  ggplot(surv_DI_fit.df, aes(time, surv, colour = Strain)) +
    geom_step(aes(linetype=Batch)) +
    geom_line(data=sreg_curves,aes(y=lnorm.sreg),color="black") +
    facet_wrap(~Strain) + 
    ggtitle("Kaplan-Meier estimates and fit to\nLognormal distribution")
```

### Cox-Proportional hazards and hazard ratios
Different approaches to survival analysis, are based on analysing the hazards. The hazard is the probability of experiencing an event at a given timepoint. Many hazard based analysis assume that hazards are proportional between treatments, meaning that they differ by a fixed factor.
Hazards were strongly influenced by Cox and a basic model is the cox proportional hazards model.
```{r coxph and PH assumption}
###Cox-Proportional hazards####
#Build model
srv_coxph <- coxph(Surv(End, Death) ~Strain + strata(Batch), data=surv_from_DI) 
###Check porportionality of hazards
cox.zph(srv_coxph, transform = "log")
````

If the proportionality of hazards is violated within the dataset can be seen from the p-value(s) returned by cox.zph. 
One can then continue with a comparison of hazard ratios, however hazard ratios are not appropriate if the assumption of proportional hazards is violated.

````{r}
####Hazard ratio
haz_rats <- hazard.ratio(x= surv_from_DI$Strain, 
                         surv.time = surv_from_DI$End, 
                         surv.event = surv_from_DI$Death, 
                         strat = surv_from_DI$Batch, 
                         method.test = "wald" ) ###Overall hazard ratios

###Pairwise hazard ratios / modified from the pairwise chisq calculation
pw_hazrats <- matrix(0., nlevels(surv_from_DI$Strain),nlevels(surv_from_DI$Strain))
for (i in 1:nlevels(surv_from_DI$Strain)) {
  for (j in (1:nlevels(surv_from_DI$Strain))[-i]) {
    datasubset <- droplevels(subset( surv_from_DI,
      surv_from_DI$Strain %in% (unique(surv_from_DI$Strain))[c(i,j)]))
    temp <- hazard.ratio(
      x= datasubset$Strain, 
      surv.time = datasubset$End, 
      surv.event = datasubset$Death, 
      strat = datasubset$Batch, 
      method.test = "likelihood.ratio" ###Define test to determine p.
      )
    pw_hazrats[i,j] <- temp$p.value
    }
}
colnames(pw_hazrats) <- levels(surv_from_DI$Strain)
rownames(pw_hazrats) <- levels(surv_from_DI$Strain)
```
```{r results = "asis"}
stargazer::stargazer(pw_hazrats,type="html",title="Pairwise hazard ratio pvalues")
```

  If the hazards are found to be non-proportional, it might be a good idea to perform survival regression analysis, or pairwise log-rank testing (see earlier) instead of hazard ratio tests. Alternatives may also come from the use of Cox mixed-effects model, similar to linear mixed-effects model, but with a different type of response variable.
  
```{r coxme}
library("coxme")
cme <-  coxme(Surv(End, Death) ~Strain + (1|Batch), data=surv_from_DI)
anova(cme)
summary(glht(cme,linfct=mcp(Strain="Tukey")))
````

# Comparison of analyis method results.

An inherent question when analyzing experimental data, is which analysis produced which result and why. 
Below, the outputs from the three parametric analysis performed are printed, so they can be compared.

```{r summaries}
summary(auc_lmer)
summary(disease_lmer)
summary(psurv_lnorm)
```

An easier comparisons might be accomplished with the compact letter display.

```{r clds}
cld(glht(auc_lmer, linfct=mcp(Strain="Tukey")))

#Compact letters for lmerTest objects are a little tricky. This solution comes from the rcompanion.
### Extract lsmeans table

lmerlsm <- difflsmeans(disease_lmer)$diffs.lsmeans.table

Comparison = str_split_fixed(rownames(lmerlsm),"Strain ",2)[,2]

### Produce compact letter display

library(rcompanion)

# cldList(comparison = Comparison,
#         p.value    = p.adjust(lmerlsm$'p-value',
#                        method =  "bonferroni")   ,
#         threshold = 0.05)


cld(glht(srv_coxph,linfct=mcp(Strain="Tukey")))
cld(glht(psurv_lnorm,linfct=mcp(Strain="Tukey")))
cld(glht(cme,linfct=mcp(Strain="Tukey")))
```

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