08_variability.rmd

---
title: Variability
bibliography: [reference/a972.bib]
---

```{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 08-variability-1}
require(tidyverse)
library(rgl)
setupKnitr()
```


# Structural variability

Here I'll explore various metrics identified in the literature for quantifying and comparing structural variability and species variability

First I'll identify trees of interest as live trees including ingrowth. I'm not going to consider initial conditions here.

Details for calculations can be found on the [Define new data](03_define.html) page.

```{r 08-variability-2}

ingrowth <- mutate(ingrowth,
  define_xy(ingrowth),
  ba = pi * (dbh / 2)^2,
  ht_p = ht,
  year = factor(year, ordered = FALSE)
)

sci_d <- d_l %>%
  mutate(year = factor(year, ordered = FALSE)) %>%
  filter(live) %>%
  bind_rows(ingrowth, .id = "cohort") %>%
  filter(year != "init")

```

Here is the graphing function and emmeans generating function I'll use for visualizing these metrics. The mixed model I'm using is of the form: `plot_level_metric = treatment + year + treatment * year` with plot as a random variable. Any metric which does not inherently produce a plot level value is averaged at the plot level. 

```{r 08-variability-3}

my_dodge <- position_dodge(0.4)

var_fig <- function(data) {
  data |>
    filter(year != "init") |>
    relevel_treatment() |>
    ggplot(aes(year, emmean, color = treatment, group = treatment)) +
      geom_line(size = 1, position = my_dodge) +
      geom_point(position = my_dodge) +
      geom_errorbar(
        aes(ymin = emmean - SE, ymax = emmean + SE),
        width = .75,
        position = my_dodge
      ) +
      theme_bw() +
      labs(y = NULL, x = "Year")
}

var_mean <- function(plot_data) {
  mod <- lmer(val ~ treatment * year + (1 | plot), data = plot_data)
  emmeans(mod,  ~ treatment + year) |> as.data.frame()
}

var_sig <- function(plot_data) {
  mod <- lmer(val ~ treatment * year + (1 | plot), data = plot_data)
  emmeans(mod, pairwise ~ treatment | year)
}

var_table <- function(mean_data, caption) {
  mean_data |>
    # filter(year != "13") |>
    group_by(treatment) |>
    arrange(treatment) |>
    color_groups(digits = 2, caption)
}

var_mod <- function(plot_data, cap) {
  mod <- lmer(val ~ treatment * year + (1 | plot), data = plot_data)
  sjPlot::tab_model(
  mod,
  show.ci = FALSE,
  show.se = TRUE,
  show.p = TRUE,
  # dv.labels = paste("Formula", to_compare),
  title = cap
  )
}

```

# Structural complexity index

SCI [@zennerNewMethodModeling2000] measures the spatial variability in a given tree attribute such as DBH or height. It is driven by differences between neighboring trees. We can calculate SCI based on residual trees only, or include ingrowth. The latter would be expected to increase the SCI, but it is unclear whether or not this kind of structural diversity is of significant importance.

It is important to keep in mind that SCI, for our study, refers to heterogeneity on the within-patch scale: at the .08 ha scale. Larger plots would be needed if we are to characterize structural variability at greater scales.

Finally, I realized that Zenner and following researchers definitions of SCI is a little confusing because they confuse bars around a vector as meaning absolute value, when in fact it means the norm of a vector (specifically the L~2~ norm, or )

# SCI summarized by treatment

```{r 08-variability-4}

sci_plot <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = calc_sci(x, y, dbh))

sci_mean <- var_mean(sci_plot)
  
var_table(sci_mean, caption = "Average SCI and confidence interval for all live trees including ingrowth")

var_fig(sci_mean)

var_mod(sci_plot, "Model for SCI")
```

Some of these trajectories look interesting, but They don't look unique, given the error bars.

# Mean directional index

Mean directional index is described slightly differently by @corral-rivasPermutationTestSpatial2010 and @pommereningReconstructingSpatialTree2008 but I believe that they are equivalent representations. I'm using the latter with four neighbors:

$$\sqrt{ \left ( \sum_{j=1}^{4} \cos{a_{ij}} \right )^2 + \left ( \sum_{j=1}^{4} \sin{a_{ij}} \right )^2 } $$

Values should range between 0 for a grided arrangement to 4. A completely random arrange results in a value of about 1.8 [@corral-rivasPermutationTestSpatial2010].

Also, with this index, I am dropping extra coincident points, keeping the first.

```{r 08-variability-5}

mdi_plot <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = mean(calc_mdi(x, y)))

mdi_mean <- var_mean(mdi_plot)
  
var_table(mdi_mean, caption = "Average mean directional index (0: grid, 4:complete aggregation")

var_fig(mdi_mean)

var_mod(mdi_plot, "model for MDI")
```

# DBH differentiation (TD)

This index is defined by @fuldnerStrukturbeschreibungBuchenEdellaubholzMischwaldern1995 as:

$$ TD_i = 1 - \frac{min(DBH_i, DBH_j)}{max(DBH_i, DBH_j)} $$

This is the single closest neighbor version, there is also a N neighbor version, Khuenhe used this one and found it was highly correlated with SCI and DBHsd

@pommereningApproachesQuantifyingForest2002 provided interpretations of this metric, in our case suggesting *average to small differentiation*, where trees are generally more than 50% (closer to 70%) the size of their neighbor. 

```{r 08-variability-6}

TD_plot <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = mean(calc_TD(x, y, dbh)))

TD_mean <- var_mean(TD_plot)
  
var_table(TD_mean, caption = "Average DBH differentiation, 0: similar tree sizes, approaching 1: more difference in neighbor tree size")

var_fig(TD_mean)

var_mod(TD_plot, "Model for TD")

```

# DBH~sd~

Standard deviation of DBH is a straightforward metric for estimating structural heterogeneity. 

```{r 08-variability-7}

dbhsd_plot <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = sd(dbh))

dbhsd_mean <- var_mean(dbhsd_plot)
  
var_table(dbhsd_mean, caption = "Standard deviation of DBH")

var_fig(dbhsd_mean)

var_mod(dbhsd_plot, "Model for dbh sd")

```

# Shannon-Weaver index (overstory)

This index calculates tree species diversity based on the number of stems or basal area [@shannonMathematicalTheoryCommunication2001]

```{r 08-variability-8}

sh_plot <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = calc_sh(spp))

sh_mean <- var_mean(sh_plot)
  
var_table(sh_mean, caption = "Average Shannon-Weaver index (stem based)")

var_fig(sh_mean)

sh_plot_ba <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = calc_sh(spp, val = ba))

sh_mean_ba <- var_mean(sh_plot_ba)
  
var_table(sh_mean_ba, caption = "Basal area based average Shannon-Weaver index")

```


# Mean DBH (DBHm)

Another metric we were going to look at is DBHm, I had previously presented this as a boxplot of trees aggregated at the treatment level. This is different than all of our other analysis which base inference off of plot-level averages, or otherwise control for the variance between plots. Here I'll look at mean dbh following the same analysis format used for the other structural variables.

```{r 08-variability-9}

dbhm_plot <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = mean(dbh))

dbhm_mean <- var_mean(dbhm_plot)
  
var_table(dbhm_mean, caption = "Mean of DBH")

var_fig(dbhm_mean)

var_mod(dbhm_plot, "Model for mean of DBH")

```

# Ducey's ratio

Here we have Ducey's ratio of the additive method of calculating SDI and the method proposed by Reineke. Potential values depend on the truncation of the underlying distribution, our min and max are 10 and 83, which, according to @duceyRatioAdditiveTraditional2009 suggests that our potential lower limit might be around 0.74 for a stand where 95% of the trees were 10 cm and the rest were 83.

Our values run from about 0.95 to about 0.98. 

```{r 08-variability-10}

ducey_plot <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = sdi_sum(dbh) / sdi_r(dbh))

ducey_mean <- var_mean(ducey_plot)

var_sig(ducey_plot)
  
var_table(ducey_mean, caption = "Ratio of SDI* to SDI")

var_fig(ducey_mean)

var_mod(ducey_plot, "Model for Ducey's ratio")
```

# Coefficient of variation

The coefficient of variation for a plot tells us about variation in tree size relative to the mean. Plots with smaller average tree size, but similar ranges of variation should result in higher CV. 

## DBH

```{r 08-variability-11}

CVdbh_plot <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = sd(dbh) / mean(dbh), DFp = sum(spp == "PSMEM") / n())

CVdbh_mean <- var_mean(CVdbh_plot)
  
var_table(CVdbh_mean, caption = "Coefficient of variation for DBH")

var_fig(CVdbh_mean) + labs(title = "Coefficient of variation for DBH")

var_sig(CVdbh_plot)$contrasts |> as.data.frame() |> kbl2()

m1 <- lmer(val ~ treatment * year + (1 | plot), data = CVdbh_plot)
m2 <- lmer(val ~ treatment * year + DFp + (1 | plot), data = CVdbh_plot)
m3 <- lmer(val ~ treatment + DFp + (1 | plot), data = CVdbh_plot)

fl <- map_chr(list(m1, m2, m3), ~ format(formula(.x)))
data.frame(formula = fl) |> kbl2(row.names = TRUE)

sjPlot::tab_model(m1, m2, m3, show.aic = TRUE, dv.labels = 1:3)

m2mean <- emmeans(m2, ~ treatment + year) |> as.data.frame()
m1mean <- emmeans(m1, ~ treatment + year) |> as.data.frame()

bind_rows(m1 = m1mean, m2 = m2mean, .id = "mod") |>
  var_fig() +
  facet_wrap(~ mod) +
  labs(title = "Compare m1 and m2 (w/o & w DFp) for CVdbh")


```

### CVdbh model

*We used the simple model with just `treatment * year` (m1) for this metric*

```{r}

var_mod(CVdbh_plot, "Model for CVdbh")

```

## Height

```{r 08-variability-12}

CVht_plot <- sci_d %>%
  group_by(treatment, year, plot) %>%
  summarize(val = sd(ht) / mean(ht), DFp = sum(spp == "PSMEM") / n())

CVht_mean <- var_mean(CVht_plot)
  
var_table(CVht_mean, caption = "Coefficient of variation for height")

var_fig(CVht_mean) + labs(title = "Coefficient of variation for height")

var_sig(CVht_plot)$contrasts |> as.data.frame() |> kbl2(digits = 3)

m1 <- lmer(val ~ treatment * year + (1 | plot), data = CVht_plot)
m2 <- lmer(val ~ treatment * year + DFp + (1 | plot), data = CVht_plot)
m3 <- lmer(val ~ treatment + DFp + (1 | plot), data = CVht_plot)

fl <- map_chr(list(m1, m2, m3), ~ format(formula(.x)))
data.frame(formula = fl) |> kbl2(row.names = TRUE)

sjPlot::tab_model(m1, m2, m3, show.aic = TRUE, dv.labels = 1:3)

```

# Combined diversity figure

```{r 08-variability-13}
structural_all <- bind_rows(
  SCI = sci_mean,
  MDI = mdi_mean,
  TD = TD_mean,
  DBHcv = CVdbh_mean,
  DBHsd = dbhsd_mean,
  `SDI*/SDI` = ducey_mean,
  .id = "metric"
) %>%
  mutate(metric = factor(metric,
    levels = c("SDI*/SDI", "DBHsd", "DBHcv", "SCI", "MDI", "TD")
  ))
```

```{r 08-variability-14}

my_dodge <- position_dodge(0.5)

structural_all %>%
  relevel_treatment() %>%
  mutate(year = paste0("20", year)) %>%
  ggplot(aes(year, emmean, group = treatment)) +
    geom_line(aes(linetype = treatment), position = my_dodge) +
    geom_point(aes(shape = treatment), size = 1.8, position = my_dodge) +
    geom_errorbar(
      aes(ymin = emmean - SE, ymax = emmean + SE),
      width = .4,
      position = my_dodge,
      alpha = 0.4
    ) +
    theme_bw() +
    labs(y = NULL, x = "Year") +
    facet_wrap(~ metric, scales = "free_y", strip.position = "left") +
    theme(
      panel.spacing.x = unit(0, "mm"),
      strip.background.y = element_blank(),
      strip.switch.pad.wrap = unit(0, "mm"),
      strip.placement = "outside",
      legend.position = "bottom",
      legend.title = element_blank(),
      legend.box.margin = margin(),
      legend.box.spacing = unit(0, "mm"),
      plot.margin = margin(1,1, unit = "mm"),
      axis.title = element_text(size = 10)
    ) +
    scale_linetype_manual(
      values = c(
        C40 = "solid",
        L40 = "dashed",
        C80 = "solid",
        L80 = "dashed",
        Control = "dotted"
      )
    ) +
    scale_shape_manual(
      values = c(
        C40 = "circle",
        L40 = "triangle",
        C80 = "circle open",
        L80 = "triangle open",
        Control = "asterisk"
      )
    )

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

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

# Understory species

## Richness

Species richness is just the count of unique species. I'll do this separately for shrub and herb species. The pattern here might suggest that crews had different skill levels in identifying herbs, or it could have to do with different seasons of observation. It is unexpected that the control has the same trend.

I used a random slope linear model here, I'm not sure this is appropriate as we are dealing with count data. I did run a glmmm with Poisson family (various links) and these all had a significantly higher BIC.

```{r 08-variability-15}

shrub_richness_plot <- veg_species %>%
  filter(covertype == "Shrub") |>
  count(name = "val", treatment, year, plot) |>
  rbind(list("L40", "08", "2L40", 0))

shrub_richness_mean <- var_mean(shrub_richness_plot)
  
var_table(shrub_richness_mean, caption = "Average understory shrub species counts")

var_fig(shrub_richness_mean)

var_mod(shrub_richness_plot, "Model for shrub richness")

```



```{r}

herb_richness_plot <- veg_species %>%
  filter(covertype == "Herb") |>
  count(name = "val", treatment, year, plot)

herb_richness_mean <- var_mean(herb_richness_plot)
  
var_table(herb_richness_mean, caption = "Average understory herb species counts")

var_fig(herb_richness_mean)

var_mod(herb_richness_plot, "Model for herb richness")

```


```{r richness-testing, eval=FALSE}
mm <- lmer(
  val ~ treatment * year + (1 | plot),
  data = herb_richness_plot
  # ,family = poisson(link = "log")
)

mm1 <- glmer(
  val ~ treatment * year + (1 | plot),
  data = herb_richness_plot
  , family = poisson(link = "log")
)

mm2 <- lm(
  val ~ treatment * year,
  data = herb_richness_plot
)

emmeans(mm2, ~ treatment + year, type = "response") |>
  as.data.frame() |>
  # rename(emmean = rate) |>
  var_fig()

emmeans(mm1, ~ treatment + year, type = "response") |>
  as.data.frame() |>
  # rename(emmean = rate) |>
  var_fig()

summary(mm2)

shrub_richness_plot |>
  group_by(treatment, year) |>
  summarise(emmean = mean(val), SE = sd(val) / sqrt(n()))
```

## Percent cover

```{r}

shrub_cover_percent <- veg_cover |>
  rename(val = cl_shrub)

shrub_cover_mean <- var_mean(shrub_cover_percent)

var_fig(shrub_cover_mean)

var_mod(shrub_cover_percent, "Model for shrub cover")


herb_cover_percent <- veg_cover |>
  rename(val = cl_herb)

herb_cover_mean <- var_mean(herb_cover_percent)

var_fig(herb_cover_mean)

var_mod(herb_cover_percent, "Model for herb cover")

```

## Understory cover

What are the dominant shrubs and herbs for each treatment/year?

It is hard to say because sword fern are treated (mostly) as herbs for the first two entries, and then as shrubs for the last, but even this is inconsistent, some third entries have sword fern listed as an herb.

```{r 08-variability-16}
spp_common_name <- function(species) {
  case_when(
    species == "POMU"     ~ "Sword fern",
    species == "IRDO"     ~ "Douglas' iris",
    species == "VISE3"    ~ "Evergreen violet",
    species == "STAJR"    ~ "Hedge nettle",
    species == "TRLA6"    ~ "Starflower",
    species == "GRASS"    ~ "Grass",
    species == "WHMO"     ~ "Modesty",
    species == "CLSI2"    ~ "Candyflower",
    species == "VAOV2"    ~ "Evergreen huckleberry",
    species == "Fireweed" ~ "Fireweed",
    species == "RUSP"     ~ "Salmonberry",
    species == "GASH"     ~ "Salal",
    species == "BENE2"    ~ "Oregon grape", 
    species == "SARA"     ~ "Red elderberry"
  )
}

veg_cover %>%
  group_by(treatment, year) %>%
  summarize(
    shrub = paste(dom_shrub, collapse = ", "),
    herb = paste(dom_herb, collapse = ", ")
  ) %>% arrange(year) %>% group_by(year) %>% color_groups(caption = "Dominant herb/shrub for four plots in each treatment/year. Evergreen huckleberey (VAOV2) and sword fern (POMU) are treated inconsistently.")

```

## Dominant understory species change

What is interesting here is that the high intensity thinnings lead consistently to heavy salmon berry shrub cover and the rest lead to heavy sword fern cover. Evergreen violet seems to be the most common herb following all treatments (except control) but the H40 and L40 appear to be more consistent/homogenous across plots.

```{r 08-variability-17}

doj <- position_dodge(width = 0.5)
veg_cover %>% 
  # filter(year != "18") %>%
  select(-c(cl_con, cl_hw, dom_con, dom_hw)) %>%
  pivot_longer(-c(year, plot, treatment, note), names_to = c(".value", "type"), names_sep = "_") %>%
  filter(type == "shrub") %>%
  mutate(
    plot2 = str_extract(plot, "\\d"),
    dom = spp_common_name(dom),
    type = if_else(type == "shrub", "Shrub", "Herb"),
    year = recode(year, `08` = "2008", `13` = "2013", `18` = "2018")
  ) %>%
  relevel_treatment() %>%
  arrange(treatment, plot, year) %>%
  ggplot(aes(year, dom, group = plot)) +
    geom_line(position = doj) +
    geom_point(aes(size = cl), position = doj, alpha = 0.5) +
    facet_grid(~ treatment, scales = "free_y") +
    theme_bw() +
    theme(
      # rect = element_rect(fill = "gray70"),
      axis.title.y = element_blank(),
      # legend.position = c(.9, .7),
      # legend.key.size = unit(2, "mm"),
      legend.box.spacing = unit(0, "mm"),
      legend.position = "bottom",
      legend.background = element_blank(),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 9)
    ) +
    scale_size_continuous(name = "Cover (%)") + # , breaks = c(15, 30, 60), limits = c(2.5, 87.5)) +
    scale_x_discrete(expand = expansion(mult = 0.25)) +
    scale_y_discrete(expand = expansion(mult = 0.2)) +
    labs(x = "Year")

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

ggsave(
  filename = "figs/understory_change.jpg",
  width = 18.2,
  height = 5.2,
  units = "cm"
)

```

change in percent cover of shrubs and herbs over time for each treatment. The fact that herb and shrub cover change so much is really confounded by the fact that sword fern (and to a lesser extent, huckleberry) are classified differently over time.

```{r 08-variability-18, include=FALSE}

# In this figure, I use mean and standard error to show variation, but in the rest of this analysis, I use emmeans of a random effects model. I'm leaving this in for now, but not including it in output.

doj <- position_dodge(width = 0.3)
veg_cover %>%
  group_by(treatment, year) %>%
  summarise(
    herb_mean = mean(cl_herb),
    herb_se = sd(cl_herb) / sqrt(n()),
    shrub_mean = mean(cl_shrub),
    shrub_se = sd(cl_shrub) / sqrt(n())
  ) %>%
  mutate(year = as.factor(year)) %>%
  pivot_longer(
    -c(treatment, year),
    names_to = c("veg_type", ".value"),
    names_sep = "_"
  ) %>%
  ggplot(aes(year, mean, color = treatment, group = treatment)) +
    geom_point(position = doj) +
    facet_wrap(~ veg_type) +
    geom_line(position = doj) +
    geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5, position = doj) +
    labs(caption = "average percent cover of shrub and herb, errorbars are 1 std. error")

```

I think what I need to do is just look at 2018 and show the total percent cover and dominant species in the same barplot. This plot shows shrub cover for each plot, colored with the dominant species for that plot.

```{r 08-variability-19}
veg_cover %>%
  filter(year == "18") %>%
  mutate(dom_shrub = spp_common_name(dom_shrub)) %>%
  ggplot(aes(treatment, fill = dom_shrub, group = fct_reorder(plot, cl_shrub))) +
    geom_col(aes(y = cl_shrub), position = "dodge") +
    labs(caption = "one bar for each plot indicating dominant species and percent cover")

veg_cover %>%
  filter(year == "18") %>%
  mutate(dom_herb = spp_common_name(dom_herb)) %>%
  ggplot(aes(treatment, fill = dom_herb, group = fct_reorder(plot, cl_herb))) +
    geom_col(aes(y = cl_herb), position = "dodge") +
    labs(caption = "one bar for each plot indicating dominant species and percent cover")
```

# Final Veg Figure

This figure displays emmeans of modeled percent cover with `treatment:year` interaction, and random intercept for `plot`. Shrubs include tanoak and cascara. Herbs include 

Error bars represent 1 SE of emmean.

```{r 08-variability-22, fig.width=fig_w(1, "in"), fig.height=9/2.54}

my_dodge <- position_dodge(0.5)

cover_all <- bind_rows(
  Shrub = shrub_cover_mean,
  Herbaceous = herb_cover_mean,
  .id = "type"
)

richness_all <- bind_rows(
  Shrub = shrub_richness_mean,
  Herbaceous = herb_richness_mean,
  .id = "type"
)
veg_all <- bind_rows(
  "Cover (%)" = cover_all,
  "Richness (count)" = richness_all,
  .id = "measure"
)

veg_all |>
  relevel_treatment() |>
  ggplot(aes(year, emmean, group = treatment)) +
    geom_line(aes(linetype = treatment), position = my_dodge) +
    geom_point(aes(shape = treatment), size = 1.8, position = my_dodge) +
    geom_errorbar(
      aes(ymin = emmean - SE, ymax = emmean + SE),
      width = .4,
      position = my_dodge,
      alpha = 0.5
    ) +
    facet_grid(measure ~ type, scales = "free_y", switch = "y") +
    theme_bw() +
    theme(
      panel.spacing.x = unit(0, "mm"),
      strip.background.y = element_blank(),
      strip.switch.pad.wrap = unit(0, "mm"),
      strip.placement = "outside",
      legend.position = "bottom",
      legend.title = element_blank(),
      legend.box.margin = margin(),
      legend.box.spacing = unit(0, "mm"),
      plot.margin = margin(1,1, unit = "mm"),
      axis.title = element_text(size = 10),
      strip.text.y = element_text(size = 10),
      legend.key.size = unit(4.8, "mm")
    ) +
    scale_linetype_manual(
      values = c(
        C40 = "solid",
        L40 = "dashed",
        C80 = "solid",
        L80 = "dashed",
        Control = "dotted"
      )
    ) +
    scale_shape_manual(
      values = c(
        C40 = "circle",
        L40 = "triangle",
        C80 = "circle open",
        L80 = "triangle open",
        Control = "asterisk"
      )
    ) +
    labs(y = NULL, x = "Year") +
    scale_x_discrete(
      labels = c("2008", "2013", "2018"),
      expand = expansion(mult = .2)
    )

ggsave(
  filename = "figs/veg_rich_cover.pdf",
  device = cairo_pdf,
  width = fig_w(1, "cm"),
  height = 9.5,
  units = "cm"
)

ggsave(
  filename = "figs/veg_rich_cover.jpg",
  width = fig_w(1, "cm"),
  height = 9.5,
  units = "cm"
)

```


This figure inlcudes percent cover and change in species richness.

```{r 08-variability-23, fig.width=3.5, fig.height=3.2}

veg_richness_plot <- veg_species %>%
  group_by(treatment, year, plot) %>%
  summarize(val = n())

veg_richness_mean <- var_mean(veg_richness_plot)

veg <- bind_rows(
  `Shrub cover (%)` = shrub_cover_mean,
  `Richness (count)` = veg_richness_mean,
  .id = "metric") |>
  mutate(
    metric = factor(metric, levels = c("Shrub cover (%)", "Richness (count)"))
  )

my_dodge <- position_dodge(0.5)

veg %>%
  relevel_treatment() %>%
  mutate(year = paste0("20", year)) %>%
ggplot(aes(year, emmean, group = treatment)) +
    geom_line(aes(linetype = treatment), position = my_dodge) +
    geom_point(aes(shape = treatment), size = 1.8, position = my_dodge) +
    geom_errorbar(
      aes(ymin = emmean - SE, ymax = emmean + SE),
      width = .4,
      position = my_dodge,
      alpha = 0.4
    ) +
    theme_bw() +
    labs(y = NULL, x = "Year") +
    facet_wrap(~ metric, scales = "free_y", strip.position = "left", ncol = 2) +
    theme(
      panel.spacing.x = unit(0, "mm"),
      strip.background.y = element_blank(),
      strip.switch.pad.wrap = unit(0, "mm"),
      strip.placement = "outside",
      strip.text = element_text(size = 10),
      axis.title = element_text(size = 10),
      # legend.position = c(0.17, 0.91),
      legend.position = "bottom",
      legend.title = element_blank(),
      legend.box.margin = margin(l = 0, unit = "mm"),
      legend.box.spacing = unit(0, "mm"),
      legend.background = element_blank(),
      legend.key.size = unit(4.8, "mm"),
      plot.margin = margin(1,1, unit = "mm")
    ) +
    scale_linetype_manual(
      values = c(
        C40 = "solid",
        L40 = "dashed",
        C80 = "solid",
        L80 = "dashed",
        Control = "dotted"
      )
    ) +
    scale_shape_manual(
      values = c(
        C40 = "circle",
        L40 = "triangle",
        C80 = "circle open",
        L80 = "triangle open",
        Control = "asterisk"
      )
    ) +
    scale_x_discrete(expand = expansion(mult = .2))

ggsave(
  filename = "figs/shrub_cover.pdf",
  device = cairo_pdf,
  width = 8.84,
  height = 8,
  units = "cm"
)

ggsave(
  filename = "figs/shrub_cover.jpg",
  width = 8.84,
  height = 8,
  units = "cm"
)
```

# Fun visual of SCI

Here I visualize the SCI_dbh for 10 random plots

```{r 08-variability-24, fig.width=8, fig.height=10, cache=TRUE}

sci_plots_d <- sci_d %>% select(x, y, z = dbh, plot, year) %>% split(~ plot + year, drop = TRUE)


visualize_plot_sdi <- function(data1) {
  data1 <- remove_coincident(data1)
  x <- data1$x - min(data1$x)
  y <- data1$y - min(data1$y)
  z <- data1$z
  col <- cm.colors(20)[1 + round(19*(z - min(z))/diff(range(z)))]
  dxy <- interp::tri.mesh(x, y, duplicate = "remove")
  persp3d(dxy, z, col = col, meshColor = "vertices")
  title3d(main = paste(data1$plot[[1]], data1$year[[1]]), edge = "x--")
  axes3d(xlen = 0, ylen = 0)
}

# visualize_plot_sdi(sci_plots_d[[30]])

open3d()
mfrow3d(nr = 5, nc = 2, sharedMouse = FALSE) 

for (i in sample(1:length(sci_plots_d), 10)) {
  next3d()
  visualize_plot_sdi(sci_plots_d[[i]])
}

rglwidget()

```

I can use the function described in `?rgl::persp3d.triSht` which uses the output of either package tripack or interp

Here are boxplots showing the variability in 2018 of some potentially important variables.

```{r 08-variability-25}
sci_d %>%
  filter(year == "18") %>%
  select(treatment, spp, dbh, ht_p, cr, ht_inc2) %>%
  pivot_longer(
    !c(treatment, spp),
    names_to = "measure",
    values_to = "value"
  ) %>%
  ggplot(aes(treatment, value)) +
    geom_boxplot() +
    facet_wrap(~ measure, scales = "free_y")
```

# References