Treatment effects on start of phenophase
dat_phenophase
## # A tibble: 4,116,485 × 16
## site canopy heat heat_name water water_name block plot species barcode
## <chr> <chr> <fct> <fct> <fct> <fct> <chr> <chr> <chr> <dbl>
## 1 cfc closed _ ambient _ ambient a a3 abiba 3148
## 2 cfc closed _ ambient _ ambient a a3 abiba 3148
## 3 cfc closed _ ambient _ ambient a a3 abiba 3148
## 4 cfc closed _ ambient _ ambient a a3 abiba 3148
## 5 cfc closed _ ambient _ ambient a a3 abiba 3148
## 6 cfc closed _ ambient _ ambient a a3 abiba 3148
## 7 cfc closed _ ambient _ ambient a a3 abiba 3148
## 8 cfc closed _ ambient _ ambient a a3 abiba 3148
## 9 cfc closed _ ambient _ ambient a a3 abiba 3148
## 10 cfc closed _ ambient _ ambient a a3 abiba 3148
## # ℹ 4,116,475 more rows
## # ℹ 6 more variables: cohort <dbl>, year <dbl>, doy <dbl>, phenophase <fct>,
## # status <dbl>, notes <chr>
dat_phenophase_time <- calc_phenophase_time(dat_phenophase)
dat_phenophase_time$phenophase %>% unique()
## [1] budbreak oneleaf mostleaf senescence leafdrop
## Levels: budbreak oneleaf mostleaf senescence leafdrop
dat_phenophase_diff <- calc_phenophase_diff(dat_phenophase_time)
p_phenophase_diff <- plot_phenophase_diff(dat_phenophase_diff %>%
filter(species %in% c("queru", "aceru")) %>%
mutate(species = factor(species,
levels = c("queru", "aceru"),
labels = c("Red oak", "Red maple")
)))
p_phenophase_diff$budbreak
ggsave(filename = "alldata/output/figures/budbreak_effect.png", p_phenophase_diff$budbreak)
p_phenophase_diff$senescence
ggsave(filename = "alldata/output/figures/senescence_effect.png", p_phenophase_diff$senescence)
Statistical test with LME.
mod <- lmerTest::lmer(diff ~ heat_name * water_name + (1 | species), data = dat_phenophase_diff %>%
filter(phenophase == "budbreak"))
summary(mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: diff ~ heat_name * water_name + (1 | species)
## Data: dat_phenophase_diff %>% filter(phenophase == "budbreak")
##
## REML criterion at convergence: 149278.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.3946 -0.3043 -0.0201 0.2236 12.1049
##
## Random effects:
## Groups Name Variance Std.Dev.
## species (Intercept) 1.725 1.314
## Residual 215.653 14.685
## Number of obs: 18175, groups: species, 21
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.1364 0.3663 24.4644 -0.372
## heat_name+1.7C -3.4117 0.2938 18145.2139 -11.611
## heat_name+3.4C -5.8146 0.2959 18146.0389 -19.648
## water_namereduced 0.9203 0.4854 17959.7990 1.896
## heat_name+1.7C:water_namereduced 0.2713 0.6846 18144.4957 0.396
## heat_name+3.4C:water_namereduced 0.6044 0.6783 18145.4162 0.891
## Pr(>|t|)
## (Intercept) 0.713
## heat_name+1.7C <2e-16 ***
## heat_name+3.4C <2e-16 ***
## water_namereduced 0.058 .
## heat_name+1.7C:water_namereduced 0.692
## heat_name+3.4C:water_namereduced 0.373
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ht_+1.7C ht_+3.4C wtr_nm h_+1.7C:
## het_nm+1.7C -0.391
## het_nm+3.4C -0.389 0.483
## watr_nmrdcd -0.233 0.294 0.291
## ht_n+1.7C:_ 0.168 -0.429 -0.207 -0.702
## ht_n+3.4C:_ 0.169 -0.211 -0.436 -0.708 0.502
mod <- lmerTest::lmer(diff ~ heat_name * water_name + (1 | species), data = dat_phenophase_diff %>%
filter(phenophase == "senescence"))
summary(mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: diff ~ heat_name * water_name + (1 | species)
## Data: dat_phenophase_diff %>% filter(phenophase == "senescence")
##
## REML criterion at convergence: 130964.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -11.4087 -0.4974 -0.0068 0.5156 5.6966
##
## Random effects:
## Groups Name Variance Std.Dev.
## species (Intercept) 5.252 2.292
## Residual 114.622 10.706
## Number of obs: 17270, groups: species, 21
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.1607 0.5305 23.0006 -0.303
## heat_name+1.7C 3.4234 0.2189 17244.7658 15.637
## heat_name+3.4C 5.0960 0.2205 17244.7765 23.113
## water_namereduced -0.3740 0.3665 17262.6229 -1.020
## heat_name+1.7C:water_namereduced 0.1076 0.5159 17244.8855 0.209
## heat_name+3.4C:water_namereduced -0.9611 0.5109 17244.7709 -1.881
## Pr(>|t|)
## (Intercept) 0.765
## heat_name+1.7C <2e-16 ***
## heat_name+3.4C <2e-16 ***
## water_namereduced 0.308
## heat_name+1.7C:water_namereduced 0.835
## heat_name+3.4C:water_namereduced 0.060 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ht_+1.7C ht_+3.4C wtr_nm h_+1.7C:
## het_nm+1.7C -0.198
## het_nm+3.4C -0.197 0.476
## watr_nmrdcd -0.118 0.287 0.284
## ht_n+1.7C:_ 0.084 -0.424 -0.202 -0.700
## ht_n+3.4C:_ 0.085 -0.206 -0.432 -0.707 0.502
Load climate data
dat_climate <- read_climate(option = "raw")
dat_climate_daily <- summ_climate_daily(dat_climate)
setwd("phenologyb4warmed/")
usethis::use_data(dat_climate_daily, overwrite = T)
setwd("..")
I saved daily climatic data. Use directly.
dat_climate_spring <- summ_climate_season(dat_climate_daily, date_start = "Mar 15", date_end = "May 15", rainfall = 1)
dat_climate_fall <- summ_climate_season(dat_climate_daily, date_start = "Jun 15", date_end = "Sep 15", rainfall = 0)
dat_climate_annual_early <- summ_climate_season(dat_climate_daily, date_start = "May 16", start_rollback = 1, date_end = "May 15", end_rollback = 0, rainfall = 0.6)
dat_climate_annual_late <- summ_climate_season(dat_climate_daily, date_start = "Sep 16", start_rollback = 1, date_end = "Sep 15", end_rollback = 0, rainfall = 0.6)
Compare experimental and temporal trends
dat_phenophase_time_summ <- summ_phenophase_time(dat_phenophase_time)
Had to remove year 2012 as the phenology looked like outliers given the temperature.
(p_trend1 <- plot_trend(
dat_phenophase = dat_phenophase_time_summ %>%
filter(phenophase == "budbreak") %>%
select(-phenophase) %>%
filter(species == "queru") %>%
filter(year != 2012),
dat_climate = dat_climate_spring,
var = "heat",
xlim = c(0, 10), ylim = c(110, 170)
))
ggsave("alldata/output/figures/trend_queru.png", p_trend1, width = 12, height = 6)
(p_trend2 <- plot_trend(
dat_phenophase = dat_phenophase_time_summ %>%
filter(phenophase == "budbreak") %>%
select(-phenophase) %>%
filter(species == "aceru") %>%
filter(year != 2012),
dat_climate = dat_climate_spring,
var = "heat",
xlim = c(0, 10), ylim = c(100, 160)
))
ggsave("alldata/output/figures/trend_aceru.png", p_trend2, width = 12, height = 6)
(p_trend3 <- plot_trend(
dat_phenophase = dat_phenophase_time_summ %>%
filter(phenophase == "senescence") %>%
select(-phenophase) %>%
filter(species == "queru") %>%
filter(year != 2012),
dat_climate = dat_climate_fall,
var = "heat",
xlim = c(15, 25), ylim = c(220, 319),
xlab = "Fall temperature (°C)", ylab = "Time of senescence (day of year)"
))
ggsave("alldata/output/figures/trend_queru_senescence.png", p_trend3, width = 12, height = 6)
(p_trend4 <- plot_trend(
dat_phenophase = dat_phenophase_time_summ %>%
filter(phenophase == "budbreak") %>%
select(-phenophase) %>%
filter(species == "queru") %>%
filter(year != 2012),
dat_climate = dat_climate_annual_early,
var = "water",
xlab = "Annual rainfall (mm)", ylab = "Time of budbreak (day of year)"
))
ggsave("alldata/output/figures/trend_queru_rainfall.png", p_trend4, width = 12, height = 6)
(p_trend5 <- plot_trend(
dat_phenophase = dat_phenophase_time_summ %>%
filter(phenophase == "senescence") %>%
select(-phenophase) %>%
filter(species == "queru") %>%
filter(year != 2012),
dat_climate = dat_climate_annual_late,
var = "water",
xlab = "Annual rainfall (mm)", ylab = "Time of senescence (day of year)"
))
ggsave("alldata/output/figures/trend_queru_rainfall_senescence.png", p_trend5, width = 12, height = 6)