---
title: "Plot and statistical analysis"
author: ""
date: "`r format(Sys.time(), '%H:%M, %d %B, %Y')`"
output:
    html_document:
      toc: true
      toc_float: true
      toc_depth: 4
      theme: cerulean
      highlight: kate
---

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

## 1. Introduction

The analysis reported here is part of the manuscript Romand et al. 2021. The script performs statistical analysis on data grouped by `Treatment` and `Plant Line` and plots a graph. Different approaches are used depending on the normality of the data as determined by the Shapiro test. If data follow a normal distribution, then the script performs Anova and post-hoc Tukey tests and draws a plot with confidence interval centered on mean. If data do not follow a normal distribution then a Kruskal-Wallis and a pot-hoc Dunn test are performed. The same plot is drawn using the medians and boostrap calculated 95% confidence intervals. The plots are also saved as pdf and vector (SVG) files.


## 2. Data

A 3 column data set in `xlsx`. Column 1, `Line`; column 2, `Treatment` and column 3 the measured quantity, in this case `G4P`. 
The 2 first columns are factors that are used to group the data.

This script can be simply modified to analyse different parameters or different treatments. It can also be modified to analyse data grouped by additional factors such as Time. 

## 3. Data analysis

### 3.1. Code

#### Set variables for personalisation

```{r}
#1. Data
#------------------------------------------------------------------------------

# Put your input file in variable DATA 
DATA <- "G4P_mutants.xlsx"


#2. Parameters for plot
#------------------------------------------------------------------------------
# Points color
COLOUR <- "#00B050"

# Select a palette for CI color
PALETTE <- "PuRd"
```

<br/>

`r {"Working directory:"}`
`r  getwd()`

<br/>

#### Setup R packages

```{r message=FALSE, results= FALSE, warning = FALSE}
if (!require(ggbeeswarm)) { install.packages("ggbeeswarm", repos = "http://cran.us.r-project.org") }
if (!require(Rmisc)) { install.packages("Rmisc", repos = "http://cran.us.r-project.org") }
if (!require(ggplot2)) { install.packages("ggplot2", repos = "http://cran.us.r-project.org") }
if (!require(RColorBrewer)) { install.packages("RColorBrewer", repos = "http://cran.us.r-project.org") }
if (!require(dplyr)) { install.packages("dplyr", repos = "http://cran.us.r-project.org") }
if (!require(rstatix)) { install.packages("rstatix", repos = "https://cloud.r-project.org") }
if (!require(purrr)) { install.packages("purrr", repos = "http://cran.us.r-project.org") }
if (!require(tidyr)) { install.packages("tidyr", repos = "http://cran.us.r-project.org") }
if (!require(magrittr)) { install.packages("magrittr", repos = "http://cran.us.r-project.org") }
if (!require(boot)) { install.packages("boot", repos = "http://cran.us.r-project.org") }
if (!require(devtools)) { install.packages("devtools", repos = "http://cran.us.r-project.org") }
if (!require(knitr)) { install.packages("knitr", repos = "http://cran.us.r-project.org") }
if (!require(svglite)) { install.packages("svglite", repos = "http://cran.us.r-project.org") }


library(ggplot2)
library(ggbeeswarm)
library(Rmisc)
library(readxl)
library(RColorBrewer)
library(rstatix)
library(boot)
library(magrittr)
library(dplyr)
library(purrr)
library(tidyr)
library(rstatix)
library(knitr)
library(svglite)

```

#### Functions used in the script

```{r}
#===============================================================================
# Draw the plot if data follow a normal distribution
#===============================================================================

plot_normal <- function(df, my_colours, my_summary) {
    p<-ggplot(data = df, aes(x=Line, y=G4P)) +
        geom_quasirandom(dodge.width=0.8,alpha = 0.6, colour = COLOUR)+
        geom_pointrange(data = my_summary, aes(ymin=G4P-ci, ymax=G4P+ci, color=Line), position=position_dodge(width=0.8))+
        scale_colour_manual(values=my_colours)+
        facet_wrap(~Treatment, strip.position = "bottom", scales = "free_x") +
        theme_classic() + 
        theme(panel.spacing = unit(0, "lines"), 
              strip.background = element_blank(),
              strip.placement = "outside")+
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
    print(p)
}

#===============================================================================
# Draw the plot if data do not follow a normal distribution
#===============================================================================

plot_not_normal <- function(df, my_colours, conf_int) {
    p<-ggplot(data=df, aes(x=Line, y=G4P)) +
        geom_quasirandom(dodge.width=0.8,alpha = 0.6, colour = COLOUR)+
        geom_linerange(data = booted_summary, aes(ymin=lower_ci_perc, ymax=upper_ci_perc, color=Line), position = position_dodge(width=0.8))+
        geom_point(data=booted_summary, aes(y=G4P, color=Line), size = 2, 
                   position=position_dodge(width=0.8))+
        scale_colour_manual(values=my_colours)+
        facet_wrap(~Treatment, strip.position = "bottom", scales = "free_x") +
        theme_classic() + 
        theme(panel.spacing = unit(0, "lines"), 
              strip.background = element_blank(),
              strip.placement = "outside")+
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
        print(p)
       }

#===============================================================================
# Check normality of the data. Exit the loop if one set of data doesn't follow a 
# normal distribution.
# Return TRUE if all data follow a normal distribution
#===============================================================================

check_normality <- function(shapiro_df) {
    # Data are supposed to follow normal distribution
    flag_normal <- TRUE
    
    for (i in 1 : nrow(shapiro_df)) {
        if(shapiro_df[i, 4] > 0.05) {
            
        } else {
            
            flag_normal <- FALSE
            break
        }
    }
    return(flag_normal)
}

#===============================================================================
#Function to calculate the median (for bootstrapping)
#===============================================================================

samplemedian <- function(d, i) {
    median(d[i])
}

#===============================================================================
#Function for Dunn test
#===============================================================================
test_dunn <- function() {
    pval <- as.data.frame(df %>% group_by(Treatment) %>% dunn_test(G4P ~ Line, p.adjust.method = "BY"))
    print(df %>% group_by(Treatment) %>% dunn_test(G4P ~ Line, p.adjust.method = "BY"))
    return(pval)
}
```

#### Main script

```{r fig.show='hide',  message=FALSE, warning=FALSE, results = 'hide', echo=TRUE, fig.retina = 1}
# Import Data
df <- read_excel(DATA, col_types = c("text", "text", "numeric"))
df$Line <- make.names(df$Line)

# Define color
my_colours = brewer.pal(n = 9, PALETTE)[9:3]

# Determining data normality status
shapiro_df <- df %>%group_by(Treatment, Line)%>%
    summarise(statistic = shapiro.test(G4P)$statistic, 
              p.value = shapiro.test(G4P)$p.value)

flag_normal <- check_normality(shapiro_df)


# Data treatement according to normality status
if(flag_normal == TRUE) {
    #print("Data follow a normal distribution")
    
    # Summary
    my_summary <- summarySE(df, measurevar="G4P", groupvars=c("Line", "Treatment"))
    
    # Plot
    plot_normal(df, my_colours, my_summary)
    
    # Stats
    anova_results <- df %>% group_by(Treatment) %>%  anova_test(G4P ~ Line)
    
    tukey_results <- df %>% group_by(Treatment) %>%  tukey_hsd(G4P ~ Line)
    
    
    # Save plot
    ggsave("my_ggplot1.pdf", width=4, height=5)
    ggsave("my_ggplot1.svg", width=4, height=5)
    
} else {
    #print("Data do not follow a normal distribution")
    
    # Calculate Median and 95%CI by bootstrap (Largely based on code by Peter Kamerman: 
    #https://www.painblogr.org/2017-10-18-purrring-through-bootstraps.html)
    
    #1. Set the seed
    set.seed(12345)
    
    #2. Create nested datafile to allow boostrapping on each subgroup.
    my_data_nested <- df %>% group_by(Line,Treatment) %>% nest()
    
    #3. Add boostrapped median values to my_nested_data
    my_data_nested %<>% dplyr::mutate(booted = purrr::map(.x = data, ~ boot::boot(data = .x$G4P,statistic = samplemedian, R = 10000, stype = "i")))
    
    #4. Plot each <S3: boot> object to check boostrapping worked!
    ## Note: Saved to an object (plots) to stop the summary being printed
    plots<-purrr::map(.x = my_data_nested$booted, ~ plot(.x))
    print(plots)
    
    #5. Add boostrapped 95%CI around median values to my_nested_data, 
    #95%CI determined based on "perc" & "bca"
    my_data_nested %<>% mutate(booted_ci = purrr::map(.x = booted, 
                                                      ~ boot::boot.ci(.x,conf = 0.95, 
                                                                      type = c("perc", "bca"))))
    
    #6. Table of each <S3: bootci> object to check
    ci_table <-purrr::map(.x = my_data_nested$booted_ci, ~ print(.x))
    
    #7. Get the stats (median and 95% CI) from the nested dataframe into a regular dataframe 
    #(my_data_booted_summary) with 95%CI determined by "bca" or "perc" method
    booted_summary <- my_data_nested %>% mutate(G4P= purrr::map(.x = booted_ci, ~ .x$t0),
                                                lower_ci_bca = purrr::map(.x = booted_ci,~ .x$bca[[4]]),	
                                                upper_ci_bca = purrr::map(.x = booted_ci,~ .x$bca[[5]]),
                                                lower_ci_perc = purrr::map(.x = booted_ci,~ .x$perc[[4]]),
                                                upper_ci_perc = purrr::map(.x = booted_ci,~ .x$perc[[5]])) %>%  
        dplyr::select(-data, -booted, -booted_ci) %>% tidyr::unnest(cols = c(G4P, lower_ci_bca, upper_ci_bca, lower_ci_perc, upper_ci_perc))
    
    # Plot
    plot_not_normal(df, my_colours, booted_summary)
    
    # Stats
     
    kruskal_pval <- (df %>% group_by(Treatment)%>%kruskal_test(G4P ~ Line)) %>% select(Treatment, p)
    pval_dunn <- test_dunn()
    
    # Save plot
    ggsave("my_ggplot1.pdf", width=4, height=5)
    ggsave("my_ggplot1.svg", width=4, height=5)
}
```


### 3.2. Results

#### a. Normality status

```{css, echo = FALSE}
caption {
  color: black;
      font-weight: bold;
      font-size: 1.0em;
    }
```


`r {kable(shapiro_df, caption = "Table 1 : Shapiro test results.", digits = 30)}`

`r if(flag_normal == TRUE) {"**Data follow a normal distribution. Now calculating means and confidence intervals for each line-treatment combination.**"} else {"**Data do not follow a normal distribution. Now calculating medians and confidence intervals by a boostrap procedure for each line-treatment combination.**"}`  

<br/>


```{r echo=FALSE, fig.retina= 1, fig.show='hold', message=FALSE, warning=FALSE, results='hide'}
if(flag_normal == FALSE) {
  plots<-purrr::map(.x = my_data_nested$booted, ~ plot(.x))
  print(plots)
}
```

`r if(flag_normal == FALSE) {"**Fig 1. : Boostrapping plots**"}`

<br/>

#### b. Final data plot  


```{r echo=FALSE,  message=FALSE, warning=FALSE, results = 'hide', fig.retina = 1}
if(flag_normal == TRUE) {
  # Plot
  plot_normal(df, my_colours, my_summary)
} else { 
  plot_not_normal(df, my_colours, booted_summary)
}  
```

`r if(flag_normal == TRUE) {"**Fig 1. : Confidence interval plot**"} else {"**Fig 2. : Confidence interval plot**"}`

<br/>

#### c. Statistical analysis

##### Summary statistics


`r if(flag_normal == TRUE) {kable(my_summary, caption = "Table 2 : Statistical summary.", digits = 3)  } else { kable(booted_summary, caption = "Table 2 : Statistical summary.", digits = 100)}`

<br/>

##### Statistical tests within each treatment condition.


`r if(flag_normal == TRUE) {kable(anova_results, caption = "Table 3 : Results of Anova test.", digits = 30) } else { kable(kruskal_pval, caption = "Table 3 : Results of Kruskal-Wallis test.", digits = 100)}`

<br/>

`r if(flag_normal == TRUE) {kable(tukey_results[,c("Treatment","group1","group2","p.adj","p.adj.signif")], caption = "Table 4 : Results of Tukey test.", digits = 30) } else { kable(pval_dunn[, c(1, 3:4, 8:10)], caption = "Table 4 : Results of Dunn test.", digits = 100)}`

<br/>

## 4. R session information

```{r}
InfoSession <- devtools::session_info()
print(InfoSession)

```

## 5. Citations
1. Bache, Stefan Milton and Wickham, Hadley (2020). magrittr: A Forward-Pipe Operator for R. R package version 2.0.1.
https://CRAN.R-project.org/package=magrittr  
  
2. Canty, Angelo and Ripley, Brian (2021). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-28.  

3. Clarke, Erik, and Scott Sherrill-Mix. 2017. Ggbeeswarm: Categorical Scatter (Violin Point) Plots. https://CRAN.R-project.org/package=ggbeeswarm.

4. Henry, Lionel and Wickham, Hadley (2020). purrr: Functional Programming Tools. R package version 0.3.4.
https://CRAN.R-project.org/package=purrr

5. Hope, Ryan M. 2013. Rmisc: Rmisc: Ryan Miscellaneous. https://CRAN.R-project.org/package=Rmisc.

6. Kassambara, Alboukadel. 2021. Rstatix: Pipe-Friendly Framework for Basic Statistical Tests. https://CRAN.R-project.org/package=rstatix.

7. Neuwirth, Erich. 2014. RColorBrewer: ColorBrewer Palettes. https://CRAN.R-project.org/package=RColorBrewer.

8. R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

9. Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

10. Wickham, Hadley, and Jennifer Bryan. 2019. Readxl: Read Excel Files. https://CRAN.R-project.org/package=readxl.

11. Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2021. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.

12. Wickham, Hadley (2021). tidyr: Tidy Messy Data. R package version 1.1.3. https://CRAN.R-project.org/package=tidyr

13. Wickham, Hadley, Hester, Jim and Chang, Winston (2021). devtools: Tools to Make Developing R Packages Easier. R package version 2.4.2. https://CRAN.R-project.org/package=devtools

