Last updated: 2022-01-28

Checks: 7 0

Knit directory: emlr_obs_analysis/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20210412) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 0a2d89a. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    data/
    Ignored:    output/other/
    Ignored:    output/publication/

Unstaged changes:
    Modified:   analysis/_site.yml
    Modified:   code/Workflowr_project_managment.R

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/vif_budgets.Rmd) and HTML (docs/vif_budgets.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
html 9753eb8 jens-daniel-mueller 2022-01-26 Build site.
html f347cd7 jens-daniel-mueller 2022-01-18 Build site.
html 513630f jens-daniel-mueller 2022-01-18 Build site.
html 1ff6966 jens-daniel-mueller 2022-01-18 Build site.
Rmd bcff7cc jens-daniel-mueller 2022-01-18 predictor analysis in child doc
html a0775c9 jens-daniel-mueller 2022-01-18 Build site.
Rmd e4051d0 jens-daniel-mueller 2022-01-18 included predictor analysis
html e2ebb58 jens-daniel-mueller 2022-01-18 Build site.
Rmd d3fbfe8 jens-daniel-mueller 2022-01-18 ordered vif_max factor levels
html d7dfc7c jens-daniel-mueller 2022-01-18 Build site.
Rmd f441bed jens-daniel-mueller 2022-01-18 rebuild without any gap filter or flagging exceptions
Rmd 53dee50 jens-daniel-mueller 2022-01-12 rebuild with correct config
html 7e0a36b jens-daniel-mueller 2021-11-21 Build site.
html e505a4b jens-daniel-mueller 2021-11-09 Build site.
html f7c3da2 jens-daniel-mueller 2021-11-03 Build site.
html e534f51 jens-daniel-mueller 2021-11-02 Build site.
html 57cfc36 jens-daniel-mueller 2021-11-01 Build site.
html 4331a22 jens-daniel-mueller 2021-10-29 Build site.
html ae5ae64 jens-daniel-mueller 2021-10-26 Build site.
html 581baa0 jens-daniel-mueller 2021-10-07 Build site.
html a7af62f jens-daniel-mueller 2021-10-06 Build site.
html f9b4f93 jens-daniel-mueller 2021-10-05 Build site.
html 960d158 jens-daniel-mueller 2021-09-29 Build site.
html 0573621 jens-daniel-mueller 2021-09-29 Build site.
html 73154c7 jens-daniel-mueller 2021-09-22 Build site.
Rmd 8f1a837 jens-daniel-mueller 2021-09-22 new vif analysis added
html ec4f702 jens-daniel-mueller 2021-09-22 Build site.
Rmd 1abd051 jens-daniel-mueller 2021-09-22 rebuildt with vif results
Rmd f474699 jens-daniel-mueller 2021-08-17 ensemble filtering
html bd9ce6b jens-daniel-mueller 2021-08-02 Build site.
Rmd 29564fb jens-daniel-mueller 2021-08-02 rebuildt with vif analysis

version_id_pattern <- "v"
config <- "vif_max"

1 Read files

# identify required version IDs

Version_IDs_1 <- list.files(path = "/nfs/kryo/work/jenmueller/emlr_cant/observations",
                            pattern = paste0("v_1", "v"))

Version_IDs_2 <- list.files(path = "/nfs/kryo/work/jenmueller/emlr_cant/observations",
                            pattern = paste0("v_2", "v"))

Version_IDs_3 <- list.files(path = "/nfs/kryo/work/jenmueller/emlr_cant/observations",
                            pattern = paste0("v_3", "v"))

Version_IDs <- c(Version_IDs_1, Version_IDs_2, Version_IDs_3)

print(Version_IDs)
 [1] "v_1v01" "v_1v02" "v_1v03" "v_1v04" "v_1v05" "v_1v06" "v_2v01" "v_2v02"
 [9] "v_2v03" "v_2v04" "v_2v05" "v_2v06" "v_3v01" "v_3v02" "v_3v03" "v_3v04"
[17] "v_3v05" "v_3v06"

1.1 Global

for (i_Version_IDs in Version_IDs) {
  # i_Version_IDs <- Version_IDs[1]
  
  print(i_Version_IDs)
  
  path_version_data     <-
    paste(path_observations,
          i_Version_IDs,
          "/data/",
          sep = "")
  
  # load and join data files
  
  dcant_budget_global <-
    read_csv(paste(path_version_data,
                   "dcant_budget_global.csv",
                   sep = ""))
  
  dcant_budget_global_mod_truth <-
    read_csv(paste(
      path_version_data,
      "dcant_budget_global_mod_truth.csv",
      sep = ""
    ))
  
  dcant_budget_global_bias <-
    read_csv(paste(path_version_data,
                   "dcant_budget_global_bias.csv",
                   sep = ""))
  
  lm_best_predictor_counts <-
    read_csv(paste(path_version_data,
                   "lm_best_predictor_counts.csv",
                   sep = ""))
  
  lm_best_dcant <-
    read_csv(paste(path_version_data,
                   "lm_best_dcant.csv",
                   sep = ""))
  
  dcant_budget_global <- bind_rows(dcant_budget_global,
                                      dcant_budget_global_mod_truth)
  
  dcant_budget_global <- dcant_budget_global %>%
    mutate(Version_ID = i_Version_IDs)
  
  dcant_budget_global_bias <- dcant_budget_global_bias %>%
    mutate(Version_ID = i_Version_IDs)
  
  lm_best_predictor_counts <- lm_best_predictor_counts %>%
    mutate(Version_ID = i_Version_IDs)
  
  lm_best_dcant <- lm_best_dcant %>%
    mutate(Version_ID = i_Version_IDs)

  params_local <-
    read_rds(paste(path_version_data,
                   "params_local.rds",
                   sep = ""))
  
  params_local <- bind_cols(
    Version_ID = i_Version_IDs,
    vif_max := str_c(params_local$vif_max, collapse = "|"),
    tref1 = params_local$tref1,
    tref2 = params_local$tref2)
  
  tref <- read_csv(paste(path_version_data,
                         "tref.csv",
                         sep = ""))
  
  params_local <- params_local %>%
    mutate(
      median_year_1 = sort(tref$median_year)[1],
      median_year_2 = sort(tref$median_year)[2],
      duration = median_year_2 - median_year_1,
      period = paste(median_year_1, "-", median_year_2)
    )
  
  if (exists("dcant_budget_global_all")) {
    dcant_budget_global_all <-
      bind_rows(dcant_budget_global_all, dcant_budget_global)
  }
  
  if (!exists("dcant_budget_global_all")) {
    dcant_budget_global_all <- dcant_budget_global
  }
  
  if (exists("dcant_budget_global_bias_all")) {
    dcant_budget_global_bias_all <-
      bind_rows(dcant_budget_global_bias_all,
                dcant_budget_global_bias)
  }

  if (!exists("dcant_budget_global_bias_all")) {
    dcant_budget_global_bias_all <- dcant_budget_global_bias
  }
  
    
  if (exists("lm_best_predictor_counts_all")) {
    lm_best_predictor_counts_all <-
      bind_rows(lm_best_predictor_counts_all, lm_best_predictor_counts)
  }
  
  if (!exists("lm_best_predictor_counts_all")) {
    lm_best_predictor_counts_all <- lm_best_predictor_counts
  }
    
  if (exists("lm_best_dcant_all")) {
    lm_best_dcant_all <-
      bind_rows(lm_best_dcant_all, lm_best_dcant)
  }
  
  if (!exists("lm_best_dcant_all")) {
    lm_best_dcant_all <- lm_best_dcant
  }
  
  if (exists("params_local_all")) {
    params_local_all <- bind_rows(params_local_all, params_local)
  }
  
  if (!exists("params_local_all")) {
    params_local_all <- params_local
  }
  
  
}
[1] "v_1v01"
[1] "v_1v02"
[1] "v_1v03"
[1] "v_1v04"
[1] "v_1v05"
[1] "v_1v06"
[1] "v_2v01"
[1] "v_2v02"
[1] "v_2v03"
[1] "v_2v04"
[1] "v_2v05"
[1] "v_2v06"
[1] "v_3v01"
[1] "v_3v02"
[1] "v_3v03"
[1] "v_3v04"
[1] "v_3v05"
[1] "v_3v06"
rm(
  dcant_budget_global,
  dcant_budget_global_bias,
  dcant_budget_global_mod_truth,
  lm_best_predictor_counts,
  lm_best_dcant,
  params_local,
  tref
)

1.2 Basins

# Version_IDs <- Version_IDs[1:length(Version_IDs)-1]

for (i_Version_IDs in Version_IDs) {
  # i_Version_IDs <- Version_IDs[1]
  
  print(i_Version_IDs)
  
  path_version_data     <-
    paste(path_observations,
          i_Version_IDs,
          "/data/",
          sep = "")
  
  # load and join data files
  
  dcant_budget_basin_AIP <-
    read_csv(paste(path_version_data,
                   "dcant_budget_basin_AIP.csv",
                   sep = ""))
  
  dcant_budget_basin_AIP_mod_truth <-
    read_csv(paste(
      path_version_data,
      "dcant_budget_basin_AIP_mod_truth.csv",
      sep = ""
    ))
  
    
  dcant_budget_basin_AIP <- bind_rows(dcant_budget_basin_AIP,
                                      dcant_budget_basin_AIP_mod_truth)
  
  dcant_budget_basin_AIP_bias <-
    read_csv(paste(path_version_data,
                   "dcant_budget_basin_AIP_bias.csv",
                   sep = ""))
  
  dcant_slab_budget_bias <-
    read_csv(paste0(path_version_data,
                    "dcant_slab_budget_bias.csv"))

  dcant_slab_budget <-
    read_csv(paste0(path_version_data,
                    "dcant_slab_budget.csv"))

  dcant_budget_basin_AIP <- dcant_budget_basin_AIP %>%
    mutate(Version_ID = i_Version_IDs)
  
  dcant_budget_basin_AIP_bias <- dcant_budget_basin_AIP_bias %>%
    mutate(Version_ID = i_Version_IDs)
  
  dcant_slab_budget <- dcant_slab_budget %>%
    mutate(Version_ID = i_Version_IDs)
  
  dcant_slab_budget_bias <- dcant_slab_budget_bias %>%
    mutate(Version_ID = i_Version_IDs)
  
  if (exists("dcant_budget_basin_AIP_all")) {
    dcant_budget_basin_AIP_all <-
      bind_rows(dcant_budget_basin_AIP_all, dcant_budget_basin_AIP)
  }
  
  if (!exists("dcant_budget_basin_AIP_all")) {
    dcant_budget_basin_AIP_all <- dcant_budget_basin_AIP
  }
  
  if (exists("dcant_budget_basin_AIP_bias_all")) {
    dcant_budget_basin_AIP_bias_all <-
      bind_rows(dcant_budget_basin_AIP_bias_all,
                dcant_budget_basin_AIP_bias)
  }
  
  if (!exists("dcant_budget_basin_AIP_bias_all")) {
    dcant_budget_basin_AIP_bias_all <- dcant_budget_basin_AIP_bias
  }
  
  if (exists("dcant_slab_budget_all")) {
    dcant_slab_budget_all <-
      bind_rows(dcant_slab_budget_all, dcant_slab_budget)
  }
  
  if (!exists("dcant_slab_budget_all")) {
    dcant_slab_budget_all <- dcant_slab_budget
  }
  
  if (exists("dcant_slab_budget_bias_all")) {
    dcant_slab_budget_bias_all <-
      bind_rows(dcant_slab_budget_bias_all,
                dcant_slab_budget_bias)
  }
  
  if (!exists("dcant_slab_budget_bias_all")) {
    dcant_slab_budget_bias_all <- dcant_slab_budget_bias
  }
  
}
[1] "v_1v01"
[1] "v_1v02"
[1] "v_1v03"
[1] "v_1v04"
[1] "v_1v05"
[1] "v_1v06"
[1] "v_2v01"
[1] "v_2v02"
[1] "v_2v03"
[1] "v_2v04"
[1] "v_2v05"
[1] "v_2v06"
[1] "v_3v01"
[1] "v_3v02"
[1] "v_3v03"
[1] "v_3v04"
[1] "v_3v05"
[1] "v_3v06"
rm(
  dcant_budget_basin_AIP,
  dcant_budget_basin_AIP_bias,
  dcant_budget_basin_AIP_mod_truth,
  dcant_slab_budget,
  dcant_slab_budget_bias
)

1.3 Basins hemisphere

# Version_IDs <- Version_IDs[1:length(Version_IDs)-1]

for (i_Version_IDs in Version_IDs) {
  # i_Version_IDs <- Version_IDs[1]
  
  print(i_Version_IDs)
  
  path_version_data     <-
    paste(path_observations,
          i_Version_IDs,
          "/data/",
          sep = "")
  
  # load and join data files
  
  dcant_budget_basin_MLR <-
    read_csv(paste(path_version_data,
                   "dcant_budget_basin_MLR.csv",
                   sep = ""))
  
  dcant_budget_basin_MLR_mod_truth <-
    read_csv(paste(
      path_version_data,
      "dcant_budget_basin_MLR_mod_truth.csv",
      sep = ""
    ))
  
    
  dcant_budget_basin_MLR <- bind_rows(dcant_budget_basin_MLR,
                                      dcant_budget_basin_MLR_mod_truth)
  
  # dcant_budget_basin_MLR_bias <-
  #   read_csv(paste(path_version_data,
  #                  "dcant_budget_basin_MLR_bias.csv",
  #                  sep = ""))
  
  dcant_budget_basin_MLR <- dcant_budget_basin_MLR %>%
    mutate(Version_ID = i_Version_IDs)
  
  # dcant_budget_basin_MLR_bias <- dcant_budget_basin_MLR_bias %>%
  #   mutate(Version_ID = i_Version_IDs)

  if (exists("dcant_budget_basin_MLR_all")) {
    dcant_budget_basin_MLR_all <-
      bind_rows(dcant_budget_basin_MLR_all, dcant_budget_basin_MLR)
  }
  
  if (!exists("dcant_budget_basin_MLR_all")) {
    dcant_budget_basin_MLR_all <- dcant_budget_basin_MLR
  }
  
  # if (exists("dcant_budget_basin_MLR_bias_all")) {
  #   dcant_budget_basin_MLR_bias_all <-
  #     bind_rows(dcant_budget_basin_MLR_bias_all,
  #               dcant_budget_basin_MLR_bias)
  # }
  # 
  # if (!exists("dcant_budget_basin_MLR_bias_all")) {
  #   dcant_budget_basin_MLR_bias_all <- dcant_budget_basin_MLR_bias
  # }
  
}
[1] "v_1v01"
[1] "v_1v02"
[1] "v_1v03"
[1] "v_1v04"
[1] "v_1v05"
[1] "v_1v06"
[1] "v_2v01"
[1] "v_2v02"
[1] "v_2v03"
[1] "v_2v04"
[1] "v_2v05"
[1] "v_2v06"
[1] "v_3v01"
[1] "v_3v02"
[1] "v_3v03"
[1] "v_3v04"
[1] "v_3v05"
[1] "v_3v06"
rm(
  dcant_budget_basin_MLR,
  # dcant_budget_basin_MLR_bias,
  dcant_budget_basin_MLR_mod_truth
)

1.4 Steady state

for (i_Version_IDs in Version_IDs) {
  # i_Version_IDs <- Version_IDs[1]
  
  print(i_Version_IDs)
  
  path_version_data     <-
    paste(path_observations,
          i_Version_IDs,
          "/data/",
          sep = "")
  
  # load and join data files
  
  dcant_obs_budget <-
    read_csv(paste0(path_version_data,
                    "anom_dcant_obs_budget.csv"))
  
  dcant_obs_budget <- dcant_obs_budget %>%
    mutate(Version_ID = i_Version_IDs)
  
  if (exists("dcant_obs_budget_all")) {
    dcant_obs_budget_all <-
      bind_rows(dcant_obs_budget_all, dcant_obs_budget)
  }
  
  if (!exists("dcant_obs_budget_all")) {
    dcant_obs_budget_all <- dcant_obs_budget
  }
  
}
[1] "v_1v01"
[1] "v_1v02"
[1] "v_1v03"
[1] "v_1v04"
[1] "v_1v05"
[1] "v_1v06"
[1] "v_2v01"
[1] "v_2v02"
[1] "v_2v03"
[1] "v_2v04"
[1] "v_2v05"
[1] "v_2v06"
[1] "v_3v01"
[1] "v_3v02"
[1] "v_3v03"
[1] "v_3v04"
[1] "v_3v05"
[1] "v_3v06"
rm(dcant_obs_budget)

1.5 Atm CO2

co2_atm <-
  read_csv(paste(path_preprocessing,
                 "co2_atm.csv",
                 sep = ""))
params_local_all <- params_local_all %>%
  mutate(vif_max = fct_inorder(vif_max))
dcant_budget_global_all <- dcant_budget_global_all %>%
  filter(estimate == "dcant", 
         method == "total") %>% 
  select(-c(estimate, method)) %>% 
  rename(dcant = value)

dcant_budget_global_all_depth <- dcant_budget_global_all

dcant_budget_global_all <- dcant_budget_global_all %>%
  filter(inv_depth == params_global$inventory_depth_standard)

dcant_budget_global_bias_all <- dcant_budget_global_bias_all %>%
  filter(estimate == "dcant") %>%
  select(-c(estimate))

dcant_budget_global_bias_all_depth <- dcant_budget_global_bias_all

dcant_budget_global_bias_all <- dcant_budget_global_bias_all %>%
  filter(inv_depth == params_global$inventory_depth_standard)
dcant_budget_basin_AIP_all <- dcant_budget_basin_AIP_all %>%
  filter(estimate == "dcant", 
         method == "total") %>% 
  select(-c(estimate, method)) %>% 
  rename(dcant = value)

dcant_budget_basin_AIP_all_depth <- dcant_budget_basin_AIP_all

dcant_budget_basin_AIP_all <- dcant_budget_basin_AIP_all %>%
  filter(inv_depth == params_global$inventory_depth_standard)

dcant_budget_basin_AIP_bias_all <- dcant_budget_basin_AIP_bias_all %>%
  filter(estimate == "dcant") %>% 
  select(-c(estimate))

dcant_budget_basin_AIP_bias_all_depth <- dcant_budget_basin_AIP_bias_all

dcant_budget_basin_AIP_bias_all <- dcant_budget_basin_AIP_bias_all %>%
  filter(inv_depth == params_global$inventory_depth_standard)
dcant_budget_basin_MLR_all <- dcant_budget_basin_MLR_all %>%
  filter(estimate == "dcant", 
         method == "total") %>% 
  select(-c(estimate, method)) %>% 
  rename(dcant = value)

# dcant_budget_basin_MLR_all_depth <- dcant_budget_basin_MLR_all

dcant_budget_basin_MLR_all <- dcant_budget_basin_MLR_all %>%
  filter(inv_depth == params_global$inventory_depth_standard)

# dcant_budget_basin_MLR_bias_all <- dcant_budget_basin_MLR_bias_all %>%
#   filter(estimate == "dcant") %>% 
#   select(-c(estimate))
# 
# dcant_budget_basin_MLR_bias_all_depth <- dcant_budget_basin_MLR_bias_all
# 
# dcant_budget_basin_MLR_bias_all <- dcant_budget_basin_MLR_bias_all %>%
#   filter(inv_depth == params_global$inventory_depth_standard)

2 Bias thresholds

global_bias_rel_max <- 12.5
global_bias_rel_max
[1] 12.5
regional_bias_rel_max <- 30
regional_bias_rel_max
[1] 30

3 Global

3.1 Individual cases

3.1.1 Absoulte values

legend_title = expression(atop(Delta * C[ant],
                               (mu * mol ~ kg ^ {
                                 -1
                               })))

dcant_budget_global_all %>%
  ggplot(aes(period, dcant, col = vif_max)) +
  geom_jitter(width = 0.05, height = 0) +
  scale_color_brewer(palette = "Dark2") +
  facet_grid(. ~ data_source) +
  ylim(0,NA) +
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        axis.title.x = element_blank())

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02

3.1.2 Biases

dcant_budget_global_bias_all %>%
  ggplot(aes(period, dcant_bias, col = vif_max)) +
  geom_hline(yintercept = 0) +
  scale_color_brewer(palette = "Dark2") +
  labs(y = expression(atop(Delta * C[ant] ~ bias,
                               (mu * mol ~ kg ^ {-1})))) +
  geom_point()

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
p_global_bias <- 
  dcant_budget_global_bias_all %>%
  ggplot() +
  geom_hline(yintercept = global_bias_rel_max * c(-1,1),
             linetype = 2) +
  geom_hline(yintercept = 0) +
  scale_color_brewer(palette = "Dark2") +
  labs(y = expression(Delta * C[ant] ~ bias ~ ("%")),
       title = "Model-based assesment") +
  theme(axis.title.x = element_blank()) +
  geom_point(aes(period, dcant_bias_rel, col = vif_max),
             alpha = 0.7) +
    theme(axis.text.x = element_text(angle = 45, hjust=1),
        axis.title.x = element_blank())

p_global_bias

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02

4 Basins

4.1 Individual cases

4.1.1 Absoulte values

dcant_budget_basin_AIP_all %>%
  ggplot(aes(period, dcant, col = vif_max)) +
  geom_jitter(width = 0.05, height = 0) +
  scale_color_brewer(palette = "Dark2") +
  facet_grid(basin_AIP ~ data_source) +
  ylim(0,NA) +
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        axis.title.x = element_blank())

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02

4.1.2 Biases

dcant_budget_basin_AIP_bias_all %>%
  ggplot(aes(period, dcant_bias, col=vif_max)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  facet_grid(basin_AIP ~ .)

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
dcant_budget_basin_AIP_bias_all %>%
  ggplot() +
  geom_tile(aes(y = 0, height = regional_bias_rel_max * 2,
                x = "2004 - 2014", width = Inf,
                fill = "bias\nthreshold"), alpha = 0.5) +
  geom_hline(yintercept = 0) +
  scale_fill_manual(values = "grey70", name = "") +
  scale_color_brewer(palette = "Dark2") +
  labs(y = expression(Delta ~ C[ant] ~ bias)) +
  theme(axis.title.x = element_blank()) +
  geom_jitter(aes(period, dcant_bias_rel, col = vif_max),
              width = 0.05, height = 0) +
  facet_grid(. ~ basin_AIP)

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
p_regional_bias <- 
  dcant_budget_basin_AIP_bias_all %>%
  ggplot() +
  geom_hline(yintercept = regional_bias_rel_max * c(-1,1),
             linetype = 2) +
  geom_hline(yintercept = 0) +
  scale_color_brewer(palette = "Dark2") +
  labs(y = expression(Delta * C[ant] ~ bias ~ ("%")),
       title = "Model-based assesment") +
  theme(axis.title.x = element_blank()) +
  geom_point(aes(period, dcant_bias_rel, col = vif_max),
             alpha = 0.7) +
    theme(axis.text.x = element_text(angle = 45, hjust=1),
        axis.title.x = element_blank()) +
  facet_grid(. ~ basin_AIP) +
  theme(
  strip.background = element_blank(),
  strip.text.x = element_blank()
)

p_regional_bias

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

4.2 Slab budgets

4.2.1 Absolute values

dcant_slab_budget_all %>%
  filter(data_source == "obs",
         period != "1994 - 2014") %>% 
  ggplot(aes(vif_max, dcant, fill = gamma_slab)) +
  geom_hline(yintercept = 0, col = "red") +
  geom_col() +
  scale_fill_scico_d(direction = -1) +
  facet_grid(basin_AIP ~ period)

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
dcant_slab_budget_all %>%
  filter(data_source == "obs",
         period != "1994 - 2014") %>%
  group_by(basin_AIP) %>%
  group_split() %>%
  map(
    ~ ggplot(data = .x,
             aes(vif_max, dcant, fill = gamma_slab)) +
      geom_hline(yintercept = 0) +
      geom_col() +
      scale_fill_scico_d(direction = -1) +
      labs(title = paste("data_source:", unique(.x$basin_AIP))) +
      facet_grid(gamma_slab ~ period)
  )
[[1]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

[[2]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

[[3]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

4.2.2 Bias

dcant_slab_budget_bias_all %>%
  filter(period != "1994 - 2014") %>%
  group_by(basin_AIP) %>%
  group_split() %>%
  # head(1) %>% 
  map(
    ~ ggplot(data = .x,
             aes(gamma_slab, dcant_bias, fill = gamma_slab)) +
      geom_col() +
      coord_flip() +
      scale_x_discrete(limits = rev) +
      scale_fill_scico_d(direction = -1) +
      facet_grid(period ~ vif_max) +
      labs(title = paste("data_source:", unique(.x$basin_AIP)))
    )
[[1]]
Warning: Removed 36 rows containing missing values (position_stack).

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

[[2]]
Warning: Removed 24 rows containing missing values (position_stack).

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

[[3]]
Warning: Removed 108 rows containing missing values (position_stack).

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

4.2.3 Spread

dcant_slab_budget_all %>%
  filter(period != "1994 - 2014",
         data_source != "mod_truth") %>%
  group_by(data_source, basin_AIP, gamma_slab, period) %>%
  summarise(dcant_range = max(dcant) - min(dcant)) %>%
  ungroup() %>%
  group_split(basin_AIP) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(gamma_slab, dcant_range, fill = gamma_slab)) +
      geom_col() +
      coord_flip() +
      scale_x_discrete(limits = rev) +
      scale_fill_scico_d(direction = -1) +
      facet_grid(period ~ data_source) +
      labs(title = paste("data_source:", unique(.x$basin_AIP)))
  )
`summarise()` has grouped output by 'data_source', 'basin_AIP', 'gamma_slab'. You can override using the `.groups` argument.
[[1]]

Version Author Date
d7dfc7c jens-daniel-mueller 2022-01-18

[[2]]

Version Author Date
d7dfc7c jens-daniel-mueller 2022-01-18

[[3]]

Version Author Date
d7dfc7c jens-daniel-mueller 2022-01-18

5 Basins hemisphere

5.1 Individual cases

5.1.1 Absoulte values

dcant_budget_basin_MLR_all %>%
  ggplot(aes(period, dcant, col = vif_max)) +
  geom_jitter(width = 0.05, height = 0) +
  scale_color_brewer(palette = "Dark2") +
  facet_grid(basin ~ data_source) +
  ylim(0,NA) +
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        axis.title.x = element_blank())

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18

6 Ensemble

6.1 Global

dcant_budget_global_all_in <- dcant_budget_global_all %>% 
  filter(data_source %in% c("mod", "obs"))

dcant_budget_global_ensemble <- dcant_budget_global_all_in %>% 
  group_by(data_source, period, tref2) %>% 
  summarise(dcant_mean = mean(dcant),
            dcant_sd = sd(dcant),
            dcant_range = max(dcant)- min(dcant)) %>% 
  ungroup()
`summarise()` has grouped output by 'data_source', 'period'. You can override using the `.groups` argument.

6.1.1 Mean

legend_title = expression(Delta * C[ant]~(PgC))

ggplot() +
  geom_col(data = dcant_budget_global_ensemble,
           aes(x = period,
               y = dcant_mean),
           fill = "darkgrey") +
  geom_errorbar(
    data = dcant_budget_global_ensemble,
    aes(
      x = period,
      y = dcant_mean,
      ymax = dcant_mean + dcant_sd,
      ymin = dcant_mean - dcant_sd
    ),
    width = 0.1
  ) +
  geom_point(
    data = dcant_budget_global_all,
    aes(period, dcant, col = vif_max),
    alpha = 0.7,
    position = position_jitter(width = 0.2, height = 0)
  ) +
  scale_y_continuous(limits = c(0,70), expand = c(0,0)) +
  scale_color_brewer(palette = "Dark2") +
  facet_grid(. ~ data_source) +
  labs(y = legend_title) +
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        axis.title.x = element_blank())

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
p_global_dcant <- ggplot() +
  geom_col(data = dcant_budget_global_ensemble %>% 
             filter(data_source == "obs"),
           aes(x = period,
               y = dcant_mean),
           fill = "darkgrey") +
    geom_point(
    data = dcant_budget_global_all %>% 
             filter(data_source == "obs"),
    aes(period, dcant, col = vif_max),
    alpha = 0.7,
    position = position_jitter(width = 0.1, height = 0)
  ) +
  geom_errorbar(
    data = dcant_budget_global_ensemble %>% 
             filter(data_source == "obs"),
    aes(
      x = period,
      y = dcant_mean,
      ymax = dcant_mean + dcant_sd,
      ymin = dcant_mean - dcant_sd
    ),
    width = 0.1
  ) +
  scale_y_continuous(limits = c(0,70), expand = c(0,0)) +
  scale_color_brewer(palette = "Dark2") +
  labs(y = legend_title,
       title = "Observation-based results") +
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank())

p_global_dcant_bias <-
p_global_dcant / p_global_bias +
  plot_layout(guides = 'collect',
              heights = c(2,1))

p_global_dcant_bias

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
# ggsave(plot = p_global_dcant_bias,
#        path = here::here("output/publication"),
#        filename = "Fig_global_dcant_budget.png",
#        height = 5,
#        width = 5)

rm(p_global_bias, p_global_dcant, p_global_dcant_bias)

6.1.2 Mean vs atm CO2

dcant_ensemble <- dcant_budget_global_ensemble %>% 
  filter(data_source == "obs",
         period != "1994 - 2014") %>% 
  select(year = tref2, dcant_mean, dcant_sd)

tcant_S04 <- bind_cols(year = 1994, dcant_mean = 118, dcant_sd = 19)

tcant_ensemble <- full_join(dcant_ensemble, tcant_S04)
Joining, by = c("year", "dcant_mean", "dcant_sd")
tcant_ensemble <- left_join(tcant_ensemble, co2_atm)
Joining, by = "year"
co2_atm_pi <- bind_cols(pCO2 = 280, dcant_mean = 0, year = 1750, dcant_sd = 0)

tcant_ensemble <- full_join(tcant_ensemble, co2_atm_pi)
Joining, by = c("year", "dcant_mean", "dcant_sd", "pCO2")
tcant_ensemble <- tcant_ensemble %>% 
  arrange(year) %>% 
  mutate(tcant = cumsum(dcant_mean),
         tcant_sd = cumsum(dcant_sd))

tcant_ensemble %>% 
  ggplot(aes(pCO2, tcant, ymin = tcant - tcant_sd, ymax = tcant + tcant_sd)) +
  geom_ribbon(fill = "grey80") +
  geom_point() +
  geom_line() +
  scale_x_continuous(breaks = seq(280, 400, 30),
                     sec.axis = dup_axis(labels =  c(1750, 1940, 1980, 2000, 2015),
                                         name = "Year")) +
  geom_text(aes(label = year), nudge_x = -5, nudge_y = 5) +
  labs(x = expression(Atmospheric~pCO[2]~(µatm)),
       y = expression(Total~oceanic~C[ant]~(PgC)))

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
# ggsave(path = "output/publication",
#        filename = "Fig_global_dcant_budget_vs_atm_pCO2.png",
#        height = 4,
#        width = 7)

6.1.3 Sum decades

dcant_budget_global_all_in_sum <-
  dcant_budget_global_all_in %>%
  filter(period != "1994 - 2014") %>%
  arrange(tref1) %>%
  group_by(data_source, vif_max) %>%
  mutate(dcant = dcant + lag(dcant)) %>% 
  ungroup() %>%
  drop_na() %>% 
  mutate(estimate = "sum")

dcant_budget_global_all_in_sum <-
  bind_rows(
    dcant_budget_global_all_in_sum,
    dcant_budget_global_all_in %>%
      filter(period == "1994 - 2014") %>%
      mutate(estimate = "direct")
  )

ggplot() +
  geom_point(
    data = dcant_budget_global_all_in_sum,
    aes(estimate, dcant, col = vif_max),
    alpha = 0.7,
    position = position_jitter(width = 0, height = 0)
  ) +
  scale_y_continuous(limits = c(0,70), expand = c(0,0)) +
  scale_color_brewer(palette = "Dark2") +
  facet_grid(. ~ data_source) +
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        axis.title.x = element_blank())

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

6.1.4 Mean bias

dcant_budget_global_ensemble_bias <- full_join(
  dcant_budget_global_ensemble %>%
    filter(data_source == "mod") %>% 
    select(period, dcant_mean, dcant_sd),
  dcant_budget_global_all %>%
    filter(data_source == "mod_truth",
           vif_max == unique(dcant_budget_global_all$vif_max)[1]) %>% 
    select(period, dcant)
)
Joining, by = "period"
dcant_budget_global_ensemble_bias <- dcant_budget_global_ensemble_bias %>% 
  mutate(dcant_mean_bias = dcant_mean - dcant,
         dcant_mean_bias_rel = 100 * dcant_mean_bias / dcant)

dcant_budget_global_ensemble_bias %>%
  ggplot(aes(period, dcant_mean_bias)) +
  geom_hline(yintercept = 0) +
  geom_point()

Version Author Date
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
dcant_budget_global_ensemble_bias %>%
  ggplot(aes(period, dcant_mean_bias_rel)) +
  geom_hline(yintercept = 0) +
  geom_point()

Version Author Date
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02

6.1.5 Vertical patterns

6.1.5.1 Absoulte values

dcant_budget_global_all_depth %>%
  filter(data_source != "mod_truth") %>% 
  group_by(data_source) %>%
  group_split() %>%
  # head(1) %>%
  map(
    ~  ggplot(data = .x,
              aes(dcant, vif_max, fill=period)) +
      geom_vline(xintercept = 0) +
      geom_col(position = "dodge") +
      scale_fill_brewer(palette = "Dark2") +
      facet_grid(inv_depth ~ .) +
      labs(title = paste("data_source:", unique(.x$data_source)))
  )
[[1]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02

[[2]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02

6.1.5.2 Biases

dcant_budget_global_bias_all_depth %>%
  ggplot(aes(dcant_bias, vif_max, fill = period)) +
  geom_vline(xintercept = 0) +
  geom_col(position = "dodge") +
  scale_fill_brewer(palette = "Dark2") +
  facet_grid(inv_depth ~ .)

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
dcant_budget_global_bias_all_depth %>%
  ggplot(aes(dcant_bias_rel, vif_max, fill = period)) +
  geom_vline(xintercept = 0) +
  geom_col(position = "dodge") +
  scale_fill_brewer(palette = "Dark2") +
  facet_grid(inv_depth ~ .)

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
rm(dcant_budget_global_all,
   dcant_budget_global_all_depth,
   dcant_budget_global_bias_all,
   dcant_budget_global_bias_all_depth,
   dcant_budget_global_ensemble,
   dcant_budget_global_ensemble_bias)

6.2 Basins

dcant_budget_basin_AIP_ensemble <- dcant_budget_basin_AIP_all %>% 
  filter(data_source %in% c("mod", "obs")) %>% 
  group_by(basin_AIP, data_source, period) %>% 
  summarise(dcant_mean = mean(dcant),
            dcant_sd = sd(dcant),
            dcant_range = max(dcant)- min(dcant)) %>% 
  ungroup()
`summarise()` has grouped output by 'basin_AIP', 'data_source'. You can override using the `.groups` argument.

6.2.1 Mean

dcant_budget_basin_AIP_ensemble %>%
  ggplot(aes(period, dcant_mean, col=basin_AIP)) +
  geom_pointrange(aes(ymax = dcant_mean + dcant_sd,
                      ymin = dcant_mean - dcant_sd),
                  shape = 21) +
  facet_grid(. ~ data_source)

Version Author Date
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
p_regional_dcant <- ggplot() +
  geom_col(
    data = dcant_budget_basin_AIP_ensemble %>%
      filter(data_source == "obs"),
    aes(x = period,
        y = dcant_mean),
    fill = "darkgrey"
  ) +
  geom_point(
    data = dcant_budget_basin_AIP_all %>%
      filter(data_source == "obs"),
    aes(period, dcant, col = vif_max),
    position = position_jitter(width = 0.1, height = 0),
    alpha = 0.7
  ) +
  geom_errorbar(
    data = dcant_budget_basin_AIP_ensemble %>%
      filter(data_source == "obs"),
    aes(
      x = period,
      y = dcant_mean,
      ymax = dcant_mean + dcant_sd,
      ymin = dcant_mean - dcant_sd
    ),
    width = 0.1
  ) +
  scale_y_continuous(limits = c(0, 35), expand = c(0, 0)) +
  scale_color_brewer(palette = "Dark2") +
  labs(y = legend_title,
       title = "Observation-based results") +
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank()) +
  facet_grid(. ~ basin_AIP)

p_regional_dcant_bias <-
p_regional_dcant / p_regional_bias +
  plot_layout(guides = 'collect',
              heights = c(2,1))

p_regional_dcant_bias
Warning: Removed 1 rows containing missing values (geom_point).

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
# ggsave(plot = p_regional_dcant_bias,
#        path = "output/publication",
#        filename = "Fig_regional_dcant_budget.png",
#        height = 5,
#        width = 10)

rm(p_regional_bias, p_regional_dcant, p_regional_dcant_bias)

6.2.2 Mean bias

dcant_budget_basin_AIP_ensemble_bias <- full_join(
  dcant_budget_basin_AIP_ensemble %>%
    filter(data_source == "mod") %>% 
    select(basin_AIP, period, dcant_mean, dcant_sd),
  dcant_budget_basin_AIP_all %>%
    filter(data_source == "mod_truth",
           vif_max == unique(dcant_budget_basin_AIP_all$vif_max)[1]) %>% 
    select(basin_AIP, period, dcant)
)
Joining, by = c("basin_AIP", "period")
dcant_budget_basin_AIP_ensemble_bias <- dcant_budget_basin_AIP_ensemble_bias %>% 
  mutate(dcant_mean_bias = dcant_mean - dcant,
         dcant_mean_bias_rel = 100 * dcant_mean_bias / dcant)


dcant_budget_basin_AIP_ensemble_bias %>%
  ggplot(aes(period, dcant_mean_bias, col = basin_AIP)) +
  geom_hline(yintercept = 0) +
  geom_point()

Version Author Date
d7dfc7c jens-daniel-mueller 2022-01-18
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
dcant_budget_basin_AIP_ensemble_bias %>%
  ggplot(aes(period, dcant_mean_bias_rel, col = basin_AIP)) +
  geom_hline(yintercept = 0) +
  geom_point()

Version Author Date
d7dfc7c jens-daniel-mueller 2022-01-18

6.2.3 Vertical patterns

6.2.3.1 Absoulte values

dcant_budget_basin_AIP_all_depth %>%
  filter(data_source != "mod_truth") %>%
  group_by(data_source) %>%
  group_split() %>%
  # head(1) %>%
  map(
    ~  ggplot(data = .x,
              aes(dcant, vif_max, fill = basin_AIP)) +
      geom_vline(xintercept = 0) +
      geom_col() +
      scale_fill_brewer(palette = "Dark2") +
      facet_grid(inv_depth ~ period) +
      labs(title = paste("data_source:", unique(.x$data_source)))
  )
[[1]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02

[[2]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

6.2.3.2 Biases

dcant_budget_basin_AIP_bias_all_depth %>%
  ggplot(aes(dcant_bias, vif_max, fill = basin_AIP)) +
  geom_vline(xintercept = 0) +
  geom_col() +
  scale_fill_brewer(palette = "Dark2") +
  facet_grid(inv_depth ~ period)

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02
dcant_budget_basin_AIP_bias_all_depth %>%
  ggplot(aes(dcant_bias_rel, vif_max, fill = basin_AIP)) +
  geom_vline(xintercept = 0) +
  geom_col(position = "dodge") +
  scale_fill_brewer(palette = "Dark2") +
  facet_grid(inv_depth ~ period)

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18
73154c7 jens-daniel-mueller 2021-09-22
ec4f702 jens-daniel-mueller 2021-09-22
bd9ce6b jens-daniel-mueller 2021-08-02

7 Steady state

dcant_obs_budget_all %>%
  group_by(inv_depth) %>%
  group_split() %>%
  # head(1) %>% 
  map(
    ~ ggplot(data = .x,
             aes(estimate, dcant_pos, fill = basin_AIP)) +
      scale_fill_brewer(palette = "Dark2") +
      geom_col() +
      facet_grid(vif_max ~ period) +
      labs(title = paste("inventory depth:",unique(.x$inv_depth)))
  )
[[1]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

[[2]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

[[3]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

[[4]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

[[5]]

Version Author Date
e2ebb58 jens-daniel-mueller 2022-01-18
d7dfc7c jens-daniel-mueller 2022-01-18

8 Predictor analysis

lm_best_predictor_counts_all <-
  full_join(lm_best_predictor_counts_all,
            params_local_all)
Joining, by = "Version_ID"
lm_best_predictor_counts_all <- lm_best_predictor_counts_all %>% 
  mutate(n_predictors_total = rowSums(across(aou:temp), na.rm = TRUE)/10)

lm_best_predictor_counts_all %>%
  ggplot(aes(x = vif_max, y = n_predictors_total)) +
  # ggdist::stat_halfeye(
  #   adjust = .5,
  #   width = .6,
  #   .width = 0,
  #   justification = -.2,
  #   point_colour = NA
  # ) +
  geom_boxplot(width = 0.5,
               outlier.shape = NA) +
  gghalves::geom_half_point(
    side = "l",
    range_scale = .4,
    alpha = .5,
    aes(col = gamma_slab)
  ) +
  scale_color_viridis_d() +
  facet_grid(basin ~ data_source)

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
1ff6966 jens-daniel-mueller 2022-01-18
lm_best_predictor_counts_all %>%
  pivot_longer(aou:temp,
               names_to = "predictor",
               values_to = "count") %>%
  group_split(predictor) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(vif_max, count, color = gamma_slab)) +
      geom_jitter(alpha = 0.5) +
      scale_color_viridis_d() +
      labs(title = paste0("predictor:", unique(.x$predictor))) +
      coord_cartesian(ylim = c(0, 10)) +
      facet_grid(basin ~ data_source)
  )
[[1]]
Warning: Removed 18 rows containing missing values (geom_point).

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
1ff6966 jens-daniel-mueller 2022-01-18

[[2]]
Warning: Removed 88 rows containing missing values (geom_point).

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
1ff6966 jens-daniel-mueller 2022-01-18

[[3]]
Warning: Removed 7 rows containing missing values (geom_point).

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
1ff6966 jens-daniel-mueller 2022-01-18

[[4]]
Warning: Removed 21 rows containing missing values (geom_point).

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
1ff6966 jens-daniel-mueller 2022-01-18

[[5]]
Warning: Removed 5 rows containing missing values (geom_point).

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
1ff6966 jens-daniel-mueller 2022-01-18

[[6]]
Warning: Removed 24 rows containing missing values (geom_point).

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
1ff6966 jens-daniel-mueller 2022-01-18

[[7]]
Warning: Removed 8 rows containing missing values (geom_point).

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
1ff6966 jens-daniel-mueller 2022-01-18
lm_best_dcant_all <-
  full_join(lm_best_dcant_all,
            params_local_all)
Joining, by = "Version_ID"
lm_best_dcant_all %>%
  count(basin, data_source, gamma_slab, vif_max, period) %>%
  ggplot(aes(vif_max, n)) +
  geom_jitter(height = 0, alpha = 0.3) +
  facet_grid(basin ~ data_source)

Version Author Date
f347cd7 jens-daniel-mueller 2022-01-18
1ff6966 jens-daniel-mueller 2022-01-18
lm_best_predictor_counts_all <-
  full_join(lm_best_predictor_counts_all,
            params_local_all)

lm_best_predictor_counts_all <- lm_best_predictor_counts_all %>% 
  mutate(n_predictors_total = rowSums(across(aou:temp), na.rm = TRUE)/10)

lm_best_predictor_counts_all %>%
  ggplot(aes(x = vif_max, y = n_predictors_total)) +
  # ggdist::stat_halfeye(
  #   adjust = .5,
  #   width = .6,
  #   .width = 0,
  #   justification = -.2,
  #   point_colour = NA
  # ) +
  geom_boxplot(width = 0.5,
               outlier.shape = NA) +
  gghalves::geom_half_point(
    side = "l",
    range_scale = .4,
    alpha = .5,
    aes(col = gamma_slab)
  ) +
  scale_color_viridis_d() +
  facet_grid(basin ~ data_source)


lm_best_predictor_counts_all %>%
  pivot_longer(aou:temp,
               names_to = "predictor",
               values_to = "count") %>%
  group_split(predictor) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(vif_max, count, color = gamma_slab)) +
      geom_jitter(alpha = 0.5) +
      scale_color_viridis_d() +
      labs(title = paste0("predictor:", unique(.x$predictor))) +
      coord_cartesian(ylim = c(0, 10)) +
      facet_grid(basin ~ data_source)
  )


lm_best_dcant_all <-
  full_join(lm_best_dcant_all,
            params_local_all)


lm_best_dcant_all %>%
  count(basin, data_source, gamma_slab, vif_max, period) %>%
  ggplot(aes(vif_max, n)) +
  geom_jitter(height = 0, alpha = 0.3) +
  facet_grid(basin ~ data_source)

sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: openSUSE Leap 15.3

Matrix products: default
BLAS:   /usr/local/R-4.1.2/lib64/R/lib/libRblas.so
LAPACK: /usr/local/R-4.1.2/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggforce_0.3.3   metR_0.11.0     scico_1.3.0     patchwork_1.1.1
 [5] collapse_1.7.0  forcats_0.5.1   stringr_1.4.0   dplyr_1.0.7    
 [9] purrr_0.3.4     readr_2.1.1     tidyr_1.1.4     tibble_3.1.6   
[13] ggplot2_3.3.5   tidyverse_1.3.1 workflowr_1.7.0

loaded via a namespace (and not attached):
 [1] fs_1.5.2           gghalves_0.1.1     lubridate_1.8.0    bit64_4.0.5       
 [5] RColorBrewer_1.1-2 httr_1.4.2         rprojroot_2.0.2    tools_4.1.2       
 [9] backports_1.4.1    bslib_0.3.1        utf8_1.2.2         R6_2.5.1          
[13] DBI_1.1.2          colorspace_2.0-2   withr_2.4.3        tidyselect_1.1.1  
[17] processx_3.5.2     bit_4.0.4          compiler_4.1.2     git2r_0.29.0      
[21] cli_3.1.1          rvest_1.0.2        xml2_1.3.3         labeling_0.4.2    
[25] sass_0.4.0         scales_1.1.1       checkmate_2.0.0    callr_3.7.0       
[29] digest_0.6.29      rmarkdown_2.11     pkgconfig_2.0.3    htmltools_0.5.2   
[33] highr_0.9          dbplyr_2.1.1       fastmap_1.1.0      rlang_0.4.12      
[37] readxl_1.3.1       rstudioapi_0.13    jquerylib_0.1.4    generics_0.1.1    
[41] farver_2.1.0       jsonlite_1.7.3     vroom_1.5.7        magrittr_2.0.1    
[45] Rcpp_1.0.8         munsell_0.5.0      fansi_1.0.2        lifecycle_1.0.1   
[49] stringi_1.7.6      whisker_0.4        yaml_2.2.1         MASS_7.3-55       
[53] grid_4.1.2         parallel_4.1.2     promises_1.2.0.1   crayon_1.4.2      
[57] haven_2.4.3        hms_1.1.1          knitr_1.37         ps_1.6.0          
[61] pillar_1.6.4       reprex_2.0.1       glue_1.6.0         evaluate_0.14     
[65] getPass_0.2-2      data.table_1.14.2  modelr_0.1.8       vctrs_0.3.8       
[69] tzdb_0.2.0         tweenr_1.0.2       httpuv_1.6.5       cellranger_1.1.0  
[73] gtable_0.3.0       polyclip_1.10-0    assertthat_0.2.1   xfun_0.29         
[77] broom_0.7.11       later_1.3.0        viridisLite_0.4.0  ellipsis_0.3.2    
[81] here_1.0.1