rsip.Rmd

---
title: "RSIP data analysis"
author: "Beni Stocker"
date: "9/19/2019"
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)))
```

Explorations with tree height
```{r}
df2 <- df %>% 
  mutate(Hs = as.numeric(Hs)) %>% 
  mutate(log_dr = log10(Dr), log_hs = log10(Hs)) %>% 
  dplyr::filter(!is.nan(log_dr) & !is.nan(log_hs) & !is.na(log_dr) & !is.na(log_hs) & !is.infinite(log_dr) & !is.infinite(log_hs))

linmod <- lm(log_dr ~ log_hs, data = df2)

plot(linmod)

df2 %>% 
  ggplot(aes(log_hs, log_dr)) + geom_point(alpha = 0.7) + geom_smooth(method = "lm")
```

Bin the data into grid points, summarise, and expand to grid
```{r}
## bin
dlon <- 5.0
dlat <- 5.0
lon_breaks <- seq(from = -180, to = 180, by = dlon)
lat_breaks <- seq(from = -90, to = 90, by = dlat)

df <- df %>%
  ungroup() %>% 
  mutate(ilon = cut(lon, 
                    breaks = lon_breaks
                    # labels = as.character(seq(length(lon_breaks)-1))
                    ),
         ilat = cut(lat, 
                    breaks = lat_breaks
                    # labels = as.character(seq(length(lat_breaks)-1))
                    )
         ) %>% 
  mutate(lon_lower = as.numeric( sub("\\((.+),.*", "\\1", ilon)),
         lon_upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", ilon) ),
         lat_lower = as.numeric( sub("\\((.+),.*", "\\1", ilat) ),
         lat_upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", ilat) )
         ) %>% 
  mutate(lon_mid = (lon_lower + lon_upper)/2,
         lat_mid = (lat_lower + lat_upper)/2)

#View(select(df, lon, lat, ilon, ilat, lon_lower, lon_upper, lat_lower, lat_upper, lon_mid, lat_mid))

df <- df %>% 
  group_by(lon_mid, lat_mid) %>% 
  drop_na(lon_mid, lat_mid)
```

Summarise
```{r}
df_agg <- df %>% 
  summarise(mean_Dr = mean(Dr, na.rm = TRUE)) %>% 
  left_join(
    summarise(df, median_Dr = median(Dr, na.rm = TRUE)),
    by = c("lon_mid", "lat_mid")
  ) %>% 
  left_join(
    summarise(df, max_Dr = max(Dr, na.rm = TRUE)),
    by = c("lon_mid", "lat_mid")
  ) %>% 
  left_join(
    summarise(df, q75 = quantile(Dr, probs = 0.75, na.rm = TRUE) ),
    by = c("lon_mid", "lat_mid")
  ) %>% 
  left_join(
    summarise(df, q90 = quantile(Dr, probs = 0.90, na.rm = TRUE) ),
    by = c("lon_mid", "lat_mid")
  ) %>% 
  left_join(
    summarise(df, q95 = quantile(Dr, probs = 0.95, na.rm = TRUE) ),
    by = c("lon_mid", "lat_mid")
  ) %>% 
  left_join(
    summarise(df, n_Dr = n()),
    by = c("lon_mid", "lat_mid")
  )
df <- df %>% ungroup()
```


## Plots

### Distribution of values by individual data points.
```{r}
df %>% 
  ggplot(aes(x = Dr)) +
  geom_histogram(binwidth = 0.5, fill = "grey70", colour = "black") +
  xlim(0, 30)
```

### Distribution of values by bins.
```{r}
df_agg %>% 
  ggplot(aes(x = max_Dr)) +
  geom_histogram(binwidth = 2, fill = "grey70", colour = "black")
```

### Global map
```{r}
gg <- plot_map3(
  df_agg %>% 
    ungroup %>% 
    dplyr::filter(n_Dr > 2) %>% 
    dplyr::rename(lon = lon_mid, lat = lat_mid) %>%
    dplyr::select(lon, lat, median_Dr),
  breaks = c(seq(0, 1000, by = 200), 1500, 2000, 3000, 5000, 7000, 10000, 15000, 20000, Inf)/1000, 
  latmin = -55, latmax = 75,
  spacing = "constant", 
  colorscale = viridis::magma,
  combine = FALSE,
  legend_title = "(m)"
  )
gg$ggmap <- gg$ggmap + labs(title = "Rooting depth", subtitle = "RSIP data binned to 5 degrees gridcells")
cowplot::plot_grid(gg$ggmap, gg$gglegend, ncol = 2, rel_widths = c(1, 0.2))
ggsave("fig/map_zroot_rsip.pdf", width = 10, height = 6)
ggsave("fig/map_zroot_rsip.png", width = 10, height = 6)
```

## Analyse by biome

This doesn't seem to work. For almost all the points given in `df` with their `lat` and `lon`, there is no biome information in the raster file. 

First, extract biome information from WWF and NCRS biome map.
```{r eval=FALSE}
## add WWF biome info
df_wwf <- ingest(
  df %>% 
    dplyr::select(sitename = id, lon, lat),
  source = "wwf",
  dir = "~/data/biomes/wwf_ecoregions/official/",
  settings = list(layer = "wwf_terr_ecos")
  )%>% 
  mutate(data = purrr::map(data, ~slice(., 1))) %>% 
  unnest(data)

save(df_wwf, file = "data/df_wwf_rsip.RData")
```

```{r}
load("data/df_wwf_rsip.RData")

df <- df %>% 
  left_join(
    df_wwf %>% 
      dplyr::select(id = sitename, biome = BIOME, biome_name = BIOME_NAME),
    by = "id"
  )


# source("R/add_biome_ncrs.R")
# df <- add_biome_ncrs(
#   df, 
#   path = "/alphadata01/bstocker/data/biomes/NCRS/global_biomes_geotiff/biomes.tif"
#   )
# 
# df_agg <- add_biome_ncrs(
#   df_agg %>% rename(lon = lon_mid, lat = lat_mid), 
#   path = "/alphadata01/bstocker/data/biomes/NCRS/global_biomes_geotiff/biomes.tif"
#   )
```

Plot data by biomes: Individual-level data
```{r}
df %>% 
  ggplot(aes(x = as.factor(biome_name), y = -Dr)) +
  geom_boxplot() +
  # geom_violin() +
  theme(axis.text.x = element_text(angle=90, hjust=1)) +
  ylim(-10,0)

ggsave("fig/zroot_biomes_rsip.pdf", height = 10, width = 10)
```

## Add model results

```{r}
df <- extract_nc(df %>% dplyr::select(id, lon, lat), "data/zroot_cwd80.nc") %>% 
  unnest(data) %>% 
  rename(zroot_cwdx80 = V1) %>% 
  mutate(zroot_cwdx80 = zroot_cwdx80 / 1000) %>% 
  right_join(df, by = c("lon", "lat", "id"))
```

## Plots

```{r}
## old color scheme
# col_mod <- "#D55E0050"
# col_obs <- "#0072B250"
#   
# col_mod_solid <- "#D55E00"
# col_obs_solid <- "#0072B2"

library(vhs)
cols <- vhs("maxell_gu")
col_mod <- add_alpha(cols[2], 0.6)
col_obs <- add_alpha(cols[3], 0.6)
  
col_mod_solid <- cols[2]
col_obs_solid <- cols[3]

## ridge plot
biome_exclude <- df %>% 
  group_by(biome_name) %>% 
  summarise(count = n()) %>% 
  dplyr::filter(count < 3) %>% 
  pull(biome_name)

gg <- df %>% 
  dplyr::filter(!(biome_name %in% biome_exclude)) %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  tidyr::pivot_longer(cols = c(mod, obs), names_to = "source", values_to = "zroot") %>% 
  dplyr::filter(!is.na(biome_name)) %>% 
  dplyr::filter(biome_name!="Mangroves") %>%   
  ggplot(aes(x = zroot, y = biome_name, color = source, point_color = source, fill = source)) +
  geom_density_ridges(
    jittered_points = TRUE, scale = .95, rel_min_height = .01,
    point_shape = "|", point_size = 1.5, size = 0.25,
    position = position_points_jitter(height = 0)
  ) +
  scale_y_discrete(expand = c(0, 0), name = "") +
  scale_x_continuous(expand = c(0, 0), name = "Rooting depth (m)", limits = c(-5, 30)) +
  scale_fill_manual(values = c(col_mod, col_obs), labels = c(expression(italic(z)[CWDX80]), "Observed"), name = "") +
  scale_color_manual(values = c(col_mod_solid, col_obs_solid), guide = "none") +
  scale_discrete_manual("point_color", values = c(col_mod_solid, col_obs_solid), guide = "none") +
  coord_cartesian(clip = "off") +
  guides(fill = guide_legend(
    override.aes = list(
      fill = c(col_mod_solid, col_obs_solid),
      color = NA, point_color = NA)
    )
  ) +
  labs(title = expression("Rooting depth, RSIP data")) +
  theme_ridges(center = TRUE)
gg

gg_fig3a <- gg
save(gg_fig3a, file = "fig/gg_fig3a.Rdata")

ggsave("fig/modobs_ridges_zroot_biome_wwf_rsip.pdf", width = 12, height = 7)
ggsave("fig/modobs_ridges_zroot_biome_wwf_rsip.png", width = 12, height = 7)
```

Scatterplot of modelled vs. observed.
```{r}
df %>% 
  analyse_modobs2("zroot_cwdx80", "Dr", type = "heat")
```

Distribution of values.
```{r}
df %>% 
  dplyr::select(observed = Dr, modelled = zroot_cwdx80) %>% 
  pivot_longer(cols = c("observed", "modelled"), names_to = "source", values_to = "zr") %>% 
  ggplot() +
  geom_histogram(aes(x = zr, y = ..density.., fill = source), 
                 binwidth = 0.5, alpha = 0.5,
                 position="identity") +
#   geom_histogram(aes(x = Dr, y = ..density..), binwidth = 0.5, fill = "grey70", colour = "black") +
#   geom_histogram(aes(x = zroot_cwdx80, y = ..density..), binwidth = 0.5, fill = "red", colour = "black", alpha = 0.5) +
  xlim(0, 30) +
  scale_fill_manual(values=c("red", "black")) +
  theme_classic()
```

Quantiles
```{r}
df_biome <- df %>% 
  group_by(biome_name) %>% 
  summarise(q10_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.1, na.rm = TRUE),
            q25_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.25, na.rm = TRUE),
            q50_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.5, na.rm = TRUE),
            q75_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.75, na.rm = TRUE),
            q90_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.9, na.rm = TRUE),
            
            q10_zroot_obs = quantile(Dr, probs = 0.1, na.rm = TRUE),
            q25_zroot_obs = quantile(Dr, probs = 0.25, na.rm = TRUE),
            q50_zroot_obs = quantile(Dr, probs = 0.5, na.rm = TRUE),
            q75_zroot_obs = quantile(Dr, probs = 0.75, na.rm = TRUE),
            q90_zroot_obs = quantile(Dr, probs = 0.9, na.rm = TRUE)
            )

out <- df_biome %>% 
  analyse_modobs2("q10_zroot_cwdx80", "q10_zroot_obs", shortsubtitle = TRUE)
gg_rsip_10 <- out$gg + labs(title = "10% quantile by biome, RSIP data",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_10
save(gg_rsip_10, file = "data/gg_rsip_10.RData")

out <- df_biome %>% 
  analyse_modobs2("q25_zroot_cwdx80", "q25_zroot_obs", shortsubtitle = TRUE)
gg_rsip_25 <- out$gg + labs(title = "25% quantile by biome, RSIP data",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_25
save(gg_rsip_25, file = "data/gg_rsip_25.RData")

out <- df_biome %>% 
  analyse_modobs2("q50_zroot_cwdx80", "q50_zroot_obs", shortsubtitle = TRUE)
gg_rsip_50 <- out$gg + labs(title = "50% quantile by biome, RSIP data",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_50
save(gg_rsip_50, file = "data/gg_rsip_50.RData")
ggsave("fig/modobs_zroot_q90_rsip.pdf", width = 6, height = 5)

out <- df_biome %>% 
  analyse_modobs2("q75_zroot_cwdx80", "q75_zroot_obs", shortsubtitle = TRUE)
gg_rsip_75 <- out$gg + labs(title = "75% quantile by biome, RSIP data",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_75
save(gg_rsip_75, file = "data/gg_rsip_75.RData")

out <- df_biome %>% 
  analyse_modobs2("q90_zroot_cwdx80", "q90_zroot_obs", shortsubtitle = TRUE)
gg_rsip_90 <- out$gg + labs(title = "90% quantile by biome, RSIP data",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_90
save(gg_rsip_90, file = "data/gg_rsip_90.RData")
```

### Analysis Fan subset

Analyse data and model-data agreement for a subset of the data where water table depth observations are available, and set modelled rooting depth to WTD if it exceeds WTD.
```{r}
df_wtd <- df %>% 
  dplyr::filter(!is.na(wtd)) %>% 

  mutate(zroot_cwdx80 = ifelse(zroot_cwdx80 > wtd, wtd, zroot_cwdx80)) %>% 
  
  ## and observed rooting depth > wtd removed
  mutate(Dr = ifelse(Dr > wtd, NA, Dr)) %>%
  
  dplyr::filter(!is.na(Dr) & !is.na(zroot_cwdx80))

print(paste("Total number of entries:", nrow(df_wtd)))
```

Plot as in Fan et al., 2017, Fig. 3F.
```{r}
df %>% 
  dplyr::filter(!is.na(wtd)) %>%   
  mutate(log_wtd = log10(wtd), log_Dr = log10(Dr)) %>% 
  ggplot(aes(log_wtd, -log_Dr)) +
  geom_point() + 
  geom_abline(intercept=0, slope=-1, linetype="dotted")
```

Plots with subset of data, where modelled values > WTD are removed.
```{r}
out <- df_wtd %>% 
  analyse_modobs2("zroot_cwdx80", "Dr")
out$gg +
  xlim(0,10) + ylim(0,10)
```

```{r}
## boxplot
df_wtd %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  pivot_longer(cols = c(obs, mod), names_to = "source", values_to = "zroot") %>% 
  ggplot(aes(x = biome_name, y = -zroot, fill = source)) +
  geom_boxplot() +
  # geom_violin() +
  # geom_bar(stat = "summary", fun.ymax = "max", position=position_dodge()) +
  theme(axis.text.x = element_text(angle=90, hjust=1)) +
  ylim(-20,0)

ggsave("fig/zroot_biomes_rsip_boxplot_modobs_wtd.pdf", height = 10, width = 10)

## barplot: quantile
df_wtd %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  pivot_longer(cols = c(obs, mod), names_to = "source", values_to = "zroot") %>% 
  ggplot(aes(x = biome_name, fill = source)) +
  stat_summary_bin(aes(y = -zroot), fun.y = "quantile", fun.args = c(probs = 0.05), geom = "bar", position=position_dodge()) +
  theme(axis.text.x = element_text(angle=90, hjust=1))

ggsave("fig/zroot_biomes_rsip_barplot_modobs_wtd.pdf", height = 10, width = 10)

## ridge plot
biome_exclude <- df_wtd %>% 
  group_by(biome_name) %>% 
  summarise(count = n()) %>% 
  dplyr::filter(count < 3) %>% 
  pull(biome_name)

gg <- df_wtd %>% 
  dplyr::filter(!(biome_name %in% biome_exclude)) %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  tidyr::pivot_longer(cols = c(mod, obs), names_to = "source", values_to = "zroot") %>% 
  dplyr::filter(!is.na(biome_name)) %>% 
  dplyr::filter(biome_name!="Mangroves") %>%   
  ggplot(aes(x = zroot, y = biome_name, color = source, point_color = source, fill = source)) +
  geom_density_ridges(
    jittered_points = TRUE, scale = .95, rel_min_height = .01,
    point_shape = "|", point_size = 1.5, size = 0.25,
    position = position_points_jitter(height = 0)
  ) +
  scale_y_discrete(expand = c(0, 0), name = "") +
  scale_x_continuous(expand = c(0, 0), name = "Rooting depth (m)", limits = c(-5, 30)) +
  scale_fill_manual(values = c(col_mod, col_obs), labels = c(expression(italic(z)[CWDX80]), "Observed"), name = "") +
  scale_color_manual(values = c(col_mod_solid, col_obs_solid), guide = "none") +
  scale_discrete_manual("point_color", values = c(col_mod_solid, col_obs_solid), guide = "none") +
  coord_cartesian(clip = "off") +
  guides(fill = guide_legend(
    override.aes = list(
      fill = c("#D55E00A0", "#0072B2A0"),
      color = NA, point_color = NA)
    )
  ) +
  labs(title = expression("Rooting depth, RSIP data, WTD subset")) +
  theme_ridges(center = TRUE)

gg

ggsave("fig/modobs_ridges_zroot_biome_wwf_rsip_wtd.pdf", width = 12, height = 7)
ggsave("fig/modobs_ridges_zroot_biome_wwf_rsip_wtd.png", width = 12, height = 7)
```

Quantiles
```{r}
df_biome_wtd <- df_wtd %>% 
  group_by(biome_name) %>% 
  summarise(q10_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.1, na.rm = TRUE),
            q25_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.25, na.rm = TRUE),
            q50_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.5, na.rm = TRUE),
            q75_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.75, na.rm = TRUE),
            q90_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.9, na.rm = TRUE),
            
            q10_zroot_obs = quantile(Dr, probs = 0.1, na.rm = TRUE),
            q25_zroot_obs = quantile(Dr, probs = 0.25, na.rm = TRUE),
            q50_zroot_obs = quantile(Dr, probs = 0.5, na.rm = TRUE),
            q75_zroot_obs = quantile(Dr, probs = 0.75, na.rm = TRUE),
            q90_zroot_obs = quantile(Dr, probs = 0.9, na.rm = TRUE)
            )

out <- df_biome_wtd %>% 
  analyse_modobs2("q10_zroot_cwdx80", "q10_zroot_obs", shortsubtitle = TRUE)
gg_rsip_wtd_10 <- out$gg + labs(title = "10% quantile by biome, RSIP data, WTD subset",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_wtd_10
save(gg_rsip_wtd_10, file = "data/gg_rsip_wtd_10.RData")

out <- df_biome_wtd %>% 
  analyse_modobs2("q25_zroot_cwdx80", "q25_zroot_obs", shortsubtitle = TRUE)
gg_rsip_wtd_25 <- out$gg + labs(title = "25% quantile by biome, RSIP data, WTD subset",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_wtd_25
save(gg_rsip_wtd_25, file = "data/gg_rsip_wtd_25.RData")

out <- df_biome_wtd %>% 
  analyse_modobs2("q50_zroot_cwdx80", "q50_zroot_obs", shortsubtitle = TRUE)
gg_rsip_wtd_50 <- out$gg + labs(title = "50% quantile by biome, RSIP data, WTD subset",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_wtd_50
save(gg_rsip_wtd_50, file = "data/gg_rsip_wtd_50.RData")
ggsave("fig/modobs_zroot_q90_rsip_wtd.pdf", width = 6, height = 5)

out <- df_biome_wtd %>% 
  analyse_modobs2("q75_zroot_cwdx80", "q75_zroot_obs", shortsubtitle = TRUE)
gg_rsip_wtd_75 <- out$gg + labs(title = "75% quantile by biome, RSIP data, WTD subset",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_wtd_75
save(gg_rsip_wtd_75, file = "data/gg_rsip_wtd_75.RData")

out <- df_biome_wtd %>% 
  analyse_modobs2("q90_zroot_cwdx80", "q90_zroot_obs", shortsubtitle = TRUE)
gg_rsip_wtd_90 <- out$gg + labs(title = "90% quantile by biome, RSIP data, WTD subset",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_wtd_90
save(gg_rsip_wtd_90, file = "data/gg_rsip_wtd_90.RData")
```

Focus on data points with large overestimation of rooting depth.
```{r}
tmp <- df %>% 
  mutate(bias = zroot_cwdx80/Dr)

tmp %>% 
  ggplot(aes(bias, ..count..)) +
  geom_histogram()

tmp <- tmp %>% 
  dplyr::filter(!is.na(wtd)) %>%
  dplyr::filter(bias > 1)
```

## Aggregated by sites

1705 distinct sites based on lon and lat info.
```{r}
# most_frequent <- function(x){x %>% table() %>% sort() %>% rev() %>% head(1) %>% names()}

df_sites <- df %>% 
  mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))) %>% 
  group_by(sitename) %>% 
  summarise(Dr = mean(Dr, na.rm = TRUE), zroot_cwdx80 = mean(zroot_cwdx80, na.rm = TRUE), wtd = mean(wtd, na.rm = TRUE))

## add WWF biome info
df_wwf_sites <- ingest(
  df %>% 
    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_sites <- df_sites %>% 
  left_join(
    df_wwf_sites %>% 
      dplyr::select(sitename, biome = BIOME, biome_name = BIOME_NAME),
    by = "sitename"
  )

df_sites %>%
  separate(sitename, into = c(NA, "lon", "lat"), sep = "_") %>% 
  mutate(lon = as.numeric(lon), lat = as.numeric(lat)) %>% 
  write_csv(file = "data/df_sites_rsip.csv")
```

Plot sites on map
```{r}
library(jcolors)
nbiomes <- df_sites %>% 
  drop_na(biome_name) %>%
  pull(biome_name) %>% 
  unique() %>% 
  length()

plot_map_simpl() +
  geom_point(data = df_sites %>% 
               drop_na(biome_name) %>% 
               separate(sitename, into = c(NA, "lon", "lat"), sep = "_") %>% 
               mutate(lon = as.numeric(lon), lat = as.numeric(lat)), 
             aes(x = lon, y = lat, color = biome_name)) +
  # scale_color_jcolors(palette = "pal8")
  # scale_color_brewer(palette = "Paired")
  scale_colour_manual(values = myjcolors(nbiomes), name = "Biome name")


ggsave("fig/sites_rsip.pdf", width = 15, height = 6)
ggsave("fig/sites_rsip.png", width = 15, height = 6)
```

Plots with data aggregated by sites.
```{r}
## boxplot
df_sites %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  pivot_longer(cols = c(obs, mod), names_to = "source", values_to = "zroot") %>% 
  ggplot(aes(x = biome_name, y = -zroot, fill = source)) +
  geom_boxplot() +
  # geom_violin() +
  # geom_bar(stat = "summary", fun.ymax = "max", position=position_dodge()) +
  theme(axis.text.x = element_text(angle=90, hjust=1)) +
  ylim(-20,0)

ggsave("fig/zroot_biomes_rsip_boxplot_modobs_sites.pdf", height = 10, width = 10)

## barplot: quantile - THIS LOOKS ACTUALLY VERY NICE FOR THE 5% QUANTILE!!!
df_sites %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  pivot_longer(cols = c(obs, mod), names_to = "source", values_to = "zroot") %>% 
  ggplot(aes(x = biome_name, fill = source)) +
  stat_summary_bin(aes(y = -zroot), fun.y = "quantile", fun.args = c(probs = 0.05), geom = "bar", position=position_dodge()) +
  theme(axis.text.x = element_text(angle=90, hjust=1))

ggsave("fig/zroot_biomes_rsip_barplot_modobs_sites.pdf", height = 10, width = 10)

## ridge plot
biome_exclude <- df_sites %>% 
  group_by(biome_name) %>% 
  summarise(count = n()) %>% 
  dplyr::filter(count < 3) %>% 
  pull(biome_name)

gg <- df_sites %>% 
  dplyr::filter(!(biome_name %in% biome_exclude)) %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  tidyr::pivot_longer(cols = c(mod, obs), names_to = "source", values_to = "zroot") %>% 
  dplyr::filter(!is.na(biome_name)) %>% 
  dplyr::filter(biome_name!="Mangroves") %>%   
  ggplot(aes(x = zroot, y = biome_name, color = source, point_color = source, fill = source)) +
  geom_density_ridges(
    jittered_points = TRUE, scale = .95, rel_min_height = .01,
    point_shape = "|", point_size = 1.5, size = 0.25,
    position = position_points_jitter(height = 0)
  ) +
  scale_y_discrete(expand = c(0, 0), name = "") +
  scale_x_continuous(expand = c(0, 0), name = "Rooting depth (m)", limits = c(-5, 30)) +
  scale_fill_manual(values = c(col_mod, col_obs), labels = c(expression(italic(z)[CWDX80]), "Observed"), name = "") +
  scale_color_manual(values = c(col_mod_solid, col_obs_solid), guide = "none") +
  scale_discrete_manual("point_color", values = c(col_mod_solid, col_obs_solid), guide = "none") +
  coord_cartesian(clip = "off") +
  guides(fill = guide_legend(
    override.aes = list(
      fill = c(col_mod_solid, col_obs_solid),
      color = NA, point_color = NA)
    )
  ) +
  theme_ridges(center = TRUE)

gg_fig3a <- gg
save(gg_fig3a, file = "data/gg_fig3a.Rdata")
gg

ggsave("fig/modobs_ridges_zroot_biome_wwf_rsip_sites.pdf", width = 12, height = 7)
ggsave("fig/modobs_ridges_zroot_biome_wwf_rsip_sites.png", width = 12, height = 7)
```
Quantiles
```{r}
df_biome_sites <- df_sites %>% 
  group_by(biome_name) %>% 
  summarise(q10_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.1, na.rm = TRUE),
            q25_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.25, na.rm = TRUE),
            q50_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.5, na.rm = TRUE),
            q75_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.75, na.rm = TRUE),
            q90_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.9, na.rm = TRUE),
            
            q10_zroot_obs = quantile(Dr, probs = 0.1, na.rm = TRUE),
            q25_zroot_obs = quantile(Dr, probs = 0.25, na.rm = TRUE),
            q50_zroot_obs = quantile(Dr, probs = 0.5, na.rm = TRUE),
            q75_zroot_obs = quantile(Dr, probs = 0.75, na.rm = TRUE),
            q90_zroot_obs = quantile(Dr, probs = 0.9, na.rm = TRUE)
            )

out <- df_biome_sites %>% 
  analyse_modobs2("q10_zroot_cwdx80", "q10_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_10 <- out$gg + labs(title = "10% quantile",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_10
save(gg_rsip_sites_10, file = "data/gg_rsip_sites_10.RData")

out <- df_biome_sites %>% 
  analyse_modobs2("q25_zroot_cwdx80", "q25_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_25 <- out$gg + labs(title = "25% quantile",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_25
save(gg_rsip_sites_25, file = "data/gg_rsip_sites_25.RData")

out <- df_biome_sites %>% 
  analyse_modobs2("q50_zroot_cwdx80", "q50_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_50 <- out$gg + labs(title = "50% quantile",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_50
save(gg_rsip_sites_50, file = "data/gg_rsip_sites_50.RData")
ggsave("fig/modobs_zroot_q90_rsip_sites.pdf", width = 6, height = 5)

out <- df_biome_sites %>% 
  analyse_modobs2("q75_zroot_cwdx80", "q75_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_75 <- out$gg + labs(title = "75% quantile",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_75
save(gg_rsip_sites_75, file = "data/gg_rsip_sites_75.RData")

out <- df_biome_sites %>% 
  analyse_modobs2("q90_zroot_cwdx80", "q90_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_90 <- out$gg + labs(title = "90% quantile",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_90
save(gg_rsip_sites_90, file = "data/gg_rsip_sites_90.RData")

df_biome_sites %>% 
  dplyr::select(biome_name, q10_zroot_obs, q10_zroot_cwdx80) %>% 
  arrange(q10_zroot_obs)

df_biome_sites %>% 
  dplyr::select(biome_name, q90_zroot_obs, q90_zroot_cwdx80) %>% 
  arrange(q90_zroot_obs)
```



### Analysis Fan subset at site-level

Analyse data and model-data agreement for a subset of the data where water table depth observations are available, and set modelled rooting depth to WTD if it exceeds WTD.
```{r}
df_wtd <- df %>% 
  dplyr::filter(!is.na(wtd)) %>% 
  mutate(wtd_limited  = ifelse(zroot_cwdx80 > wtd, TRUE, FALSE)) %>% 
  mutate(zroot_cwdx80 = ifelse(zroot_cwdx80 > wtd, wtd, zroot_cwdx80))

print(paste("Fraction of entries WTD-limited:", sum(df_wtd$wtd_limited, na.rm = TRUE)/nrow(df_wtd)))

df_wtd <- df_wtd %>% 
  dplyr::filter(!is.na(Dr) & !is.na(zroot_cwdx80))

df_sites_wtd <- df_wtd %>% 
  mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))) %>% 
  group_by(sitename) %>% 
  summarise(Dr = mean(Dr, na.rm = TRUE), zroot_cwdx80 = mean(zroot_cwdx80, na.rm = TRUE), wtd = mean(wtd, na.rm = TRUE)) %>% 
  left_join(
    df_wwf_sites %>% 
      dplyr::select(sitename, biome = BIOME, biome_name = BIOME_NAME),
    by = "sitename"
  )

print(paste("Total number of entries:", nrow(df_wtd)))
print(paste("Total number of sites:", nrow(df_sites_wtd)))
```

Plots
```{r}
## boxplot
df_sites_wtd %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  pivot_longer(cols = c(obs, mod), names_to = "source", values_to = "zroot") %>% 
  ggplot(aes(x = biome_name, y = -zroot, fill = source)) +
  geom_boxplot() +
  # geom_violin() +
  # geom_bar(stat = "summary", fun.ymax = "max", position=position_dodge()) +
  theme(axis.text.x = element_text(angle=90, hjust=1)) +
  ylim(-20,0)

ggsave("fig/zroot_biomes_rsip_boxplot_modobs_sites_wtd.pdf", height = 10, width = 10)

## barplot: quantile - THIS LOOKS ACTUALLY VERY NICE FOR THE 5% QUANTILE!!!
df_sites_wtd %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  pivot_longer(cols = c(obs, mod), names_to = "source", values_to = "zroot") %>% 
  ggplot(aes(x = biome_name, fill = source)) +
  stat_summary_bin(aes(y = -zroot), fun.y = "quantile", fun.args = c(probs = 0.05), geom = "bar", position=position_dodge()) +
  theme(axis.text.x = element_text(angle=90, hjust=1))

ggsave("fig/zroot_biomes_rsip_barplot_modobs_sites_wtd.pdf", height = 10, width = 10)

## ridge plot
biome_exclude <- df_sites_wtd %>% 
  group_by(biome_name) %>% 
  summarise(count = n()) %>% 
  dplyr::filter(count < 3) %>% 
  pull(biome_name)

df_sites_wtd %>% 
  dplyr::filter(!(biome_name %in% biome_exclude)) %>% 
  rename(obs = Dr, mod = zroot_cwdx80) %>% 
  tidyr::pivot_longer(cols = c(mod, obs), names_to = "source", values_to = "zroot") %>% 
  dplyr::filter(!is.na(biome_name)) %>% 
  dplyr::filter(biome_name!="Mangroves") %>%   
  ggplot(aes(x = zroot, y = biome_name, color = source, point_color = source, fill = source)) +
  geom_density_ridges(
    jittered_points = TRUE, scale = .95, rel_min_height = .01,
    point_shape = "|", point_size = 1.5, size = 0.25,
    position = position_points_jitter(height = 0)
  ) +
  scale_y_discrete(expand = c(0, 0), name = "") +
  scale_x_continuous(expand = c(0, 0), name = "Rooting depth (m)", limits = c(-5, 30)) +
  scale_fill_manual(values = c(col_mod, col_obs), labels = c(expression(italic(z)[CWDX80]), "Observed"), name = "") +
  scale_color_manual(values = c(col_mod_solid, col_obs_solid), guide = "none") +
  scale_discrete_manual("point_color", values = c(col_mod_solid, col_obs_solid), guide = "none") +
  coord_cartesian(clip = "off") +
  guides(fill = guide_legend(
    override.aes = list(
      fill = c("#D55E00A0", "#0072B2A0"),
      color = NA, point_color = NA)
    )
  ) +
  labs(title = expression("Rooting depth, RSIP data - WTD limit")) +
  theme_ridges(center = TRUE)

ggsave("fig/modobs_ridges_zroot_biome_wwf_rsip_sites_wtd.pdf", width = 12, height = 7)
ggsave("fig/modobs_ridges_zroot_biome_wwf_rsip_sites_wtd.png", width = 12, height = 7)
```

Quantiles
```{r}
df_biome_sites_wtd <- df_sites_wtd %>% 
  group_by(biome_name) %>% 
  summarise(q10_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.1, na.rm = TRUE),
            q25_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.25, na.rm = TRUE),
            q50_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.5, na.rm = TRUE),
            q75_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.75, na.rm = TRUE),
            q90_zroot_cwdx80 = quantile(zroot_cwdx80, probs = 0.9, na.rm = TRUE),
            
            q10_zroot_obs = quantile(Dr, probs = 0.1, na.rm = TRUE),
            q25_zroot_obs = quantile(Dr, probs = 0.25, na.rm = TRUE),
            q50_zroot_obs = quantile(Dr, probs = 0.5, na.rm = TRUE),
            q75_zroot_obs = quantile(Dr, probs = 0.75, na.rm = TRUE),
            q90_zroot_obs = quantile(Dr, probs = 0.9, na.rm = TRUE)
            )

out <- df_biome_sites_wtd %>% 
  analyse_modobs2("q10_zroot_cwdx80", "q10_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_wtd_10 <- out$gg + labs(title = "10% quantile, WTD limited",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_wtd_10
save(gg_rsip_sites_wtd_10, file = "data/gg_rsip_sites_wtd_10.RData")

out <- df_biome_sites_wtd %>% 
  analyse_modobs2("q25_zroot_cwdx80", "q25_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_wtd_25 <- out$gg + labs(title = "25% quantile, WTD limited",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_wtd_25
save(gg_rsip_sites_wtd_25, file = "data/gg_rsip_sites_wtd_25.RData")

out <- df_biome_sites_wtd %>% 
  analyse_modobs2("q50_zroot_cwdx80", "q50_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_wtd_50 <- out$gg + labs(title = "50% quantile, WTD limited",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_wtd_50
save(gg_rsip_sites_wtd_50, file = "data/gg_rsip_sites_wtd_50.RData")
ggsave("fig/modobs_zroot_q90_rsip_sites_wtd.pdf", width = 6, height = 5)

out <- df_biome_sites_wtd %>% 
  analyse_modobs2("q75_zroot_cwdx80", "q75_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_wtd_75 <- out$gg + labs(title = "75% quantile, WTD limited",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_wtd_75
save(gg_rsip_sites_wtd_75, file = "data/gg_rsip_sites_wtd_75.RData")

out <- df_biome_sites_wtd %>% 
  analyse_modobs2("q90_zroot_cwdx80", "q90_zroot_obs", 
                  shortsubtitle = TRUE, rsquared = FALSE)
gg_rsip_sites_wtd_90 <- out$gg + labs(title = "90% quantile, WTD limited",
              x = expression(italic(z)[CWDX80] ~ " (m)"),
              y = expression("Observed" ~ italic(z)[r] ~ " (m)"))
gg_rsip_sites_wtd_90
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_sites2 %>%  -->
<!--   dplyr::select(sitename, Dr_sitemean = Dr, Dr_res_sitemean = Dr_res) %>%  -->
<!--   right_join( -->
<!--     df2 %>%  -->
<!--       mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))),  -->
<!--     by = "sitename" -->
<!--     ) -->
<!-- ``` -->

<!-- Use data only for species that appear in at least 3 sites => 72 species. -->
<!-- ```{r} -->
<!-- use_species <- df2 %>% -->
<!--   filter(Species != "NA") %>%  -->
<!--   mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))) %>%  -->
<!--   dplyr::select(sitename, Species) %>%  -->
<!--   distinct() %>%  -->
<!--   group_by(Species) %>%  -->
<!--   summarise(n = n()) %>%  -->
<!--   dplyr::filter(n >= 3) %>%  -->
<!--   pull(Species) -->
<!-- ``` -->

<!-- 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 >= 3) %>%  -->
<!--   pull(sitename) -->
<!-- ``` -->

<!-- Apply these two criteria. This leaves 286 data points. -->
<!-- ```{r} -->
<!-- df3 <- df2 %>%  -->
<!--   dplyr::filter(Species %in% use_species & sitename %in% use_sites) -->
<!--   # dplyr::filter(Species %in% use_species2) -->

<!-- saveRDS(df3, file = "data/df_tga.rds") -->
<!-- ``` -->


<!-- Plot. -->
<!-- ```{r eval=FALSE} -->
<!-- df3 %>%  -->
<!--   ggplot(aes(x = Dr_sitemean, y = Dr)) +  # , color = Species -->
<!--   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 -->
<!-- ```{r} -->
<!-- df3 %>%  -->
<!--   group_by(Species) %>%  -->
<!--   ggplot(aes(x = Dr_sitemean, y = Dr, 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) -->

<!-- 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 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 = 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(),  -->
<!--             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 <- 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 = 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(gg_rsip_sites_10, gg_rsip_sites_90, gg_rsip_sites_wtd_10, gg_rsip_sites_wtd_90, ncol = 4, labels = c('b', 'c', 'd', 'e'))
plot_grid(gg_fig3a, bottom_row, ncol = 1, labels = c('a', ''), rel_heights = c(0.7, 0.4))
ggsave("fig/fig4.pdf", width = 12, height = 9)
ggsave("fig/fig4.png", width = 12, height = 9)
```