03_Prep_data.Rmd

# Prepare data for meta-analysis
## Set environment and import data  

```{r load_libraries, include=FALSE}
if (!require("pacman"))
   install.packages("pacman")
pacman::p_load(tidyverse, kableExtra, agricolae
)
```

```{r import_data}
source("R/import_data.R")
PM_MB_means <- import_data()

# Calculate AUDPC using final severity and first sign of disease - 7 days
source("R/AUDPC_for_single_obs.R") 
```


## Standardising sample variance  
Meta-analyses use variance from the data to calculate the effect sizes. The type of variance needs to be the same across all the data incorporated in the meta-analysis.  

We are narrowing the meta-analysis to only include demethylation inhibitors, tebuconazole and propiconazole. Within studies which have used these fungicides, what type of error was reported?  

```{r variance_types}
PM_MB_means %>%
   filter(fungicide_ai == "tebuconazole" |
             fungicide_ai == "propiconazole",
          row_spacing == 0.75) %>%
   group_by(trial_ref, location, Y_error_type, fungicide_ai) %>%
   summarise(DMI_trts = length(Y_error_type)) %>%
   arrange(Y_error_type) %>%
   kable(caption = "Number of treatments for each trial and the error type reported") %>%
   kable_styling(fixed_thead = T) %>% 
   footnote(general = "DMI - demethylation inhibitors (tebuconazole and propiconazole).")
```

The experiment summaries report either standard deviation, least square differences or no error term (`NA`). Below we will convert LSD and standard deviation to sample variance. In addition, we need to impute variances for the trials which did not report any form of variance (`NA`).  



### Converting least-squares to sample variance

```{r all_LSD_trials}
Trials <- PM_MB_means %>%
   filter(Y_error_type == "lsd (P=0.05)") %>%
   select(trial_ref) %>%
   distinct()
```


The following code cnverts 'least square differences' (LSD) to sampling variance as per the method in [Nugugi et.al (2011)](https://apsjournals.apsnet.org/doi/10.1094/PHYTO-08-10-0221).
Ideas for using a T-critical value of 1.697 came from reading a [statisics-how-to website](https://www.statisticshowto.datasciencecentral.com/how-to-calculate-the-least-significant-difference-lsd/).
A [table of T-critical values](https://www.statisticshowto.datasciencecentral.com/tables/t-distribution-table/) was consulted where 30 degrees of freedom within experiment was used to find the 'T-crit' value of 1.697.  

![Calculating sampling variance from mean squared error Paul et al. (2008)](Paul_etal_2008.PNG)  

```{r LSD_2_Sampling Variance}

for (Trial2012 in Trials$trial_ref) {
   Tcrit <-
      1.697  # 0.05 and DFw = 30 # T critical value for which there is significant difference between two groups; This value is the same for both trials "mung1112/01" and "mung1112/02"
   DFw <-
      sum(PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$replicates) -
      length(PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$replicates) # degrees of freedom within groups
   
   LSD <-
      PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$yield_error[1]
   
   
   V_yield <-
      (PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$replicates[1] * (LSD /
                                                                             Tcrit) ^ 2) / 2 # method - Nugugi et.al (2011) DOI: 10.1094/phyto-08-10-0221
   
   
   for (j in seq_along(PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$replicates)) {
      if (j == 1) {
         Sv <-
            vector(length = length(PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$replicates))
      } # reset Sv vector if this is the first loop
      
      
      
      
      ###############################################3
      # ___ Adjusting formula to suit our data_______
      #
      #  Sv[j]^2 = MSE/(PM_MB_means[PM_MB_means$trial_ref == Trial2012,]$replicates[j] *
      #                   PM_MB_means[PM_MB_means$trial_ref == Trial2012,]$grain_yield.t.ha[j])
      # Balance equation
      # Sv[j] = V_yield /
      #    (PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$replicates[j] *
      #       as.numeric(PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$grain_yield.t.ha[j])^2)
      # #################################################################
      
      
      
      
      Sv[j] = V_yield /
         PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$replicates[j]
      
      # Sv[j] = MSE /
      #    (PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$replicates[j] *
      #       as.numeric(PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$grain_yield.t.ha[j])^2)
      #
      
      ### -> This formula should be V = (n * (LSD/1.96)^2) / 2
      
   }
   
   # and replacing the values in the dataset
   
   PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$yield_error <-
      Sv
   PM_MB_means[PM_MB_means$trial_ref == Trial2012, ]$Y_error_type <-
      "SamplingVariance"
   
}
```


### Converting standard deviation to sample variance
Now to estimate the variance from the standard deviation

```{r stdev2SampVariance}
# subset data by those which describe the Yield error with standard deviation
Trials <- PM_MB_means %>%
   filter(Y_error_type == "stdev") %>%
   select(trial_ref) %>%
   distinct()


# calculate the Sampling variance from the standard deviation
# standard error of the mean equals the standard deviation divided by the square root of the number of samples
# Sv <- sem^2 <- stdev/sqrt(n)
# Sampling variance equals the square of the standard error.
# Sv <- (stdev/sqrt(n))^2

# Therefore

for (j in seq_along(PM_MB_means[PM_MB_means$trial_ref %in% Trials$trial_ref,]$yield_error)) {
   if (j == 1) {
      Sv <-
         vector(length = length(PM_MB_means[PM_MB_means$trial_ref %in% Trials$trial_ref,]$yield_error))
   }# create empty vector at first iteration with length of j
   
   Sv[j] = (PM_MB_means[PM_MB_means$trial_ref %in% Trials$trial_ref,]$yield_error[j]) ^
      2
}

# and replacing the values in the dataset
PM_MB_means[PM_MB_means$trial_ref %in% Trials$trial_ref,]$yield_error <-
   Sv
PM_MB_means[PM_MB_means$trial_ref %in% Trials$trial_ref,]$Y_error_type <-
   "SamplingVariance"
```


## Yield saving  
Yield saving and Yield saving proportional to the control.  

Data from `r length(PM_MB_means$trial_ref)` trials were collected and the means of each treatment calculated (`R/Cleaning historical data.R`). These data were combined into one file `1902 powdery mildew-Mungbean - Collated means.csv`.

```{r Calculate_yield-gains}
PM_MB_means$grain_yield.t.ha. <-
   as.numeric(PM_MB_means$grain_yield.t.ha.) # NAs produced due to some cells having text description to why there is no specific data

for (i in unique(PM_MB_means$trial_ref)) {
   dat1 <- PM_MB_means[PM_MB_means$trial_ref == i, ]
   
   if (any(is.na(unique(dat1$row_spacing)))) {
      warning(
         paste(
            unique(dat1$trial_ref),
            "at",
            unique(dat1$location),
            "in year",
            unique(dat1$year),
            "has unknown row.spacing (NA)\n 'yield_gain' and 'proportional yield gain' not calculated\n"
         )
      )
      next()
   }
   
   for (j in unique(dat1$row_spacing)) {
      dat2 <- dat1[dat1$row_spacing == j, ]
      controlY <-
         mean(dat2[dat2$fungicide_ai == "control", "grain_yield.t.ha."])
      
      dat2$yield_gain <- dat2$grain_yield.t.ha. - controlY
      dat2$prop_YG <- dat2$yield_gain  /  controlY
      
      dat2[dat2$fungicide_ai == "control", c("yield_gain", "prop_YG")] <-
         NA
      
      dat1[dat1$row_spacing == j, ] <- dat2
   }
   PM_MB_means[PM_MB_means$trial_ref == i, ] <- dat1
}

# ____________________________________________
# calculate yield gain and prop_yield gain for experiments with no row spacing data

dat1 <- PM_MB_means[PM_MB_means$trial_ref == "AM1303", ]

controlY <-
   mean(dat1[dat1$fungicide_ai == "control", "grain_yield.t.ha."])

dat1$yield_gain <- dat1$grain_yield.t.ha. - controlY
dat1$prop_YG <- dat1$yield_gain / controlY
dat1[dat1$fungicide_ai == "control", c("yield_gain", "prop_YG")] <-
   c(NA, NA)
PM_MB_means[PM_MB_means$trial_ref == "AM1303", ] <- dat1

write.csv(PM_MB_means,
          "cache/1911_PM_MB_means&Ygains.csv")
```

## Reduce data set to required variables  
Let's tidy up the data set into a data frame of only the variables which are required for the analysis.  

First remove the irrelevant columns/variables. Then calculate fungicide application timing variables relative to the first sign of disease. Entries with NA in the grain yield or powdery mildew severity are also removed from the data.    

```{r clustered_fungicide_applications}
source("R/slimming_PM_dat.R") # see R script for the code which reduces the data.

slim_PM_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) %>%
   filter(!is.na(grain_yield.t.ha)) %>%
   filter(!is.na(PM_final_severity))

```

### Cluster fungicide applications  
We are going to define 'fungicide application timing' as a categorical variable, relative to the first sign of powdery mildew.  

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.  

Follow up fungicide applications (cluster_2 and cluster_3) were also categorised as either, early, recommended or late.  

```{r cluster}
slim_PM_dat %<>%
   mutate(
      cluster_1 = case_when(
         fungicide_timing_1 < 0 ~ "Early",
         fungicide_timing_1 >= 0 &
            fungicide_timing_1 < 4 ~ "Recommended",
         TRUE ~ "Late"
      )
   ) %>%
   mutate(
      cluster_2 = case_when(
         is.na(fungicide_timing_2) ~ NA_character_,
         fungicide_timing_2 < 13 ~ "Early",
         fungicide_timing_2 >= 13 &
            fungicide_timing_2 < 17 ~ "Recommended",
         TRUE ~ "Late"
      )
   ) %>%
   mutate(
      cluster_3 = case_when(
         fungicide_timing_3 < 14 ~ "Early",
         is.na(fungicide_timing_3) ~ NA_character_,
         TRUE ~ "Recommended"
      )
   )
```

Let's look at these clusters to see where we have the most degrees of freedom.

```{r cluster_clusters}
clusters <-
   unique(slim_PM_dat[c("cluster_1", "cluster_2", "cluster_3")])

for (i in 1:nrow(clusters)) {
   if (i == 1) {
      clusters$n <- NA
      clusters$n_x_trial <- NA
   }
   
   clusters[i, "n"] <-
      nrow(slim_PM_dat[slim_PM_dat$cluster_1 == clusters[i, "cluster_1"] &
                          slim_PM_dat$cluster_2 == clusters[i, "cluster_2"] &
                          slim_PM_dat$cluster_3 == clusters[i, "cluster_3"],])
   
   clusters[i, "n_x_trial"] <-
      length(unique(slim_PM_dat[slim_PM_dat$cluster_1 == clusters[i, "cluster_1"] &
                                   slim_PM_dat$cluster_2 == clusters[i, "cluster_2"] &
                                   slim_PM_dat$cluster_3 == clusters[i, "cluster_3"],
                                "trial_ref"]))
   
}

knitr::kable(clusters)
```

From this we can see that there are 16 different cluster combinations, which would make the analysis too difficult. I will reduce this to 6 levels within one variable.  

Single sprays either early, recommended timing or late.

Two sprays or more where the first spray was either early, recommended timing, or late.

```{r simple_clusters}
slim_PM_dat <- slim_PM_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 ~ "Late_plus",
         TRUE ~ "Other"
      )
   )

slim_PM_dat[slim_PM_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

```{r tableSprayManagement}
table(slim_PM_dat$spray_management)

# For each treatment, when were the first fungicide applications made 
# relative to the first sign of disease
table(slim_PM_dat$spray_management,slim_PM_dat$fungicide_timing_1)
```

There are few 'Early_plus' and 'late' treatments, these treatments therefore will have too few comparisons with other treatments in the meta-analysis to be provide accurate results. Therefore we will remove 'Early_plus' and 'Late' from the analysis.  


```{r simpler_clusters_REMOVE}
slimmer_PM_dat <-
   slim_PM_dat[slim_PM_dat$spray_management != "Late" &
                  slim_PM_dat$spray_management != "Early_plus", ]

write.csv(slimmer_PM_dat, file = "cache/slimmer_PM_clusterdat.csv", row.names = FALSE)
```

Our data currently uses sample variance to define the yield error. We should expect that when yield increases so does the variance. Let's quickly plot this to check.

```{r Var_mean_sqaures}
# redefine sample variance calculated from mean square
slimmer_PM_dat$vi <-
   slimmer_PM_dat$Y_Msquare /
   (slimmer_PM_dat$n * slimmer_PM_dat$grain_yield.t.ha ^ 2)

slimmer_PM_dat %>%
   ggplot(aes(x = grain_yield.t.ha, y = vi)) +
   geom_point()
```

We see here that the sample variance generally decreases or has no association with the grain yield. This might be due to the unpredictable nature of mungbeans, however going forward we will use mean squares which defines between trial variation and is a the more generic error term.  

```{r Var_MeanSquarePlot}
slimmer_PM_dat %>%
   ggplot(aes(x = grain_yield.t.ha, y = Y_Msquare)) +
   geom_point()
```

&nbsp;  
&nbsp;  
&nbsp;  

## AUDPC visualise raw data  

>How well does 'area under the disease progress curve' (AUDPC) predict grain yield?  

and _Does spray_management scenario mitigate yield loss through reducing AUDPC?_  

We would expect to see grain yield decline with increasing AUDPC. Let's just plot yield vs AUDPC.  

```{r Yield_vs_AUDPC}

plot(grain_yield.t.ha. ~ AUDPC_m, data = PM_MB_means)
abline(lm(grain_yield.t.ha. ~ AUDPC_m, data = PM_MB_means))

```

There appears to be no relationship between the two continuous variables.  
Does a linear model confirm this?  

```{r AUDPClm}
summary(lm(grain_yield.t.ha. ~ AUDPC_m, data = PM_MB_means))
```

There is no relationship as per the linear model.  

Perhaps the variation due to trial location, season, row spacing, and other seasonal variables are preventing us from seeing a relationship. Let's do the plot again for each trial where there are no missing values.  
Note: The plots are ordered by mean yield of the trial.

```{r AUDPCxTrial}
sm2lg <- PM_MB_means %>%
   group_by(trial_ref) %>%
   summarise(TY = mean(grain_yield.t.ha., na.rm = TRUE)) %>%
   arrange(TY)

for (i in sm2lg$trial_ref) {
   dat1 <- PM_MB_means[PM_MB_means$trial_ref == i, ]
   if (any(is.na(dat1$grain_yield.t.ha.))) {
      next
   }
   if (any(is.na(dat1$AUDPC_m))) {
      next
   }
   
   plot(
      grain_yield.t.ha. ~ AUDPC_m,
      data = dat1,
      main = paste(dat1$location[1], (dat1$year[1])),
      ylim = c(0, max(PM_MB_means$grain_yield.t.ha., na.rm = TRUE)),
      xlim = c(0, max(PM_MB_means$AUDPC_m, na.rm = TRUE))
   )
   abline(lm(grain_yield.t.ha. ~ AUDPC_m, data = dat1))
   points(grain_yield.t.ha. ~ AUDPC_m,
          data = dat1[dat1$fungicide_ai == "control",],
          col = "red")
   
}
```

This sequence of graphs show that AUDPC and disease incidence does not impact on yields if the season is a low yielding, ie lower than 1 tonne per hectare.  
The higher the yield the greater greater the influence on the disease on the yield especially if the crop is high yielding with AUDPC also being large enough to reduce yield.  

Let's have a look again at the same sequence of graphs and analyses for area under the disease progress _stairs_.  

```{r Yield_vs_AUDPS}
plot(grain_yield.t.ha. ~ AUDPS_m, data = PM_MB_means)
abline(lm(grain_yield.t.ha. ~ AUDPS_m, data = PM_MB_means))

summary(lm(grain_yield.t.ha. ~ AUDPS_m, data = PM_MB_means))

sm2lg <- PM_MB_means %>%
   group_by(trial_ref) %>%
   summarise(TY = mean(grain_yield.t.ha., na.rm = TRUE)) %>%
   arrange(TY)



for (i in sm2lg$trial_ref) {
   dat1 <- PM_MB_means[PM_MB_means$trial_ref == i,]
   if (any(is.na(dat1$grain_yield.t.ha.))) {
      next
   }
   if (any(is.na(dat1$AUDPS_m))) {
      next
   }
   
   
   plot(
      grain_yield.t.ha. ~ AUDPS_m,
      data = dat1,
      main = paste(dat1$location[1], (dat1$year[1])),
      ylim = c(0, max(
         PM_MB_means$grain_yield.t.ha., na.rm = TRUE
      )),
      xlim = c(0, max(PM_MB_means$AUDPS_m, na.rm = TRUE))
   )
   abline(lm(grain_yield.t.ha. ~ AUDPS_m, data = dat1))
   points(grain_yield.t.ha. ~ AUDPS_m,
          data = dat1[dat1$fungicide_ai == "control", ],
          col = "red")
   
}

```


```{r delete_this_chunk, eval=FALSE, include=FALSE}
PM_MB_means$infected_days <- as.numeric(PM_MB_means$harvest_date - PM_MB_means$first_sign_disease)

plot(AUDPC_m ~ infected_days, data = PM_MB_means)
plot(grain_yield.t.ha. ~ infected_days, data = PM_MB_means)

PM_MB_means$harvest_date - PM_MB_means$first_sign_disease
PM_MB_means[0 > as.Date(PM_MB_means$harvest_date, format = "%Y-%m-%d")-
as.Date(PM_MB_means$first_sign_disease, format = "%Y-%m-%d"),]
```