01_TrialsForInclusionInAnalysis.Rmd

# Trials considered for inclusion in meta-analysis  

Trials testing fungicide efficacy on powdery mildew in mungbean in the Grains Research and Development Corporation (GRDC) northern grains region were sourced for a meta-analysis.

Data were collected from a number of researchers from the Department of Agriculture and Fisheries (DAF), the University of Southern Queensland (USQ) and the Northern Growers Alliance (NGA).
We would like to acknowledge in particular Professor Malcolm Riley's and Dr. Sue Thompson's effort in establishing and coordinating early field trials between 2010 and 2014 and GRDC for funding many of the trials which were included in this meta-analysis.

<br>  

### Load required R packages  

```{r MEsummary_Libraries, echo=TRUE, results='hide'}
if (!require("pacman"))
   install.packages("pacman")

pacman::p_load(tidyverse, kableExtra, leaflet, htmltools, lubridate, here)
```

## Criteria for inclusion in meta-analysis  

   1. A field trial testing the efficacy of fungicide on powdery mildew afflicted mungbean plants in eastern Australia.  
   2. Trial data needed to include:  
      i) Fungicide active ingredients.  
      ii) A demethylase inhibitor (DMI) fungicide.  
      iii) The date at which powdery mildew first appeared in the trial.  
      iv) Disease severity at the end of the growing seasons.  
      v) Fungicide application dates.  
      vi) Fungicide dose.  
      vii) Crop yield.  
      viii) Treatment means and accompanying variance.

<br>  
<br>  

## Import data

**Import data from trials which the raw data was available**

```{r import_dataRAW}
PM_MB_raw <-
   rbind(
      read.csv(
         here("cache/2010 PMmung Hermitage means.csv"),
         stringsAsFactors = FALSE
      ),
      read.csv(here("cache/2011 PMmung Herm means.csv"),
               stringsAsFactors = FALSE),
      read.csv(
         here("cache/2011 PMmung Kingaroy means.csv"),
         stringsAsFactors = FALSE
      ),
      read.csv(here("cache/AM1305_Goolhi_means.csv"),
               stringsAsFactors = FALSE),
      read.csv(
         here("cache/AM1304-MarysMount_means.csv"),
         stringsAsFactors = FALSE
      ),
      read.csv(
         here("cache/AM1303-Premer-Disease_means.csv"),
         stringsAsFactors = FALSE
      ),
      read.csv(
         here("cache/BB1305_Millmerran_means.csv"),
         stringsAsFactors = FALSE
      ),
      read.csv(here("cache/Herm_16_means.csv"),
               stringsAsFactors = FALSE),
      read.csv(here("cache/King_16_means.csv"),
               stringsAsFactors = FALSE),
      read.csv(here("cache/Fogerty_17_mean.csv"),
               stringsAsFactors = FALSE),
      read.csv(here("cache/Hermitage_17_mean.csv"),
               stringsAsFactors = FALSE),
      read.csv(here("cache/Wellcamp_18_mean.csv"),
               stringsAsFactors = FALSE),
      read.csv(here("cache/Hermitage_19_mean.csv"),
               stringsAsFactors = FALSE)
   )
write.csv(
   PM_MB_raw,
   file = here("cache/PM_Mungbean_SummaryOfTrialsWithRawData"),
   row.names = FALSE
)
```

<br>

Import dataset of manual entries  
**Import dataset which were manually entered from trial reports**  
```{r import_dataMan}
PM_MB_man <-
   read.csv(here("cache/PM_Mungbean_DataManuallyEntered.csv"),
            stringsAsFactors = FALSE)
```

Bind raw and transcribed datasets.  
```{r bind_data}
PM_MB_dat <-
   rbind(PM_MB_raw,
         PM_MB_man)
```


## Summary of all trials  

```{r data_summary}
PM_MB_dat %>%
   arrange(year, location) %>%
   group_by(trial_ref) %>%
   summarise(
      Year = unique(year),
      Location = unique(location),
      `Replicates per treatment` = if (all(replicates[1] == replicates)) {
         as.character(replicates[1])
      } else{
         as.character(paste(min(replicates), "-", max(replicates)))
      },
      `Trial design` = unique(trial_design),
      `Planting date` = unique(planting_date),
      `First sign of disease` = unique(first_sign_disease),
      `Fungicide treatments` = length(n_treatment),
   ) %>%
   arrange(Year) %>%
   kable(
      caption = "Description of Experiments",
      align = "c",
      col.names = c(
         "Unique Trial\nReference",
         "Year",
         "Location",
         "Replicates\nper treatments",
         "Trial design",
         "Planting date",
         "First sign\nof disease",
         "Fungicide\n treatments"
      )
   )  %>%
   kable_styling(fixed_thead = TRUE, full_width = TRUE) %>%
   column_spec(c(3, 5:6), width = "3cm") %>%
   scroll_box(height = "500px") 
```

### Experiment locations

Trial locations occurring in the GRDC Northern Grains Region.

```{r Experiment_locations, echo=FALSE}
experiment_sites <- data.frame(
   location = c(
      "Hermitage Research Station",
      "Red Vale & Kingaroy Research Station",
      "Gatton",
      "Emerald Agricultural College",
      "Bongeen 1 & 2",
      "NGA - Premer",
      "NGA - Millmerran",
      "NGA - Marys Mount",
      "NGA - Goolhi",
      "Missen Flats",
      "Wellcamp"
   ),
   lat = c(
      -28.215 ,
      -26.58,
      -27.56,
      -23.54,
      -27.51666,
      -31.4514,
      -27.87653,
      -31.0245,
      -31,
      -27.905313,
      -27.564277
   ),
   lon = c(
      152.1003,
      151.83,
      152.32,
      148.19,
      151.450,
      149.9,
      151.2567,
      150.0618,
      149.755,
      151.9772,
      151.867879
   )
)

write.csv(experiment_sites,
          "cache/Mungbean_experiment_sites.csv",
          row.names = FALSE)

leaflet(experiment_sites) %>%
   addProviderTiles("Esri.WorldTopoMap") %>%
   setView(lat = -27.31,
           lng = 150.456,
           zoom = 5) %>%
   addMarkers(
      lng = ~ lon,
      lat = ~ lat,
      popup = ~ htmlEscape(location)
   )
```


## Subset data to selection criteria  

Let's apply the selection criteria our dataset.  
All trials identified for this meta-analysis reported: 

* fungicide active ingredient,
* dose,
* first sign of disease. 

Therefore no trials need to be removed to satisfy these criteria.  


### Retain trials including: Yield

Three trials omitted due to not reporting yields  

```{r TrialWithoutYield}
PM_MB_dat %>%
   filter(is.na(grain_yield.t.ha.)) %>%
   distinct(trial_ref, year, location)

PM_MB_dat <- PM_MB_dat %>%
   filter(!is.na(grain_yield.t.ha.))
```

###Retain trials including: Disease severity

One more was trial removed as it did not report disease severity.  

```{r TrialsWithotDisease}
PM_MB_dat %>%
   filter(is.na(PM_final_severity)) %>%
   distinct(trial_ref, year, location)

PM_MB_dat <- PM_MB_dat %>%
   filter(!is.na(PM_final_severity))
```
### Retain trials including: fungicide application dates

No trials from our subset need to be removed for not reporting fungicide application dates.

```{r FungicideApplicationDates}
PM_MB_dat %>%
   group_by(trial_ref) %>%
   summarise(No_Record_Of_Fungicide_Application_Dates = all(is.na(fungicide_application_1)),
             .groups = 'drop') %>%
   filter(No_Record_Of_Fungicide_Application_Dates)
```

### Retain trials including: fungicide dose

No trials from our subset need to be removed for not reporting fungicide dose.

```{r FungicideDoseTrials}
PM_MB_dat %>%
   group_by(trial_ref) %>%
   summarise(No_Record_Of_Fungicide_Dose = all(is.na(dose_ai.ha)),
             .groups = 'drop') %>%
   filter(No_Record_Of_Fungicide_Dose)
```

### Exclude fungicides tested in too few trials

The meta-analysis should be focused on fungicides with the same mode of action.
Fungicides from the best represented group will be retained.

```{r Fungicides}
PM_MB_dat %>%
   group_by(fungicide_ai, trial_ref) %>%
   summarise() %>%
   count(sort = TRUE) %>%
   rename(Trials = n) %>%
   ggplot(aes(x = reorder(fungicide_ai, Trials), y = Trials)) +
   xlab("Fungicide active ingredient") +
   ylab("N Trials") +
   geom_col() +
   ggtitle(label = "Number of trials in which the\nspecified fungicide was used") +
   coord_flip() +
   theme_classic()
```

The demethylation inhibitors (DMI), tebuconazole and propiconazole, are used in the highest frequencies.
The DMIs have the same fungicide mode of action and are good candidates to be pooled in the meta-analysis.

Amistar Xtra and Custodia both contain strobilurin and triazole, however, because they contain differing dose ratios (inverted) pooling may not be appropriate therefore they will be removed.  


### Retain only fungicides with DMI action

```{r filterDMI_trials}
DMI_Trials <-
   PM_MB_dat %>%
   filter(fungicide_ai == "tebuconazole" |
             fungicide_ai == "propiconazole") %>%
   distinct(trial_ref) %>%
   pull()

PM_MB_dat <-
   PM_MB_dat[PM_MB_dat$trial_ref %in% DMI_Trials, ]
```

Now remove any non-DMI treatments or controls.

```{r KeepOnlyDMI}
PM_MB_dat <-
   PM_MB_dat %>%
   filter(
      fungicide_ai == "control" |
         fungicide_ai == "tebuconazole" |
         fungicide_ai == "propiconazole"
   )
```


### Classify dataset variables

Importantly, the class of each variable in the data should be defined, retaining only the variables relevant to the analysis. 

```{r classifyData}
source("R/classify_PM_data.R") # use a custom code

# classify dataset with Trials in both 
PM_MB_dat <-
   classify_PM_data(PM_MB_dat)

```

## Defining spray schedule variable  

First are going to calculate the time between the first sign of disease and the fungicide applications.

```{r clustered_fungicide_applications}
# collective data
PM_MB_dat %<>%
   mutate(fungicide_timing_1 = fungicide_application_1 - first_sign_disease) %>%
   mutate(fungicide_timing_2 = fungicide_application_2 - fungicide_application_1) %>%
   mutate(fungicide_timing_3 = fungicide_application_3 - fungicide_application_2)

```

To ensure sufficient number of replicates, treatments were binned by the first fungicide application date into three categorical variables relating to when the first fungicide application was made, relative to the first sign of disease.  

These categorical variables are named:  
   - **Early**: First fungicide application was prior to first sign of disease.  
   - **Recommended**: First fungicide application was applied on the day powdery mildew was observed, or within three days of first sign.  
   - **Late**: First fungicide application was four or more days after first sign of disease being observed.  

The number of follow-up sprays need also be defined.

```{r TreatmentTable}
data.frame(
   TreatmentName = c(
      "Early",
      "Recommended",
      "Late",
      "EarlyPlus",
      "RecommendedPlus",
      "LatePlus"
   ),
   n_sprays = rep(c("Single", "Two - Three"), each = 3),
   DaysRelativeToFirstSign = c(
      "Prior to first sign of Powdery Mildew",
      "1 - 3 days after first sign of Powdery Mildew",
      "7 - 8 days after first sign of Powdery Mildew"
   )
) %>%
   kable()
```

Simplify the clusters for the data-set
```{r simple_Y_clusters}
PM_MB_dat <- PM_MB_dat %>%
   mutate(
      spray_management = case_when(
         fungicide_timing_1 < 0 &
            is.na(fungicide_application_2) &
            is.na(fungicide_application_3) ~ "Early",
         fungicide_timing_1 >= 0 &
            fungicide_timing_1 < 4 &
            is.na(fungicide_application_2) &
            is.na(fungicide_application_3) ~ "Recommended",
         fungicide_timing_1 >= 4 &
            is.na(fungicide_application_2) &
            is.na(fungicide_application_3) ~ "Late",
         fungicide_timing_1 < 0 &
            !is.na(fungicide_application_2) ~ "Early_plus",
         fungicide_timing_1 >= 0 &
            fungicide_timing_1 < 4 &
            !is.na(fungicide_application_2) ~ "Recommended_plus",
         fungicide_timing_1 >= 4 &
            !is.na(fungicide_application_2) ~ "Late_plus",
         TRUE ~ "Other"
      )
   )

PM_MB_dat[PM_MB_dat$fungicide_ai == "control",
          c(
             "fungicide_timing_1",
             "fungicide_timing_2",
             "fungicide_timing_3",
             "spray_management"
          )] <- "control"
```


Now to view the number break-down of the `spray_management` treatments in the *yield* and *disease severity* data-sets

```{r tableSprayManagement_Y}
table(PM_MB_dat$spray_management)
```

'Early_plus' treatments are few in number, these treatments will have too few comparisons with other treatments in the meta-analysis to provide accurate results. Therefore we will remove 'Early_plus' from the analysis.

```{r simpler_clusters_remove}
PM_MB_dat <-
   PM_MB_dat %>%
   filter(spray_management != "Early_plus") %>%
   dplyr::select(
      -c(
         fungicide_timing_1,
         # remove the following columns which no longer have use
         fungicide_timing_2,
         fungicide_timing_3,
         fungicide_application_1,
         fungicide_application_2,
         fungicide_application_3,
         fungicide_application_4,
         fungicide_application_5,
         fungicide_application_6,
         fungicide_application_7,
      )
   )
```

## Identify variance

We need variance accompanying the mean observations for each treatment to successfully undertake a meta-analysis.
Let's investigate which of the included trials included variance for both response variables of interest, yield and disease severity.

```{r Trials_For_Inclusion}
PM_MB_dat %>%
   mutate(MA_analysis = case_when(
      is.na(yield_error & disease_error) == FALSE ~ "Yield and Disease",
      is.na(yield_error) == FALSE & 
         is.na(disease_error) ~ "Yield",
      is.na(disease_error) == FALSE & 
         is.na(yield_error) ~ "Disease",
      TRUE ~ "Nil"
   )) %>%
   distinct(trial_ref, location, year, MA_analysis) %>%
   arrange(year) %>%
   kable()
```

Remove one trial in Emerald which did not report disease variance or yield variance.

```{r remove_nil_var}
# Remove trials without variance for both Disease and yield for data_set `PM_MB_dat
PM_MB_dat <- 
   PM_MB_dat %>%
   mutate(contains_var = case_when(
      is.na(yield_error) &
         is.na(disease_error) ~ NA_character_,
      TRUE ~ trial_ref)) %>%
   filter(trial_ref %in% contains_var) %>%
   dplyr::select(-contains_var) 
```

Following an inspection of the yield and disease variance we determined imputation was not suitable and therefore we must exclude trials that did not report variance.
For each analysis, yield and disease severity a separate data-set will be created, including only trials that report variance for the respective response variable

One trial will be removed from the yield data set `PM_dat_Y` for not reporting yield variance.

```{r TrialsWithVariance}
PM_MB_dat %>%
   mutate(No_Var = is.na(yield_error)) %>%
   filter(No_Var == TRUE) %>%
   distinct(trial_ref, year, location)

PM_dat_Y <- 
   PM_MB_dat %>%
   mutate(contains_var = case_when(
      is.na(yield_error) ~ NA_character_,
      TRUE ~ trial_ref)) %>%
   filter(trial_ref %in% contains_var) %>%
   dplyr::select(-contains_var)
```

Two more trials will be removed for our disease severity data set `PM_dat_D` for not reporting yield variance.

```{r TrialsWithVariance_D}
PM_MB_dat %>%
   mutate(No_Var = is.na(disease_error)) %>%
   filter(No_Var == TRUE) %>%
   distinct(trial_ref, year, location)

PM_dat_D <- 
   PM_MB_dat %>%
   mutate(contains_var = case_when(
      is.na(disease_error) ~ NA_character_,
      TRUE ~ trial_ref)) %>%
   filter(trial_ref %in% contains_var) %>%
   dplyr::select(-contains_var)
```


## Standardise Variance

All the variance types need to be converted to sample variance.
First we will start by standardising yield variance according to the method reported by [Ngugi et.al (2011)](https://apsjournals.apsnet.org/doi/10.1094/PHYTO-08-10-0221).

## Yield variance standardisation

```{r StandardiseYieldVariance}
PM_dat_Y <-
   PM_dat_Y %>%
   mutate(
      vi =
         case_when(
            Y_error_type == "stdev" ~ yield_error ^ 2,
            Y_error_type == "lsd (P=0.05)" ~
               (replicates * ((yield_error / 1.96) ^ 2) / 2)
         ),
      id = row_number(),
      spray_management = fct_relevel(spray_management, sort)
   ) %>%
   dplyr::select(-c(yield_error,
                    Y_error_type))

```

## Disease severity variance standardisation

```{r StandardiseSeverityVariance}
PM_dat_D <-
   PM_dat_D %>%
   mutate(
      vi =
         case_when(
            D_error_type == "stdev" ~ disease_error ^ 2,
            D_error_type == "lsd (P=0.05)" ~
               (replicates * ((disease_error / 1.96) ^ 2) / 2)
         ),
      id = row_number(),
      spray_management = fct_relevel(spray_management, sort)
   ) %>%
   dplyr::select(-c(disease_error,
                    D_error_type))

```

## Save the cleaned data

```{r saveData01}
write.csv(PM_dat_Y, file = "cache/PM_yield_clean_data.csv", row.names = FALSE)
write.csv(PM_dat_D, file = "cache/PM_disease_clean_data.csv", row.names = FALSE)
write.csv(PM_MB_dat, file = "cache/PM_MB_clean_data.csv", row.names = FALSE)
save(PM_MB_dat,
     PM_dat_D,
     PM_dat_Y,
     file = here("cache/ImportDataAndSelectTrials01.Rdata"))
```