Preparing R:

  • Loading the required packages:

  • Ensure the proper working directory: (consider deleting)

getwd()
## [1] "/Users/jroweibekwe/Library/Mobile Documents/com~apple~CloudDocs/MSc/Bio 9915/Bio_9915_R_Project/Bio_9915_Final_Project"

Import the spring wheat data and climate data dataset:

historical <- read.csv("/Users/jroweibekwe/Library/Mobile Documents/com~apple~CloudDocs/MSc/Bio 9915/Bio_9915_R_Project/Bio_9915_Final_Project/historical_yield_temp.csv")

View(historical)
  • Inspect the data structure:
str(historical)
## 'data.frame':    1406 obs. of  11 variables:
##  $ CAR.ID       : int  4820 4820 4820 4820 4820 4820 4820 4820 4820 4820 ...
##  $ province     : chr  "Alberta" "Alberta" "Alberta" "Alberta" ...
##  $ year         : int  1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 ...
##  $ yield        : num  34.2 30 36.1 34 39.4 33.8 42.1 35.8 48.1 44.5 ...
##  $ april..4.    : num  8.78 6.82 5.62 6.38 6.83 ...
##  $ may..5.      : num  13.1 14.5 11.3 11.4 11 ...
##  $ june..6.     : num  17.9 19.3 16.8 15.9 14.9 ...
##  $ july..7.     : num  17.9 19 20.6 18.8 18.2 ...
##  $ august..8.   : num  14.6 17.2 17.9 19.2 20.1 ...
##  $ september..9.: num  13.5 12.6 12.9 15.4 12.7 ...
##  $ X            : logi  NA NA NA NA NA NA ...
head(historical)
##   CAR.ID province year yield april..4. may..5. june..6. july..7. august..8.
## 1   4820  Alberta 1987  34.2     8.775  13.129   17.900   17.927     14.645
## 2   4820  Alberta 1988  30.0     6.820  14.490   19.295   19.011     17.239
## 3   4820  Alberta 1989  36.1     5.623  11.295   16.787   20.589     17.895
## 4   4820  Alberta 1990  34.0     6.375  11.437   15.937   18.793     19.226
## 5   4820  Alberta 1991  39.4     6.825  11.010   14.940   18.175     20.060
## 6   4820  Alberta 1992  33.8     7.585  11.430   16.560   15.695     15.925
##   september..9.  X
## 1        13.491 NA
## 2        12.588 NA
## 3        12.889 NA
## 4        15.370 NA
## 5        12.730 NA
## 6        10.585 NA

Descriptive Analysis:

Average and group the yield by province to get the yearly yield for each province

historical_yield_per_province <-historical %>%
  group_by(year, province) %>%
  summarize(MeanYieldProvince = mean(yield))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by year and province.
## ℹ Output is grouped by year.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(year, province))` for per-operation grouping
##   (`?dplyr::dplyr_by`) instead.

Yield by province over time:

Some datsets missing variables, maybe average the years or focus on Saskatchewan data up to 2019

ggplot(historical_yield_per_province, aes(x = year, y = MeanYieldProvince)) +
  geom_line() +
  geom_point(size = 1) +
  facet_wrap(~province, ncol = 1) +
  labs(x = "Year", y = "Mean Spring Wheat Yield") +
  theme_bw()
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 22 rows containing missing values or values outside the scale range
## (`geom_point()`).

Create Average temperatures per growing season each year variable:

historical <- mutate(historical, MeanTempYearly = (april..4. + may..5. + june..6. +
    july..7. + august..8. + september..9.
  ) / 6)

Average the yearly temperatures and group them by province to get the yearly temperatures per province

historical_temp_per_province <-historical %>%
  group_by(year, province) %>%
  summarize(
    MeanTempProvince = mean(MeanTempYearly))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by year and province.
## ℹ Output is grouped by year.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(year, province))` for per-operation grouping
##   (`?dplyr::dplyr_by`) instead.

Temperature by province over time:

ggplot(historical_temp_per_province, aes(x = year, y = MeanTempProvince)) +
  geom_line() +
  geom_point(size = 1) +
  facet_wrap(~province, ncol = 1) +
  labs( x = "Year",
    y = "Average Temperature per Growing Season (April-September) (°C)") +
  theme_bw()
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).

Temperature and Yield by province over time:

#Modelling the Data: Central question: Is the yield of spring wheat impacted by the temperature over time? -linear regression

model_yield<- lm(yield ~ MeanTempYearly, data = historical)

Checking modelling assumptions: -data seems normal but this model may not be the best fit -maybe create a model for each province or use mean temp province against mean yield province

autoplot(model_yield, smooth.colour = NA)
## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 1339 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1339 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1339 rows containing missing values or values outside the scale range
## (`geom_line()`).

#Modelling Interpetations -large F -very small p-value -signif?

anova(model_yield)
## Analysis of Variance Table
## 
## Response: yield
##                  Df Sum Sq Mean Sq F value    Pr(>F)    
## MeanTempYearly    1   5901  5900.7  41.424 1.704e-10 ***
## Residuals      1337 190452   142.4                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_yield)
## 
## Call:
## lm(formula = yield ~ MeanTempYearly, data = historical)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.353  -8.283  -2.165   7.446  37.919 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     59.1540     3.1996  18.488  < 2e-16 ***
## MeanTempYearly  -1.5730     0.2444  -6.436  1.7e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.94 on 1337 degrees of freedom
##   (67 observations deleted due to missingness)
## Multiple R-squared:  0.03005,    Adjusted R-squared:  0.02933 
## F-statistic: 41.42 on 1 and 1337 DF,  p-value: 1.704e-10

Replot the data and the model:

Next stepss….consider a model for each province or use only one province data