---
title: "The `forestecology` R package for fitting and assessing neighborhood models of the effect of interspecific competition on the growth of trees"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{The `forestecology` R package for fitting and assessing neighborhood models of the effect of interspecific competition on the growth of trees}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  fig.width = 16 / 2,
  fig.height = 9 / 2,
  message = FALSE,
  warning = FALSE,
  cache = FALSE
)
set.seed(76)
```

# Purpose

This package is designed to work for spatially mapped, repeat censused forests plots. The package has commands to fit models of tree growth based on neighborhood competition which can be used to estimate species-specific competition coefficients. The model fits can then be evaluated using a spatial cross-validation scheme to detect possible overfitting. Additionally, these models can test whether the species identity of competitors matters using a permutation test-style shuffling of competitor identity (under the null hypothesis) and subsequently evaluating if model performance changes. See Allen and Kim (2020) [A permutation test and spatial cross-validation approach to assess models of interspecific competition between trees](https://doi.org/10.1371/journal.pone.0229930) for a full description.  


# Example analysis

We present an example analysis using toy data pre-loaded into the package where we will:

1. Compute growth of trees based on census data
1. Add spatial information
1. Identify all focal and corresponding competitor trees
1. Fit model and make predictions
1. Run spatial cross-validation

```{r}
library(ggplot2)
library(dplyr)
library(tidyr)
library(stringr)
library(forestecology)
library(patchwork)
library(blockCV)

# Resolve conflicting functions
filter <- dplyr::filter
```

## Compute growth of trees based on census data

The starting point of our analysis are data from two repeat censuses `census_1_ex` and `census_2_ex`. For example, consider the forest census data in `census_1_ex`.

```{r}
census_1_ex
```

We convert the `census_1_ex` data frame to an object of type `sf` and then plot using `geom_sf()`. 

```{r}
ggplot() +
  geom_sf(
    data = census_1_ex %>% sf::st_as_sf(coords = c("gx", "gy")),
    aes(col = sp, size = dbh)
  )
```

We first combine data from two repeat censuses into a single `growth` data frame that has the average annual growth of all trees alive at both censuses that aren't resprouts at the second census per Allen and Kim (2020).

```{r}
growth_ex <-
  compute_growth(
    census_1 = census_1_ex %>% 
      mutate(sp = to_any_case(sp) %>% factor()),
    census_2 = census_2_ex %>% 
      filter(!str_detect(codes, "R")) %>% 
      mutate(sp = to_any_case(sp) %>% factor()),
    id = "ID"
  ) %>% 
  # Compute basal area:
  mutate(basal_area = 0.0001 * pi * (dbh1 / 2)^2)
```

## Add spatial information

Our growth model assumes that two individual trees compete if they are less than a pre-specified distance `comp_dist` apart. Furthermore, we define a buffer region of size `comp_dist` from the boundary of the study region.

```{r}
# Set competitor distance
comp_dist <- 1

# Add buffer variable to growth data frame
growth_ex <- growth_ex %>%
  add_buffer_variable(size = comp_dist, region = study_region_ex)

# Optional: Create sf representation of buffer region
buffer_region <- study_region_ex %>% 
  compute_buffer_region(size = comp_dist)
```

In the visualization below, the solid line represents the boundary of the study region while the dashed line delimits the buffer region within. All trees outside this buffer region (in red) will be our "focal" trees of interest in our model since we have complete competitor information on all of them. All trees inside this buffer region (in blue) will only be considered as "competitor" trees to "focal" trees.

```{r}
base_plot <- ggplot() +
  geom_sf(data = study_region_ex, fill = "transparent") +
  geom_sf(data = buffer_region, fill = "transparent", linetype = "dashed")

base_plot + 
  geom_sf(data = growth_ex, aes(col = buffer), size = 2)
```

Next we add information pertaining to our spatial cross-validation scheme. We first manually define the spatial blocks that will act as our cross-validation folds and convert them to an `sf` object using the `sf_polygon()` function from the `sfheaders` package.

```{r}
fold1 <- rbind(c(0, 0), c(5, 0), c(5, 5), c(0, 5), c(0, 0))
fold2 <- rbind(c(5, 0), c(10, 0), c(10, 5), c(5, 5), c(5, 0))

blocks_ex <- bind_rows(
  sf_polygon(fold1),
  sf_polygon(fold2)
) %>%
  mutate(folds = c(1, 2) %>% factor())
```

Next we assign each tree to the correct folds using the `foldID` variable of the output returned by the `spatialBlock()` function from the [`blockCV`](https://github.com/rvalavi/blockCV) package.

```{r}
SpatialBlock_ex <- blockCV::spatialBlock(
  speciesData = growth_ex, k = 2, selection = "systematic", blocks = blocks_ex,
  showBlocks = FALSE, verbose = FALSE
)

growth_ex <- growth_ex %>%
  mutate(foldID = SpatialBlock_ex$foldID %>% factor())
```

In the visualization below, the spatial blocks that act as our cross-validation folds are delineated in orange. The shape of each point indicates which fold each tree has been assigned to. 

```{r}
base_plot + 
  geom_sf(data = growth_ex, aes(col = buffer, shape = foldID), size = 2) +
  geom_sf(data = blocks_ex, fill = "transparent", col = "orange")
```

## Compute focal versus competitor tree information

Based on our `growth` data frame, we now explicitly define all "focal" trees and their respective "competitor" trees in a `focal_vs_comp` data frame. This data frame has rows corresponding to each focal tree, and all information about its competitors are saved in the list-column variable `comp`. We implemented this nested format using `nest()` in order to minimize redundancy, given that the same tree can act as a competitor multiple times. 

```{r}
focal_vs_comp_ex <- growth_ex %>%
  create_focal_vs_comp(comp_dist, blocks = blocks_ex, id = "ID", comp_x_var = "basal_area")
focal_vs_comp_ex
```

Using `unnest()` we can fully expand the competitor information saved in the `focal_vs_comp` data frame. For example, the tree with `focal_ID` equal to 2 located at (1.5, 2.5) has two competitors within `comp_dist` distance from it.  

```{r}
focal_vs_comp_ex %>% 
  unnest(cols = "comp")
```

## Fit model and make predictions

We then fit our competitor growth model as specified in Allen and Kim (2020).

```{r}
comp_bayes_lm_ex <- focal_vs_comp_ex %>%
  comp_bayes_lm(prior_param = NULL)
```

The resulting output is an `comp_bayes_lm` object containing the posterior distribution of all linear regression parameters, the intercept, the slope for dbh for each species, and a matrix of all species pairs competitive effects on growth. The S3 object class is associated with several methods.

```{r}
# Print
comp_bayes_lm_ex

# Posterior distributions (plots combined with patchwork pkg)
p1 <- autoplot(comp_bayes_lm_ex, type = "intercepts")
p2 <- autoplot(comp_bayes_lm_ex, type = "dbh_slopes")
p3 <- autoplot(comp_bayes_lm_ex, type = "competition")
(p1 | p2) / p3
```

Furthermore, we can apply a `predict()` method to the resulting `comp_bayes_lm` object to obtain fitted/predicted values of this model. We append these `growth_hat` values to our `focal_vs_comp` data frame.

```{r}
focal_vs_comp_ex <- focal_vs_comp_ex %>%
  mutate(growth_hat = predict(comp_bayes_lm_ex, newdata = focal_vs_comp_ex))
focal_vs_comp_ex
```

We then compute the root mean squared error (RMSE) of the observed versus fitted growths as a measure of our model's fit.

```{r}
focal_vs_comp_ex %>%
  rmse(truth = growth, estimate = growth_hat) %>%
  pull(.estimate)
```

## Run spatial cross-validation

Whereas in our example above we fit our model to the entirety of the data and then generate fitted/predicted growths on this same data, we now apply the same model with spatial cross-validation. All the trees in a given fold will be given a turn as the "test" data while the trees in all remaining folds will be the "training" data. We then fit the model to the training data, but compute fitted/predicted growths for the separate and independent data. 

```{r}
focal_vs_comp_ex <- focal_vs_comp_ex %>%
  run_cv(comp_dist = comp_dist, blocks = blocks_ex)
```

Note the increase in RMSE, reflecting the fact that our original estimate of model error was overly optimistic as it did not account for spatial autocorrelation.

```{r}
focal_vs_comp_ex %>%
  rmse(truth = growth, estimate = growth_hat) %>%
  pull(.estimate)
```