diff --git a/create_suppl_info.Rmd b/create_suppl_info.Rmd
index 1e21e81..e26f4a1 100644
--- a/create_suppl_info.Rmd
+++ b/create_suppl_info.Rmd
@@ -96,13 +96,128 @@ if (!file.exists(filn)){
 
 Visualise
 ```{r}
-## s_def diagnostic
-df$out[[1]]$gg_fet
+# myplot <- function(n, df){
+#   print(n)  
+#   df$out[[n]]$gg_fet
+# }
+# purrr::map(as.list(seq(36)), ~myplot(., df))
+
+## example for no decline
+gg1 <- df$out[[2]]$gg_fet[[1]] +
+  labs(title = NULL, y = "EF (unitless)") +
+  theme_classic()
+
+## example for S0 diagnosed (9, 34, 21)
+gg2 <- df$out[[34]]$gg_fet[[1]] +
+  labs(title = NULL, y = "EF (unitless)") +
+  theme_classic()
+
+gg3 <- df$out[[21]]$gg_fet[[1]] +
+  labs(title = NULL, y = "EF (unitless)") +
+  theme_classic()
+
+## example for flattening
+gg4 <- df$out[[31]]$gg_fet[[1]] +
+  labs(title = NULL, y = "EF (unitless)") +
+  theme_classic()
+
+cowplot::plot_grid(gg1, gg4, gg2, gg3, nrow = 2, labels = c('a', 'b', 'c', 'd'))
+ggsave("fig/plot_test_s0_diag.pdf", width = 12, height = 8)
+ggsave("fig/plot_test_s0_diag.png", width = 12, height = 8)
+```
+
+Visualise
+```{r}
+tmp <- df %>% 
+  select(-lon, -lat) %>% 
+  unnest(out)
+  mutate(gg = purrr::map(out, ~pull(slice(., 1), gg_fet)),
+         cwd_lue0_fet = purrr::map_dbl(out, ~pull(., cwd_lue0_fet)))
+  
+## example for no decline
+gg1 <- tmp$gg_fet[[2]] +
+  labs(title = NULL, y = "EF (unitless)") +
+  theme_classic()
+
+## example for S0 diagnosed (9, 34, 21)
+gg2 <- tmp$gg_fet[[32]] +
+  labs(title = NULL, y = "EF (unitless)") +
+  theme_classic()
+
+gg3 <- tmp$gg_fet[[21]] +
+  labs(title = NULL, y = "EF (unitless)") +
+  theme_classic()  
+  
+## example for flattening
+gg4 <- tmp$gg_fet[[29]] +
+  labs(title = NULL, y = "EF (unitless)") +
+  theme_classic()  
+
+cowplot::plot_grid(gg1, gg4, gg2, gg3, nrow = 2, labels = c('a', 'b', 'c', 'd'))
+
+## retain selected
+tmp2 <- tmp %>% 
+  slice(c(2, 32, 21, 29))
+
+# ## s_def diagnostic
+# df$out[[1]]$gg_fet
+# 
+# ## water balance time series
+# df$out[[1]]$data[[1]] %>% 
+#   ggplot(aes(time, bal)) + 
+#   geom_line()
+# 
+# df$out[[1]]$data[[1]] %>% 
+#   ggplot(aes(deficit, fet)) + 
+#   geom_point()
+# 
+# df$out[[1]]$data[[1]] %>% 
+#   ggplot(aes(NR, et)) + 
+#   geom_point() +
+#   geom_smooth(method = "lm")
+```
+
+Plot the same from the model outputs.
+```{r}
+# ## determine chunk used for the four sites in tmp2
+# df_sites_ichunk <- read_csv("data/df_sites_rsip.csv") %>%
+#   mutate(idx = 1:n()) %>%
+#   mutate(chunk = rep(1:as.integer(50),
+#                      each = (nrow(.)/as.integer(50)), len = nrow(.)))
+# 
+# tmp2 <- tmp2 %>% 
+#   left_join(
+#     df_sites_ichunk %>% 
+#       select(lon, lat, Dr, chunk)
+#   )
+
+read_myfile <- function(ichunk, use_whc){
+  path <- "./data/out_rsofun_cwdx/"
+  filename <- file.path(path, paste0("out_rsofun_cwdx_whc_", as.character(use_whc), "_ichunk_", as.character(ichunk), ".rds"))
+  if (file.exists(filename)){
+    df <- readRDS(filename)
+    return(df)
+  } else {
+    return(tibble())
+  }
+}
+
+set.seed(1982)
+df_rsofun <- purrr::map_dfr(
+      as.list(sample(seq(50), 3)),
+      ~read_myfile(., "200")
+      ) %>% 
+      mutate(setup = "whc_200")
+
+# ## example of S0 detected
+# df_rsofun$out_cwd_lue0[[20]]$gg +
+#   geom_vline(xintercept = df_rsofun$whc[20])
+
+```
+
+```{r}
+purrr::map(seq(1:30), ~{print(df_rsofun$out_cwd_lue0[[.]]$gg)})
 
-## water balance time series
-df$out[[1]]$data[[1]] %>% 
-  ggplot(aes(time, bal)) + 
-  geom_line()
 ```
 
 ## Magnitudes of EF
diff --git a/rsip.Rmd b/rsip.Rmd
index da7981f..9a6b9b8 100644
--- a/rsip.Rmd
+++ b/rsip.Rmd
@@ -897,18 +897,97 @@ save(gg_rsip_sites_wtd_90, file = "data/gg_rsip_sites_wtd_90.RData")
 
 ## Trait gradient analysis
 
+First, keep only trees, shrubs, and semi-shrubs, and remove the significant effect of shoot height. Note that 3230 data points have missing Hs. The filtering here reduces the dataset to 2287 entries.
+```{r}
+df2 <- df %>% 
+  mutate(Hs = as.numeric(Hs)) %>% 
+  filter(Hs > 0 & Dr > 0) %>% 
+  filter(Growth_form %in% c("Tree", "Semi-shrub", "Shrub"))
+
+df2 %>% 
+  ggplot(aes(Hs, Dr)) +
+  geom_point(alpha = 0.5) +
+  scale_x_log10() +
+  scale_y_log10() +
+  geom_smooth(method = "lm")
+```
+
+```{r}
+linmod <- lm(log(Dr) ~ log(Hs), 
+             data = df2, 
+             na.action = "na.omit")
+```
+
+```{r}
+plot(linmod)
+hist(linmod$residuals)
+```
+Consider residuals for the TGA.
+```{r}
+df2 <- df2 %>% 
+  mutate(Dr_res = linmod$residuals)
+```
+
+
+<!-- Consider Hs-normalised values -->
+<!-- ```{r} -->
+<!-- df %>%  -->
+<!--   mutate(Dr_norm = Dr / Hs^linmod$coefficients["log(Hs)"]) %>%  -->
+<!--   ggplot(aes(Hs, Dr_norm)) + -->
+<!--   geom_point(alpha = 0.5) + -->
+<!--   geom_smooth(method = "lm") -->
+
+<!-- df %>%  -->
+<!--   ggplot(aes(Hs, Dr)) + -->
+<!--   geom_point(alpha = 0.5) + -->
+<!--   geom_smooth(method = "lm") -->
+<!-- ``` -->
+
+Aggregate to sites.
+```{r}
+df_sites2 <- df2 %>% 
+  mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))) %>% 
+  group_by(sitename) %>% 
+  summarise(Dr = mean(Dr, na.rm = TRUE),
+            Dr_res = mean(Dr_res, na.rm = TRUE))
+
+## add WWF biome info
+df_wwf_sites2 <- ingest(
+  df2 %>% 
+    distinct(lon, lat) %>% 
+    mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))),
+  source = "wwf",
+  dir = "~/data/biomes/wwf_ecoregions/official/",
+  settings = list(layer = "wwf_terr_ecos")
+  )%>% 
+  mutate(data = purrr::map(data, ~slice(., 1))) %>% 
+  unnest(data)
+
+df_sites2 <- df_sites2 %>% 
+  left_join(
+    df_wwf_sites2 %>% 
+      dplyr::select(sitename, biome = BIOME, biome_name = BIOME_NAME),
+    by = "sitename"
+  )
+
+df_sites2 %>%
+  separate(sitename, into = c(NA, "lon", "lat"), sep = "_") %>% 
+  mutate(lon = as.numeric(lon), lat = as.numeric(lat)) %>% 
+  write_csv(file = "data/df_sites_rsip_tga.csv")
+```
+
 Add site mean to full data
 ```{r}
-df2 <- df_sites %>% 
-  dplyr::select(sitename, Dr_sitemean = Dr, zroot_cwdx80_sitemean = zroot_cwdx80) %>% 
+df2 <- df_sites2 %>% 
+  dplyr::select(sitename, Dr_sitemean = Dr, Dr_res_sitemean = Dr_res) %>% 
   right_join(
-    df %>% 
+    df2 %>% 
       mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))), 
     by = "sitename"
     )
 ```
 
-Use data only for species that appear in at least N sites.
+Use data only for species that appear in at least 3 sites => 72 species.
 ```{r}
 use_species <- df2 %>%
   filter(Species != "NA") %>% 
@@ -917,32 +996,28 @@ use_species <- df2 %>%
   distinct() %>% 
   group_by(Species) %>% 
   summarise(n = n()) %>% 
-  dplyr::filter(n >= 5) %>% 
+  dplyr::filter(n >= 3) %>% 
   pull(Species)
 ```
 
-Use data only for sites where at least M species are present.
+Use data only for sites where at least 3 species are present => 148 sites.
 ```{r}
 use_sites <- df2 %>% 
   dplyr::select(sitename, Species) %>% 
   distinct() %>% 
   group_by(sitename) %>% 
   summarise(n = n()) %>% 
-  dplyr::filter(n >= 5) %>% 
+  dplyr::filter(n >= 3) %>% 
   pull(sitename)
 ```
 
-subset data to keep only those filtered above and for which at least five points per species remain. 51 species at the end, 402 data points.
+Apply these two criteria. This leaves 286 data points.
 ```{r}
-use_species2 <- df2 %>% 
-  dplyr::filter(Species %in% use_species & sitename %in% use_sites) %>% 
-  group_by(Species) %>% 
-  summarise(n = n()) %>% 
-  dplyr::filter(n >= 5) %>% 
-  pull(Species)
-
 df3 <- df2 %>% 
-  dplyr::filter(Species %in% use_species & sitename %in% use_sites & Species %in% use_species2)
+  dplyr::filter(Species %in% use_species & sitename %in% use_sites)
+  # dplyr::filter(Species %in% use_species2)
+
+saveRDS(df3, file = "data/df_tga.rds")
 ```
 
 
@@ -950,7 +1025,13 @@ Plot.
 ```{r eval=FALSE}
 df3 %>% 
   ggplot(aes(x = Dr_sitemean, y = Dr)) +  # , color = Species
-  geom_point()
+  geom_point() +
+  geom_abline(slope = 1, intercept = 0, linetype = "dotted")
+
+df3 %>% 
+  ggplot(aes(x = Dr_res_sitemean, y = Dr_res)) +  # , color = Species
+  geom_point() +
+  geom_abline(slope = 1, intercept = 0, linetype = "dotted")
 ```
 
 Plot just the lines
@@ -963,18 +1044,31 @@ df3 %>%
   theme_classic()
   # geom_point(alpha = 0.3)
 
+df3 %>% 
+  group_by(Species) %>% 
+  ggplot(aes(x = Dr_res_sitemean, y = Dr_res, group = Species, color = Species)) +
+  geom_smooth(method = "lm", se = FALSE, size = 0.5, alpha = 0.2) +
+  geom_abline(intercept=0, slope=1, linetype="dotted") +
+  theme_classic() +
+  geom_point(alpha = 0.3)
+
 ggsave("fig/tga_rsip.png", width = 8, height = 4)
 ```
 
 
-Fit linear regressions by species.
+Fit linear regressions by species for Dr (not residual).
 ```{r}
+get_width <- function(df){
+  df %>% pull(Dr) %>% range() %>% diff()
+}
+
 df_tga <- df3 %>% 
   dplyr::filter(Species %in% use_species & sitename %in% use_sites) %>% 
   group_by(Species) %>% 
   nest() %>% 
-  mutate(linmod_n = purrr::map(data, ~lm(Dr ~ Dr_sitemean, data = .))) %>% 
-  mutate(slope_n = purrr::map_dbl(linmod_n, ~coef(.)[2])) %>% 
+  mutate(linmod = purrr::map(data, ~lm(Dr ~ Dr_sitemean, data = .)),
+         width = purrr::map_dbl(data, ~get_width(.))) %>% 
+  mutate(slope = purrr::map_dbl(linmod, ~coef(.)[2])) %>% 
   left_join(df3 %>% 
               dplyr::select(Species, Family) %>% 
               distinct(), 
@@ -982,17 +1076,90 @@ df_tga <- df3 %>%
 
 df_tga %>% 
   ggplot() +
-  geom_histogram(aes(slope_n, ..density..), fill = "royalblue", binwidth = 0.4, alpha = 0.5) +
-  geom_density(aes(slope_n, ..density..), color = "royalblue") +
+  geom_histogram(aes(slope, ..density..), fill = "royalblue", binwidth = 0.4, alpha = 0.5) +
+  geom_density(aes(slope, ..density..), color = "royalblue") +
   xlim(-2,3)
 
 df_tga %>% 
+  ggplot(aes(width, slope)) +
+  geom_point() +
+  geom_smooth(method = "lm") +
+  geom_hline(yintercept = 1, linetype = "dotted")
+
+linmod2 <- lm(slope ~ width, data = df_tga)
+summary(linmod2)
+
+df_tga %>% 
+  group_by(Family) %>% 
+  summarise(slope = mean(slope)) %>% 
+  mutate(Family = forcats::fct_reorder(Family, slope)) %>% 
+  drop_na() %>% 
+  ggplot(aes(Family, slope)) +
+  geom_bar(stat = "identity") +
+  coord_flip()
+```
+
+Fit linear regressions by species for `Dr_res` (residual of Dr from model log(Dr) ~ log(Hs)).
+```{r}
+get_width_res <- function(df){
+  df %>% pull(Dr_res) %>% range() %>% diff()
+}
+
+df_tga_res <- df3 %>% 
+  dplyr::filter(Species %in% use_species & sitename %in% use_sites) %>% 
+  group_by(Species) %>% 
+  nest() %>% 
+  mutate(linmod = purrr::map(data, ~lm(Dr_res ~ Dr_res_sitemean, data = .)),
+         width = purrr::map_dbl(data, ~get_width_res(.))) %>% 
+  mutate(slope = purrr::map_dbl(linmod, ~coef(.)[2])) %>% 
+  left_join(df3 %>% 
+              dplyr::select(Species, Family) %>% 
+              distinct(), 
+            by = "Species") %>% 
+  
+  ## remove outlier slope
+  filter(slope < 30)
+
+df_tga_res %>% 
+  ggplot() +
+  geom_histogram(aes(slope, ..density..), fill = "royalblue", binwidth = 0.4, alpha = 0.5) +
+  geom_density(aes(slope, ..density..), color = "royalblue") +
+  xlim(-2,3)
+
+
+df_tga_res %>% 
+  ggplot(aes(width, slope)) +
+  geom_point() +
+  geom_smooth(method = "lm") +
+  geom_hline(yintercept = 1, linetype = "dotted")
+
+df_tga_res %>% 
+  ggplot(aes(width, abs(slope))) +
+  geom_point() +
+  geom_smooth(method = "lm") +
+  geom_hline(yintercept = 1, linetype = "dotted")
+
+linmod2 <- lm(slope ~ width, data = df_tga_res)
+summary(linmod2)
+
+## species-level
+df_tga_res %>% 
+  drop_na() %>% 
+  ggplot(aes(forcats::fct_reorder(Species, slope), slope)) +
+  geom_bar(stat = "identity") +
+  geom_hline(yintercept = 1.0, linetype = "dotted") +
+  coord_flip()
+
+## family level
+df_tga_res %>% 
   group_by(Family) %>% 
-  summarise(slope_n = mean(slope_n)) %>% 
-  mutate(Family = forcats::fct_reorder(Family, slope_n)) %>% 
+  summarise(slope = mean(slope), n = n()) %>% 
+  mutate(family_n = paste0(Family, " (", as.character(n), ")")) %>% 
+  mutate(family_n = forcats::fct_reorder(family_n, slope)) %>% 
   drop_na() %>% 
-  ggplot(aes(Family, slope_n)) +
+  ggplot(aes(family_n, slope)) +
   geom_bar(stat = "identity") +
+  geom_hline(yintercept = 1.0, linetype = "dotted") +
   coord_flip()
 ```
 
diff --git a/rsip_tga.Rmd b/rsip_tga.Rmd
new file mode 100644
index 0000000..3d80ffe
--- /dev/null
+++ b/rsip_tga.Rmd
@@ -0,0 +1,374 @@
+---
+title: "Trait gradient analysis of rooting depth"
+author: "Beni Stocker"
+output:
+  html_document:
+    toc: true
+    toc_depth: 3
+    toc_float: true
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+library(readr)
+library(dplyr)
+library(tidyr)
+library(rbeni)
+library(raster)
+library(tibble)
+library(ingestr)
+library(ggplot2)
+library(ggridges)
+library(ggrepel)
+library(cowplot)
+```
+
+The data is downloaded from the google spreadsheet [RSIP Working Copy], tab 'Analysis sheet', into a CSV (19.9.2019) and saved as `data/RSIP_Analysis_sheet.csv`.
+
+Use `gsheet::gsheet2tbl`
+
+Read the data.
+```{r}
+df <- read_csv("~/data/rootingdepth/rsip/RSIP_Analysis_sheet_210409.csv") %>%      # previously done with "data/RSIP_Analysis_sheet.csv"
+# df <- read_csv("~/data/rootingdepth/rsip/RSIP_Analysis_sheet_210721.csv") %>% 
+  rename(lon = Long, lat = Lat) %>% 
+  rowid_to_column(var = "id") %>% 
+  
+  ## problem: some have a reference error
+  dplyr::filter(lon != "#REF!") %>% 
+  mutate(lon = as.numeric(lon), lat = as.numeric(lat), 
+         Dr = as.numeric(Dr),
+         wtd = as.numeric(Water_Table_Depth_Fan))
+  
+print(paste("Total number of entries:", nrow(df)))
+```
+
+
+## Trait gradient analysis
+
+First, remove data from herbaceous plants that occurr at sites together with woody plants. Group grasses and forbs into "herbaceous" and trees, shrubs, and semi-shrubs into "woody".
+```{r}
+df2 <- df %>% 
+  mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))) %>% 
+  mutate(type = ifelse(Growth_form %in% c("Forb", "Grass"), "herb",
+                       ifelse(Growth_form %in% c("Tree", "Shrub", "Semi-shrub"), "woody", NA))) %>% 
+  filter(!is.na(type))
+
+find_coexisting <- function(vec){
+  vec <- unique(vec)
+  if (length(vec) > 1){
+    return(TRUE)
+  } else {
+    return(FALSE)
+  }
+}
+
+tmp <- df2 %>% 
+  group_by(sitename) %>% 
+  summarise(coexisting = find_coexisting(type))
+
+df3 <- df2 %>% 
+  left_join(tmp, by = "sitename") %>% 
+  filter(!(type == "herb" & coexisting))
+```
+
+## Plot depth heigh relationship
+
+Remove the significant effect of shoot height. Note that 3230 data points have missing Hs. The filtering here reduces the dataset to 2287 entries. 
+```{r}
+df4 <- df3 %>% 
+  mutate(Hs = as.numeric(Hs)) %>% 
+  filter(Hs > 0 & Dr > 0)
+
+df4 %>% 
+  ggplot(aes(Hs, Dr)) +
+  geom_point(alpha = 0.5) +
+  scale_x_log10() +
+  scale_y_log10() +
+  geom_smooth(method = "lm") +
+  theme_classic() +
+  labs(x = expression(italic("H")[s] ~ "(m)"),
+       y = expression(italic("z")[r] ~ "(m)"))
+
+ggsave("fig/height_depth.pdf", width = 6, height = 5)
+ggsave("fig/height_depth.png", width = 6, height = 5)
+```
+
+```{r}
+linmod <- lm(log(Dr) ~ log(Hs), 
+             data = df4, 
+             na.action = "na.omit")
+```
+
+```{r}
+plot(linmod)
+hist(linmod$residuals)
+```
+Consider residuals for the TGA.
+```{r}
+df4 <- df4 %>% 
+  mutate(Dr_res = linmod$residuals)
+```
+
+
+<!-- Consider Hs-normalised values -->
+<!-- ```{r} -->
+<!-- df %>%  -->
+<!--   mutate(Dr_norm = Dr / Hs^linmod$coefficients["log(Hs)"]) %>%  -->
+<!--   ggplot(aes(Hs, Dr_norm)) + -->
+<!--   geom_point(alpha = 0.5) + -->
+<!--   geom_smooth(method = "lm") -->
+
+<!-- df %>%  -->
+<!--   ggplot(aes(Hs, Dr)) + -->
+<!--   geom_point(alpha = 0.5) + -->
+<!--   geom_smooth(method = "lm") -->
+<!-- ``` -->
+
+## Aggregate to sites
+
+```{r}
+df_sites2 <- df4 %>% 
+  mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))) %>% 
+  group_by(sitename) %>% 
+  summarise(Dr = mean(Dr, na.rm = TRUE),
+            Dr_res = mean(Dr_res, na.rm = TRUE))
+
+## add WWF biome info
+df_wwf_sites2 <- ingest(
+  df4 %>% 
+    distinct(lon, lat) %>% 
+    mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))),
+  source = "wwf",
+  dir = "~/data/biomes/wwf_ecoregions/official/",
+  settings = list(layer = "wwf_terr_ecos")
+  )%>% 
+  mutate(data = purrr::map(data, ~slice(., 1))) %>% 
+  unnest(data)
+
+df_sites2 <- df_sites2 %>% 
+  left_join(
+    df_wwf_sites2 %>% 
+      dplyr::select(sitename, biome = BIOME, biome_name = BIOME_NAME),
+    by = "sitename"
+  )
+
+df_sites2 %>%
+  separate(sitename, into = c(NA, "lon", "lat"), sep = "_") %>% 
+  mutate(lon = as.numeric(lon), lat = as.numeric(lat)) %>% 
+  write_csv(file = "data/df_sites_rsip_tga.csv")
+```
+
+Add site mean to full data
+```{r}
+df4 <- df_sites2 %>% 
+  dplyr::select(sitename, Dr_sitemean = Dr, Dr_res_sitemean = Dr_res) %>% 
+  right_join(
+    df4 %>% 
+      mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))), 
+    by = "sitename"
+    )
+```
+
+
+## Filter data
+
+Use data only for sites where at least 3 data points are available. Reduces data from 1,497 to 1,197 points.
+```{r}
+use_sites <- df4 %>% 
+  dplyr::select(sitename, Species) %>% 
+  group_by(sitename) %>% 
+  summarise(n = n()) %>% 
+  dplyr::filter(n >= 3) %>% 
+  pull(sitename)
+
+df5 <- df4 %>% 
+  dplyr::filter(sitename %in% use_sites)
+```
+
+Use data only for species that appear in at least 3 sites => 35 species and 166 data points.
+```{r}
+use_species <- df5 %>%
+  filter(Species != "NA") %>% 
+  dplyr::select(sitename, Species) %>% 
+  distinct() %>% 
+  group_by(Species) %>% 
+  summarise(n = n()) %>% 
+  dplyr::filter(n >= 3) %>% 
+  pull(Species)
+
+df6 <- df5 %>% 
+  dplyr::filter(Species %in% use_species)
+```
+
+```{r}
+# test : number of species per site
+df6 %>% 
+  dplyr::select(sitename, Species) %>% 
+  distinct() %>% 
+  group_by(sitename) %>% 
+  summarise(n = n())
+
+# test : number of sites per species
+df6 %>% 
+  dplyr::select(sitename, Species) %>% 
+  distinct() %>% 
+  group_by(Species) %>% 
+  summarise(n = n()) %>% 
+  arrange(n)
+
+saveRDS(df6, file = "data/df_tga.rds")
+```
+
+## TGA
+
+Plot.
+```{r eval=FALSE}
+df6 %>% 
+  ggplot(aes(x = Dr_res_sitemean, y = Dr_res)) +  # , color = Species
+  geom_point() +
+  geom_abline(slope = 1, intercept = 0, linetype = "dotted")
+```
+
+Plot just the lines
+```{r}
+gg1 <- df6 %>% 
+  group_by(Family) %>% 
+  ggplot(aes(x = Dr_res_sitemean, y = Dr_res, group = Species, color = Family)) +
+  geom_smooth(method = "lm", se = FALSE, size = 0.5, alpha = 0.2) +
+  geom_abline(intercept=0, slope=1, linetype="dotted") +
+  theme_classic() +
+  geom_point(alpha = 0.3) +
+  labs(x = expression("Site mean" ~ italic("z")[r] ~ "(m)"),
+       y = expression(italic("z")[r] ~ "(m)")) +
+  theme(legend.title = element_blank())
+
+gg1
+ggsave("fig/tga_rsip.pdf", width = 6, height = 4)
+ggsave("fig/tga_rsip.png", width = 6, height = 4)
+```
+
+
+<!-- Fit linear regressions by species for Dr (not residual). -->
+<!-- ```{r} -->
+<!-- get_width <- function(df){ -->
+<!--   df %>% pull(Dr) %>% range() %>% diff() -->
+<!-- } -->
+
+<!-- df_tga <- df5 %>%  -->
+<!--   dplyr::filter(Species %in% use_species & sitename %in% use_sites) %>%  -->
+<!--   group_by(Species) %>%  -->
+<!--   nest() %>%  -->
+<!--   mutate(linmod = purrr::map(data, ~lm(Dr ~ Dr_sitemean, data = .)), -->
+<!--          width = purrr::map_dbl(data, ~get_width(.))) %>%  -->
+<!--   mutate(slope = purrr::map_dbl(linmod, ~coef(.)[2])) %>%  -->
+<!--   left_join(df5 %>%  -->
+<!--               dplyr::select(Species, Family) %>%  -->
+<!--               distinct(),  -->
+<!--             by = "Species") -->
+
+<!-- df_tga %>%  -->
+<!--   ggplot() + -->
+<!--   geom_histogram(aes(slope, ..density..), fill = "royalblue", binwidth = 0.4, alpha = 0.5) + -->
+<!--   geom_density(aes(slope, ..density..), color = "royalblue") + -->
+<!--   xlim(-2,3) -->
+
+<!-- df_tga %>%  -->
+<!--   ggplot(aes(width, slope)) + -->
+<!--   geom_point() + -->
+<!--   geom_smooth(method = "lm") + -->
+<!--   geom_hline(yintercept = 1, linetype = "dotted") -->
+
+<!-- linmod2 <- lm(slope ~ width, data = df_tga) -->
+<!-- summary(linmod2) -->
+
+<!-- df_tga %>%  -->
+<!--   group_by(Family) %>%  -->
+<!--   summarise(slope = mean(slope)) %>%  -->
+<!--   mutate(Family = forcats::fct_reorder(Family, slope)) %>%  -->
+<!--   drop_na() %>%  -->
+<!--   ggplot(aes(Family, slope)) + -->
+<!--   geom_bar(stat = "identity") + -->
+<!--   coord_flip() -->
+<!-- ``` -->
+
+Fit linear regressions by species for `Dr_res` (residual of Dr from model log(Dr) ~ log(Hs)).
+```{r}
+get_width_res <- function(df){
+  df %>% pull(Dr_res) %>% range() %>% diff()
+}
+
+df_tga_res <- df5 %>% 
+  dplyr::filter(Species %in% use_species & sitename %in% use_sites) %>% 
+  group_by(Species) %>% 
+  nest() %>% 
+  mutate(linmod = purrr::map(data, ~lm(Dr_res ~ Dr_res_sitemean, data = .)),
+         width = purrr::map_dbl(data, ~get_width_res(.))) %>% 
+  mutate(slope = purrr::map_dbl(linmod, ~coef(.)[2])) %>% 
+  left_join(df5 %>% 
+              dplyr::select(Species, Family) %>% 
+              distinct(), 
+            by = "Species") %>% 
+  
+  ## remove outlier slope
+  filter(slope < 30) %>% 
+  mutate(n = purrr::map_int(data, ~nrow(.)))
+
+gg2 <- df_tga_res %>% 
+  ggplot() +
+  geom_histogram(aes(slope, ..count..), binwidth = 0.4, alpha = 0.7) +
+  # geom_density(aes(slope, ..density..), color = "red") +
+  geom_vline(xintercept = 1.0, linetype = "dotted") +
+  theme_classic() +
+  labs(x = "Slope", y = "Count")
+
+
+df_tga_res %>% 
+  ggplot(aes(width, slope)) +
+  geom_point() +
+  geom_smooth(method = "lm") +
+  geom_hline(yintercept = 1, linetype = "dotted")
+
+df_tga_res %>% 
+  ggplot(aes(width, abs(slope))) +
+  geom_point() +
+  geom_smooth(method = "lm") +
+  geom_hline(yintercept = 1, linetype = "dotted")
+
+linmod2 <- lm(slope ~ width, data = df_tga_res)
+summary(linmod2)
+
+## species-level
+gg3 <- df_tga_res %>% 
+  drop_na() %>% 
+  ggplot(aes(forcats::fct_reorder(Species, slope), slope)) +
+  geom_bar(stat = "identity") +
+  geom_hline(yintercept = 1.0, linetype = "dotted") +
+  coord_flip() +
+  labs(y = "Slope", x = "")
+
+## family level
+df_tga_res %>% 
+  group_by(Family) %>% 
+  summarise(slope = mean(slope), n = n()) %>% 
+  mutate(family_n = paste0(Family, " (", as.character(n), ")")) %>% 
+  mutate(family_n = forcats::fct_reorder(family_n, slope)) %>% 
+  drop_na() %>% 
+  ggplot(aes(family_n, slope)) +
+  geom_bar(stat = "identity") +
+  geom_hline(yintercept = 1.0, linetype = "dotted") +
+  coord_flip()
+```
+
+## Publication figures
+
+```{r}
+bottom_row <- plot_grid(gg2, gg3, ncol = 4, labels = c('b', 'c'), rel_widths = c(0.3, 0.7))
+plot_grid(gg1, bottom_row, ncol = 1, labels = c('a', ''))
+
+top_row <- plot_grid(gg1, gg2, ncol = 2, labels = c('a', 'b'), rel_widths = c(0.7, 0.3))
+plot_grid(top_row, gg3, nrow = 2, labels = c('', 'c'), rel_heights = c(1, 1.5) )
+
+ggsave("fig/tga.pdf", width = 12, height = 9)
+ggsave("fig/tga.png", width = 12, height = 9)
+```
diff --git a/rsofun_rsip_mct.Rmd b/rsofun_rsip_mct.Rmd
index c272961..efdb483 100644
--- a/rsofun_rsip_mct.Rmd
+++ b/rsofun_rsip_mct.Rmd
@@ -15,6 +15,7 @@ library(rsofun)
 library(ingestr)
 library(rbeni)
 library(cowplot)
+library(ggridges)
 ```
 
 ## Get sites
@@ -109,11 +110,11 @@ write_csv(df_sites, file = "data/df_sites_rsip.csv")
 
 ## Get forcing data
 
-Use `submit_forcing_rsip.sh` to run `rscript_forcing_rsip.R`. This saves files in subdirectory `"./data/forcing_rsip/"`.
+Use `submit_forcing_rsip.sh` to run `rscript_forcing_rsip.R`. Thid saves files in subdirectory `"./data/forcing_rsip/"`.
 
 ## Test
 
-Use `submit_rsofun_cwdx.sh` to run `rscript_rsofun_cwdx.R` for different S_0 (WHC) and diagnose it following the S_dEF method. This writes files into the directory `"./data/out_rsofun_cwdx/"`. 
+Use `submit_rsofun_cwdx.sh` to run `rscript_rsofun_cwdx.R` for different S_0 (WHC) and diagnose it following the S_dEF method. This runs `run_rsofun_cwdx_by_chunk.R` and saves files in subdirectory `"./data/out_rsofun_cwdx/"`. 
 
 ### Get data
 
@@ -212,9 +213,8 @@ ggsave("fig/plot_rsofun_rsip.png", width = 10, height = 8)
 
 
 ### ridge plot
-```{r}
-library(ggridges)
 
+```{r}
 gg4 <- df %>% 
   dplyr::filter(!(biome_name %in% biome_exclude)) %>%
   filter(setup == "whc_NA") %>%