06_summary.rmd

---
title: "Summary statistics"
---

```{r setup, include=FALSE}
knitr::knit_hooks$set(time_it = time_it)
knitr::opts_chunk$set(echo = TRUE, warning = TRUE, message = FALSE, time_it = TRUE, cache.extra = "wombat")
```

```{r 06-summary-1}
library(tidyverse)
library(dplyr)
```

# Summary measures

I'm making the assumption now that all trees recorded as cut were alive in the initial conditions. For all of the following I use trees with status 1, 2, 3, 15, 16. Fifteen and 16 are cut trees, but in the long data format, they only occur on entries where year = "init" so they do not affect other periods. Status 1, 2 and 3 represent second-growth, advanced regen, and old-growth, respectively.

Generally, the `summarise()` command in each section defines the output of a given measure. The `mutate()` command is also used for calculating new columns. Calculations are done by grouping levels (generally by treatment, by year).

```{r 06-summary-2}
# I'll store all my metrics in this list
measures <- list()
```

## Density

I'll look at a quick summary of quartiles of the number of tree observations by plot for each year. Initial plots had between 83 and 164 observations and subsequent years are roughly parallel. Density is being calculated for all live trees, including those broken above dbh.

```{r 06-summary-3}

d_l %>%
  filter(status %in% c(1, 2, 3, 15, 16, 31)) %>%
  group_by(year, plot) %>%
  summarize(n = n()) %>%
  summarize(as_tibble_row(quantile(n)))

measures$density <- function(data, species) {
  if(is.null(species)) species <- unique(data$spp)
  data %>%
    filter(spp %in% species, live | status %in% c(15, 16)) %>%
    group_by(treatment, year, plot) %>%
    # expansion factor 12.5 (12.5 * .08 = 1 ha)
    summarise(density_plot = n() * 12.5) %>%
    summarise(
      mean = mean(density_plot),
      n = n(),
      std = sd(density_plot)
    )
}
```

## Mortality

New mortality here is ascribed to observations of the current period which were alive in the previous period and are not alive in the current period and reported in tph. Currently, this is undefined for the first period ("init") and is zero for "08" (because "init" and "08" are really only one observation). The values for 2013 and 2018 record the new mortality in the previous period, i.e.: 08-13 and 13-18, respectively.

```{r 06-summary-4}
measures$mort <- function(data, species) {
  if(is.null(species)) species <- unique(data$spp)
  data %>%
    filter(spp %in% species) %>%
    group_by(tree_id) %>%
    mutate(mort = ((live == FALSE) * (lag(live) == TRUE)) > 0) %>%
    group_by(year, treatment, plot) %>%
    summarize(mort_plot = sum(mort) * 12.5) %>%
    summarize(
      mean = mean(mort_plot),
      n = n(),
      std = sd(mort_plot)
    )
}
```

## Qmd

DBH is actually QMD here. It will be calculated based on all live trees, including those broken above dbh and with broken and dead tops (most of these are SESE and assumed to be alive).

```{r 06-summary-5}
measures$dbh <- function(data, species) {
  if(is.null(species)) species <- unique(data$spp)
  data %>%
    filter(spp %in% species, live | status %in% c(15, 16)) %>%
    group_by(treatment, year, plot) %>%
    summarise(plot_mean = sqrt(mean(dbh^2))) %>%
    summarise(
      mean = mean(plot_mean),
      n = n(),
      std = sd(plot_mean)
    )
}

measures$dom_dbh <- function(data, species) {
  if(is.null(species)) species <- unique(data$spp)
  data %>%
    filter(live | status %in% c(15, 16)) %>%
    group_by(plot, year) %>%
    slice_max(order_by = dbh, n = 8, with_ties = FALSE) %>%
    filter(spp %in% species) %>%
    group_by(treatment, year) %>%
    summarise(
      mean = mean(dbh),
      n = n(),
      std = sd(dbh))
}
```

## Heights

Heights will be determined using only trees with unbroken, live crowns.

Because there are about 14 observations in 2008 and 2013 that do not have height measurements, I am omitting these from analysis for now. Also, I am only looking at trees that have height measurements in 2018. This means that of the 3,538 total height observations, I am only looking at 1,929 over the the three measurement periods.

For dominant heights, I am using dominant dbh trees *not including broken and dead tops*. Not all of these have heights. For this reason, there are fewer than 4 trees per plot.

```{r 06-summary-6}

# get ids of trees that have heights in 2018
ht_id_2018 <- filter(d_l, year == "18", !is.na(ht)) %>% pull(tree_id)

# average height I AM OMMITING MISSING HEIGHT MEASUREMENTS.
measures$ht <- function(data, species) {
  if(is.null(species)) species <- unique(data$spp)
  data %>%
    filter(
      spp %in% species,
      live | status %in% c(15, 16),
      !get_cond(c(2, 3)),
      tree_id %in% ht_id_2018,
      !is.na(ht)
    ) %>%
  group_by(treatment, year, plot) %>%
  summarise(plot_mean = mean(ht)) %>%
  summarise(
    mean = mean(plot_mean),
    n = n(),
    std = sd(plot_mean)
  )
}

# Dominant height
# I AM CALCULATING DOM HEIGHT ONLY FROM THOSE AVAILABLE
measures$dom_ht <- function(data, species) {
  if(is.null(species)) species <- unique(data$spp)
  data %>%
    filter(
      live | status %in% c(15, 16),
      !get_cond(2, 3),
      tree_id %in% ht_id_2018
    ) %>%
    group_by(plot, year) %>%
    slice_max(order_by = dbh, n = 8, with_ties = FALSE) %>%
    filter(spp %in% species) %>%
    group_by(treatment, year) %>%
    summarise(
      mean = mean(ht, na.rm = TRUE),
      n = sum(!is.na(ht)),
      std = sd(ht, na.rm = TRUE)
    )
}

```

## Basal Area

Basal area, similar to density will be calculated for all live trees.

```{r 06-summary-7}

measures$ba <- function(data, species) {
  if(is.null(species)) species <- unique(data$spp)
  data %>%
    filter(spp %in% species, live | status %in% c(15, 16)) %>%
    group_by(treatment, year, plot) %>%
    summarise(
      # I'm calculating BA (m^2/ha) by plot
      ba_plot = sum(dbh^2 * pi / 40000, na.rm = TRUE) * 12.5,
      .groups = "drop_last"
    ) %>%
    summarize(
      # And averagging over plots within a treatment
      mean = mean(ba_plot),
      n = n(),
      std = sd(ba_plot),
      .groups = "drop_last"
    )
}
```
### Compare to Teraoka (2016)

I want to create another function for calculating basal area. This should return basal area in m^2^/ha:

1. by species
2. for all species
3. for a subset of species
4. it should include all live + harvested trees

It turns out that this function is almost identical to my previous one, except that it allows to calculate by species.

```{r 06-summary-8}
library(tidyverse)

get_ba <- function(data, by_species = FALSE, species_subset = NULL) {
  if(is.null(species_subset)) species_subset <- unique(data$spp)
  if(by_species) {
    grp_vars <- rlang::exprs(treatment, spp, year, plot)
  } else {
    grp_vars <- rlang::exprs(treatment, year, plot)
  }
  data %>%
    filter(spp %in% species_subset, live | status %in% c(15, 16)) %>%
    group_by(!!!grp_vars) %>%
    summarise(
      # I'm calculating BA (m^2/ha) by plot
      ba_plot = sum(dbh^2 * pi / 40000, na.rm = TRUE) * 12.5,
      n_ba = n(),
      .groups = "drop_last"
    ) %>%
    summarize(
      # And averagging over plots within a treatment
      mean = mean(ba_plot),
      n = sum(n_ba),
      std = sd(ba_plot),
      .groups = "drop_last"
    )
}
```

Here I'll plot BA by species and year, similar to Teraoka (2016)

```{r 06-summary-9}
teraoka_colors <- c("#04315c", "#EB4332", "#86E981", "#557AF0", "#B3AFAE")
get_ba(d_l, by_species = TRUE) %>%
  filter(year %in% c("08", "13")) %>%
  mutate(
    treatment = factor(treatment, levels = c("C", "H40", "L40", "H80", "L80")),
    spp = factor(spp, levels = c("TSHE", "SESE3", "PSMEM", "PISI", "ALRU2"))
  ) %>%
  ggplot(aes(x = treatment, y = mean)) +
    geom_col(aes(fill = spp)) +
    facet_wrap(vars(year)) +
    labs(x = "Treatment", fill = "Species", y = "BA (m^2 / ha)") +
    scale_fill_manual(values = teraoka_colors)
```

It looks pretty similar, but there are discrepancies. Next I'll print out the percent changes in basal area, again similar to Teraoka et al. (2016), but for both 5-year growth periods.

I define a percent change function which takes output from one of my metric functions (expecting the variable "mean_[variable_name]")

```{r 06-summary-10}
pct_change <- function(data, var) {
  mutate(data, pct_change = ({{var}} - lag({{var}})) / lag({{var}})) %>%
  select(-starts_with(c("n", "sd")))
}

get_ba(d_l) %>%
  pct_change(mean)
```

## SDI

For SDI, I use the same trees that were used for basal area and density: all standing trees (broken or not) that are taller than breast height.

```{r 06-summary-11}
# do any target trees have missing dbh?
any(nrow(filter(d_l, status %in% c(1, 2, 3, 15, 16), is.na(dbh))))

measures$sdi <- function(data, species) {
  if(is.null(species)) species <- unique(data$spp)
  data %>%
    filter(spp %in% species, live | status %in% c(15, 16)) %>%
    group_by(treatment, year, plot) %>%
    summarise(
      sdi_plot = sum(12.5 * (dbh * 0.04)^1.605, na.rm = TRUE)
    ) %>%
    summarize(
      mean = mean(sdi_plot),
      n = n(),
      std = sd(sdi_plot)
    )
}

```

## Percent redwood by BA

```{r}

measures$p_redwood <- function(data, species) {
  if(is.null(species)) species <- unique(data$spp)
  data %>%
    filter(spp %in% species, live | status %in% c(15, 16)) %>%
    group_by(treatment, year, plot) %>%
    mutate(redwood_ba = if_else(spp == "SESE3", ba, 0)) %>%
    summarise(
      p_redwood_plot = 100 * sum(redwood_ba) / sum(ba)
    ) %>%
    summarize(
      mean = mean(p_redwood_plot),
      n = n(),
      std = sd(p_redwood_plot)
    )
}

```

## Density of dominant

```{r 06-summary-12}

```


# Initial summary tables

## All species

```{r 06-summary-13, echo=FALSE}
# Call all summary functions with a given dataset and speces list
# and join the results into a wide dataframe.

get_measures <- function(data = d_l, species = NULL) {
  measures %>%
  map(exec, data, species) %>%
  bind_rows(.id = "measure")
}

# takes a long dataset of summary values

make_summary <- function(data, measures = c("dbh", "ht", "density", "ba", "sdi")) {
  data %>%
  filter(measure %in% measures) %>%
  relevel_treatment() %>%
  relabel_measure() %>%
  relocate(measure) %>%
  mutate(across(where(is.double), ~ sprintf("%.1f", .x))) %>%
  mutate(mean = paste0(mean, " (", std, ")")) %>%
  select(measure, treatment, year, mean) %>%
  pivot_wider(names_from = treatment, values_from = mean, names_sort = TRUE) 
}

```

```{r 06-summary-14}

get_measures(d_l) %>%
  make_summary() %>%
  kableExtra::kbl() %>%
  kableExtra::kable_classic(full_width = FALSE, html_font = "Garamond") %>%
  kableExtra::add_header_above(c(" " = 2, "Treatment" = 5)) %>%
  kableExtra::collapse_rows(columns = 1, valign = "top")

```

## Redwood

```{r 06-summary-15}

get_measures(species = "SESE3") %>%
  make_summary() %>%
  kableExtra::kbl() %>%
  kableExtra::kable_classic(full_width = FALSE, html_font = "Garamond") %>%
  kableExtra::add_header_above(c(" " = 2, "Treatment" = 5)) %>%
  kableExtra::collapse_rows(columns = 1, valign = "top")

```

## Douglas-fir

```{r 06-summary-16}

get_measures(species = "PSMEM") %>%
  make_summary() %>%
  kableExtra::kbl() %>%
  kableExtra::kable_classic(full_width = FALSE, html_font = "Garamond") %>%
  kableExtra::add_header_above(c(" " = 2, "Treatment" = 5)) %>%
  kableExtra::collapse_rows(columns = 1, valign = "top")

```

# Figures

```{r sum17, cache=TRUE, fig.height = 5.91, fig.width = 7.16, echo = FALSE, message=FALSE, warning=FALSE}

fix_year <- function(year) {
  case_when(year == "init" ~2007.5,
            year == "08" ~2008,
            year == "13" ~2013,
            year == "18" ~2018)
}

fig_data <- bind_rows(
  tibble(spp = "SESE3", get_measures(species = "SESE3")),
  tibble(spp = "PSMEM", get_measures(species = "PSMEM"))
)

# choose which stand summary variables to include
fig_include_vars <- c("density", "dbh", "ba", "ht")

make_fig <- function(data = fig_data) {
  yr_breaks <- c(2008, 2013, 2018)
  data %>%
    filter(measure %in% fig_include_vars, year != "init") %>%
    relevel_treatment() %>%
    relabel_measure(fig = TRUE) %>%
    mutate(
      year = fix_year(year),
      spp = recode(spp, SESE3 = "Redwood", PSMEM = "Douglas-fir")) %>%
    ggplot(aes(year, mean, color = spp, shape = spp)) +
      geom_line(position = position_dodge(width = .5)) +
      geom_point(position = position_dodge(width = .5)) +
      geom_errorbar(
        aes(
          ymin = mean - qt(0.975, 3) * std / sqrt(n),
          ymax = mean + qt(0.975, 3) * std / sqrt(n)),
        width = .75,
        position = position_dodge(.5)
      ) +
      facet_grid(
        vars(measure),
        vars(treatment),
        scales = "free_y",
        switch = "y",
        labeller = labeller(measure = label_parsed)
      ) +
      theme_bw() +
      theme(
        panel.spacing.x = unit(1, "lines"),
        strip.background.y = element_blank(),
        strip.placement = "outside",
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.box.spacing = element_blank()
      ) +
      scale_x_continuous(breaks = yr_breaks) +
      # scale_y_continuous(breaks = scales::breaks_extended(4)) +
      scale_color_manual("spp",
        values = c(Redwood = "black", "Douglas-fir" = "#969696")
      ) +
      scale_shape_manual("spp",
        values = c(Redwood = 16, "Douglas-fir" = 17)
      ) +
      labs(y = NULL, x = "Year")
}

make_fig()

ggsave(
  filename = "figs/summary_fig.pdf",
  device = cairo_pdf,
  width = 18.2,
  height = 15,
  units = "cm"
)

ggsave(
  filename = "figs/summary_fig.jpg",
  width = 18.2,
  height = 15,
  units = "cm"
)
```