Merge pull request #83 from ronnyhdez/T82

Process to create residuals distribution plots

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

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))
 
@@ -275,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) %>% 
@@ -307,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) %>% 
@@ -330,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) %>% 
@@ -353,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) %>% 
@@ -376,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) %>% 
@@ -394,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) %>%
@@ -408,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) %>%
@@ -436,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))
 
@@ -487,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)))) %>% 
@@ -514,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) %>% 
@@ -537,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) %>% 
@@ -560,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) %>% 
@@ -583,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) %>% 
@@ -601,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) %>%
@@ -615,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) %>%
@@ -643,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))