Visualize differences
dat_phenophase_diff <- calc_phenophase_diff(dat_phenophase_time_clean)
plot_phenophase_diff(dat_phenophase_diff %>% filter(species == "queru", phenophase %in% c("senescence", "leafdrop")) %>% tidy_phenophase_name(), simple = F)
plot_phenophase_diff(dat_phenophase_diff %>% filter(phenophase == "senescence") %>% tidy_species_name() %>% tidy_phenophase_name(), simple = T)
plot_phenophase_diff(dat_phenophase_diff %>% filter(phenophase == "leafdrop") %>% tidy_species_name() %>% tidy_phenophase_name(), simple = T)
Bayesian LME
dat_phenophase_time <- dat_phenophase_time_clean %>%
group_by(species) %>%
mutate(group = str_c(site, year, sep = "_") %>% factor() %>% as.integer()) %>%
ungroup() %>%
tidy_treatment_code()
\begin{aligned} y_{i,t,s} &\sim \mathcal{N}(\mu_{i,t,s}, \sigma^2) \newline \mu_{i,t,s} &= \mu + \beta_1 W_i + \beta_2 D_i + \beta_3 (W_i \times D_i) + \beta_4 C_i + \beta_5 (W_i \times C_i) + \alpha_{t,s} \newline \alpha_{t,s} &\sim \mathcal{N}(0, \tau^2) \newline \newline \mu &\sim \mathcal{U}(181, 365) \newline \beta_{1:5} &\sim \mathcal{N}(0, 10) \newline \sigma &\sim \mathcal{N}(0, 4) \newline \tau &\sim \mathcal{N}(0, 4) \end{aligned}
test_phenophase_lme(data = dat_phenophase_time, path = "alldata/intermediate/phenophase/uni/", season = "fall", chains = 1, iter = 4000, version = 2, num_cores = 5)
tidy_stanfit_all(path = "alldata/intermediate/phenophase/uni/", season = "fall", num_cores = 20)
calc_bayes_derived(path = "alldata/intermediate/phenophase/uni/", season = "fall", type = "phenophase", random = F, num_cores = 20)
df_lme_all <- read_bayes_all(path = "alldata/intermediate/phenophase/uni/", season = "fall", full_factorial = T, derived = F) %>%
tidy_phenophase_name(season = "fall")
p_lme_diagnostics <- plot_bayes_diagnostics(df_MCMC = df_lme_all %>% filter(species == "queru"), plot_corr = F)
p_lme_diagnostics$p_MCMC
p_lme_diagnostics$p_posterior
df_lme_all_derived <- read_bayes_all(path = "alldata/intermediate/phenophase/uni/", season = "fall", full_factorial = T, derived = T, tidy_mcmc = T) %>%
tidy_species_name() %>%
tidy_phenophase_name(season = "fall")
p_lme_summ <- plot_bayes_summary(df_lme_all_derived, plot_corr = F)
p_lme_summ$p_coef_line
Fall vs. spring phenophases
dat_phenophase_time_clean %>%
tidy_phenophase_name() %>%
tidy_species_name() %>%
plot_spring_fall()
## $`leaf coloring`
##
## $`leaf drop`
