---
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 import

The script will analyse a data file placed in the working directory with the following layout. 

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

In this example we have 2 levels "N" and "NoN" for Treatment and 3 levels "WT2", "QM", "rsh1-1" for Line.

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 <- "QY.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(gridExtra)) { install.packages("gridExtra", repos = "http://cran.us.r-project.org") }
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(gridExtra)
library(ggplot2)
library(ggbeeswarm)
library(Rmisc)
library(readxl)
library(RColorBrewer)
library(rstatix)
library(gridExtra)
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=FvFm)) +
        geom_quasirandom(dodge.width=0.8,alpha = 0.6, colour = COLOUR)+
        geom_pointrange(data = my_summary, aes(ymin=FvFm-ci, ymax=FvFm+ci, color=Line), position=position_dodge(width=0.8))+
        scale_colour_manual(values=my_colours)+
        scale_y_continuous(breaks=seq(0,0.8,0.2), limits=c(0,0.95),expand = c(0, 0))+
        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=FvFm)) +
        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=FvFm, color=Line), size = 2, 
                   position=position_dodge(width=0.8))+
        scale_colour_manual(values=my_colours)+
        scale_y_continuous(breaks=seq(0,0.8,0.2), limits=c(0,0.95),expand = c(0, 0))+
        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(FvFm ~ Line, p.adjust.method = "BY"))
    print(df %>% group_by(Treatment) %>% dunn_test(FvFm ~ 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"))

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

# Determine data normality using the Shapiro test
shapiro_df <- df %>%group_by(Treatment, Line)%>%
    summarise(statistic = shapiro.test(FvFm)$statistic, 
              p.value = shapiro.test(FvFm)$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="FvFm", groupvars=c("Line", "Treatment"))
    
    # Plot
    plot_normal(df, my_colours, my_summary)
    
    # Stats
    anova_results <- df %>% group_by(Treatment) %>%  anova_test(FvFm ~ Line)
    
    tukey_results <- df %>% group_by(Treatment) %>%  tukey_hsd(FvFm ~ 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$FvFm,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(FvFm= 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(FvFm, 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(FvFm ~ 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 test and summary statistics

```{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. : Plots of bootstrapping process.**"}`

<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

