Reproducible report:
Silicon, biological seed treatment and cutting reduce the intensity of leaf spot diseases affecting Lolium multiflorum
Created: 2021-10-06
Last update: 2022-01-27
Essential packages
library(tidyverse)
library(gsheet)
library(cowplot)
library(patchwork)
library(ggthemes)
theme_set(theme_half_open())dat <- read_csv("data/dat_treatments.csv")Curves 2018 check
p2018 <- dat %>%
filter(treat == "CHK") %>%
filter(cut == "0C") %>%
filter(year == "2018") %>%
group_by(time, treat, cut) %>%
summarize(
`Brown spot` = mean(bip),
`Parrot's eye` = mean(cer),
`Gray leaf spot` = mean(pyr),
mean_all = mean(all)
) %>%
# select(-mean_all) %>%
ungroup() %>%
pivot_longer(4:6, names_to = "disease", values_to = "sev") %>%
ggplot(aes(time, sev, fill = disease, group = disease)) +
geom_area(alpha = 0.7) +
theme_clean() +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white"),
legend.text = element_text(face = c(rep("italic", 5), rep("plain", 5)))
) +
scale_fill_canva() +
ylim(0, 25) +
labs(
fill = "Disease", x = "Day after disease onset (17 August)",
y = "Percent area affected (%)",
title = "2018"
)## `summarise()` has grouped output by 'time', 'treat'. You can override using the `.groups` argument.## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.p2018
Curves 2019 check
p2019 <- dat %>%
filter(treat == "CHK") %>%
filter(cut == "0C") %>%
filter(year == "2019") %>%
group_by(time, treat, cut) %>%
summarize(
`Brown spot` = mean(bip),
`Parrot's eye` = mean(cer),
`Gray leaf spot` = mean(pyr),
mean_all = mean(all)
) %>%
# select(-mean_all) %>%
ungroup() %>%
pivot_longer(4:6, names_to = "disease", values_to = "sev") %>%
ggplot(aes(time, sev, fill = disease, group = disease)) +
geom_area(alpha = 0.7) +
theme_clean() +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white"),
legend.text = element_text(face = c(rep("italic", 5), rep("plain", 5)))
) +
scale_fill_canva() +
ylim(0, 25) +
labs(
fill = "Disease", x = "Day after disease onset (10 May)",
y = "Percent area affected (%)",
title = "2019"
)## `summarise()` has grouped output by 'time', 'treat'. You can override using the `.groups` argument.## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.p2019
Figure 1
(p2018 | p2019) +
plot_layout(guides = "collect")
ggsave("figs/p1.png", width = 9, height = 4, bg = "white")All curves
dat2 <- dat %>%
filter(year == "2018") %>%
group_by(time, treat, cut) %>%
summarize(
`Brown spot` = mean(bip),
`Parrot's eye` = mean(cer),
`Gray leaf spot` = mean(pyr),
mean_all = mean(all)
) %>%
# select(-mean_all) %>%
ungroup() %>%
pivot_longer(4:6, names_to = "disease", values_to = "sev")## `summarise()` has grouped output by 'time', 'treat'. You can override using the `.groups` argument.All curves 2018 (Fig. S1)
p20180cut <- dat2 %>%
filter(cut == "0C") %>%
ggplot(aes(time, sev, color = treat, groups = treat)) +
geom_line() +
geom_point() +
scale_color_colorblind() +
facet_wrap(~disease, ncol = 1) +
ylim(0, 18) +
theme_clean() +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
labs(
title = "2018",
subtitle = "Uncut",
x = "Day after 17 August",
y = "Percent area affected (%)"
)
p20181cut <- dat2 %>%
filter(cut == "1C") %>%
ggplot(aes(time, sev, color = treat, groups = treat)) +
geom_line() +
geom_point() +
annotate("segment",
linetype = "dotted", x = 12, xend = 12, y = 12, yend = 0,
arrow = arrow()
) +
annotate("text", x = 12, y = 14, label = "29 Aug") +
scale_color_colorblind() +
facet_wrap(~disease, ncol = 1) +
theme_clean() +
ylim(0, 18) +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
labs(
title = "2018",
subtitle = "One cut",
x = "Day after 17 August",
y = "Percent area affected (%)"
)
p20182cut <- dat2 %>%
filter(cut == "2C") %>%
ggplot(aes(time, sev, color = treat, groups = treat)) +
geom_line() +
geom_point() +
annotate("segment",
linetype = "dotted", x = 12, xend = 12, y = 12, yend = 0,
arrow = arrow()
) +
annotate("text", x = 12, y = 14, label = "29 Aug") +
annotate("segment",
linetype = "dotted", x = 49, xend = 49, y = 12, yend = 0,
arrow = arrow()
) +
annotate("text", x = 49, y = 14, label = "5 Oct") +
scale_color_colorblind() +
facet_wrap(~disease, ncol = 1) +
theme_clean() +
ylim(0, 18) +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
labs(
title = "2018",
subtitle = "Two cuts",
x = "Day after 17 August",
y = "Percent area affected (%)"
)
(p20180cut + p20181cut + p20182cut) +
plot_layout(guides = "collect")
ggsave("figs/p9.png", width = 12, height = 8, bg = "white")All curves 2019 (Fig. S2)
dat3 <- dat %>%
filter(year == "2019") %>%
group_by(time, treat, cut) %>%
summarize(
Bipolaris = mean(bip),
Cercospora = mean(cer),
Pyricularia = mean(pyr),
mean_all = mean(all)
) %>%
ungroup() %>%
pivot_longer(4:7, names_to = "disease", values_to = "sev") %>%
filter(disease == "mean_all")## `summarise()` has grouped output by 'time', 'treat'. You can override using the `.groups` argument.p20190cut <- dat3 %>%
filter(cut == "0C") %>%
ggplot(aes(time, sev, color = treat, groups = treat)) +
geom_line() +
geom_point() +
scale_color_colorblind() +
theme_clean() +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
labs(
title = "2019",
subtitle = "Uncut",
x = "Day after 10 May",
y = "Percent area affected (%)"
)
p20191cut <- dat3 %>%
filter(cut == "1C") %>%
ggplot(aes(time, sev, color = treat, groups = treat)) +
geom_line() +
geom_point() +
scale_color_colorblind() +
theme_clean() +
annotate("segment",
linetype = "dotted", x = 17, xend = 17, y = 20, yend = 0,
arrow = arrow()
) +
annotate("text", x = 17, y = 20, label = "27 May") +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
labs(
title = "2019",
subtitle = "One cut",
x = "Day after 10 May",
y = "Percent area affected (%)"
)
p20192cut <- dat3 %>%
filter(cut == "2C") %>%
ggplot(aes(time, sev, color = treat, groups = treat)) +
geom_line() +
geom_point() +
scale_color_colorblind() +
theme_clean() +
annotate("segment",
linetype = "dotted", x = 17, xend = 17, y = 20, yend = 0,
arrow = arrow()
) +
annotate("text", x = 17, y = 20, label = "27 May") +
annotate("segment",
linetype = "dotted", x = 38, xend = 38, y = 20, yend = 0,
arrow = arrow()
) +
annotate("text", x = 38, y = 20, label = "8 Jul") +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
labs(
title = "2019",
subtitle = "Two cuts",
x = "Day after 10 May",
y = "Percent area affected (%)"
)
(p20190cut + p20191cut + p20192cut) +
plot_layout(guides = "collect")
ggsave("figs/p11.png", width = 12, height = 3, bg = "white")AUDPC calculation
library(epifitter)dat1 <- dat %>%
group_by(year, treat, cut, rep) %>%
summarize(audpc = AUDPC(time, all))## `summarise()` has grouped output by 'year', 'treat', 'cut'. You can override using the `.groups` argument.dat1_2018 <- dat1 %>%
filter(year == 2018) %>%
mutate(audpc2 = audpc / 72)
dat1_2019 <- dat1 %>%
filter(year == 2019) %>%
mutate(audpc2 = audpc / 119)
dat2 <- rbind(dat1_2018, dat1_2019)Box plots AUDPC
p2 <- dat2 %>%
ggplot(aes(treat, audpc2, color = cut)) +
geom_boxplot(outlier.colour = NA) +
geom_jitter(color = "grey80", width = 0.1) +
facet_grid(~cut) +
ylim(0, 17) +
scale_color_colorblind() +
theme_clean() +
theme(
legend.position = "none",
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
labs(
x = "Treatment",
y = "Standardized AUDPC",
color = "No. of cuts"
)Mixed model
library(lme4)
library(DHARMa)
library(performance)m1 <- lmer(audpc2 ~ treat * cut + (1 | year),
data = dat2
)
plot(simulateResiduals(m1))
check_model(m1)
car::Anova(m1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: audpc2
## Chisq Df Pr(>Chisq)
## treat 226.908 3 < 2.2e-16 ***
## cut 109.278 2 < 2.2e-16 ***
## treat:cut 19.959 6 0.002816 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(report)
report(m1)## 'interpret_d()' is now deprecated. Please use 'interpret_cohens_d()'.
## 'interpret_d()' is now deprecated. Please use 'interpret_cohens_d()'.## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer) to predict audpc2 with treat and cut (formula: audpc2 ~ treat * cut). The model included year as random effect (formula: ~1 | year). The model's total explanatory power is substantial (conditional R2 = 0.80) and the part related to the fixed effects alone (marginal R2) is of 0.74. The model's intercept, corresponding to treat = CHK and cut = 0C, is at 11.07 (95% CI [9.40, 12.74], t(82) = 13.16, p < .001). Within this model:
##
## - The effect of treat [Si] is statistically significant and negative (beta = -5.04, 95% CI [-6.59, -3.49], t(82) = -6.48, p < .001; Std. beta = -1.48, 95% CI [-1.93, -1.02])
## - The effect of treat [T] is statistically significant and negative (beta = -4.43, 95% CI [-5.98, -2.88], t(82) = -5.70, p < .001; Std. beta = -1.30, 95% CI [-1.75, -0.85])
## - The effect of treat [TSi] is statistically significant and negative (beta = -5.01, 95% CI [-6.56, -3.47], t(82) = -6.44, p < .001; Std. beta = -1.47, 95% CI [-1.92, -1.02])
## - The effect of cut [1C] is statistically significant and positive (beta = 1.59, 95% CI [0.05, 3.14], t(82) = 2.05, p = 0.044; Std. beta = 0.47, 95% CI [0.01, 0.92])
## - The effect of cut [2C] is statistically significant and negative (beta = -3.84, 95% CI [-5.39, -2.29], t(82) = -4.93, p < .001; Std. beta = -1.13, 95% CI [-1.58, -0.67])
## - The interaction effect of cut [1C] on treat [Si] is statistically non-significant and negative (beta = -1.98, 95% CI [-4.17, 0.21], t(82) = -1.80, p = 0.075; Std. beta = -0.58, 95% CI [-1.22, 0.06])
## - The interaction effect of cut [1C] on treat [T] is statistically significant and negative (beta = -3.98, 95% CI [-6.17, -1.79], t(82) = -3.62, p < .001; Std. beta = -1.17, 95% CI [-1.81, -0.53])
## - The interaction effect of cut [1C] on treat [TSi] is statistically non-significant and negative (beta = -1.14, 95% CI [-3.32, 1.05], t(82) = -1.03, p = 0.305; Std. beta = -0.33, 95% CI [-0.97, 0.31])
## - The interaction effect of cut [2C] on treat [Si] is statistically non-significant and positive (beta = 0.78, 95% CI [-1.41, 2.96], t(82) = 0.70, p = 0.483; Std. beta = 0.23, 95% CI [-0.41, 0.87])
## - The interaction effect of cut [2C] on treat [T] is statistically non-significant and positive (beta = 0.04, 95% CI [-2.15, 2.23], t(82) = 0.04, p = 0.971; Std. beta = 0.01, 95% CI [-0.63, 0.65])
## - The interaction effect of cut [2C] on treat [TSi] is statistically non-significant and positive (beta = 0.10, 95% CI [-2.09, 2.29], t(82) = 0.09, p = 0.925; Std. beta = 0.03, 95% CI [-0.61, 0.67])
##
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed usinglibrary(emmeans)
library(multcompView)
library(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:patchwork':
##
## area## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geysercld(emmeans(m1, ~ cut | treat))## treat = CHK:
## cut emmean SE df lower.CL upper.CL .group
## 2C 7.23 0.842 2.69 4.3706 10.09 1
## 0C 11.07 0.842 2.69 8.2098 13.93 2
## 1C 12.67 0.842 2.69 9.8046 15.53 2
##
## treat = Si:
## cut emmean SE df lower.CL upper.CL .group
## 2C 2.97 0.842 2.69 0.1052 5.83 1
## 1C 5.64 0.842 2.69 2.7821 8.50 2
## 0C 6.03 0.842 2.69 3.1694 8.89 2
##
## treat = T:
## cut emmean SE df lower.CL upper.CL .group
## 2C 2.84 0.842 2.69 -0.0216 5.70 1
## 1C 4.25 0.842 2.69 1.3883 7.11 1
## 0C 6.64 0.842 2.69 3.7770 9.50 2
##
## treat = TSi:
## cut emmean SE df lower.CL upper.CL .group
## 2C 2.32 0.842 2.69 -0.5400 5.18 1
## 0C 6.06 0.842 2.69 3.1955 8.92 2
## 1C 6.52 0.842 2.69 3.6551 9.38 2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.pairs(emmeans(m1, ~ cut | treat))## treat = CHK:
## contrast estimate SE df t.ratio p.value
## 0C - 1C -1.595 0.778 83 -2.049 0.1070
## 0C - 2C 3.839 0.778 83 4.933 <.0001
## 1C - 2C 5.434 0.778 83 6.983 <.0001
##
## treat = Si:
## contrast estimate SE df t.ratio p.value
## 0C - 1C 0.387 0.778 83 0.498 0.8726
## 0C - 2C 3.064 0.778 83 3.937 0.0005
## 1C - 2C 2.677 0.778 83 3.440 0.0026
##
## treat = T:
## contrast estimate SE df t.ratio p.value
## 0C - 1C 2.389 0.778 83 3.069 0.0081
## 0C - 2C 3.799 0.778 83 4.881 <.0001
## 1C - 2C 1.410 0.778 83 1.812 0.1720
##
## treat = TSi:
## contrast estimate SE df t.ratio p.value
## 0C - 1C -0.460 0.778 83 -0.591 0.8256
## 0C - 2C 3.735 0.778 83 4.800 <.0001
## 1C - 2C 4.195 0.778 83 5.391 <.0001
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesEfficacy plots
dat_control <- read_csv("data/dat_efficacy.csv")## Rows: 11 Columns: 4## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Treat
## dbl (3): audpc, audpc2, Efficacy##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.p3 <- dat_control %>%
ggplot(aes(reorder(Treat, Efficacy), Efficacy, color = Efficacy)) +
# geom_col(size = 0.4)+
geom_point(size = 4) +
geom_errorbar(aes(ymin = 0, ymax = Efficacy), width = 0) +
scale_color_gradient_tableau() +
coord_flip() +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
y = "Reduction relative to the uncut check treatment (%)",
x = "",
fill = ""
) +
theme_clean() +
theme(
legend.position = "none",
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
)
ggsave("figs/control.png",
width = 6, height = 4,
bg = "white"
)
p3
## Figure 2
(p2 | p3) +
plot_annotation(tag_levels = "A") +
plot_layout(widths = c(1.5, 1))
ggsave("figs/p2.png", width = 12, height = 5, bg = "white")Foliar Si
dat_si <- read_csv("data/Si_leaf.csv")## Rows: 72 Columns: 4## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Treat, Cut
## dbl (2): Si, year##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Plots
psi1 <- dat_si %>%
filter(year == "2018") %>%
ggplot(aes(Cut, Si, color = Treat)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 0.3) +
facet_wrap(~Treat, ncol = 4) +
theme_clean() +
scale_color_colorblind() +
theme(
legend.position = "none",
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
ylim(20, 65) +
stat_summary(fun.data = "mean_cl_boot", size = 0.3) +
labs(
title = "2018",
x = "Treatment / Cut",
y = "Foliar Silicon (mg Si/kg)"
)
psi11 <- dat_si %>%
filter(year == "2019") %>%
ggplot(aes(Cut, Si, color = Treat)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 0.3) +
facet_wrap(~Treat, ncol = 4) +
theme_clean() +
scale_color_colorblind() +
theme(
legend.position = "none",
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
ylim(20, 65) +
stat_summary(fun.data = "mean_cl_boot", size = 0.3) +
labs(
title = "2019",
x = "Treatment / Cut", y = "Foliar Silicon (mg Si/kg)"
)Predict Si
m2 <- aov(Si ~ Treat, data = dat_si)
anova(m2)## Analysis of Variance Table
##
## Response: Si
## Df Sum Sq Mean Sq F value Pr(>F)
## Treat 3 3545.3 1181.76 14.929 1.446e-07 ***
## Residuals 68 5382.9 79.16
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m2)## Anova Table (Type II tests)
##
## Response: Si
## Sum Sq Df F value Pr(>F)
## Treat 3545.3 3 14.929 1.446e-07 ***
## Residuals 5382.9 68
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(simulateResiduals(m2))## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
library(ggthemes)
dat_si2 <- read_csv("data/Si_soil.csv")## Rows: 72 Columns: 4## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Treat, Cut
## dbl (2): Si, year##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.psi2 <- dat_si2 %>%
filter(year == "2018") %>%
ggplot(aes(Cut, Si, color = Treat)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 0.3) +
facet_wrap(~Treat, ncol = 4) +
theme_clean() +
scale_color_colorblind() +
ylim(0, 10) +
theme(
legend.position = "none",
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
stat_summary(fun.data = "mean_cl_boot", size = 0.3) +
labs(title = "2018", x = "Treatment / Cut", y = "Soil Silicon (mg Si/kg)")
psi3 <- dat_si2 %>%
filter(year == "2019") %>%
ggplot(aes(Cut, Si, color = Treat)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 0.3) +
facet_wrap(~Treat, ncol = 4) +
theme_clean() +
ylim(0, 10) +
scale_color_colorblind() +
theme(
legend.position = "none",
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
stat_summary(fun.data = "mean_cl_boot", size = 0.3) +
labs(
title = "2019",
x = "Treatment / Cut",
y = "Soil Silicon (mg Si/kg)"
)Figure 3
(psi2 + psi3) / (psi1 + psi11) +
plot_annotation(tag_levels = "A")## Warning: Removed 1 rows containing non-finite values (stat_summary).## Warning: Removed 1 rows containing missing values (geom_point).
ggsave("figs/p5.png", width = 8, height = 7)## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Warning: Removed 1 rows containing missing values (geom_point).Weather 2018
library(lubridate)##
## Attaching package: 'lubridate'## The following object is masked from 'package:cowplot':
##
## stamp## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, uniondat_wth <- read_csv("data/wth_2018.csv")## Rows: 275 Columns: 5## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (4): Tmax, Tmin, Rain, RH##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dat_wth <- dat_wth %>%
separate(Date, into = c("Month", "Day")) %>%
mutate(Year = "2018") %>%
unite(date, Year, Month, Day, sep = "-") %>%
mutate(date = as.Date(date))met_2018 <- dat_wth %>%
ggplot(aes(date, RH)) +
scale_x_date(
date_labels = "%b %d",
date_breaks = "2 week"
) +
geom_line(linetype = "dotted") +
theme_clean() +
labs(
y = "Value",
x = "Date",
title = "2018"
) +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
geom_col(aes(x = date, y = Rain), color = "steelblue", fill = "steelblue") +
geom_line(aes(x = date, y = Tmin), color = "navy") +
geom_line(aes(x = date, y = Tmax), color = "darkorange") +
annotate(geom = "text", x = as.Date("2018-12-05"), y = 85, label = "RH") +
annotate(geom = "text", x = as.Date("2018-12-05"), y = 30, label = "TMax") +
annotate(geom = "text", x = as.Date("2018-12-05"), y = 20, label = "TMin") +
annotate(geom = "text", x = as.Date("2018-12-05"), y = 5, label = "Rain")
ggsave("figs/p6.png", width = 12, height = 5)Weather 2019
library(lubridate)
dat_wth2 <- read_csv("data/wth_2019.csv")## Rows: 275 Columns: 5## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (4): Tmax, Tmin, Rain, RH##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dat_wth2 <- dat_wth2 %>%
separate(Date, into = c("Month", "Day")) %>%
mutate(Year = "2019") %>%
unite(date, Year, Month, Day, sep = "-") %>%
mutate(date = as.Date(date))library(ggthemes)
met_2019 <- dat_wth2 %>%
ggplot(aes(date, RH)) +
scale_x_date(
date_labels = "%b %d",
date_breaks = "2 week"
) +
geom_line(linetype = "dotted") +
theme_clean() +
labs(
y = "Value",
x = "Date",
title = "2019"
) +
theme(
plot.background = element_rect(colour = "white"),
legend.background = element_rect(colour = "white")
) +
geom_col(data = dat_wth2, aes(x = date, y = Rain), color = "steelblue", fill = "steelblue") +
geom_line(data = dat_wth2, aes(x = date, y = Tmin), color = "navy") +
geom_line(data = dat_wth2, aes(x = date, y = Tmax), color = "darkorange") +
annotate(geom = "text", x = as.Date("2019-12-05"), y = 85, label = "RH") +
annotate(geom = "text", x = as.Date("2019-12-05"), y = 30, label = "TMax") +
annotate(geom = "text", x = as.Date("2019-12-05"), y = 20, label = "TMin") +
annotate(geom = "text", x = as.Date("2019-12-05"), y = 5, label = "Rain")
ggsave("figs/p7.png", width = 12, height = 5)Figure S3
(met_2018 / met_2019) +
plot_annotation(tag_levels = "A")
ggsave("figs/p8.png", width = 12, height = 7)