---
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 a high quality reference data set that were both laser scanner and destructively sampled. The laser scanning was done in leaf-off conditions with a Riegl VZ-400. The destructive sampling contains total branch and main stem dry mass, and also basic density for both the main stem and the branches. We used different versions of TreeQSM with the same input parameters per tree to test our model.

## TreeQSM v2.4.1

Below are the mass estimates and statistics using TreeQSM v2.4.1 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.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

It is important to note that 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"
  )
```