complete.analyses.Rmd

---
title: "R Notebook"
output:
  html_document:
    df_print: paged
  pdf_document: default
---
# About 

This is an example of fungicide sensitivity test in vitro for Sclerotinia sclerotiorum by discriminatory concetrations approach:

Serial dilution method with 6 concentrations for each fungicide from 22 baseline isolates (never exposed to fungicides) + 22 from between farmer fields and fungicide field trials

EC50 estimated from dose-response curve by ezec-package R (Kamvar, 2016)

Identify a single concentration that can be a discriminatory concentration

Use DC to assess sensitivity of large collection

Data analyses using reproducible methods

The 23 own created functions were compelled in a package FastconcR in https://github.com/EdgarNietoFungi/FastconcR



#  Background   
Knowledge of baseline sensitivity of a pathogen  population to fungicides is important for detection of resistance and for achieving effective disease management strategies (Brent and Holloman 2007). QoI, DMI, MBC and SDHI fungicides are used to control _S. sclerotiorum_ in an IPM approach in soybean and dry bean in the U.S. DMI fungicides are considered to medium risk of resistance emergence, whereas QoI, SDHI and MBC are high risk (Fungicide Resistance Action Committee 2016). MBC and dicarboximide fungicides are used for white mold in Brazil, wherein resistance to MBC fungicide, Thiophanate methyl, has been reported (Lehner et al. 2015). Fungicide resistance has also been reported in China to the MBC, Carbendazim and dicarboximide, Dimethachlon (Ma et al. 2009).

**Commonly Used Fungicides For WM**
------------------------------------
|Fungicide| Efficacy of Control| Launched|
|:---:|:----:|:-----:|:-----:|
|Thiophanate methyl (MBC)| fair control| 1970|
|Tetraconazole (DMI)|  fair control| 1988|
|Picoxystrobin (QoI)| good-very good control| 2000|
|Boscalid (SDHI)| very good control| 2003|

Examples of these include: **Topsin, Domark, Aproach & Endura**. Several other products are mixtures

MBC: Methyl Benzimidazole Carbamates 
DMI: DeMethylation Inhibitors 
QoI: Quinone outside Inhibitors 
SDHI: Succinate dehydrogenase inhibitors 

Field applications of fungicides in the United States on soybean show control of _S. sclerotiorum_ is low to fair for QoI, DMI and MBC fungicides, and good for SDHI fungicides (The North Central Regional Committee on Soybean Diseases 2015). However, it is unknown whether efficacy is related to reduce sensitivity within populations under study. 

 
```{r  Packages, warning= FALSE, include=FALSE}
#
using<-function(...) {
    libs<-unlist(list(...))
    req<-unlist(lapply(libs,require,character.only=TRUE))
    need<-libs[req==FALSE]
    if(length(need)>0){ 
        install.packages(need)
        lapply(need,require,character.only=TRUE)
    }
}

using("tidyverse","broom","ezec","DescTools", "agricolae", "multcompView","ggpubr", "patchwork")

# library(tidyverse)
# library(broom)
# library(ezec)
# library(DescTools)


#library(rstatix)
#library(ggpubr)

# library(agricolae)
# library(multcompView)
# library(ggpubr)
# library(patchwork)

#library(plyr)

###
```
# Hypotheses  

 
1. Isolates _S. sclerotiorum_ collected from producer fields and fungicide research trials, which already have been in contact with **MBC (FRAC group 1), DMI (FRAC group 3), SdhI (FRAC group 7), and QoI (FRAC group 11) fungicides** will have less sensitivity **(except FRAC group 3)** compared to the baseline, because these fungicides have a high risk of resistance and resistance to mostly these groups have been reported previously in _S. sclerotiorum_;
2. Isolates _S. sclerotiorum_ collected from dry bean will have lower sensitivity compared to the isolates from soybean for the three group fungicides **(except FRAC group 3)**, which have fewer fungicide applications;
3. Isolates  _S. sclerotiorum_ from Brazil and Mexico will have less sensitivity compared to the USA isolates for the three group fungicides **(except FRAC group 3)** because of different populations and management practices used in tropical and sub-tropical environment.


```{r Objects, warning=FALSE}
baseline.isolates <- c(1, 8, 12, 20, 21, 74, 87, 118, 123, 129, 449, 461, 467, 475, 558, 564, 568, 581, 645, 667, 800)

survey.isolates <-  c( 2385, 2386, 2098, 2099, 2100, 2139, 2140, 2143, 2220, 2222, 2223, 2320, 2362, 2388, 2390 )

Farmer.fields <- c(136:442, 455, 456, 466, 468, 470, 471, 478:554, 602:610, 612:635, 671, 672, 682:690, 695:797, 810:823, 834, 835, 848:1024, 1255: 1326, 1491: 1500, 1661:1670, 1831: 2246, 2303: 2342,2362: 2381, 2393: 2573)

baseline.isolates.2 <- c(1: 136, 444: 454, 457: 465, 467, 469, 472:477, 555:601, 611, 636:670, 673: 680, 691:694 , 798: 809,  824: 833, 836: 847)

drybean <- c(1, 5, 12,13:118, 123:128, 133:135, 145, 146, 152, 155:160, 182: 185, 194:200, 205, 220, 223, 248, 253:255, 274, 279, 280, 290, 294, 304: 309, 323, 358, 393, 395, 396, 397, 400:402, 405, 408, 409, 434, 443:465, 467, 469: 493, 495:505, 555:613, 615: 760, 762: 819, 824:833, 835:855, 858:921, 966: 971, 980, 981, 985: 996, 998: 1009, 1220:1225 , 1857: 1940, 2442: 2569)

soybean <- c(143, 147, 181, 187, 188, 189,202, 257: 259, 264:268, 276, 281, 289, 293, 295, 310 , 399, 412:417, 419, 425:427, 439, 440, 494, 506:554, 834, 972:979, 982, 997, 1010: 1022, 1025:1219, 1229: 1856, 1941:2441 )

SI.Production.Field <- c( 698:741, 743:744, 746:760, 762: 778, 
786: 797, 810:819, 848: 855, 858: 914)

SI.Screening.Nursery.Field <- c(444:454, 457:465, 467, 469, 472:477, 555: 582, 584: 601, 611, 636:670, 673:681,691:694, 798:809, 824: 833 , 836:847)# also taking into account WM Monitor 
mexican.isolates <- c(242, 248, 615:632, 779, 1857: 1940)
brazilian.isolates <-  c(400:402, 972:979, 1010:1022, 1227:1254, 2417:2569)


whitemold.inventory <- read_csv("data/whitemold.inventory.csv") %>% mutate(Field_Year= paste (Field, Year,sep = "_year_"))

for (variable in unique(whitemold.inventory$State)) {
    assign(variable, whitemold.inventory %>% filter (State == variable) %>% pull(ID), envir = .GlobalEnv)
}

for (variable in unique(whitemold.inventory$Field_Year)) {
    assign(variable, whitemold.inventory %>% filter (Field_Year == variable)%>% pull(ID), envir = .GlobalEnv)
}
IDA <- ID
#MAKE THE SAME BY COUNTY
for (variable in unique(whitemold.inventory$County)) {
    assign(variable, whitemold.inventory %>% filter (County == variable)%>% pull(ID), envir = .GlobalEnv)
}


```
# Functions
```{r Functions}
#Get range function
get_range <- function(mynumber) {
  bb <- boxplot.stats(mynumber)
  cc <- bb$stats
  dd <- max(cc)
  ee <- min(cc)
  return(data.frame(upper = dd, lower = ee))
}

##Function.Reading data for serial.dilution
reading_data_serial_dilution <- function(filename){
  data <-
  read.csv(filename)
 # data <- subset(dareading_data("data/Benomyl.Colletotrichum.csv")ta, select = -X)
  data$repeats <- rep_len(1:4, length.out = nrow(data))
  
   # renaming ID for the current names
 data$ID[ data$ID =="12B"]  <- "12" 
 #Although 74SS1 belong to # 74, the one what I used was specifically this:74SS1
 data$ID[ data$ID =="74SS1"]  <- "74" 
 data$ID[ data$ID =="62-02"]  <- "2098" 
 data$ID[ data$ID =="62-03"]  <- "2099" 
 data$ID[ data$ID =="62-04"]  <- "2100" 
 data$ID[ data$ID =="78-01"]  <- "2139" 
 data$ID[ data$ID =="78-02"]  <- "2140" 
 data$ID[ data$ID =="78-05"]  <- "2143" 
 data$ID[ data$ID =="S-01"]  <- "2320" 
 data$ID[ data$ID =="H-01"]  <- "2220" 
 data$ID[ data$ID =="H-03"]  <- "2222" 
 data$ID[ data$ID =="H-04"]  <- "2223"
 data$ID[ data$ID =="I-20"]  <- "2362"
 data$ID[ data$ID =="419"]  <- "2385"
 data$ID[ data$ID =="413"]  <- "2386"
 data$ID[ data$ID =="W212"]  <- "2388"
 data$ID[ data$ID =="318"]  <- "2390"
 data$ID[data$ID == "54C"]  <- "2453"
 data$ID[data$ID == "65B"]  <- "2528"
 data$ID[data$ID == "51C"]  <- "2449"
 data$ID[data$ID == "71B"]  <- "2505"
 data$ID[data$ID == "64D"]  <- "2518"
 data$ID[data$ID == "53B"]  <- "2536"
 data$ID[data$ID == "60A"]  <- "2559"
 data$ID[data$ID == "H-01"]  <- "2220"
 data$ID[data$ID == "H-03"]  <- "2222"
 data$ID[data$ID == "419"]  <- "2385"
 data$ID[data$ID == "78-02"]  <- "2140"

 
 data$ID <- as.numeric(data$ID)
  data <- data %>%
         mutate(polar= replace(polar, polar == 0, 0.6)) %>%  #replacing 0 cm growth for the size of plug that is 0.6
    mutate(ecuatorial= replace(ecuatorial, ecuatorial == 0, 0.6)) %>% #replacing 0 cm growth for the size of plug that is 0.6
    # group_by(ID, experimental_replicate, dose, ecuatorial, polar, repeats) %>%
    mutate(growth = ((ecuatorial + polar) / 2)) %>% 
  group_by(ID, dose) %>%
  mutate(growth_range = list(get_range(growth))) %>%
  unnest(cols = c(growth_range)) %>% 
  filter(growth <= upper & growth >= lower) %>%
  dplyr::rename(response = growth) %>%
  ungroup() %>%
  select(c(ID, experimental_replicate, repeats, dose, response))
    }


##Function.Reading data for Discriminatory concentration
reading_data_discriminatory_concentration <- function(filename) {
  data <- filename
  # data <- subset(data, select = -X)
  #data$repeats <- rep_len(1:3, length.out = nrow(data))
  # renaming ID for the current names
  data$ID[data$ID == "12B"]  <- "12"
  #Although 74SS1 belong to # 74, the one what I used was specifically this:74SS1
  data$ID[data$ID == "74SS1"]  <- "74"
  data$ID[data$ID == "62-02"]  <- "2098"
  data$ID[data$ID == "62-03"]  <- "2099"
  data$ID[data$ID == "62-04"]  <- "2100"
  data$ID[data$ID == "78-01"]  <- "2139"
  data$ID[data$ID == "78-02"]  <- "2140"
  data$ID[data$ID == "78-05"]  <- "2143"
  data$ID[data$ID == "S-01"]  <- "2320"
  data$ID[data$ID == "H-01"]  <- "2220"
  data$ID[data$ID == "H-03"]  <- "2222"
  data$ID[data$ID == "H-04"]  <- "2223"
  data$ID[data$ID == "I-20"]  <- "2362"
  data$ID[data$ID == "419"]  <- "2385"
  data$ID[data$ID == "413"]  <- "2386"
  data$ID[data$ID == "W212"]  <- "2388"
  data$ID[data$ID == "318"]  <- "2390"
  data$ID[data$ID == "54C"]  <- "2453"
  data$ID[data$ID == "65B"]  <- "2528"
  data$ID[data$ID == "51C"]  <- "2449"
  data$ID[data$ID == "71B"]  <- "2505"
  data$ID[data$ID == "64D"]  <- "2518"
  data$ID[data$ID == "53B"]  <- "2536"
  data$ID[data$ID == "60A"]  <- "2559"
  data$ID[data$ID == "H-01"]  <- "2220"
  data$ID[data$ID == "H-03"]  <- "2222"
  data$ID[data$ID == "419"]  <- "2385"
  data$ID[data$ID == "78-02"]  <- "2140"
  data$ID <- as.numeric(data$ID)
  data <- data %>%
    mutate(response = replace(response, response == 0, 0.6)) %>%
    #replacing 0 cm growth for the size of plug that is 0.6
    # group_by(ID, experimental_replicate, concentration, growth, repeats) %>%
    group_by(ID) %>%
    mutate(response_range = list(get_range(response))) %>%
    unnest() %>%
    mutate(control_range = list(get_range(control))) %>%
    unnest() %>%
    filter(response <= upper &
             response >= lower,
           control <= upper1 & control >= lower1) %>%
    ungroup() 
}

#getting EC50
getting_EC50 <- function(filename){
  
  getting.EC50 <- EC_table(filename, form = response ~ dose)
getting.EC50 <- getting.EC50 %>% 
  dplyr::rename(ID = sample ) %>% # renaming
 # mutate(ID = as.factor(ID)) %>% 
  dplyr::rename(EC50 = Estimate.50) %>%
   group_by(ID, EC50) %>%
  mutate(source = ifelse(
    ID %in% baseline.isolates.2,# if they are in this group
    "baseline",
    ifelse(
      ID %in% survey.isolates, # if they are in this group, that is named farmer fiels also ahead
      "survey_isolates",
      "treatmentyear2016_isolates" # if they are in this group, that is named fungicide field trials also ahead
      )), Host = ifelse(
    ID %in% drybean,
    "Drybean", ifelse(
      ID %in% soybean,
      "Soybean","different_host"))) %>% ungroup() %>% 
  mutate(
  source = as.factor(source),
  Host = as.factor(Host)) %>%
#taking out "different host" becuse does not belong to drybean nor soybean
  filter(!ID == 8, !ID == 129) %>% 
      mutate(source = recode(source, survey_isolates = "Farmer fields", treatmentyear2016_isolates = "Fungicide field trials", baseline = "Baseline")) %>% mutate(source = as.factor(source)) %>% 
  ungroup()
}
##
getting_DC <- function(filename){
   data <-filename
  
  data1<- data %>% group_by(ID, dose) %>% summarise(mean_response = mean(response, na.rm = TRUE)) %>%
    ungroup() %>% 
  spread(dose, mean_response) %>% 
    dplyr::rename(control = `0`) %>% #taking out "different host" becuse does not belong to drybean nor soybean
  filter(!ID == 8, !ID == 129) 
    
    
  data2 <- as_data_frame( lapply(data1,function(x) {x[1] ; x}))
    colnames(data2)[-c(1:2)] <- paste("RG_", colnames(data2[,-c(1:2)]), sep = "")
  data3 <- data2 %>% mutate_each(funs((./control)*100), starts_with("RG_")) %>% 
    select(-c(control)) 
} 
##
  getting_DC_2 <- function(filename){
   data <-filename
  
data1 <-data %>% #joining
  # normalizing data
  mutate(logEC50 = (log(EC50))) #%>% 
  }

  lineal_model <- function (vardep, varindep, data) {
    lm(paste(vardep, "~", varindep), data = data)
  }

  wrangling_DC<- function(filename){
   data <-filename

data1 <-data %>%
 mutate(source = ifelse(
    ID %in% baseline.isolates.2,
    "Baseline",
    ifelse(
      ID %in% Farmer.fields,
      "Producer Fields","Fungicide Field Trials"
    )), Host = ifelse(
    ID %in% drybean,
    "Drybean", ifelse(
      ID %in% soybean,
      "Soybean","different_host")))  %>%
  mutate(
  source = as.factor(source),
  Host = as.factor(Host)) %>%
#taking out "different host" becuse does not belong to drybean nor soybean
  filter(!ID == 8, !ID == 129) %>%
  mutate(
    country = ifelse(
      ID %in% mexican.isolates,
      "Mexico",
      ifelse(ID %in% brazilian.isolates,
             "Brazil", "USA")))%>% 
 mutate(Field_Year= case_when(
ID %in% NA_year_1973~"NA_year_1973",
ID %in% NA_year_1974~"NA_year_1974",
ID %in% NA_year_1963~"NA_year_1963",
ID %in% NA_year_1975~"NA_year_1975",
ID %in% NA_year_1977~"NA_year_1977",
ID %in% NA_year_1976~"NA_year_1976",
ID %in% NA_year_1978~"NA_year_1978",
ID %in% NA_year_1979~"NA_year_1979",
ID %in% NA_year_1980~"NA_year_1980",
ID %in% NA_year_2003~"NA_year_2003",
ID %in% NA_year_1982~"NA_year_1982",
ID %in% NA_year_1983~"NA_year_1983",
ID %in% NA_year_1984~"NA_year_1984",
ID %in% NA_year_1985~"NA_year_1985",
ID %in% NA_year_1987~"NA_year_1987",
ID %in% NA_year_1990~"NA_year_1990",
ID %in% NA_year_1965~"NA_year_1965",
ID %in% NA_year_1989~"NA_year_1989",
ID %in% NA_year_1991~"NA_year_1991",
ID %in% NA_year_1992~"NA_year_1992",
ID %in% NA_year_1969~"NA_year_1969",
ID %in% NA_year_1981~"NA_year_1981",
ID %in% NA_year_1993~"NA_year_1993",
ID %in% NA_year_1994~"NA_year_1994",
ID %in% NA_year_1995~"NA_year_1995",
ID %in% NA_year_1996~"NA_year_1996",
ID %in% NA_year_1997~"NA_year_1997",
ID %in% WM_production_31_year_1997~"WM_production_31_year_1997",
ID %in% WM_production_32_year_1997~"WM_production_32_year_1997",
ID %in% NA_year_1988~"NA_year_1988",
ID %in% WM_production_33_year_1997~"WM_production_33_year_1997",
ID %in% NA_year_NA~"NA_year_NA",
ID %in% NA_year_1998~"NA_year_1998",
ID %in% NA_year_1999~"NA_year_1999",
ID %in% NA_year_2000~"NA_year_2000",
ID %in% NA_year_2002~"NA_year_2002",
ID %in% WM_nursery_1_year_2003~"WM_nursery_1_year_2003",
ID %in% WM_nursery_2_year_2003~"WM_nursery_2_year_2003",
ID %in% WM_nursery_3_year_2003~"WM_nursery_3_year_2003",
ID %in% WM_nursery_4_year_2003~"WM_nursery_4_year_2003",
ID %in% WM_nursery_5_year_2003~"WM_nursery_5_year_2003",
ID %in% WM_production_35_year_2003~"WM_production_35_year_2003",
ID %in% NA_year_2004~"NA_year_2004",
ID %in% NA_year_2005~"NA_year_2005",
ID %in% NA_year_2006~"NA_year_2006",
ID %in% WM_production_36_year_2004~"WM_production_36_year_2004",
ID %in% WM_production_37_year_2004~"WM_production_37_year_2004",
ID %in% WM_production_38_year_2004~"WM_production_38_year_2004",
ID %in% WM_production_39_year_2004~"WM_production_39_year_2004",
ID %in% WM_production_40_year_2004~"WM_production_40_year_2004",
ID %in% WM_production_41_year_2004~"WM_production_41_year_2004",
ID %in% WM_production_42_year_2004~"WM_production_42_year_2004",
ID %in% WM_production_43_year_2004~"WM_production_43_year_2004",
ID %in% WM_production_44_year_2004~"WM_production_44_year_2004",
ID %in% WM_nursery_6_year_2004~"WM_nursery_6_year_2004",
ID %in% WM_nursery_7_year_2004~"WM_nursery_7_year_2004",
ID %in% WM_nursery_8_year_2004~"WM_nursery_8_year_2004",
ID %in% WM_nursery_9_year_2004~"WM_nursery_9_year_2004",
ID %in% WM_nursery_10_year_2004~"WM_nursery_10_year_2004",
ID %in% WM_nursery_11_year_2004~"WM_nursery_11_year_2004",
ID %in% WM_nursery_12_year_2004~"WM_nursery_12_year_2004",
ID %in% WM_nursery_24_year_2005~"WM_nursery_24_year_2005",
ID %in% WM_nursery_13_year_2005~"WM_nursery_13_year_2005",
ID %in% WM_nursery_14_year_2005~"WM_nursery_14_year_2005",
ID %in% WM_nursery_15_year_2005~"WM_nursery_15_year_2005",
ID %in% WM_nursery_16_year_2005~"WM_nursery_16_year_2005",
ID %in% WM_nursery_17_year_2005~"WM_nursery_17_year_2005",
ID %in% WM_nursery_18_year_2005~"WM_nursery_18_year_2005",
ID %in% WM_nursery_19_year_2006~"WM_nursery_19_year_2006",
ID %in% NA_year_2007~"NA_year_2007",
ID %in% WM_production_1_year_2007~"WM_production_1_year_2007",
ID %in% WM_production_2_year_2007~"WM_production_2_year_2007",
ID %in% WM_production_3_year_2007~"WM_production_3_year_2007",
ID %in% WM_production_4_year_2007~"WM_production_4_year_2007",
ID %in% WM_production_5_year_2007~"WM_production_5_year_2007",
ID %in% WM_production_6_year_2007~"WM_production_6_year_2007",
ID %in% WM_production_7_year_2007~"WM_production_7_year_2007",
ID %in% WM_production_8_year_2007~"WM_production_8_year_2007",
ID %in% WM_production_9_year_2007~"WM_production_9_year_2007",
ID %in% NA_year_2008~"NA_year_2008",
ID %in% WM_production_10_year_2008~"WM_production_10_year_2008",
ID %in% WM_nursery_20_year_2008~"WM_nursery_20_year_2008",
ID %in% WM_production_11_year_2008~"WM_production_11_year_2008",
ID %in% WM_production_12_year_2008~"WM_production_12_year_2008",
ID %in% WM_production_13_year_2008~"WM_production_13_year_2008",
ID %in% NA_year_2009~"NA_year_2009",
ID %in% WM_nursery_21_year_2008~"WM_nursery_21_year_2008",
ID %in% WM_production_38_year_2009~"WM_production_38_year_2009",
ID %in% WM_nursery_22_year_2009~"WM_nursery_22_year_2009",
ID %in% WM_nursery_23_year_2009~"WM_nursery_23_year_2009",
ID %in% WM_production_14_year_2009~"WM_production_14_year_2009",
ID %in% WM_production_15_year_2009~"WM_production_15_year_2009",
ID %in% WM_production_16_year_2009~"WM_production_16_year_2009",
ID %in% WM_production_17_year_2009~"WM_production_17_year_2009",
ID %in% WM_production_18_year_2009~"WM_production_18_year_2009",
ID %in% WM_production_19_year_2009~"WM_production_19_year_2009",
ID %in% WM_production_20_year_2009~"WM_production_20_year_2009",
ID %in% WM_production_21_year_2010~"WM_production_21_year_2010",
ID %in% WM_production_22_year_2010~"WM_production_22_year_2010",
ID %in% WM_production_23_year_2010~"WM_production_23_year_2010",
ID %in% WM_production_24_year_2010~"WM_production_24_year_2010",
ID %in% WM_production_25_year_2010~"WM_production_25_year_2010",
ID %in% WM_production_26_year_2010~"WM_production_26_year_2010",
ID %in% WM_production_27_year_2010~"WM_production_27_year_2010",
ID %in% WM_production_28_year_2010~"WM_production_28_year_2010",
ID %in% WM_production_29_year_2010~"WM_production_29_year_2010",
ID %in% WM_production_30_year_2010~"WM_production_30_year_2010",
ID %in% NA_year_2011~"NA_year_2011",
ID %in% NA_year_2012~"NA_year_2012",
ID %in% NA_year_2013~"NA_year_2013",
ID %in% NA_year_2014~"NA_year_2014",
ID %in% NA_year_2016~"NA_year_2014",
ID %in% po_year_2016~"NA_year_2016",
ID %in% flo_year_2016~"flo_year_2016",
ID %in% lan_year_2016~"lan_year_2016",
ID %in% ho_year_2017~"ho_year_2017",
ID %in% an_year_2017~"an_year_2017",
ID %in% dod_year_2017~"dod_year_2017",
ID %in% va_year_2017~"va_year_2017",
ID %in% na_year_2017~"na_year_2017",
ID %in% sa_year_2017~"sa_year_2017",
ID %in% han_year_2017~"han_year_2017",
ID %in% cu_year_2017~"cu_year_2017",
ID %in% mon_year_2017~"mon_year_2017",
ID %in% ant_year_2017~"ant_year_2017",
ID %in% holt_year_2017~"holt_year_2017",
ID %in% arr_year_2017~"arr_year_2017",
ID %in% ini_year_2017~"ini_year_2017",
ID %in% ley_year_2017~"ley_year_2017",
ID %in% pre_year_2017~"pre_year_2017",
ID %in% ru_year_2017~"ru_year_2017",
ID %in% ca_year_2017~"ca_year_2017",
ID %in% st_year_2017~"st_year_2017",
ID %in% Al_year_2015~"Al_year_2015",
ID %in% Hi_year_2015~"Hi_year_2015",
ID %in% In_year_2015~"In_year_2015",
ID %in% mont_year_2015~"mont_year_2015",
ID %in% Sand_year_2015~"Sand_year_2015",
ID %in% Sani_year_2015~"Sani_year_2015",
ID %in% St_year_2015~"St_year_2015",
ID %in% wau_year_2015~"wau_year_2015",
ID %in% la_year_2015~"la_year_2015",
ID %in% co_year_2015~"co_year_2015",
ID %in% cern_year_2015~"cern_year_2015",
ID %in% cer_year_2015~"cer_year_2015",
ID %in% NA_year_2010~"NA_year_2010",
ID %in% NA_year_2017~"NA_year_2017",
TRUE ~ "no specific field"),
  State= case_when(
ID %in% NE~"NE",
ID %in% CO~"CO",
ID %in% IDA~"ID",
ID %in% CA~"CA",
ID %in% NY~"NY",
ID %in% NA~"NA",
ID %in% WA~"WA",
ID %in% NC~"NC",
ID %in% AZ~"AZ",
ID %in% MT~"MT",
ID %in% WI~"WI",
ID %in% MI~"MI",
ID %in% ON~"ON",
ID %in% AB~"AB",
ID %in% KS~"KS",
ID %in% ND~"ND",
ID %in% MD~"MD",
ID %in% Arusha~"Arusha",
ID %in% SK~"SK",
ID %in% OK~"OK",
ID %in% DF~"DF",
ID %in% OH~"OH",
ID %in% SIN~"SIN",
ID %in% IA~"IA",
ID %in% PA~"PA",
ID %in% Herts~"Herts",
ID %in% Unknown~"Unknown",
ID %in% FL~"FL",
ID %in% GA~"GA",
ID %in% MS~"MS",
ID %in% LA~"LA",
ID %in% KY~"KY",
ID %in% OR~"OR",
ID %in% MG~"MG",
ID %in% GO~"GO",
ID %in% Herverlee~"Herverlee",
# ID %in% Mooi Rive~"Mooi Rive",
ID %in% Marathon~"Marathon",
ID %in% SC~"SC",
ID %in% MN~"MN",
ID %in% Blenheim~"Blenheim",
ID %in% Barranquita~"Barranquita",
ID %in% Victoria~"Victoria",
ID %in% Tasmania~"Tasmania",
ID %in% DE~"DE",
ID %in% Guerbigny~"Guerbigny",
ID %in% QUE~"QUE",
ID %in% Kigali~"Kigali",
ID %in% birdseed~"birdseed",
ID %in% Kessenich~"Kessenich",
ID %in% Elen~"Elen",
ID %in% Moshi~"Moshi",
ID %in% Peronne~"Peronne",
ID %in% Loudéac~"Loudéac",
# ID %in% Bordeaux-Hourtin~"Bordeaux-Hourtin",
ID %in% BA~"BA",
ID %in% RS~"RS",
ID %in% MO~"MO",
ID %in% PR~"PR",
TRUE ~ "no specific state"
), 
County= case_when(
ID %in% Lambari~"Lambari",
ID %in% `Porto Firme`~"Porto Firme",
ID %in% `Unaí`~"Unai",
TRUE ~ "no specific county"
)) %>% 
  mutate(
    Field_Year =  ifelse(
      County == "Lambari",
      "Lambari",
      ifelse(
        County == "Porto Firme",
        "Porto Firme",
        ifelse(County == "Unai", "Unai", Field_Year
        )))) %>% mutate(
Field_Year = recode(
      Field_Year,
      WM_production_21_year_2010 = "NE.10.fld21.DB",
      WM_production_1_year_2007 = "ND.07.fld01.DB",
      WM_production_9_year_2007 = "WA.07.fld09.DB",
      lan_year_2016 = "MI.16.East.Lansing.SB",
      WM_production_10_year_2008 = "WA.08.fld10.DB",
      WM_production_25_year_2010 = "ND.10.fld25.DB",
      mont_year_2015= "MI.15.Montcalm.SB",
      na_year_2017 = "IA.17.Nashua.SB",
      mon_year_2017 = "MI.17.Montcalm.SB",
      Baseline = "Baseline",
      han_year_2017 = "WI.17.Hancock.SB",
      NA_year_2010 = "BR.2010.fld01.DB",
      pre_year_2017 = "MX.2017.pre.DB",
      ca_year_2017 = "MX.2017.ca.DB",
      ley_year_2017 = "MX.2017.ley.DB",
      cer_year_2015 = "WI.2015.cer.SB",
      `Porto Firme` = "BR.2010.por.DB",
      `Unai` = "BR.2010.unai.DB",
      `Lambari` = "BR.2010.lam.DB"
        )) 
  
  }

  #NEED TO MODIFI HERE WHEN BRAZIL, MAKE BT COUNTY NOR EVEN STATE
##wrangling survey to get the EC50DC
  wrangling_survey <- function(filename) {
    data <- filename
    data1 <- data %>%
      group_by(ID) %>% 
      summarise(
        mean_response = mean(response, na.rm = TRUE),
        mean_control = mean(control, na.rm = TRUE)
      ) %>% ungroup() %>%
      mutate(RG = (mean_response / mean_control) * 100)
  }

   wrangling_DC_2<- function(filename){
   data <-filename
   data1 <-data %>%
    group_by(ID) %>% #keeping the highest value of each ID value
  top_n(1, EC50DC) %>% 
  ungroup()
   }
  
  wrangling_DC_3<- function(filename){
   data <-filename
   data1 <-data %>%
   mutate(Field_Year =  ifelse(source == "Baseline", "Baseline", Field_Year),
    Field_Year = as.factor(Field_Year)) %>% 
   
  group_by(Field_Year) %>% filter(n()>=9) %>% 
     ungroup() %>% 
     filter(!Field_Year =="no specific field" ) %>% 
     mutate(Field_Year = fct_reorder(Field_Year, EC50DC, mean))
}

   #test source
  
  test_kruskal_source <- function(filename, vardep, varindep) {
    data <- filename
    vardep <- data$EC50DC
    varindep <- data$source
    data1 <-  data %>%
      do(tidy(shapiro.test(vardep)))
    if (data1$p.value > 0.05)  {
      print ("normal")# I have to add the ANOVA
    } else {
      kruskal.test(as.numeric(vardep) ~ as.factor(varindep), data = data1)
         }}
  
  test_Dunn_source <- function(filename, vardep, varindep) {
    data <- filename
    vardep <- data$EC50DC
    varindep <- data$source
    data1 <-  data %>%
      do(tidy(shapiro.test(vardep)))
    if (data1$p.value > 0.05)  {
      print ("normal")# I have to add the ANOVA
    } else {
      kruskal.test(as.numeric(vardep) ~ as.factor(varindep), data = data1)
     if (data1$p.value > 0.05) {
        print ("No stats difference")
      } else{
        DunnTest(as.numeric(vardep) ~ as.factor(varindep), data = data1, method="bonferroni")}}}
  
   #test host

   test_kruskal_host <- function(filename, vardep, varindep) {
    data <- filename
    vardep <- data$EC50DC
    varindep <- data$Host
    data1 <-  data %>%
      do(tidy(shapiro.test(vardep))) 
      if (data1$p.value > 0.05)  {
      print ("normal")# I have to add the ANOVA
    } else {
      kruskal.test(as.numeric(vardep) ~ as.factor(varindep), data = data1)
      }}    
  
  test_Dunn_host <- function(filename, vardep, varindep) {
    data <- filename
    vardep <- data$EC50DC
    varindep <- data$Host
    data1 <-  data %>%
      do(tidy(shapiro.test(vardep)))
    if (data1$p.value > 0.05)  {
      print ("normal")# I have to add the ANOVA
    } else {
      kruskal.test(as.numeric(vardep) ~ as.factor(varindep), data = data1)
     if (data1$p.value > 0.05) {
        print ("No stats difference")
      } else{
        DunnTest(as.numeric(vardep) ~ as.factor(varindep), data = data1, method="bonferroni")}}}
  
  test_field <- function(filename, vardep, varindep) {
    data <- filename
    vardep <- data$EC50DC
    varindep <- data$Field_Year
    data1 <-  data %>%
      do(tidy(shapiro.test(vardep)))
    if (data1$p.value > 0.05)  {
      print ("normal")# I have to add the ANOVA
    } else {
      kruskal.test(as.numeric(vardep) ~ as.factor(varindep), data = data1)
     if (data1$p.value > 0.05) {
        print ("No stats difference")
     } 
      
      else{
      DunnTest(as.numeric(vardep) ~ as.factor(varindep), data = data1, method="bonferroni")}
    
    
    }}
  ###
      myplot_model_1_boscalid <- function(filename, RG, logEC50, lm)  { 
    data <- filename
    RG <- data$RG
    logEC50 <- data$logEC50
    Host <- data$Host
    source <- data$source
ggplot(data = filename,  aes(x = RG, y = logEC50))  +
  geom_point(aes(
    shape = Host,
    fill= source
  ),
  size = 2, stroke= 1, colour= "black")  + scale_shape_manual(values = c(21, 24))+ scale_fill_grey(start = 0, end = 1) +
      geom_smooth(
        method = "lm",
        se = FALSE,
        colour = "black",
        size = 0.3
      ) +
      theme(
        plot.title = element_text(
          size = 18,
          face = "bold",
          hjust = 0.5,
          family = "Arial"
        ),
        axis.title = element_text(
          size = 18,
          face = "bold",
          hjust = 0.5
        ),
        axis.text = element_text(
          face = "bold",
          size = 18,
          family = "Arial"
        ),
        panel.background = element_rect(fill = "white", colour = "grey50")
      ) + labs(x =
                 "Relative growth (%)", y = expression(bold(Log  (EC[bold("50")])))) + geom_abline(
                   intercept = 0,
                   slope = 1,
                   color = "black",
                   lty = "dashed",
                   size = 0.3
                  ) + expand_limits(x = c(15, 65)) + scale_x_continuous(name = waiver(), breaks = c(15, 25, 35, 45, 55, 65))  + geom_label(aes(x = 53,
                                                                                                                                                                                        y = -2.75)
                                                                                                                                                                                    ,
                                                                                                                                                                                    label = c(paste  (
                                                                                                                                                                                      " Y =",
                                                                                                                                                                                      paste0 (round(summary(lm)[[4]][2], 4), "x"),
                                                                                                                                                                                      "-",
                                                                                                                                                                                      round(summary(lm)[[4]][1], 3) *
                                                                                                                                                                                        -1,
                                                                                                                                                                                      paste  ("\n R =",
                                                                                                                                                                                              paste0 (round(
                                                                                                                                                                                                summary(lm)[[9]][1], 4
                                                                                                                                                                                              ), "\n p < 0.001"))
                                                                                                                                                                                    )),size = 3.5,
                                                                                                                                                                                    fontface = "bold")

 }
    ### tetraconazole
  
   myplot_model_1_tetraconazole <- function(filename, RG, logEC50, lm)  { 
    data <- filename
    RG <- data$RG
    logEC50 <- data$logEC50
    Host <- data$Host
    source <- data$source
ggplot(data = filename,  aes(x = RG, y = logEC50))  +
  geom_point(aes(
    shape = Host,
    fill= source
  ),
  size = 2, stroke= 1, colour= "black")  + scale_shape_manual(values = c(21, 24))+ scale_fill_grey(start = 0, end = 1) +
      geom_smooth(
        method = "lm",
        se = FALSE,
        colour = "black",
        size = 0.3
      ) +
      theme(
        plot.title = element_text(
          size = 18,
          face = "bold",
          hjust = 0.5,
          family = "Arial"
        ),
        axis.title = element_text(
          size = 18,
          face = "bold",
          hjust = 0.5
        ),
        axis.text = element_text(
          face = "bold",
          size = 18,
          family = "Arial"
        ),
        panel.background = element_rect(fill = "white", colour = "grey50")
      ) + labs(x =
                 "Relative growth (%)", y = expression(bold(Log  (EC[bold("50")])))) + geom_abline(
                   intercept = 0,
                   slope = 1,
                   color = "black",
                   lty = "dashed",
                   size = 0.3
                  ) + expand_limits(x = c(15, 65)) + scale_x_continuous(name = waiver(), breaks = c(15, 25, 35, 45, 55, 65))  + geom_label(aes(x = 53,
                                                                                                                                                                                        y = -0.75)
                                                                                                                                                                                    ,
                                                                                                                                                                                    label = c(paste  (
                                                                                                                                                                                      " Y =",
                                                                                                                                                                                      paste0 (round(summary(lm)[[4]][2], 4), "x"),
                                                                                                                                                                                      "-",
                                                                                                                                                                                      round(summary(lm)[[4]][1], 3) *
                                                                                                                                                                                        -1,
                                                                                                                                                                                      paste  ("\n R =",
                                                                                                                                                                                              paste0 (round(
                                                                                                                                                                                                summary(lm)[[9]][1], 4
                                                                                                                                                                                              ), "\n p < 0.001"))
                                                                                                                                                                                    )),size = 3.5,
                                                                                                                                                                                    fontface = "bold")

 }
  
   #### picoxystrobin
  
   myplot_model_1_picoxystrobin <- function(filename, RG, logEC50, lm)  { 
    data <- filename
    RG <- data$RG
    logEC50 <- data$logEC50
    Host <- data$Host
    source <- data$source
ggplot(data = filename,  aes(x = RG, y = logEC50))  +
  geom_point(aes(
    shape = Host,
    fill= source
  ),
  size = 2, stroke= 1, colour= "black")  + scale_shape_manual(values = c(21, 24))+ scale_fill_grey(start = 0, end = 1) +
      geom_smooth(
        method = "lm",
        se = FALSE,
        colour = "black",
        size = 0.3
      ) +
      theme(
        plot.title = element_text(
          size = 18,
          face = "bold",
          hjust = 0.5,
          family = "Arial"
        ),
        axis.title = element_text(
          size = 18,
          face = "bold",
          hjust = 0.5
        ),
        axis.text = element_text(
          face = "bold",
          size = 18,
          family = "Arial"
        ),
        axis.text.y = element_text(
                                                                                                                                                                                                                                angle = 20,
                                                                                                                                       hjust = 1),
        panel.background = element_rect(fill = "white", colour = "grey50")
      ) + labs(x =
                 "Relative growth (%)", y = expression(bold(Log  (EC[bold("50")])))) + geom_abline(
                   intercept = 0,
                   slope = 1,
                   color = "black",
                   lty = "dashed",
                   size = 0.3
                  ) + expand_limits(x = c(15, 65)) + scale_x_continuous(name = waiver(), breaks = c(15, 25, 35, 45, 55, 65))  + geom_label(aes(x = 28,
                                                                                                                                                                                        y = -4.30)
                                                                                                                                                                                    ,
                                                                                                                                                                                    label = c(paste  (
                                                                                                                                                                                      " Y =",
                                                                                                                                                                                      paste0 (round(summary(lm)[[4]][2], 4), "x"),
                                                                                                                                                                                      "-",
                                                                                                                                                                                      round(summary(lm)[[4]][1], 3) *
                                                                                                                                                                                        -1,
                                                                                                                                                                                      paste  ("\n R =",
                                                                                                                                                                                              paste0 (round(
                                                                                                                                                                                                summary(lm)[[9]][1], 4
                                                                                                                                                                                              ), "\n p < 0.001"))
                                                                                                                                                                                    )),size = 3.5,
                                                                                                                                                                                    fontface = "bold")

 }
  
  
    ####
   myplot_model_2 <- function(filename, EC50, Estimate.50DC, Host, source) {
     data <- filename
    EC50 <- data$EC50
    Estimate.50DC <- data$Estimate.50DC
    Host <- data$Host
    source <- data$source
    ggplot(data = filename,  aes(x = EC50, y = Estimate.50DC))  +
  geom_point(aes(
    shape = Host,
    fill= source
  ),
  size = 2, stroke= 1, colour= "black")  + scale_shape_manual(values = c(21, 24))+ scale_fill_grey(start = 0, end = 1)  +
   geom_smooth(method = "lm",
               se = FALSE,
               colour = "black",
               size = 0.3) +
   theme(
     plot.title = element_text(
       size = 18,
       face = "bold",
       hjust = 0.5,
       family = "Arial"
     ),
     axis.title = element_text(size = 18, face = "bold", hjust = 0.5),
     axis.text = element_text(
       face = "bold",
       size = 12,
       family = "Arial"
     ),
     panel.background = element_rect(fill = "white", colour = "grey50")
   ) +
   labs(x =
          expression(bold(EC[bold("50")]) ~ (bold(mg/L ~ bold(
            "a.i."
          )))), y = expression(bold(EC[bold("50") ~ (bold("D"))]) ~ (bold(mg/L ~ bold(
            "a.i."
          ))))) + geom_abline(
            intercept = 0,
            slope = 1,
            color = "black",
            lty = "dashed",
            size = 0.3)
   }
#label 
  generate_label <- function(HSD){
    # Extract labels and factor levels from Tukey post-hoc 
    Dunn.levels <- HSD[[1]][,2]
    Dunn.labels <- multcompLetters(Dunn.levels)['Letters']
    plot.labels <- names(Dunn.labels[['Letters']])
    plot.levels <- data.frame(plot.labels, Dunn.levels, labels = Dunn.labels[['Letters']],
                              stringsAsFactors = FALSE)
        return(plot.levels)
}
   ###
generate_label_2 <- function(HSD, flev, surveyy, highlabel){
  # Extract labels and factor levels from Tukey post-hoc 
  Dunn.levels <- HSD[[1]][,2]
  Dunn.labels <- multcompLetters(Dunn.levels)['Letters']
  plot.labels <- names(Dunn.labels[['Letters']])
  # Get highest quantile for Dunn's 5 number summary and add a bit of space to buffer between    
  # upper quantile and label placement
  boxplot.df <- plyr:: ddply(surveyy, flev, function (x) max(fivenum(x$EC50DC)) + highlabel)
    # Create a data frame out of the factor levels and Dunn's homogenous group letters
  plot.levels <- data.frame(plot.labels, labels = Dunn.labels[['Letters']],
                            stringsAsFactors = FALSE)
  
  # Merge it with the labels
  labels.df <- merge(plot.levels, boxplot.df, by.x = 'plot.labels', by.y = flev, sort = FALSE)
  
  return(labels.df)
}
###
myplot_comparison <- function(filename, Field_Year,EC50DC, highlabel) {
  data <- filename
  Field_Year <- data$Field_Year
  EC50DC<- data$EC50DC
  
  ggplot(data = filename,  aes(x = Field_Year, y = EC50DC)) + geom_boxplot() +
    geom_text(data = generate_label_2(test_field(filename), 'Field_Year', filename,highlabel),
              aes(x = plot.labels, y = V1, label = labels), size= 3) +     labs(x = "Field", y = expression(bold(EC[bold("50") ~ (bold("D"))]) ~
                                                                                                     (bold(mg/L ~ bold(
                                                                                                       "a.i."
                                                                                                     ))))) + theme(
                                                                                                       panel.border = element_rect(
                                                                                                         colour = "black",
                                                                                                         fill = NA,
                                                                                                         size = 0.5
                                                                                                       ),
                                                                                                       axis.title = element_text(size = 11, face = "bold", hjust = 0.5),
                                                                                                       axis.text.x = element_text(
                                                                                                         face = "bold",
                                                                                                         size = 8,
                                                                                                         family = "Arial",
                                                                                                         angle = 30,
                                                                                                         hjust = 1
                                                                                                       ),
                                                                                                       axis.text.y = element_text(
                                                                                                         face = "bold",
                                                                                                         size = 8,
                                                                                                         family = "Arial"
                                                                                                       ),
                                                                                                       panel.background = element_rect(fill = "white", colour = "grey50")
                                                                                                     ) 
  
  
}

 myplot_baseline <- function(filename, EC50DC, source) {
  data <- filename
  EC50DC <- data$EC50DC
  source <- data$source
  
  data %>% filter(source == "Baseline") %>% mutate(
    Baseline = ifelse(
      ID %in% SI.Screening.Nursery.Field,
      # if they are in this group
      "Contemporary isolates*",
      "Historic isolates**"
    )
  ) %>%   ggplot(aes(x = Baseline, y = EC50DC, )) + geom_violin() + geom_jitter(aes(shape = Baseline), size = 6) +     labs(x = "Baseline", y = expression(bold(EC[bold("50") ~ (bold("D"))]) ~
                                                                                                                                                             (bold(mg/L ~ bold(
                                                                                                                                                               "a.i."
                                                                                                                                                             ))))) + theme(
                                                                                                                                                               panel.border = element_rect(
                                                                                                                                                                 colour = "black",
                                                                                                                                                                 fill = NA,
                                                                                                                                                                 size = 1
                                                                                                                                                               ),
                                                                                                                                                               axis.title = element_text(
                                                                                                                                                                 size = 18,
                                                                                                                                                                 face = "bold",
                                                                                                                                                                 hjust = 0.5
                                                                                                                                                               ),
                                                                                                                                                               axis.text.x = element_text(
                                                                                                                                                                 face = "bold",
                                                                                                                                                                 size = 10,
                                                                                                                                                                 family = "Arial",
                                                                                                                                                                 hjust = 0.5
                                                                                                                                                               ),
                                                                                                                                                               axis.text.y = element_text(
                                                                                                                                                                 face = "bold",
                                                                                                                                                                 size = 10,
                                                                                                                                                                 family = "Arial"
                                                                                                                                                               ),
                                                                                                                                                               panel.background = element_rect(fill = "white", colour = "grey50")
                                                                                                                                                             )
}
```
# Fungicides:
**1. Thiophanate methyl**  
**2. Tetraconazole**  
**3. Boscalid**  
**4. Picoxystrobin**  


```{r Reading Data, echo = TRUE, warning=FALSE}
#TM
hey <-
  reading_data_serial_dilution("data/TM_serialdi_BL-WMN_survey-WMNCSRP(validationdata).csv") %>%
  filter(!experimental_replicate == 1 |
           !dose == 1.5 | !repeats == 3)  #something weird happened with the repeats 3 so Iam revoming those
#Check here because there are more repeats for control, in that case create a new variable or column
## First try with postdoc Thomas
survey.TM.thomas <- read.csv("data/DD_survey-WMNCSRP.csv") %>%
  select(
    ID,
    repeats,
    experimental_replicate,
    TM.ecuatorial,
    TM.polar,
    control.ecuatorial,
    control.polar
  ) %>%
  mutate(TM.growth = ((TM.ecuatorial + TM.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = TM.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()


## Second try with Edgar
survey.TM.USA <- 
      read.csv("data/DD_survey-WMNCSRP_TM.csv") %>% 
mutate(TM.growth = ((TM.ecuatorial + TM.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = TM.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()

## Third try with Cristian
survey.TM.more.USA <- 
 read.csv("data/TM.Cristian.csv") %>% 
mutate(TM.growth = ((TM.ecuatorial + TM.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = TM.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()

#Tetraconazole
hey2 <- reading_data_serial_dilution("data/tetraconazole_serialdi_BL-WMN_survey-WMNCSRP(validationdata).csv")
## First try with postdoc Thomas
survey.tetraconazole.thomas <- read.csv("data/DD_survey-WMNCSRP.csv") %>%
  select(
    ID,
    repeats,
    experimental_replicate,
    tetraconazole.ecuatorial,
    tetraconazole.polar,
    control.ecuatorial,
    control.polar
  ) %>%
  mutate(tetraconazole.growth = ((tetraconazole.ecuatorial + tetraconazole.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = tetraconazole.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()

## Second try with Edgar
survey.tetraconazole.USA <- 
      read.csv("data/DD_survey-WMNCSRP_tetraconazole.csv") %>% 
mutate(tetraconazole.growth = ((tetraconazole.ecuatorial + tetraconazole.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = tetraconazole.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()


## Third try with Cristian
survey.tetraconazole.more.USA <- 
 read.csv("data/tetraconazole.Cristian.csv") %>% 
mutate(tetraconazole.growth = ((tetraconazole.ecuatorial + tetraconazole.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = tetraconazole.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>%  
reading_data_discriminatory_concentration()

## Fourth try with Edgar 
survey.tetraconazole.MX.BZ <- 
 read.csv("data/DD_survey-WMNCSRP_tetraconazole.Mexican.Brazilian.csv") %>% mutate(tetraconazole.growth = ((tetraconazole.ecuatorial + tetraconazole.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = tetraconazole.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()


#Boscalid
hey3 <- reading_data_serial_dilution("data/boscalid_serialdilution_baseline(validationdata).csv") %>% filter(!dose ==  0.4 & !ID == 800)
## First try with postdoc Thomas
survey.boscalid.thomas <- read.csv("data/DD_survey-WMNCSRP.csv") %>%
  select(
    ID,
    repeats,
    experimental_replicate,
    boscalid.ecuatorial,
    boscalid.polar,
    control.ecuatorial,
    control.polar
  ) %>%
  mutate(boscalid.growth = ((boscalid.ecuatorial + boscalid.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = boscalid.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()


   ## Second try with Edgar
survey.boscalid.USA <- 
      read.csv("data/DD_survey-WMNCSRP_boscalid.csv") %>% 
mutate(boscalid.growth = ((boscalid.ecuatorial + boscalid.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = boscalid.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()


 ## Third try with Cristian
survey.boscalid.more.USA <- 
 read.csv("data/Boscalid.Cristian.csv") %>% 
mutate(boscalid.growth = ((boscalid.ecuatorial + boscalid.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = boscalid.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()


## Fourth try with Edgar 
survey.boscalid.MX.BZ <- 
 read.csv("data/DD_survey-WMNCSRP_boscalid.Brazilian.Mexican.csv") %>% mutate(boscalid.growth = ((boscalid.ecuatorial + boscalid.polar) / 2),
         control.growth = ((control.ecuatorial + control.polar) / 2)) %>% 
dplyr::rename(response = boscalid.growth, control = control.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()


# #Picoxystrobin
hey4 <- reading_data_serial_dilution("data/picoxystrobin_serialdi_BL-WMN_survey-WMNCSRP(validationdata).csv")
## First try with postdoc Thomas
survey.picoxystrobin.thomas <- read.csv("data/DD_survey-WMNCSRP.csv") %>%
  select(
    ID,
    repeats,
    experimental_replicate,
    picoxystrobin.ecuatorial,
    picoxystrobin.polar,
     sham.ecuatorial,
    sham.polar
  ) %>%
  mutate(picoxystrobin.growth = ((picoxystrobin.ecuatorial + picoxystrobin.polar) / 2),
         sham.growth = ((sham.ecuatorial + sham.polar) / 2)) %>% 
dplyr::rename(response = picoxystrobin.growth, control = sham.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
reading_data_discriminatory_concentration()
 
## Second try with Edgar
survey.picoxystrobin.USA <- 
      read.csv("data/DD_survey-WMNCSRP_picoxystrobin.csv") %>% 
mutate(picoxystrobin.growth = ((picoxystrobin.ecuatorial + picoxystrobin.polar) / 2),
         sham.growth = ((sham.ecuatorial + sham.polar) / 2)) %>% 
dplyr::rename(response = picoxystrobin.growth, control = sham.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
  reading_data_discriminatory_concentration()


## Third try with Cristian
survey.picoxystrobin.more.USA <- 
 read.csv("data/picoxystrobin.Cristian.csv") %>% 
mutate(picoxystrobin.growth = ((picoxystrobin.ecuatorial + picoxystrobin.polar) / 2),
         sham.growth = ((sham.ecuatorial + sham.polar) / 2)) %>% 
dplyr::rename(response = picoxystrobin.growth, control = sham.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
 reading_data_discriminatory_concentration()

## Fourth try with Edgar 
survey.picoxystrobin.MX.BZ <- 
 read.csv("data/DD_survey-WMNCSRP_picoxystrobin.Mexican.Brazilian.csv") %>% mutate(picoxystrobin.growth = ((picoxystrobin.ecuatorial + picoxystrobin.polar) / 2),
         sham.growth = ((sham.ecuatorial + sham.polar) / 2)) %>% 
dplyr::rename(response = picoxystrobin.growth, control = sham.growth) %>% 
  select(c(ID, experimental_replicate, repeats, response, control)) %>% 
 reading_data_discriminatory_concentration()
```
#  Fungicide sensitivity by fungicide
## Serial dilution  
Fungicide sensitivity was determined from 22 baseline isolates plus 21 additional isolates from FF and FFT randomly selected. Serial dilution  by fungicide eas carried out from 5-6 concentrations with 4 repetitions per isolate and 2 experimental replications. Plates were incubated in darkness at 23 ± 2 °C for 42 h and diameter measured with digital calipers. A dose-response curve was fit to estimate the EC50 for each fungicide.  


```{r EC50 serial dilution, echo=TRUE, warning=FALSE}
# TM
TM.EC50 <-  getting_EC50(hey)
# Tetraconazole
tetraconazole.EC50 <- getting_EC50(hey2)
# Boscalid
boscalid.EC50 <- getting_EC50(hey3)
#Picoxystrobin
picoxystrobin.EC50 <- getting_EC50(hey4)
```


## Procedure to get the Discriminatory Concentration (DC)  
We identified the concentration with the best prediction of EC50 for each fungicide, known as discriminatory concentration (DC), by linear regression of **% mycelial growth vs. log(EC50)**. The concentration yielding the **highest coefficient of determination (r2)** was selected as the DC for each fungicide. For the remaining isolates among the 207, growth on the DC was used to estimate EC50, which is termed EC50(D).

```{r Getting Discriminatory dose, warning=FALSE, echo=TRUE}
#TM
  #getting_EC50 releases ID with factor
TM <- getting_DC(hey) 
  TM.joined <-  TM %>% left_join(TM.EC50)
  
TM.joined.2 <-  getting_DC_2(TM.joined) 
  
#Boscalid
boscalid <- getting_DC(hey3) 
  boscalid.joined <-  boscalid%>% left_join(boscalid.EC50)
  
 boscalid.joined.2 <-  getting_DC_2(boscalid.joined) 
 ###Dose chosen as DC. Although there is one dose with higher r2, there is no difference between and this was selected
# Linear model of log EC50 and relative growth at 0.2 ppm, check normality and homogeneity of variances
finalRG_0.2 <- lm(logEC50 ~ RG_0.2,  boscalid.joined.2)
summary(finalRG_0.2)
#plot(finalRG0.2)

##taking out outliers from Q-Q PLOT $ LEVERAGE TEST, ISOLATES # 1032 from treatment, 2222 & 2223 from Farmer fields##
boscalid.joined.2.clean <-
  boscalid.joined.2 %>% filter(!ID == "1032", !ID == "2222",!ID == "2223") %>% 
  select(ID, logEC50, RG_0.2, EC50) %>% 
  dplyr::rename(RG = RG_0.2) %>% 
   wrangling_DC()

summary(boscalid.lm <- lineal_model("logEC50","RG",data =  boscalid.joined.2.clean))
   # Getting the EC50DC according to the model
boscalid.EC50DC <- boscalid.joined.2.clean %>% mutate(Estimate.50DC =
                                     #using the model the intercept and the coefficient of the model
                                    exp(boscalid.lm[[1]][[1]] + boscalid.lm[[1]][[2]] * RG)) %>%# here I need to do this dynamic
  select(ID, Estimate.50DC, EC50)  %>%  
  wrangling_DC()

# Linear Regression graph of log EC50 and relative growth at 0.2 ppm
###

myplot_model_1_boscalid(boscalid.joined.2.clean, lm = boscalid.lm) 
 
######### Linear model of EC50  and EC50DC just to corroborate is a GOOD RELATIONSHIP, check normality and homogeneity of variances

summary(boscalid.lm.2 <- lineal_model("Estimate.50DC","EC50",data =  boscalid.EC50DC))

round(boscalid.lm.2[[1]][[2]],4)
round(boscalid.lm.2[[1]][[1]],3) 

myplot_model_2(boscalid.EC50DC)+ expand_limits(x = c(0.04, 0.25), y = c(0.04, 0.25))  + scale_x_continuous(name = waiver(),
                                                                             breaks = c(0.05, 0.1, 0.15, 0.20, 0.25)) + scale_y_continuous(breaks = c(0.05, 0.1, 0.15, 0.20, 0.25)) + geom_label(aes(x = 0.20,
                                                    y = 0.10)
                                                ,
                                                label = c(paste  (
                                                  " Y =",
                                                  paste0 (round(summary(boscalid.lm.2)[[4]][2], 4), "x"),
                                                  
                                                  round(summary(boscalid.lm.2)[[4]][1], 3) *
                                                    -1,
                                                  paste  ("\n R =",
                                                          paste0 (round(
                                                           summary(boscalid.lm.2)[[9]][1], 4
                                                          ), "\n p < 0.001"))
                                                )),size = 5,
                                                fontface = "bold")
                                               
                                               
                                               



round(boscalid.lm[[1]][[2]],4)
round(boscalid.lm[[1]][[1]],3) 
#Tetraconazole
  
tetraconazole <- getting_DC(hey2) 
  tetraconazole.joined <-  tetraconazole%>% left_join(tetraconazole.EC50)
  
tetraconazole.joined.2 <-  getting_DC_2(tetraconazole.joined) 
 ###Dose chosen as DC. Although there is one dose with higher r2, there is no difference between and this was selected
# Linear model of log EC50 and relative growth at 2 ppm, check normality and homogeneity of variances

finalRG_2 <- lm(logEC50 ~ RG_2, tetraconazole.joined.2)
summary(finalRG_2)
#plot(finalRG0.2)

##taking out outliers from Q-Q PLOT $ LEVERAGE TEST, ISOLATES # 1032 from treatment, 2222 & 2223 from Farmer fields##
tetraconazole.joined.2.clean <-
  tetraconazole.joined.2 %>% filter(!ID == "21", !ID == "558",!ID == "667") %>% 
   select(ID, logEC50, RG_2, EC50) %>% 
  dplyr::rename(RG = RG_2) %>% 
   wrangling_DC()

summary(tetraconazole.lm <- lineal_model("logEC50","RG",data =  tetraconazole.joined.2.clean))
   # Getting the EC50DC according to the model
tetraconazole.EC50DC <- tetraconazole.joined.2.clean %>% mutate(Estimate.50DC =
                                     #using the model the intercept and the coefficient of the model
                                    exp(tetraconazole.lm[[1]][[1]] + tetraconazole.lm[[1]][[2]] * RG)) %>%# here I need to do this dynamic
  select(ID, Estimate.50DC, EC50)  %>%  wrangling_DC()

# Linear Regression graph of log EC50 and relative growth at 0.2 ppm
###

myplot_model_1_tetraconazole(tetraconazole.joined.2.clean, lm = tetraconazole.lm) 
 
 

######### Linear model of EC50  and EC50DC just to corroborate is a GOOD RELATIONSHIP, check normality and homogeneity of variances

summary(tetraconazole.lm.2 <- lineal_model("Estimate.50DC","EC50",data =  tetraconazole.EC50DC))

round(tetraconazole.lm.2[[1]][[2]],4)
round(tetraconazole.lm.2[[1]][[1]],3)

myplot_model_2(tetraconazole.EC50DC)+ expand_limits(x = c(0.27, 1.8), y = c(0.27, 1.8))  + scale_x_continuous(name = waiver(),
                                                                             breaks = c(0.3, 0.6, 0.9, 1.2, 1.5,1.8)) + scale_y_continuous(breaks = c(0.3, 0.6, 0.9, 1.2, 1.5,1.8)) + geom_label(aes(x = 1.5,
                  y = 0.6)
              ,
              label = c(paste  (
                " Y =",
                paste0 (round(summary(tetraconazole.lm.2)[[4]][2], 4), "x"),
                round(summary(tetraconazole.lm.2)[[4]][1], 3) *
                  -1,
                paste  ("\n R =",
                        paste0 (round(
                          summary(tetraconazole.lm.2)[[9]][1], 4
                        ), "\n p < 0.001"))
              )),size = 5,
              fontface = "bold")
 ## Picoxystrobin

picoxystrobin<- getting_DC(hey4) 
  picoxystrobin.joined <-  picoxystrobin%>% left_join(picoxystrobin.EC50)
  
 picoxystrobin.joined.2 <-  getting_DC_2(picoxystrobin.joined) 
 ###Dose chosen as DC. Although there is one dose with higher r2, there is no difference between and this was selected
# Linear model of log EC50 and relative growth at 0.01 ppm, check normality and homogeneity of variances
finalRG_0.01 <- lm(logEC50 ~ RG_0.01,  picoxystrobin.joined.2)
summary(finalRG_0.01)

##taking out outliers from Q-Q PLOT $ LEVERAGE TEST, ISOLATES # 1032 from treatment, 2222 & 2223 from Farmer fields##
picoxystrobin.joined.2.clean <-
  picoxystrobin.joined.2 %>% filter(!ID == "1032", !ID == "558",!ID == "581") %>% 
  select(ID, logEC50, RG_0.01, EC50) %>% 
  dplyr::rename(RG = RG_0.01) %>% 
   wrangling_DC()


summary(picoxystrobin.lm <- lineal_model("logEC50","RG",data =  picoxystrobin.joined.2.clean))
   # Getting the EC50DC according to the model
picoxystrobin.EC50DC <- picoxystrobin.joined.2.clean %>% mutate(Estimate.50DC =
                                     #using the model the intercept and the coefficient of the model
                                    exp(picoxystrobin.lm[[1]][[1]] + picoxystrobin.lm[[1]][[2]] * RG)) %>%# here I need to do this dynamic
  select(ID, Estimate.50DC, EC50)  %>%  wrangling_DC()

# Linear Regression graph of log EC50 and relative growth at 0.2 ppm
###
summary(picoxystrobin.lm.2 <- lineal_model("Estimate.50DC","EC50",data =  picoxystrobin.EC50DC))

round(picoxystrobin.lm.2[[1]][[2]],4)
round(picoxystrobin.lm.2[[1]][[1]],3)

myplot_model_1_picoxystrobin(picoxystrobin.joined.2.clean, lm = picoxystrobin.lm) 


######### Linear model of EC50  and EC50DC just to corroborate is a GOOD RELATIONSHIP, check normality and homogeneity of variances
myplot_model_2(picoxystrobin.EC50DC)+ expand_limits(x = c(0.005, 0.018), y = c(0.005, 0.018))  + scale_x_continuous(name = waiver(),
                                                                             breaks = c(0.0075, 0.01, 0.0125, 0.0150, 0.018)) + scale_y_continuous(breaks = c(0.0075, 0.01, 0.0125, 0.0150, 0.018)) + geom_label(aes(x = 0.015,
                  y = 0.009)
              ,
              label = c(paste  (
                " Y =",
                paste0 (round(summary(picoxystrobin.lm.2)[[4]][2], 4), "x"),
                round(summary(picoxystrobin.lm.2)[[4]][1], 3) *
                  -1,
                paste  ("\n R =",
                        paste0 (round(
                          summary(picoxystrobin.lm.2)[[9]][1], 4
                        ), "\n p < 0.001"))
              )),size = 5,
              fontface = "bold")
 
#Function to clean
picoxystrobin.joined.2 %>% filter(!EC50 %in% Outlier(EC50))
```



### Reading from White Mold Inventory; needed downstream
```{r White Mold Inventory}
#based on object table1 which has all the complete isolates used in this study,  

 # I have to keep Field adn Year separate tahts why reading again the main table from data
 
 whitemold.inventory.2 <- read_csv("data/whitemold.inventory.csv")
 
 Colletotrichum <-
    whitemold.inventory.2 %>% 
   # filter(ID %in% lasiodiplodia) %>% 
   mutate(country = ifelse(
      ID %in% mexican.isolates,
      "Mexico",
      ifelse(ID %in% brazilian.isolates,
             "Brazil", "USA")),
      source = ifelse(
    ID %in% baseline.isolates.2,
    "Baseline",
    ifelse(
      ID %in% Farmer.fields,
      "Producer Fields","Fungicide Field Trials"
    )), Host = ifelse(
    ID %in% drybean,
    "Drybean", ifelse(
      ID %in% soybean,
      "Soybean","different_host")),
 source = as.factor(source),
  Host = as.factor(Host),
 country = as.factor(country),
 Field = as.factor(Field),
 Year = as.factor(Year),
 State = as.factor(State),
 County = as.factor(County),
 Fungicide_current_season = as.factor(Fungicide_current_season),
 molecule.s._current_season = as.factor(molecule.s._current_season),
 Group_MOA_current_season = as.factor(Group_MOA_current_season),
 Applications = as.factor(Applications),
 Rate = as.factor(Rate)
      ) %>%    select(
     ID,
     source,#we will change later to the real variable "source"
     Host,
     country, 
     Field,
     Year,
     State,
     County,
     Fungicide_current_season,
     molecule.s._current_season,
     Group_MOA_current_season,
     Applications,
     Rate
   ) 
 
 controls.in.the.plots <- Colletotrichum %>% filter( Fungicide_current_season == "control") %>% select(ID) %>% pull()

 controls.in.the.plots.2 <- c(1175, 1195, 1327, 1443, 1446, 1461, 1541, 1710)
```
##  Use of the Discriminatory Concentration DC  
We randomly selected 154 isolates: 86 from dry beans (from 9 states with 1-17 isolates/state, and 2 states from Brazil) and 66 from soybeans (6 NE, 18 IA, 13 WI & 18 MI and 4 from Mexico)  from the collection (baseline, FF, and FF) plus 5 same isolates used in the serial dilution approach


```{r Using the DC and the model}
TM.joined.2 <- TM.joined.2 %>% select(ID, EC50)

survey.TM <- bind_rows(survey.TM.thomas,
survey.TM.USA,
survey.TM.more.USA,
#survey.TM.MX.BZ, there are the same the other fungicides
TM.joined.2) 

#Have to analyze the controls

#wranglind_DC_2

survey.TM.2 <- survey.TM%>% 
  wrangling_DC() %>% 
  mutate(Field_Year =  ifelse(source == "Baseline", "Baseline", Field_Year),
    Field_Year = as.factor(Field_Year)) %>% 
  dplyr::rename(EC50DC = EC50) %>% #change here to be able to join to other datframes to get the total used
  select(ID, EC50DC) %>% 
   # is here where now I have to remove the isolates from controls from fields and potato host isolates provided by Dr.Chilvers at the end
  dplyr::filter(! ID %in% controls.in.the.plots,  ! ID %in% 2570:2574)
  
survey.boscalid.thomas.2 <- wrangling_survey(survey.boscalid.thomas) %>%
  mutate(EC50DC = exp (boscalid.lm [[1]][[1]] + boscalid.lm[[1]][[2]] * RG))# exponential is the opposite of log
        #using the model the intercept and the coefficient of the model

survey.boscalid.USA.2 <- wrangling_survey(survey.boscalid.USA) %>%
  mutate(EC50DC = exp (boscalid.lm [[1]][[1]] + boscalid.lm[[1]][[2]] * RG))#

##
survey.boscalid.more.USA <- survey.boscalid.more.USA %>% mutate(response = (response/10), control = (control/10))




survey.boscalid.more.USA.2 <- wrangling_survey(survey.boscalid.more.USA) %>%
  mutate(EC50DC = exp (boscalid.lm [[1]][[1]] + boscalid.lm[[1]][[2]] * RG))#
##

survey.boscalid.MX.BZ.2 <- wrangling_survey(survey.boscalid.MX.BZ) %>%
  mutate(EC50DC = exp (boscalid.lm [[1]][[1]] + boscalid.lm[[1]][[2]] * RG))#
       
#Gathering 
boscalid.EC50.2 <- boscalid.EC50 %>% select(ID, EC50) %>% 
  dplyr::rename(EC50DC = EC50)#just renaming in purpose to join the following dataframe

survey.boscalid <- bind_rows(survey.boscalid.thomas.2,
                             survey.boscalid.USA.2,
                             survey.boscalid.more.USA.2,
                             survey.boscalid.MX.BZ.2, 
                             boscalid.EC50.2) %>% select(ID,EC50DC)%>% 
   # is here where now I have to remove the isolates from controls from fields and potato host isolates provided by Dr. Chilvers at the end
  dplyr::filter(! ID %in% controls.in.the.plots, ! ID %in% 2570:2574 )

#Have to analyze the controls
#fix the function here
#wranglind_DC_2
#Check this
IDA <- ID
survey.boscalid.complete <- survey.boscalid%>% 
  wrangling_DC() %>%  
 wrangling_DC_2() 

survey.boscalid.fields <- survey.boscalid%>% 
  wrangling_DC() %>% 
   wrangling_DC_2() %>% 
  wrangling_DC_3()
  
###For comparisons between hosts I had to take out the ones coming from baseline otherwise would not be fair

survey.boscalid.complete.host <- survey.boscalid.complete %>% 
  filter(!source=="Baseline")

# Tetraconazole
#any particular case check so far by this way
survey.tetraconazole.thomas.2 <- wrangling_survey(survey.tetraconazole.thomas) %>%
  mutate(EC50DC = exp (tetraconazole.lm [[1]][[1]] + tetraconazole.lm[[1]][[2]] * RG))# exponential is the opposite of log
        #using the model the intercept and the coefficient of the model
##
survey.tetraconazole.USA.2 <- wrangling_survey(survey.tetraconazole.USA) %>%
  mutate(EC50DC = exp (tetraconazole.lm [[1]][[1]] + tetraconazole.lm[[1]][[2]] * RG))#

##
survey.tetraconazole.more.USA <- survey.tetraconazole.more.USA %>% mutate(response = (response/10), control = (control/10))

survey.tetraconazole.more.USA.2 <- wrangling_survey(survey.tetraconazole.more.USA) %>%
  mutate(EC50DC = exp (tetraconazole.lm [[1]][[1]] + tetraconazole.lm[[1]][[2]] * RG))#

survey.tetraconazole.MX.BZ.2 <- wrangling_survey(survey.tetraconazole.MX.BZ) %>%
  mutate(EC50DC = exp (tetraconazole.lm [[1]][[1]] + tetraconazole.lm[[1]][[2]] * RG))#
       
#Gathering 
tetraconazole.EC50.2 <- tetraconazole.EC50 %>% select(ID, EC50) %>% 
  dplyr::rename(EC50DC = EC50)#just rebaming in purpose to join the following dataframe

survey.tetraconazole <- bind_rows(survey.tetraconazole.thomas.2,
                             survey.tetraconazole.USA.2,
                             survey.tetraconazole.more.USA.2,
                             survey.tetraconazole.MX.BZ.2, 
                             tetraconazole.EC50.2) %>% select(ID,  EC50DC) %>% 
# is here where now I have to remove the isolates from controls from fields and potato host isolates provided by Dr. Chilvers at the end
  dplyr::filter(! ID %in% controls.in.the.plots, ! ID %in% 2570:2574 )
#Have to analyze the controls

#wranglind_DC_2

survey.tetraconazole.complete <- survey.tetraconazole%>% 
  wrangling_DC() %>% 
  wrangling_DC_2()

survey.tetraconazole.fields <- survey.tetraconazole%>% 
  wrangling_DC() %>%
  wrangling_DC_2() %>% 
  wrangling_DC_3()
  # count(Field_Year) %>% 
  #  filter(n >= 9)




###For comparisons between hosts I had to take out the ones coming from baseline otherwise would not be fair

survey.tetraconazole.complete.host <- survey.tetraconazole.complete %>% 
  filter(!source=="Baseline")



# picoxystrobin

#any particular case check so far by this way
survey.picoxystrobin.thomas.2 <- wrangling_survey(survey.picoxystrobin.thomas) %>%
  mutate(EC50DC = exp (picoxystrobin.lm [[1]][[1]] + picoxystrobin.lm[[1]][[2]] * RG))# exponential is the opposite of log
        #using the model the intercept and the coefficient of the model
##
# 
survey.picoxystrobin.USA.2 <- wrangling_survey(survey.picoxystrobin.USA) %>%
  mutate(EC50DC = exp (picoxystrobin.lm [[1]][[1]] + picoxystrobin.lm[[1]][[2]] * RG))#

##
survey.picoxystrobin.more.USA <- survey.picoxystrobin.more.USA %>% mutate(response = (response/10), control = (control/10))

survey.picoxystrobin.more.USA.2 <- wrangling_survey(survey.picoxystrobin.more.USA) %>%
  mutate(EC50DC = exp (picoxystrobin.lm [[1]][[1]] + picoxystrobin.lm[[1]][[2]] * RG))#
##

survey.picoxystrobin.MX.BZ.2 <- wrangling_survey(survey.picoxystrobin.MX.BZ) %>%
  mutate(EC50DC = exp (picoxystrobin.lm [[1]][[1]] + picoxystrobin.lm[[1]][[2]] * RG))#
       
#Gathering 
picoxystrobin.EC50.2 <- picoxystrobin.EC50 %>% select(ID, EC50) %>% 
  dplyr::rename(EC50DC = EC50)#just rebaming in purpose to join the following dataframe

survey.picoxystrobin <- bind_rows(survey.picoxystrobin.thomas.2,
                             survey.picoxystrobin.USA.2,
                             survey.picoxystrobin.more.USA.2,
                             survey.picoxystrobin.MX.BZ.2, 
                             picoxystrobin.EC50.2) %>% select(ID,  EC50DC) %>% 
# is here where now I have to remove the isolates from controls from fields and potato host isolates provided by Dr. Chilvers at the end
  dplyr::filter(! ID %in% controls.in.the.plots, ! ID %in% 2570:2574 )
#Have to analyze the controls

survey.picoxystrobin.complete <- survey.picoxystrobin%>% 
  wrangling_DC() %>% 
  wrangling_DC_2()

survey.picoxystrobin.fields <- survey.picoxystrobin%>% 
  wrangling_DC() %>%
  wrangling_DC_2() %>% 
  wrangling_DC_3()
  
###For comparisons between hosts I had to take out the ones coming from baseline otherwise would not be fair

survey.picoxystrobin.complete.host <- survey.picoxystrobin.complete %>% 
  filter(!source=="Baseline")
```
#  Statistic Analysis 
## Boscalid 
### By Source 
```{r Stat Anal Boscalid By Source}
test_kruskal_source(survey.boscalid.complete %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() )

test_Dunn_source(survey.boscalid.complete %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() )

generate_label( test_Dunn_source(survey.boscalid.complete %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() ))

survey.boscalid.complete %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup()  %>% 
  group_by(source) %>% summarize(round(range(EC50DC),3), average= round(mean(EC50DC), 2), N= n(),  sd = sd(EC50DC), se = round(sd/sqrt(N),3))

```

### By Host  

```{r Stat Anal Boscalid By Host}
#generate_label(
  test_kruskal_host(survey.boscalid.complete.host %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() )
 
#generate_label
  test_Dunn_host(survey.boscalid.complete.host %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() )
 
survey.boscalid.complete.host %>%
  group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() %>% 
  group_by(Host) %>% summarize(round(range(EC50DC),3), average= round(mean(EC50DC),2), N= n(),  sd = sd(EC50DC), se = round(sd/sqrt(N),3))
```

### By Field

```{r Stat Anal Boscalid By Field}
generate_label(  test_field (survey.boscalid.fields))

survey.boscalid.fields.2 <- survey.boscalid.fields %>% 
mutate(Field_difference = ifelse(
    Field_Year == "Baseline"| Field_Year == "NE.10.fld21.DB"|
      Field_Year == "ND.07.fld01.DB"|Field_Year == "WA.07.fld09.DB",
    "Yes","No")) %>% 
  mutate(Country_difference = ifelse(
    country == "USA",
    "Yes","No"))
  
####
myplot_comparison(survey.boscalid.fields.2,highlabel = 0.01) + expand_limits(y = c(0.06, 0.18)) +
  scale_y_continuous(breaks = c(0.06, 0.08, 0.1, 0.12, 0.14, 0.16, 0.18)) + geom_hline(yintercept =
                                                                                         0.1606226, linetype = "dashed")# make dynamic the value of intercept, create a new fucntion

```
 

## Tetraconazole 

### By Source
```{r Stat Anal Tetraconazole By Source}
test_kruskal_source(survey.tetraconazole.complete  %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() )

test_Dunn_source(survey.tetraconazole.complete  %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() )

generate_label(test_Dunn_source(survey.tetraconazole.complete  %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() ))

survey.tetraconazole.complete %>% 
 group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() %>% 
  group_by(source) %>% summarize(range(EC50DC), average= mean(EC50DC), N= n(),  sd = sd(EC50DC), se = round(sd/sqrt(N),3))
```
 
### By Host
```{r Stat Anal Tetraconazole By Host}
#generate_label
  test_kruskal_host(survey.tetraconazole.complete.host %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup())
  
#generate_label
  test_Dunn_host(survey.tetraconazole.complete.host %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup())
  #)

survey.tetraconazole.complete.host %>%
group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() %>% 
  group_by(Host) %>% summarize(range(EC50DC), average= mean(EC50DC), N= n(),  sd = sd(EC50DC), se = round(sd/sqrt(N),3))
```


### By Field
```{r Stat Anal Tetraconazole By Field}

generate_label( test_field (survey.tetraconazole.fields))

survey.tetraconazole.fields.2 <- survey.tetraconazole.fields %>% 
mutate(Field_difference = ifelse(
    Field_Year == "Baseline"| Field_Year == "NE.10.fld21.DB"|
      Field_Year == "BR.2010.unai.DB",
    "Yes","No")) %>% 
  mutate(Country_difference = ifelse(
    country == "USA",
    "Yes","No"))
####
myplot_comparison(survey.tetraconazole.fields.2, highlabel = 0.1)  + expand_limits(y = c(0.5, 1.6)) + scale_y_continuous(breaks = c(0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6)) 

```


## Picoxystrobin  
### By Source  
```{r Stat Anal Picoxystrobin By Source}
 test_kruskal_source(survey.picoxystrobin.complete %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup())
 ##
test_Dunn_source(survey.picoxystrobin.complete %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup())
# generate the label
generate_label( test_Dunn_source(survey.picoxystrobin.complete %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup()))

survey.picoxystrobin.complete %>% 
  group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup()%>% 
  group_by(source) %>% summarize(round(range(EC50DC),3), average= round(mean(EC50DC), 2), N= n(),  sd = sd(EC50DC), se = round(sd/sqrt(N),3))

```


### By Host  
```{r Stat Anal Picoxystrobin By Host}
#generate_label
  test_kruskal_host(survey.picoxystrobin.complete.host  %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup())
  
#generate_label
  test_Dunn_host(survey.picoxystrobin.complete.host  %>% group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup())

survey.picoxystrobin.complete.host %>%
   group_by(ID) %>%
   top_n(1, EC50DC) %>% #keeping the highest value of each ID value
   ungroup() %>% 
  group_by(Host) %>% summarize(round(range(EC50DC),3), average= round(mean(EC50DC),2), N= n(),  sd = sd(EC50DC), se = round(sd/sqrt(N),3))

```

### By Field

```{r Stat Anal Picoxystrobin By Field}
test_field (survey.picoxystrobin.fields)

survey.picoxystrobin.fields.2  <- survey.picoxystrobin.fields %>% 
mutate(Field_difference = ifelse(
    Field_Year == "Baseline"| Field_Year == "ND.07.fld01.DB"|
      Field_Year == "BR.2010.unai.DB"|Field_Year == "MX.2017.ca.DB"|Field_Year == "MX.2017.ley.DB" ,
    "Yes","No"))%>% 
  mutate(Country_difference = ifelse(
    country == "USA",
    "Yes","No"))
####
myplot_comparison(survey.picoxystrobin.fields.2,highlabel = 0.005)  + scale_y_continuous(breaks = c(0.008, 0.01, 0.012, 0.014, 0.016, 0.018, 0.020, 0.022, 0.024)) + geom_hline(yintercept =
                                                                                                                                                       0.02148932, linetype = "dashed")
  
```


# Table Publication
```{r Table Publication}
##
table1 <- 
survey.boscalid.complete %>% 
 bind_rows(wrangling_DC(survey.TM.2)) %>% 
  bind_rows(survey.tetraconazole.complete) %>% 
  bind_rows(survey.picoxystrobin.complete) %>% 
select(-c(EC50DC)) %>% 
  distinct(ID, .keep_all = TRUE) %>% 
group_by(source, Host, country, Field_Year, State ) %>%
  summarize(N= n()) %>%
    ungroup() %>%
        add_row(source = "Total", N  = sum(.$N))


hola.2 <- survey.TM.2 %>% select(-c(EC50DC))  %>% 
  pull(ID) %>% unique()
## Corroborating
survey.boscalid.complete %>% 
 bind_rows(survey.TM.2) %>% 
  bind_rows(survey.tetraconazole.complete) %>% 
  bind_rows(survey.picoxystrobin.complete) %>% 
select(-c(EC50DC)) %>% 
  pull(ID) %>% 
  unique()
```


## Less sensive isolates from boscalid and picoxystrobin
```{r Less sensive isolates}
boscalid.less.sensive <-survey.boscalid.fields.2 %>% arrange(desc(EC50DC))

boscalid.less.sensive[1:11,] %>% arrange(ID)

boscalid.less.sensive.2 <-boscalid.less.sensive %>% select(ID) %>% pull()

boscalid.less.sensive.3 <- sort(boscalid.less.sensive.2[1:10])
boscalid.less.sensive.3

picoxystrobin.less.sensive <-survey.picoxystrobin.fields.2 %>% arrange(desc(EC50DC))

picoxystrobin.less.sensive[1:11,] %>% arrange(ID)

picoxystrobin.less.sensive.2 <-picoxystrobin.less.sensive %>% select(ID) %>% pull()
picoxystrobin.less.sensive.3 <- sort(picoxystrobin.less.sensive.2[1:10])
picoxystrobin.less.sensive.3
```

## Corroboration complete isolates  
```{r Corroboration complete isolates}
survey.picoxystrobin.complete %>% group_by( Host) %>% summarise(ave=mean(EC50DC))

 survey.boscalid.fields %>% select(-c(EC50DC)) %>% 
 distinct(ID, .keep_all = TRUE) %>% 
group_by(Field_Year ) %>%
summarize(N= n()) 
```

# Supplementary Table: Fungicides applied in the field
```{r Supplementary Table}
#based on object table1 which is correct, just 
lasiodiplodia <- survey.boscalid.complete %>% 
 bind_rows(wrangling_DC(survey.TM.2)) %>% 
  bind_rows(survey.tetraconazole.complete) %>% 
  bind_rows(survey.picoxystrobin.complete) %>% 
select(-c(EC50DC)) %>% 
  distinct(ID, .keep_all = TRUE) %>% 
  select(ID) %>% 
  pull()

 colnames(whitemold.inventory)
 # I have to keep Field and Year separate tahts why reading again the main table from data
 
tableS1 <-  wrangling_DC(filename = whitemold.inventory) %>% mutate(source = as.factor(source),
                        Host = as.factor(Host),
                        country = as.factor(country),
                        Field = as.factor(Field),
                        Year = as.factor(Year),
                        State = as.factor(State),
                        County = as.factor(County),
                        Fungicide_current_season = as.factor(Fungicide_current_season),
                        molecule.s._current_season = as.factor(molecule.s._current_season),
                        Group_MOA_current_season = as.factor(Group_MOA_current_season),
                        Applications = as.factor(Applications),
                        Rate = as.factor(Rate))%>% select(
                          ID,
                          source,#we will change later to the real variable "source"
                          Host,
                          country, 
                          Field,
                          Year,
                          State,
                          County,
                          Fungicide_current_season,
                          molecule.s._current_season,
                          Group_MOA_current_season,
                          Applications,
                          Rate, Field_Year) %>% filter(ID %in% lasiodiplodia )
```

# Supplementary Figure Baseline Distribution 
```{r Supplementary Figure Baseline}
a.tet.baseline <-
  myplot_baseline(survey.tetraconazole.complete)  + expand_limits(y = c(0.2, 1.8)) + scale_y_continuous(breaks = c(0.2,0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,1.8 )) +   theme(
      legend.position = "none"
    ) + labs(tag = "A")  +     theme(plot.tag = element_text(face = "bold")) 
  ##  
b.bos.baseline <- myplot_baseline(survey.boscalid.complete)+ expand_limits(y = c(0.04, 0.24)) +  scale_y_continuous(breaks = c(0.04,0.06, 0.08, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.22, 0.24))+   theme(
      axis.text.x = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      legend.position = "none"
    ) + labs(tag = "B")  +     theme(plot.tag = element_text(face = "bold"))
####
c.pico.baseline <- myplot_baseline(survey.picoxystrobin.complete)  + expand_limits(y = c(0.008, 0.024)) +  scale_y_continuous(breaks = c(0.008, 0.01, 0.012, 0.014, 0.016, 0.018, 0.020, 0.022, 0.024 ))+   theme(
      axis.title.y = element_blank(),
      legend.position = "none"
    ) + labs(tag = "C")  +     theme(plot.tag = element_text(face = "bold"))
#Graph together
a.tet.baseline| (b.bos.baseline/c.pico.baseline)
```

# Figures 1 and 2 article  
**A. Tetraconazole**  
**B. Boscalid**  
**C. Picoxystrobin**
```{r Figures 1 and 2 Article}
a.tet <-
  myplot_model_1_tetraconazole (tetraconazole.joined.2.clean, lm = tetraconazole.lm)+   theme(
     legend.position = "none"
  ) + labs(tag = "A")  +     theme(  plot.tag = element_text(face = "bold")) 
#
b.bos <-
  myplot_model_1_boscalid (boscalid.joined.2.clean, lm = boscalid.lm)  +   theme(
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    legend.position = "none"
  ) + labs(tag = "B")  +     theme(  plot.tag = element_text(face = "bold")) 
###
c.pico <-
  myplot_model_1_picoxystrobin (picoxystrobin.joined.2.clean, lm = picoxystrobin.lm)  +   theme(axis.title.y = element_blank(), legend.position = "none")+ labs(tag = "C")  +     theme(  plot.tag = element_text(face = "bold")) 

#Graph together
a.tet | (b.bos / c.pico)
###
aa.tet <-
  myplot_model_2(tetraconazole.EC50DC) + expand_limits(x = c(0.27, 1.8), y = c(0.27, 1.8))  + scale_x_continuous(name = waiver(),
                                                                                                                 breaks = c(0.3, 0.6, 0.9, 1.2, 1.5, 1.8)) + scale_y_continuous(breaks = c(0.3, 0.6, 0.9, 1.2, 1.5, 1.8)) +   theme(
                                                                                                                                                                                                       legend.position = "none"
                                                                                                                 ) + geom_label(
                                                                                                                   aes(x = 0.6,
                                                                                                                       y = 1.5)
                                                                                                                   ,
                                                                                                                   label = c(paste  (
                                                                                                                     " Y =",
                                                                                                                     paste0 (round(summary(tetraconazole.lm.2)[[4]][2], 4), "x"),
                                                                                                                     round(summary(tetraconazole.lm.2)[[4]][1], 3) *
                                                                                                                       -1,
                                                                                                                     paste  ("\n R =",
                                                                                                                             paste0 (round(
                                                                                                                               summary(tetraconazole.lm.2)[[9]][1], 4
                                                                                                                             ), "\n p < 0.001"))
                                                                                                                   )),
                                                                                                                   size = 3,
                                                                                                                   fontface = "bold"
                                                                                                                 ) + labs(tag = "A")  +     theme(  plot.tag = element_text(face = "bold")) 

bb.bos <-
  myplot_model_2(boscalid.EC50DC) + expand_limits(x = c(0.04, 0.25), y = c(0.04, 0.25))  + scale_x_continuous(name = waiver(),
                                                                                                              breaks = c(0.05, 0.1, 0.15, 0.20, 0.25)) + scale_y_continuous(breaks = c(0.05, 0.1, 0.15, 0.20, 0.25)) +   theme(axis.title.x   = element_blank(), axis.title.y = element_blank(), legend.position = "none") + geom_label(
                                                                                                                aes(x = 0.08,
                                                                                                                    y = 0.2)
                                                                                                                ,
                                                                                                                label = c(paste  (
                                                                                                                  " Y =",
                                                                                                                  paste0 (round(summary(boscalid.lm.2)[[4]][2], 4), "x"),
                                                                                                                  
                                                                                                                  round(summary(boscalid.lm.2)[[4]][1], 3) *
                                                                                                                    -1,
                                                                                                                  paste  ("\n R =",
                                                                                                                          paste0 (round(
                                                                                                                            summary(boscalid.lm.2)[[9]][1], 4
                                                                                                                          ), "\n p < 0.001"))
                                                                                                                )),
                                                                                                                size = 3,
                                                                                                                fontface = "bold"
                                                                                                              ) + labs(tag = "B")  +     theme(  plot.tag = element_text(face = "bold")) 
###
 cc.pico <-
   myplot_model_2(picoxystrobin.EC50DC) + expand_limits(x = c(0.005, 0.018), y = c(0.005, 0.018))  + scale_x_continuous(name = waiver(),
                                                                                                                        breaks = c(0.0075, 0.01, 0.0125, 0.0150, 0.018)) + scale_y_continuous(breaks = c(0.0075, 0.01, 0.0125, 0.0150, 0.018))   + theme(
                                                                                                                          axis.title.y = element_blank(),
                                                                                                                          legend.position = "none",
                                                                                                                          axis.text.y = element_text(angle = 20,
                                                                                                                                                     hjust = 1), axis.text.x = element_text(angle = 20,
                                                                                                                                                     hjust = 1)
                                                                                                                        ) + geom_label(
                                                                                                                          aes(x = 0.0075,
                                                                                                                              y = 0.0148)
                                                                                                                          ,
                                                                                                                          label = c(paste  (
                                                                                                                            " Y =",
                                                                                                                            paste0 (round(summary(picoxystrobin.lm.2)[[4]][2], 4), "x"),
                                                                                                                            round(summary(picoxystrobin.lm.2)[[4]][1], 3) *
                                                                                                                              -1,
                                                                                                                            paste  ("\n R =",
                                                                                                                                    paste0 (round(
                                                                                                                                      summary(picoxystrobin.lm.2)[[9]][1], 4
                                                                                                                                    ), "\n p < 0.001"))
                                                                                                                          )),
                                                                                                                          size = 3,
                                                                                                                          fontface = "bold"
                                                                                                                        )  + labs(tag = "C")  +     theme(plot.tag = element_text(face = "bold"))
 ## together
aa.tet | (bb.bos / cc.pico)

```

### Figure 3 By Fields Article  
**A. Tetraconazole**  
**B. Boscalid**  
**C. Picoxystrobin**
```{r Figure 3 Article}

aaa.tet <-
  myplot_comparison(survey.tetraconazole.fields.2, highlabel = 0.1)  + expand_limits(y = c(0.4, 2.5)) + scale_y_continuous(breaks = c(0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4)) + theme(plot.margin = unit(c(0.5, 0.5, 0.5, 0.25), "cm")) +   
   labs(tag = "A")  +     theme(plot.tag = element_text(face = "bold"))

bbb.bos <-
  myplot_comparison(survey.boscalid.fields.2, highlabel = 0.017) + expand_limits(y = c(0.05, 0.24)) +
  scale_y_continuous(breaks = c(0.06, 0.1, 0.14,  0.18,  0.22)) + geom_hline(yintercept =  boscalid.less.sensive$EC50DC[9]
                                                                                                  , linetype = "dashed") +  theme(plot.margin = unit(c(0, 0, 0, 1), "cm")) +
  theme(
    axis.title.x   = element_blank(),
    axis.title.y = element_blank(),
    legend.position = "none"
  ) +  #theme(plot.margin = unit(c(-1, 1, 1, 1), "cm"))
  labs(tag = "B")  +     theme(plot.tag = element_text(face = "bold")) # make dynamic the value of intercept, create a new fucntion

ccc.pico <-
  myplot_comparison(survey.picoxystrobin.fields.2, highlabel = 0.007)  + expand_limits(y = c(0.008, 0.032)) +
  scale_y_continuous(
    breaks = c(
      0.008,
      0.012,
      0.016,
      0.020,
      0.024,
      0.028,
      0.032
    )
  ) + geom_hline(yintercept =
                   picoxystrobin.less.sensive$EC50DC[10] , linetype = "dashed") +# +   theme(plot.margin = unit(c(1, 1, 1, 1), "cm"))  
theme(axis.title.y = element_blank(), legend.position = "none") +   theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))  + labs(tag = "C")  +     theme(plot.tag = element_text(face = "bold"))

#Together
aaa.tet| (bbb.bos/ccc.pico)

```