---
title: "Approaches to analyze experimental *R. solanacearum* infections:\n*Supplementary Material 3:*\nDatasets from the Allen Lab\n Part I: Data import, and formatting" 
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, datasets from Lowe, et al., MPMI, 2015 and Lowe-Power, et al., mBio, 2016 are prepared for analysis. Each dataset is analyzed seperately in part II and III, respectively.

# Data Type
The disease index score is based the severity of symptoms. A symptom read-out (usually "% wilted leaves") is turned into a numeric score. This score ranges from 0-4, and each number is used to reflect a certain interval of disease severity. The value of the disease index can never be below 0, or above 4. Therefore, the disease index is not a continuous response variable, but is ordinal. 

# Data Import

## Specifying the data

Specify directory (if running from console..)
```{r workingdirectory, eval=F}
###Install all required packages:
#install.packages(c("tidyr","dplyr","ggplot2","MESS","lme4","lmerTest","multcomp","survival","psm","coxme","stargazer","survcomp","rms","modelr","party","stringr"))


###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)
````

Specify tables.
```{r input definition}

###Name of the file to be read
mpmi2015 <- c("S3_LOWE_MPMI_2015.csv") 
mBio2016 <- c("S3_Lowe-Power_mBio_2016.csv")
````

These datasets from the research group of Caitilyn Allen at University of Wisconsin have previously been published. 

## Reading data and formatting

The data can be read into *R* and stored in a data frame. I am calling that data frame disease.index
````{r read data, eval=T}
di_mpmi <- as.data.frame(read.table(mpmi2015, header=T,
                                          sep=";" , ###Sets the seperator of the csv file
                                          dec=","), ###Sets the decimal operator of the csv file
                                          stringsAsFactors=F) ###
di_mbio <- as.data.frame(read.table(mBio2016, header=T,
                                          sep=";" , ###Sets the seperator of the csv file
                                          dec=","), ###Sets the decimal operator of the csv file
                                          stringsAsFactors=F) ###
````
  Using str() one can see if the tables were properly imported.
  
````{r inspect structure}
str(di_mpmi)
str(di_mbio)
````

These tables need to be made long.

````{r load tidyr, message=F}
library("tidyr")
```
````{r generate long table}
mpmi2015_long <- gather(di_mpmi, key=Plant, 
                     value=DI,
                     c(A, B, C, D, E, F ,G ,H ,I ,J, K, L, M ,N, O, P), 
                     na.rm=F)
mbio2016_long <- gather(di_mbio, key=Plant, 
                     value=DI,
                     c(A, B, C, D, E, F ,G ,H ,I ,J), 
                     na.rm=F)
```

Now, it is time to generate the "subject" variable, unique for each subject. A unique identifier can be generated by combining information from Strain, Plant and Batch. 

```{r}
mpmi2015_long$subject <- interaction(mpmi2015_long$Strain,
                                     mpmi2015_long$Plant,
                                     mpmi2015_long$Batch)
mbio2016_long$subject <- interaction(mbio2016_long$Strain,
                                     mbio2016_long$Plant,
                                     mbio2016_long$Batch)
###Factorize batch information
mpmi2015_long$Batch <- as.factor(mpmi2015_long$Batch)
mbio2016_long$Batch <- as.factor(mbio2016_long$Batch)
````

##Defining the contrasts

  A crucial step that will influence all outputs of the statistical analysis is setting the "contrasts". Contrasts specifies the "baseline" for each of the variables. "Treatment" contrasts specify that the first alphabetical level will be used as a reference for all others (see Strain below), while a "sum" constrast means that the reference value is the mean across all levels of that variable.
  
```{r contrasts}
####Specify what should be "appropriate" contrasts
contrasts(mpmi2015_long$Strain) <- "contr.Treatment"###First alphabetical strain will be the baseline
contrasts(mbio2016_long$Strain) <- "contr.Treatment"###First alphabetical strain will be the baseline

contrasts(mpmi2015_long$Batch) <- "contr.sum" ###Batches will be averaged to generate the baseline!
contrasts(mbio2016_long$Batch) <- "contr.sum"

```
  
##Adding a censoring variable

The below adds the variable that defines observations that are relevant for disease development.

```{r censor dead plants}
##Censoring for both datasets.

di_long <- na.omit(mpmi2015_long)
###Slightly modified with as.numeric
oldw <- getOption("warn")
options(warn = -1)
###Interesting R-related observation, the below does not work when subsetting is done with filter(), 
###because filter does not retain rownames!
for (i in 1:max(as.numeric(di_long$subject))) { ###Go by subject
  dummy1 <- di_long[as.numeric(di_long$subject)==i,] ##Create a first dummy object, that is a subset of the full data containing the current subject 
  if(min(dummy1$DI) == 0){
    #remove those observations that are before disease onset, except the one directly before disease onset.
    dummy1 <- dummy1[dummy1$DPI %in% (max(dummy1$DPI[dummy1$DI==0]):max(dummy1$DPI)),]
  } 
  if (max(dummy1$DI) == 4) { ###If this subject dies at some point
    dummy2 <- dummy1[dummy1$DI==4,] ###Make a new dummy object, that only contains those recordings where DI=4
    NEW <- di_long[as.numeric(di_long$subject)==i & (di_long$DPI %in% min(dummy1$DPI):min(dummy2$DPI)),] ###Generate data subset "NEW", which contains those observations for a subject, that are between (including) the last recording where DI=0 and the first recording where disease index is 4.
  } else { ###If this dubject does not die
    NEW <- dummy1 ###New is the same as dummy1
  }
  di_long$Useful[rownames(di_long) %in% rownames(NEW)] <- c("Yes") ###All of those row(names) that are part of the "NEW" object are useful. Therefore these receive status "Yes" in column "Useful"
}
di_long$Useful[which(is.na(di_long$Useful))] <- c("No") ###Those that are not yes (and therefore are NA) become No
rm(dummy1,dummy2)

mpmi2015_long <- di_long


##Same for mBio dataset
di_long <- na.omit(mbio2016_long)

###Interesting R-related observation, the below does not work when subsetting is done with filter(), 
###because filter does not retain rownames!
for (i in 1:max(as.numeric(di_long$subject))) { ###Go by subject
  dummy1 <- di_long[as.numeric(di_long$subject)==i,] ##Create a first dummy object, that is a subset of the full data containing the current subject 
  if(min(dummy1$DI) == 0){
    #remove those observations that are before disease onset, except the one directly before disease onset.
    dummy1 <- dummy1[dummy1$DPI %in% (max(dummy1$DPI[dummy1$DI==0]):max(dummy1$DPI)),]
  } 
  if (max(dummy1$DI) == 4) { ###If this subject dies at some point
    dummy2 <- dummy1[dummy1$DI==4,] ###Make a new dummy object, that only contains those recordings where DI=4
    NEW <- di_long[as.numeric(di_long$subject)==i & (di_long$DPI %in% min(dummy1$DPI):min(dummy2$DPI)),] ###Generate data subset "NEW", which contains those observations for a subject, that are between (including) the last recording where DI=0 and the first recording where disease index is 4.
  } else { ###If this dubject does not die
    NEW <- dummy1 ###New is the same as dummy1
  }
  di_long$Useful[rownames(di_long) %in% rownames(NEW)] <- c("Yes") ###All of those row(names) that are part of the "NEW" object are useful. Therefore these receive status "Yes" in column "Useful"
}
di_long$Useful[which(is.na(di_long$Useful))] <- c("No") ###Those that are not yes (and therefore are NA) become No
rm(dummy1,dummy2)
mbio2016_long <- di_long
options(warn = oldw)
```

## Table quality check

```{r load broom, message=F}
library("broom")
```
```{r inspect table}
str(mpmi2015_long)
tidy(mpmi2015_long) ###This can be used to assess the descriptive statistics of the data.

str(mbio2016_long)
tidy(mbio2016_long) ###This can be used to assess the descriptive statistics of the data. 

```

#Table export

```{r table saving, eval = T}
write.csv(mbio2016_long, file = "S3_mBio2016_long.csv")
write.csv(mpmi2015_long, file = "S3_mpmi2015_long.csv")
```

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