Yield loss data

We will import a dataset we obtained from collaborators and prepared for the analysis. Each row is a simulation run for one of the ten planting days in the month of June for each of the 27 years (1980 to 2007).

library(tidyverse)
simul_loss <- read.csv("data/yield_loss_simulations.csv")

Decadal estimates of yield loss

Let’s split the dataset into three periods and get the summary statistics.

## Yield loss 1980s
loss_80 <- simul_loss %>%
  filter(
    class_year != "1990-1999",
    class_year != "2000-2007"
  )
sml_80 <- summary(loss_80$relative_loss)

## Yield loss 1990s
loss_90 <- simul_loss %>%
  filter(
    class_year != "1980-1989",
    class_year != "2000-2007"
  )
sml_90 <- summary(loss_90$relative_loss)

## Yield loss 2000s
loss_20 <- simul_loss %>%
  filter(
    class_year != "1980-1989",
    class_year != "1990-1999"
  )
sml_20 <- summary(loss_20$relative_loss)

library(knitr)
table_loss <- frame_data(~"Decade", ~"Mean", ~"Median", ~"Min", ~"Max", "1980s", sml_80[4], sml_80[3], sml_80[1], sml_80[6], "1990s", sml_90[4], sml_90[3], sml_90[1], sml_90[6], "2000s", sml_20[4], sml_20[3], sml_20[1], sml_20[6])

table_loss

GAM analysis

A generalized additive model (GAM) was fitted to the data using the gam function of mgcv package. We were intersted in testing whether there was significant upward trend in yield loss in the time series, as well as to test the effect of two sowing periods (sowing dates before or after June 15th).

library(mgcv)

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:lme4':
## 
##     lmList

## The following object is masked from 'package:dplyr':
## 
##     collapse

## This is mgcv 1.8-24. For overview type 'help("mgcv-package")'.

# Check the distribution of the variable - trend
simul_loss$class_sow2 <- as.ordered(as.factor(simul_loss$class_sow))


mod_gam2 <- gam(relative_loss ~ class_sow2 + s(year) + s(year, by = class_sow2, bs = "cr"), data = simul_loss)

summary(mod_gam2)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## relative_loss ~ class_sow2 + s(year) + s(year, by = class_sow2, 
##     bs = "cr")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.6776     0.2367  36.663  < 2e-16 ***
## class_sow2.L  -1.6209     0.3347  -4.842 2.18e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf Ref.df     F  p-value    
## s(year)                      8.784  8.982 9.374 1.48e-12 ***
## s(year):class_sow2June 15-30 3.011  3.738 0.761    0.483    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.298   Deviance explained =   33%
## GCV = 16.499  Scale est. = 15.686    n = 280