---
title: "Validation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Validation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r, echo=FALSE, message=FALSE}
# Plotting Packages
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)
library(ggpmisc)
library(yardstick)
library(gt)
```

How do we know if Real Twig returns accurate volume estimates? We rigorously tested our method against multiple high quality reference data sets that were both laser scanner and destructively sampled. The laser scanning was done in leaf-off conditions with a Riegl VZ-400. We used different versions of TreeQSM with the same input parameters per tree and different data sets to test our method. A detailed discussion of the results and implications can be found in *Morales & MacFarlane (2024)*: <https://doi.org/10.1093/forestry/cpae046>

The graph and table below contain destructive sampling with total mass and basic density (wood + bark) across four data sets, three of which are publicly available and are linked below. The mass estimates and statistics use TreeQSM v2.4.1 with its built in tapering compared to Real Twig on the same QSMs.

-   Burt et al. (2021): <https://zenodo.org/records/4056899>, <https://zenodo.org/records/4056903>

-   Demol et al. (2021): <https://zenodo.org/records/4557401>

-   Hackenberg et al. (2015b): <https://www.simpleforest.org/pages/dataQuercusPetraea.html>

```{r, echo=FALSE, results='hide', message=FALSE, warning=FALSE}
# Import all data
file <- system.file("extdata/validation2.csv", package = "rTwig")
data_comp <- read.csv(file)

p_all <- ggplot(data = data_comp, aes(x = dbh.cm, y = tperr, color = version, fill = version)) +
  geom_point(aes(shape = Dataset, color = version, fill = version)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs")) +
  geom_hline(yintercept = 0) +
  labs(
    x = "DBH (cm)",
    y = "Total Mass Error (%)",
    color = "Method",
    fill = "Method",
    shape = "Dataset"
  ) +
  theme_classic() +
  scale_color_manual(values = c("#D41159", "#1A85FF")) + 
  scale_fill_manual(values = c("#D41159", "#1A85FF")) +
  geom_hline(
    yintercept = 10,
    linetype = "dashed",
    color = "black",
    linewidth = 0.25
  ) +
  geom_hline(
    yintercept = -10,
    linetype = "dashed",
    color = "black",
    linewidth = 0.25,
    show.legend = TRUE
  )
```

```{r, echo=FALSE, fig.width=7, fig.height=4, fig.align='center', warning=FALSE}
p_all
```

```{r, echo=FALSE, warning=FALSE}
### ALL TREES ##################################################################
total_stats <- data_comp %>%
  ungroup() %>%
  group_by(version) %>%
  summarize(
    MRE.pct = mean(tperr, na.rm = TRUE),
    RMSE.kg = sqrt(mean(tdiff^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Mt.DS, na.rm = TRUE) * 100
  )

CCC_total <- data_comp %>%
  ungroup() %>%
  group_by(version) %>%
  yardstick::ccc(Mt.TLS, Mt.DS) %>%
  select(.estimate) %>%
  rename(CCC = 1)

total_stats <- bind_cols(total_stats, CCC_total) %>%
  mutate(type = "total") %>%
  relocate(type, .before = MRE.pct) %>%
  select(-type)

data <- total_stats %>%
  mutate_if(is.numeric, round, 3) %>%
  rename(
    "Mean Relative Error (%)" = MRE.pct,
    "RMSE (kg)" = RMSE.kg,
    "Relative RMSE (%)" = RRMSE.pct
  )

data %>%
  gt() %>%
  # tab_header(title = unique(.$`_data`$version)) %>%
  # cols_hide(version) %>%
  cols_label(
    version = "Method",
  ) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  tab_options(
    table.border.top.style = "hidden",
    table.border.bottom.style = "hidden",
    table_body.hlines.color = "white",
    table.font.color = "black",
    heading.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    table_body.border.bottom.color = "black"
  )
```

## TreeQSM v2.4.1

The following figures look at the Harvard Forest data set using different versions of TreeQSM. The destructive sampling data contains total branch and main stem dry mass, and also basic density for both the main stem and the branches. This allows us to test for compensating errors, and look at the true differences between TreeQSM versions. We used the same input parameters per tree and TreeQSM version to test our model.

```{r, echo=FALSE, results='hide', message=FALSE, warning=FALSE}
# Import Data
file <- system.file("extdata/validation.csv", package = "rTwig")
all_biomass <- read.csv(file) %>%
  rename(Version = version)

# Distinct colors
factors <- distinct(all_biomass, Version)
cbPalette <- c("#009E73", "#CC79A7", "#56B4E9", "#999999", "#E69F00", "#009E73", "#F0E442", "#0072B2", "#D55E00")
cbPalette <- cbPalette[1:nrow(factors)]
factors <- bind_cols(factors, cbPalette) %>%
  rename(colors = 2)

# Filter Data
all_biomass <- all_biomass %>%
  left_join(factors, by = join_by(Version)) %>%
  filter(Version %in% c("TreeQSM v2.4.1", "Real Twig (TreeQSM v2.4.1)")) %>%
  mutate(id = paste0(SpCode, Tree))

# Total Biomass
p1 <- ggplot(data = all_biomass, aes(x = Mt.DS, y = Mt.TLS, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name), ) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_abline() +
  stat_poly_eq(aes(
    label = paste(after_stat(eq.label), after_stat(rr.label), sep = "~~~")
  )) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "TLS Biomass (kg)",
    title = "Total",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Main Stem Biomass
p2 <- ggplot(data = all_biomass, aes(x = Md.DS, y = Md.TLS, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_abline() +
  stat_poly_eq(aes(
    label = paste(after_stat(eq.label), after_stat(rr.label), sep = "~~~")
  )) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "", # "TLS Biomass (kg)",
    title = "Main Stem",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Branch Biomass
p3 <- ggplot(data = all_biomass, aes(x = Mb.DS, y = Mb.TLS, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_abline() +
  stat_poly_eq(aes(
    label = paste(after_stat(eq.label), after_stat(rr.label), sep = "~~~")
  )) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "", # "TLS Biomass (kg)",
    title = "Branch",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Total Biomass % Error
p4 <- ggplot(data = all_biomass, aes(x = Mt.DS, y = tperr, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "% TLS Biomass Error",
    title = "Total % Error",
    shape = "Speces",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Main Stem Biomass % Error
p5 <- ggplot(data = all_biomass, aes(x = Md.DS, y = sperr, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "Reference Biomass (kg)",
    y = "", # "% TLS Biomass Error",
    title = "Main Stem % Error",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Branch Biomass % Error
p6 <- ggplot(data = all_biomass, aes(x = Mb.DS, y = tperr, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "", # "% TLS Biomass Error",
    title = "Branch % Error",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))
```

```{r, echo=FALSE, fig.width=7, fig.height=4, fig.align='center', warning=FALSE}
ggarrange(p1, p4, common.legend = TRUE, legend = "bottom")
ggarrange(p2, p5, common.legend = TRUE, legend = "bottom")
ggarrange(p3, p6, common.legend = TRUE, legend = "bottom")
```

```{r, echo=FALSE, warning=FALSE}
### TOTALS #####################################################################
total_stats <- all_biomass %>%
  group_by(Version) %>%
  summarise(
    MRE.pct = mean(tperr),
    RMSE.kg = sqrt(mean(tdiff^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Mt.DS) * 100
  )

CCC_total <- all_biomass %>%
  group_by(Version) %>%
  ccc(Mt.TLS, Mt.DS) %>%
  select(Version, .estimate) %>%
  rename(CCC = 2)

total_stats <- left_join(total_stats, CCC_total, by = "Version") %>%
  mutate(type = "total")

### MAIN STEM ##################################################################
stem_stats <- all_biomass %>%
  group_by(Version) %>%
  summarise(
    MRE.pct = mean(sperr),
    RMSE.kg = sqrt(mean(sdif^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Md.DS) * 100
  )

CCC_stem <- all_biomass %>%
  group_by(Version) %>%
  ccc(Md.TLS, Md.DS) %>%
  select(Version, .estimate) %>%
  rename(CCC = 2)

stem_stats <- left_join(stem_stats, CCC_stem, by = "Version") %>%
  mutate(type = "stem")

### BRANCHES ###################################################################
branch_stats <- all_biomass %>%
  group_by(Version) %>%
  summarise(
    MRE.pct = mean(bperr),
    RMSE.kg = sqrt(mean(bdiff^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Mb.DS) * 100
  )

CCC_branch <- all_biomass %>%
  group_by(Version) %>%
  ccc(Mb.TLS, Mb.DS) %>%
  select(Version, .estimate) %>%
  rename(CCC = 2)

branch_stats <- left_join(branch_stats, CCC_branch, by = "Version") %>%
  mutate(type = "branch")

all_stats <- bind_rows(total_stats, stem_stats, branch_stats)

### TABLES ###################################################################
data <- all_stats %>%
  mutate_if(is.numeric, round, 3) %>%
  rename(
    "Mean Relative Error (%)" = MRE.pct,
    "RMSE (kg)" = RMSE.kg,
    "Relative RMSE (%)" = RRMSE.pct
  ) %>%
  pivot_longer(cols = 2:5, names_to = "metric", values_to = "value") %>%
  pivot_wider(names_from = type, values_from = value) %>%
  group_by(Version) %>%
  group_split()

tables <- vector(mode = "list", length(data))

data[[2]] %>%
  gt() %>%
  tab_header(title = "TreeQSM v2.4.1") %>%
  cols_hide(Version) %>%
  cols_label(
    metric = "Metric",
    total = "Total Woody AGB",
    stem = "Main Stem Biomass",
    branch = "Branch Biomass"
  ) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  tab_options(
    table.border.top.style = "hidden",
    table.border.bottom.style = "hidden",
    table_body.hlines.color = "white",
    table.font.color = "black",
    heading.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    table_body.border.bottom.color = "black"
  )

data[[1]] %>%
  gt() %>%
  tab_header(title = "Real Twig (TreeQSM v2.4.1)") %>%
  cols_hide(Version) %>%
  cols_label(
    metric = "Metric",
    total = "Total Woody AGB",
    stem = "Main Stem Biomass",
    branch = "Branch Biomass"
  ) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  tab_options(
    table.border.top.style = "hidden",
    table.border.bottom.style = "hidden",
    table_body.hlines.color = "white",
    table.font.color = "black",
    heading.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    table_body.border.bottom.color = "black"
  )
```

## TreeQSM v2.3.0

Below are the mass estimates and statistics using TreeQSM v2.3.0 with its built in tapering and Real Twig applied to the same QSMs.

```{r, echo=FALSE, results='hide', message=FALSE, warning=FALSE}
# Import Data
file <- system.file("extdata/validation.csv", package = "rTwig")
all_biomass <- read.csv(file) %>%
  rename(Version = version)

# Distinct colors
factors <- distinct(all_biomass, Version)
cbPalette <- c("#009E73", "#CC79A7", "#56B4E9", "#999999", "#E69F00", "#009E73", "#F0E442", "#0072B2", "#D55E00")
cbPalette <- cbPalette[1:nrow(factors)]
factors <- bind_cols(factors, cbPalette) %>%
  rename(colors = 2)

# Filter Data
all_biomass <- all_biomass %>%
  left_join(factors, by = join_by(Version)) %>%
  filter(Version %in% c("TreeQSM v2.3.0", "Real Twig (TreeQSM v2.3.0)")) %>%
  mutate(id = paste0(SpCode, Tree))

# Total Biomass
p1 <- ggplot(data = all_biomass, aes(x = Mt.DS, y = Mt.TLS, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name), ) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_abline() +
  stat_poly_eq(aes(
    label = paste(after_stat(eq.label), after_stat(rr.label), sep = "~~~")
  )) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "TLS Biomass (kg)",
    title = "Total",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Main Stem Biomass
p2 <- ggplot(data = all_biomass, aes(x = Md.DS, y = Md.TLS, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_abline() +
  stat_poly_eq(aes(
    label = paste(after_stat(eq.label), after_stat(rr.label), sep = "~~~")
  )) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "", # "TLS Biomass (kg)",
    title = "Main Stem",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Branch Biomass
p3 <- ggplot(data = all_biomass, aes(x = Mb.DS, y = Mb.TLS, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_abline() +
  stat_poly_eq(aes(
    label = paste(after_stat(eq.label), after_stat(rr.label), sep = "~~~")
  )) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "", # "TLS Biomass (kg)",
    title = "Branch",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Total Biomass % Error
p4 <- ggplot(data = all_biomass, aes(x = Mt.DS, y = tperr, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "% TLS Biomass Error",
    title = "Total % Error",
    shape = "Speces",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Main Stem Biomass % Error
p5 <- ggplot(data = all_biomass, aes(x = Md.DS, y = sperr, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "Reference Biomass (kg)",
    y = "", # "% TLS Biomass Error",
    title = "Main Stem % Error",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Branch Biomass % Error
p6 <- ggplot(data = all_biomass, aes(x = Mb.DS, y = tperr, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "", # "% TLS Biomass Error",
    title = "Branch % Error",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))
```

```{r, echo=FALSE, fig.width=7, fig.height=4, fig.align='center', warning=FALSE, message=FALSE}
ggarrange(p1, p4, common.legend = TRUE, legend = "bottom")
ggarrange(p2, p5, common.legend = TRUE, legend = "bottom")
ggarrange(p3, p6, common.legend = TRUE, legend = "bottom")
```

```{r, echo=FALSE, message=FALSE, warning=FALSE}
### TOTALS #####################################################################
total_stats <- all_biomass %>%
  group_by(Version) %>%
  summarise(
    MRE.pct = mean(tperr),
    RMSE.kg = sqrt(mean(tdiff^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Mt.DS) * 100
  )

CCC_total <- all_biomass %>%
  group_by(Version) %>%
  ccc(Mt.TLS, Mt.DS) %>%
  select(Version, .estimate) %>%
  rename(CCC = 2)

total_stats <- left_join(total_stats, CCC_total, by = "Version") %>%
  mutate(type = "total")

### MAIN STEM ##################################################################
stem_stats <- all_biomass %>%
  group_by(Version) %>%
  summarise(
    MRE.pct = mean(sperr),
    RMSE.kg = sqrt(mean(sdif^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Md.DS) * 100
  )

CCC_stem <- all_biomass %>%
  group_by(Version) %>%
  ccc(Md.TLS, Md.DS) %>%
  select(Version, .estimate) %>%
  rename(CCC = 2)

stem_stats <- left_join(stem_stats, CCC_stem, by = "Version") %>%
  mutate(type = "stem")

### BRANCHES ###################################################################
branch_stats <- all_biomass %>%
  group_by(Version) %>%
  summarise(
    MRE.pct = mean(bperr),
    RMSE.kg = sqrt(mean(bdiff^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Mb.DS) * 100
  )

CCC_branch <- all_biomass %>%
  group_by(Version) %>%
  ccc(Mb.TLS, Mb.DS) %>%
  select(Version, .estimate) %>%
  rename(CCC = 2)

branch_stats <- left_join(branch_stats, CCC_branch, by = "Version") %>%
  mutate(type = "branch")

all_stats <- bind_rows(total_stats, stem_stats, branch_stats)

### TABLES ###################################################################
data <- all_stats %>%
  mutate_if(is.numeric, round, 3) %>%
  rename(
    "Mean Relative Error (%)" = MRE.pct,
    "RMSE (kg)" = RMSE.kg,
    "Relative RMSE (%)" = RRMSE.pct
  ) %>%
  pivot_longer(cols = 2:5, names_to = "metric", values_to = "value") %>%
  pivot_wider(names_from = type, values_from = value) %>%
  group_by(Version) %>%
  group_split()

tables <- vector(mode = "list", length(data))

data[[2]] %>%
  gt() %>%
  tab_header(title = "TreeQSM v2.3.0") %>%
  cols_hide(Version) %>%
  cols_label(
    metric = "Metric",
    total = "Total Woody AGB",
    stem = "Main Stem Biomass",
    branch = "Branch Biomass"
  ) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  tab_options(
    table.border.top.style = "hidden",
    table.border.bottom.style = "hidden",
    table_body.hlines.color = "white",
    table.font.color = "black",
    heading.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    table_body.border.bottom.color = "black"
  )

data[[1]] %>%
  gt() %>%
  tab_header(title = "Real Twig (TreeQSM v2.3.0)") %>%
  cols_hide(Version) %>%
  cols_label(
    metric = "Metric",
    total = "Total Woody AGB",
    stem = "Main Stem Biomass",
    branch = "Branch Biomass"
  ) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  tab_options(
    table.border.top.style = "hidden",
    table.border.bottom.style = "hidden",
    table_body.hlines.color = "white",
    table.font.color = "black",
    heading.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    table_body.border.bottom.color = "black"
  )
```

## SimpleForest

Real Twig was not tested with SimpleForest during its development. While Real Twig does improve volume estimates for SimpleForest versus its built in allometric corrections, there are still improvements to be made, as SimpleForest QSM cylinders are generally much more overestimated than TreeQSM cylinders, making the identification of "good" cylinders difficult.

Below are the mass estimates and statistics using SimpleForest v5.3.2 with its built in vessel volume correction, and Real Twig applied to the same QSMs.

```{r, echo=FALSE, results='hide', message=FALSE, warning=FALSE}
# Import Data
file <- system.file("extdata/validation.csv", package = "rTwig")
all_biomass <- read.csv(file) %>%
  rename(Version = version)

# Distinct colors
factors <- distinct(all_biomass, Version)
cbPalette <- c("#009E73", "#CC79A7", "#56B4E9", "#999999", "#E69F00", "#009E73", "#F0E442", "#0072B2", "#D55E00")
cbPalette <- cbPalette[1:nrow(factors)]
factors <- bind_cols(factors, cbPalette) %>%
  rename(colors = 2)

# Filter Data
all_biomass <- all_biomass %>%
  left_join(factors, by = join_by(Version)) %>%
  filter(Version %in% c("SimpleForest", "Real Twig (SimpleForest)")) %>%
  mutate(id = paste0(SpCode, Tree))

# Total Biomass
p1 <- ggplot(data = all_biomass, aes(x = Mt.DS, y = Mt.TLS, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name), ) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_abline() +
  stat_poly_eq(aes(
    label = paste(after_stat(eq.label), after_stat(rr.label), sep = "~~~")
  )) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "TLS Biomass (kg)",
    title = "Total",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Main Stem Biomass
p2 <- ggplot(data = all_biomass, aes(x = Md.DS, y = Md.TLS, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_abline() +
  stat_poly_eq(aes(
    label = paste(after_stat(eq.label), after_stat(rr.label), sep = "~~~")
  )) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "", # "TLS Biomass (kg)",
    title = "Main Stem",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Branch Biomass
p3 <- ggplot(data = all_biomass, aes(x = Mb.DS, y = Mb.TLS, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_abline() +
  stat_poly_eq(aes(
    label = paste(after_stat(eq.label), after_stat(rr.label), sep = "~~~")
  )) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "", # "TLS Biomass (kg)",
    title = "Branch",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Total Biomass % Error
p4 <- ggplot(data = all_biomass, aes(x = Mt.DS, y = tperr, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "% TLS Biomass Error",
    title = "Total % Error",
    shape = "Speces",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Main Stem Biomass % Error
p5 <- ggplot(data = all_biomass, aes(x = Md.DS, y = sperr, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "Reference Biomass (kg)",
    y = "", # "% TLS Biomass Error",
    title = "Main Stem % Error",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))

# Branch Biomass % Error
p6 <- ggplot(data = all_biomass, aes(x = Mb.DS, y = tperr, color = Version, fill = Version)) +
  geom_point(aes(shape = scientific.name)) +
  geom_smooth(method = "lm", formula = y ~ x) +
  geom_hline(yintercept = 0) +
  labs(
    x = "", # "Reference Biomass (kg)",
    y = "", # "% TLS Biomass Error",
    title = "Branch % Error",
    shape = "Species",
    color = "Method",
    fill = "Method"
  ) +
  theme_classic() +
  scale_color_manual(values = unique(all_biomass$colors)) +
  scale_fill_manual(values = unique(all_biomass$colors))
```

```{r, echo=FALSE, fig.width=7, fig.height=4, fig.align='center', warning=FALSE, message=FALSE}
ggarrange(p1, p4, common.legend = TRUE, legend = "bottom")
ggarrange(p2, p5, common.legend = TRUE, legend = "bottom")
ggarrange(p3, p6, common.legend = TRUE, legend = "bottom")
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}
### TOTALS #####################################################################
total_stats <- all_biomass %>%
  group_by(Version) %>%
  summarise(
    MRE.pct = mean(tperr),
    RMSE.kg = sqrt(mean(tdiff^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Mt.DS) * 100
  )

CCC_total <- all_biomass %>%
  group_by(Version) %>%
  ccc(Mt.TLS, Mt.DS) %>%
  select(Version, .estimate) %>%
  rename(CCC = 2)

total_stats <- left_join(total_stats, CCC_total, by = "Version") %>%
  mutate(type = "total")

### MAIN STEM ##################################################################
stem_stats <- all_biomass %>%
  group_by(Version) %>%
  summarise(
    MRE.pct = mean(sperr),
    RMSE.kg = sqrt(mean(sdif^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Md.DS) * 100
  )

CCC_stem <- all_biomass %>%
  group_by(Version) %>%
  ccc(Md.TLS, Md.DS) %>%
  select(Version, .estimate) %>%
  rename(CCC = 2)

stem_stats <- left_join(stem_stats, CCC_stem, by = "Version") %>%
  mutate(type = "stem")

### BRANCHES ###################################################################
branch_stats <- all_biomass %>%
  group_by(Version) %>%
  summarise(
    MRE.pct = mean(bperr),
    RMSE.kg = sqrt(mean(bdiff^2, na.rm = TRUE)),
    RRMSE.pct = RMSE.kg / mean(Mb.DS) * 100
  )

CCC_branch <- all_biomass %>%
  group_by(Version) %>%
  ccc(Mb.TLS, Mb.DS) %>%
  select(Version, .estimate) %>%
  rename(CCC = 2)

branch_stats <- left_join(branch_stats, CCC_branch, by = "Version") %>%
  mutate(type = "branch")

all_stats <- bind_rows(total_stats, stem_stats, branch_stats)

### TABLES ###################################################################
data <- all_stats %>%
  mutate_if(is.numeric, round, 3) %>%
  rename(
    "Mean Relative Error (%)" = MRE.pct,
    "RMSE (kg)" = RMSE.kg,
    "Relative RMSE (%)" = RRMSE.pct
  ) %>%
  pivot_longer(cols = 2:5, names_to = "metric", values_to = "value") %>%
  pivot_wider(names_from = type, values_from = value) %>%
  group_by(Version) %>%
  group_split()

tables <- vector(mode = "list", length(data))

data[[2]] %>%
  gt() %>%
  tab_header(title = "SimpleForest") %>%
  cols_hide(Version) %>%
  cols_label(
    metric = "Metric",
    total = "Total Woody AGB",
    stem = "Main Stem Biomass",
    branch = "Branch Biomass"
  ) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  tab_options(
    table.border.top.style = "hidden",
    table.border.bottom.style = "hidden",
    table_body.hlines.color = "white",
    table.font.color = "black",
    heading.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    table_body.border.bottom.color = "black"
  )

data[[1]] %>%
  gt() %>%
  tab_header(title = "Real Twig (SimpleForest)") %>%
  cols_hide(Version) %>%
  cols_label(
    metric = "Metric",
    total = "Total Woody AGB",
    stem = "Main Stem Biomass",
    branch = "Branch Biomass"
  ) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  tab_options(
    table.border.top.style = "hidden",
    table.border.bottom.style = "hidden",
    table_body.hlines.color = "white",
    table.font.color = "black",
    heading.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    table_body.border.bottom.color = "black"
  )
```