From 2dbf64c4ca4240d43b67b771f487ee85acf93488 Mon Sep 17 00:00:00 2001 From: ronnyhdez Date: Wed, 23 Aug 2023 17:54:57 -0600 Subject: [PATCH 1/4] Ref #82 process to create residuals distribution plots --- scripts/lm_preparation.R | 33 ++++++++++++++ test_plot_residuals.R | 96 +++++++++++++++++++++++++++++++++------- 2 files changed, 114 insertions(+), 15 deletions(-) diff --git a/scripts/lm_preparation.R b/scripts/lm_preparation.R index 6afd7fa..0ed8ad4 100644 --- a/scripts/lm_preparation.R +++ b/scripts/lm_preparation.R @@ -33,6 +33,14 @@ ind_mich <- michigan_monthly_500 %>% ind_sites <- bind_rows(ind_bar, ind_bor, ind_mich) +# Function to select variables for residuals plots +select_augmented <- function(data, variable) { + data %>% + select(variable, augmented) %>% + unnest(cols = c(augmented)) %>% + select(variable, gpp_dt_vut_ref, .resid) +} + # Linear model for all sites [A | Diff VIs] all_sites_lm <- ind_sites %>% select(-kndvi_mean) %>% @@ -45,6 +53,9 @@ all_sites_lm <- ind_sites %>% augmented = map(fit, augment) ) +all_sites_augmented_monthly <- select_augmented(all_sites_lm, "index") %>% + mutate(site = "All") + all_fit <- all_sites_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -77,6 +88,9 @@ evi_lm <- ind_sites %>% augmented = map(fit, augment) ) +evi_augmented_monthly <- select_augmented(evi_lm, "site") %>% + mutate(index = "EVI") + evi_fit <- evi_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -100,6 +114,9 @@ ndvi_lm <- ind_sites %>% augmented = map(fit, augment) ) +ndvi_augmented_monthly <- select_augmented(ndvi_lm, "site") %>% + mutate(index = "NDVI") + ndvi_fit <- ndvi_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -123,6 +140,9 @@ nirv_lm <- ind_sites %>% augmented = map(fit, augment) ) +nirv_augmented_monthly <- select_augmented(nirv_lm, "site") %>% + mutate(index = "NIRv") + nirv_fit <- nirv_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -146,6 +166,9 @@ cci_lm <- ind_sites %>% augmented = map(fit, augment) ) +cci_augmented_monthly <- select_augmented(cci_lm, "site") %>% + mutate(index = "CCI") + cci_fit <- cci_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -164,6 +187,10 @@ vis_site_glance_monthly <- bind_rows(evi_glance_monthly, nirv_glance_monthly, cci_glance_monthly) +vis_site_augmented_monthly <- bind_rows(evi_augmented_monthly, + ndvi_augmented_monthly, + nirv_augmented_monthly, + cci_augmented_monthly) # # Linear model for all sites kNDVI # kndvi_lm <- ind_sites %>% @@ -202,6 +229,9 @@ all_vis_lm <- ind_sites %>% augmented = map(fit, augment) ) +all_vis_augmented_monthly <- select_augmented(all_vis_lm, "site") %>% + mutate(index = "All") + all_fit <- all_vis_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -233,6 +263,9 @@ all_sites_all_indices <- lm(gpp_dt_vut_ref ~ evi_mean + metrics <- augment(all_sites_all_indices) %>% select(gpp_dt_vut_ref, .resid) +all_sites_all_vis_augmented_monthly <- metrics %>% + mutate(index = "All", site = "All") + rmse <- sqrt(mean((metrics$.resid)^2)) mae <- mean(abs(metrics$.resid)) diff --git a/test_plot_residuals.R b/test_plot_residuals.R index 8f538c4..3178cb7 100644 --- a/test_plot_residuals.R +++ b/test_plot_residuals.R @@ -6,19 +6,18 @@ all_sites_lm library(cowplot) library(ggridges) -all_sites_lm %>% - # mutate(rmse = map_dbl(augmented, ~sqrt(mean((.x$.resid)^2)))) %>% - select(index, augmented) %>% - unnest(cols = c(augmented)) %>% - select(index, gpp_dt_vut_ref, .fitted, .resid) %>% - ggplot(aes(x = .resid, fill = index)) + - geom_density(alpha = .7) + - scale_fill_viridis_d() + - scale_y_continuous(expand = c(0, 0)) + - theme_half_open(12) - +# all_sites_lm %>% +# # mutate(rmse = map_dbl(augmented, ~sqrt(mean((.x$.resid)^2)))) %>% +# select(index, augmented) %>% +# unnest(cols = c(augmented)) %>% +# select(index, gpp_dt_vut_ref, .fitted, .resid) %>% +# ggplot(aes(x = .resid, fill = index)) + +# geom_density(alpha = .7) + +# scale_fill_viridis_d() + +# scale_y_continuous(expand = c(0, 0)) + +# theme_half_open(12) -all_sites_lm %>% +all_sites_lm %>% # mutate(rmse = map_dbl(augmented, ~sqrt(mean((.x$.resid)^2)))) %>% select(index, augmented) %>% unnest(cols = c(augmented)) %>% @@ -32,9 +31,6 @@ all_sites_lm %>% theme_ridges() + theme(legend.position = "none") -# all_sites_lm %>% -# unnest(augmented) %>% glimpse() - all_sites_lm %>% # mutate(rmse = map_dbl(augmented, ~sqrt(mean((.x$.resid)^2)))) %>% select(index, augmented) %>% @@ -45,6 +41,76 @@ all_sites_lm %>% scale_color_viridis_d() + theme_ridges() +evi_lm %>% + # mutate(rmse = map_dbl(augmented, ~sqrt(mean((.x$.resid)^2)))) %>% + select(site, augmented) %>% + unnest(cols = c(augmented)) %>% + select(site, gpp_dt_vut_ref, .fitted, .resid) %>% + rename("residuals" = ".resid") %>% + ggplot(aes(x = residuals, y = site, fill = site)) + + geom_density_ridges() + + geom_vline(xintercept = 0, colour = "#FF5500", + linewidth = 0.7, linetype = "dashed") + + scale_fill_viridis_d() + + theme_ridges() + + theme(legend.position = "none") + +# ------------------------------------------------------------------------------ +# Probando la union de todos los datos: +vis_site_augmented_monthly %>% + bind_rows(all_sites_augmented_monthly, + all_sites_all_vis_augmented_monthly, + all_vis_augmented_monthly) %>% + mutate(index = case_when( + index == "evi_mean" ~ "EVI", + index == "ndvi_mean" ~ "NDVI", + index == "nirv_mean" ~ "NIRv", + index == "cci_mean" ~ "CCI", + .default = index + )) %>% + rename("residuals" = ".resid") %>% + group_by(site, index) %>% + nest() %>% + arrange(site) + + + + + + +residuals_data <- vis_site_augmented_monthly %>% + bind_rows(all_sites_augmented_monthly, + all_sites_all_vis_augmented_monthly, + all_vis_augmented_monthly) %>% + mutate(index = case_when( + index == "evi_mean" ~ "EVI", + index == "ndvi_mean" ~ "NDVI", + index == "nirv_mean" ~ "NIRv", + index == "cci_mean" ~ "CCI", + .default = index + )) %>% + rename("residuals" = ".resid") + +create_residuals_plot <- function(data, site) { + data %>% + filter(site == site) %>% + ggplot(aes(x = residuals, y = index, fill = index)) + + geom_density_ridges() + + geom_vline(xintercept = 0, colour = "#FF5500", + linewidth = 0.7, linetype = "dashed") + + scale_fill_viridis_d() + + theme_ridges() + + theme(legend.position = "none") +} + +# Borden +residuals_data %>% + create_residuals_plot("Borden") + +# Bartlett + + + # Si tengo tiempo, revisar este proceso From c87bf6b3c55859cfec687742f2417f111dd75749 Mon Sep 17 00:00:00 2001 From: ronnyhdez Date: Wed, 23 Aug 2023 19:00:45 -0600 Subject: [PATCH 2/4] Ref #82 function and map to create monthly residuals plots --- test_plot_residuals.R | 54 ++++++++++++++++++++++++++----------------- 1 file changed, 33 insertions(+), 21 deletions(-) diff --git a/test_plot_residuals.R b/test_plot_residuals.R index 3178cb7..69551b9 100644 --- a/test_plot_residuals.R +++ b/test_plot_residuals.R @@ -57,28 +57,28 @@ evi_lm %>% # ------------------------------------------------------------------------------ # Probando la union de todos los datos: -vis_site_augmented_monthly %>% - bind_rows(all_sites_augmented_monthly, - all_sites_all_vis_augmented_monthly, - all_vis_augmented_monthly) %>% - mutate(index = case_when( - index == "evi_mean" ~ "EVI", - index == "ndvi_mean" ~ "NDVI", - index == "nirv_mean" ~ "NIRv", - index == "cci_mean" ~ "CCI", - .default = index - )) %>% - rename("residuals" = ".resid") %>% - group_by(site, index) %>% - nest() %>% - arrange(site) +# vis_site_augmented_monthly %>% +# bind_rows(all_sites_augmented_monthly, +# all_sites_all_vis_augmented_monthly, +# all_vis_augmented_monthly) %>% +# mutate(index = case_when( +# index == "evi_mean" ~ "EVI", +# index == "ndvi_mean" ~ "NDVI", +# index == "nirv_mean" ~ "NIRv", +# index == "cci_mean" ~ "CCI", +# .default = index +# )) %>% +# rename("residuals" = ".resid") %>% +# group_by(site, index) %>% +# nest() %>% +# arrange(site) -residuals_data <- vis_site_augmented_monthly %>% +monthly_residuals_data <- vis_site_augmented_monthly %>% bind_rows(all_sites_augmented_monthly, all_sites_all_vis_augmented_monthly, all_vis_augmented_monthly) %>% @@ -93,21 +93,33 @@ residuals_data <- vis_site_augmented_monthly %>% create_residuals_plot <- function(data, site) { data %>% - filter(site == site) %>% + filter(site == {{site}}) %>% ggplot(aes(x = residuals, y = index, fill = index)) + geom_density_ridges() + geom_vline(xintercept = 0, colour = "#FF5500", linewidth = 0.7, linetype = "dashed") + scale_fill_viridis_d() + + scale_x_continuous(n.breaks = 10) + theme_ridges() + theme(legend.position = "none") } -# Borden -residuals_data %>% - create_residuals_plot("Borden") +monthly_residuals_data %>% + create_residuals_plot("All") + +# Loop over all the sites +sites <- c("Bartlett", "Borden", "Michigan", "All") + +residuals_plots <- map(sites, ~ create_residuals_plot(residuals_data, .x)) + + +plot_grid(residuals_plots[[1]], + residuals_plots[[2]], + residuals_plots[[3]], + residuals_plots[[4]], + nrow = 1) + -# Bartlett From 0f252bc479385ffb6a689ed201a74a102921fe5f Mon Sep 17 00:00:00 2001 From: ronnyhdez Date: Wed, 23 Aug 2023 22:16:52 -0600 Subject: [PATCH 3/4] Ref #82 lm residuals plots --- scripts/lm_preparation.R | 54 ++++++++++++++++++++++- test_plot_residuals.R | 93 +++++++++++++++++++++++++++++----------- 2 files changed, 120 insertions(+), 27 deletions(-) diff --git a/scripts/lm_preparation.R b/scripts/lm_preparation.R index 0ed8ad4..5d42b68 100644 --- a/scripts/lm_preparation.R +++ b/scripts/lm_preparation.R @@ -308,6 +308,9 @@ all_sites_lm <- ind_sites %>% augmented = map(fit, augment) ) +all_sites_augmented_weekly <- select_augmented(all_sites_lm, "index") %>% + mutate(site = "All") + all_fit <- all_sites_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -340,6 +343,9 @@ evi_lm <- ind_sites %>% augmented = map(fit, augment) ) +evi_augmented_weekly <- select_augmented(evi_lm, "site") %>% + mutate(index = "EVI") + evi_fit <- evi_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -363,6 +369,9 @@ ndvi_lm <- ind_sites %>% augmented = map(fit, augment) ) +ndvi_augmented_weekly <- select_augmented(ndvi_lm, "site") %>% + mutate(index = "NDVI") + ndvi_fit <- ndvi_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -386,6 +395,9 @@ nirv_lm <- ind_sites %>% augmented = map(fit, augment) ) +nirv_augmented_weekly <- select_augmented(nirv_lm, "site") %>% + mutate(index = "NIRv") + nirv_fit <- nirv_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -409,6 +421,9 @@ cci_lm <- ind_sites %>% augmented = map(fit, augment) ) +cci_augmented_weekly <- select_augmented(cci_lm, "site") %>% + mutate(index = "CCI") + cci_fit <- cci_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -427,6 +442,11 @@ vis_site_glance_weekly <- bind_rows(evi_glance_weekly, nirv_glance_weekly, cci_glance_weekly) +vis_site_augmented_weekly <- bind_rows(evi_augmented_weekly, + ndvi_augmented_weekly, + nirv_augmented_weekly, + cci_augmented_weekly) + # Linear model for all VI's [C | All VIs] all_vis_lm <- ind_sites %>% select(-kndvi_mean) %>% @@ -441,6 +461,9 @@ all_vis_lm <- ind_sites %>% augmented = map(fit, augment) ) +all_vis_augmented_weekly <- select_augmented(all_vis_lm, "site") %>% + mutate(index = "All") + all_fit <- all_vis_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -469,10 +492,12 @@ all_sites_all_indices <- lm(gpp_dt_vut_ref ~ evi_mean + cci_mean, data = ind_sites) # summary(all_sites_all_indices) - metrics <- augment(all_sites_all_indices) %>% select(gpp_dt_vut_ref, .resid) +all_sites_all_vis_augmented_weekly <- metrics %>% + mutate(index = "All", site = "All") + rmse <- sqrt(mean((metrics$.resid)^2)) mae <- mean(abs(metrics$.resid)) @@ -520,6 +545,9 @@ all_fit <- all_sites_lm %>% select(-data, -fit, -glanced) %>% mutate(site = "All") +all_sites_augmented_daily <- select_augmented(all_sites_lm, "index") %>% + mutate(site = "All") + all_sites_glance_daily <- all_sites_lm %>% mutate(rmse = map_dbl(augmented, ~sqrt(mean((.x$.resid)^2))), mae = map_dbl(augmented, ~mean(abs(.x$.resid)))) %>% @@ -547,6 +575,9 @@ evi_lm <- ind_sites %>% augmented = map(fit, augment) ) +evi_augmented_daily <- select_augmented(evi_lm, "site") %>% + mutate(index = "EVI") + evi_fit <- evi_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -570,6 +601,9 @@ ndvi_lm <- ind_sites %>% augmented = map(fit, augment) ) +ndvi_augmented_daily <- select_augmented(ndvi_lm, "site") %>% + mutate(index = "NDVI") + ndvi_fit <- ndvi_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -593,6 +627,9 @@ nirv_lm <- ind_sites %>% augmented = map(fit, augment) ) +nirv_augmented_daily <- select_augmented(nirv_lm, "site") %>% + mutate(index = "NIRv") + nirv_fit <- nirv_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -616,6 +653,9 @@ cci_lm <- ind_sites %>% augmented = map(fit, augment) ) +cci_augmented_daily <- select_augmented(cci_lm, "site") %>% + mutate(index = "CCI") + cci_fit <- cci_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -634,6 +674,11 @@ vis_site_glance_daily <- bind_rows(evi_glance_daily, nirv_glance_daily, cci_glance_daily) +vis_site_augmented_daily <- bind_rows(evi_augmented_daily, + ndvi_augmented_daily, + nirv_augmented_daily, + cci_augmented_daily) + # Linear model for all VI's [C | All VIs] all_vis_lm <- ind_sites %>% select(-kndvi_mean) %>% @@ -648,6 +693,9 @@ all_vis_lm <- ind_sites %>% augmented = map(fit, augment) ) +all_vis_augmented_daily <- select_augmented(all_vis_lm, "site") %>% + mutate(index = "All") + all_fit <- all_vis_lm %>% unnest(tidied) %>% select(-data, -fit, -glanced) %>% @@ -676,10 +724,12 @@ all_sites_all_indices <- lm(gpp_dt_vut_ref ~ evi_mean + cci_mean, data = ind_sites) # summary(all_sites_all_indices) - metrics <- augment(all_sites_all_indices) %>% select(gpp_dt_vut_ref, .resid) +all_sites_all_vis_augmented_daily <- metrics %>% + mutate(index = "All", site = "All") + rmse <- sqrt(mean((metrics$.resid)^2)) mae <- mean(abs(metrics$.resid)) diff --git a/test_plot_residuals.R b/test_plot_residuals.R index 69551b9..e7fb2f0 100644 --- a/test_plot_residuals.R +++ b/test_plot_residuals.R @@ -76,8 +76,24 @@ evi_lm %>% - +# Function to create residuals plot +create_residuals_plot <- function(data, site) { + data %>% + filter(site == {{site}}) %>% + ggplot(aes(x = residuals, y = index, fill = index)) + + geom_density_ridges() + + geom_vline(xintercept = 0, colour = "#FF5500", + linewidth = 0.7, linetype = "dashed") + + scale_fill_viridis_d() + + # scale_x_continuous(limits = c(-30, 20), n.breaks = 10) + + theme_ridges() + + theme(legend.position = "none") +} +# Vector with sites for looping in the plots +sites <- c("Bartlett", "Borden", "Michigan", "All") + +# Monthly monthly_residuals_data <- vis_site_augmented_monthly %>% bind_rows(all_sites_augmented_monthly, all_sites_all_vis_augmented_monthly, @@ -91,33 +107,60 @@ monthly_residuals_data <- vis_site_augmented_monthly %>% )) %>% rename("residuals" = ".resid") -create_residuals_plot <- function(data, site) { - data %>% - filter(site == {{site}}) %>% - ggplot(aes(x = residuals, y = index, fill = index)) + - geom_density_ridges() + - geom_vline(xintercept = 0, colour = "#FF5500", - linewidth = 0.7, linetype = "dashed") + - scale_fill_viridis_d() + - scale_x_continuous(n.breaks = 10) + - theme_ridges() + - theme(legend.position = "none") -} - -monthly_residuals_data %>% - create_residuals_plot("All") - # Loop over all the sites -sites <- c("Bartlett", "Borden", "Michigan", "All") - -residuals_plots <- map(sites, ~ create_residuals_plot(residuals_data, .x)) +monthly_residuals_plots <- + map(sites, ~ create_residuals_plot(monthly_residuals_data, .x)) + +# Weekly +weekly_residuals_data <- vis_site_augmented_weekly %>% + bind_rows(all_sites_augmented_weekly, + all_sites_all_vis_augmented_weekly, + all_vis_augmented_weekly) %>% + mutate(index = case_when( + index == "evi_mean" ~ "EVI", + index == "ndvi_mean" ~ "NDVI", + index == "nirv_mean" ~ "NIRv", + index == "cci_mean" ~ "CCI", + .default = index + )) %>% + rename("residuals" = ".resid") +# Loop over all the sites +weekly_residuals_plots <- + map(sites, ~ create_residuals_plot(weekly_residuals_data, .x)) + +# Daily +daily_residuals_data <- vis_site_augmented_daily %>% + bind_rows(all_sites_augmented_daily, + all_sites_all_vis_augmented_daily, + all_vis_augmented_daily) %>% + mutate(index = case_when( + index == "evi_mean" ~ "EVI", + index == "ndvi_mean" ~ "NDVI", + index == "nirv_mean" ~ "NIRv", + index == "cci_mean" ~ "CCI", + .default = index + )) %>% + rename("residuals" = ".resid") -plot_grid(residuals_plots[[1]], - residuals_plots[[2]], - residuals_plots[[3]], - residuals_plots[[4]], - nrow = 1) +# Loop over all the sites +daily_residuals_plots <- + map(sites, ~ create_residuals_plot(daily_residuals_data, .x)) + +plot_grid(monthly_residuals_plots[[1]] + labs(title = "monthly 1"), + monthly_residuals_plots[[2]], + monthly_residuals_plots[[3]], + monthly_residuals_plots[[4]] + labs(title = "monthly 4"), + weekly_residuals_plots[[1]], + weekly_residuals_plots[[2]], + weekly_residuals_plots[[3]], + weekly_residuals_plots[[4]], + daily_residuals_plots[[1]], + daily_residuals_plots[[2]], + daily_residuals_plots[[3]], + daily_residuals_plots[[4]], + nrow = 3, + ncol = 4) From 7177264bb7c3527fd9417a3a5f18d2cb3ac3db25 Mon Sep 17 00:00:00 2001 From: ronnyhdez Date: Sat, 26 Aug 2023 16:00:36 -0600 Subject: [PATCH 4/4] Ref #82 final fix on residual plot creation --- chapter_2_lm.qmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chapter_2_lm.qmd b/chapter_2_lm.qmd index 7feda7f..e3de57e 100644 --- a/chapter_2_lm.qmd +++ b/chapter_2_lm.qmd @@ -7,7 +7,7 @@ editor_options: ## Introduction -Gross Primary Production (GPP) is the total amount of carbon fixation by plants +Vegetation Gross Primary Production (GPP) is the total amount of carbon fixation by plants through photosynthesis [@badgley_terrestrial_2019]. Quantifying GPP is essential for understanding land-atmosphere carbon exchange [@kohler_global_2018], ecosystem function, and ecosystem responses to climate change