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

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 7,
  dpi = 100,
  out.width = "95%"
)
```

## Summary

- [Description](#description)
- [Getting ready](#getting-ready)
- [Loading example data](#loading-example-data)
- [Simulating a random community](#simulating-a-random-community)
  - [Effect of background density](#effect-of-background-density)
  - [Effect of proportion arguments](#effect-of-proportion-arguments)
- [Simulating nested communities](#simulating-nested-communities)
  - [Effect of proportion argument](#effect-of-proportion-argument)
  - [Effect of bias argument](#effect-of-bias-argument)
- [Simulating niche conservatism in communities](#simulating-niche-conservatism-in-communities)
  - [Effect of background density](#effect-of-background-density)
  - [Effect of proportion arguments](#effect-of-proportion-arguments)
- [Predictions for communities](#predictions-for-communities)
  - [Predict to data frames](#predictions-to-data-frames)
  - [Predict to SpatRaster](#predictions-to-spatraster)
  - [Truncating predictions](#truncating-predictions)
  - [Simple community outcomes](#simple-community-outcomes)
- [Save and import](#save-and-import)

<hr>

## Description

This vignette shows how to simulate virtual communities using functions in the nicheR package. For our purposes, we will consider a community a set of species niches (ellipses) that are distributed in a given environmental space (E-space). The simulation of virtual communities is useful to generate hypothetical scenarios of community assembly and explore patterns derived from the way species are distributed in environmental and geographic space (G-space). 

The main functions that automate community simulations in `nicheR` are:

- `random_ellipses()`: generates a set of random ellipses in the environmental space provided. A scenario in which niche conservatism is not present.
- `nested_ellipses()`: generates a set of nested ellipses taking the reference niche as the biggest possible ellipse. A scenario in which all species niches are contained within the reference niche, and each is contained within the previous one.
- `conserved_ellipses()`: generates a set of ellipses in the environmental space provided aiming for a set of ellipses that are similar to the reference niche. A scenario of niche conservatism.

<br>

## Getting ready

If `nicheR`  has not been installed yet, please do so. See the [Main guide](https://castanedam.github.io/nicheR/index.html) for installation instructions.

Use the following lines of code to load `nicheR` and other packages needed for this vignette, and to set a working directory (if necessary).

<u>Note</u>: We will display functions from other packages as `package::function()`. 

```{r get_ready, results='hide', message=FALSE, warning=FALSE}
# Load packages
library(nicheR)
library(terra)

# Current directory
getwd()

# Define new directory
#setwd("YOUR/DIRECTORY")  # modify if setting a new directory

# Saving original plotting parameters
original_par <- par(no.readonly = TRUE)
```

<br>

## Loading example data

The lines of code below help us load the data to run our example community simulations. The data are included in the `nicheR` package and consist of an `nicheR_ellipsoid` object with a reference niche defined by two environmental variables (bio1 and bio12) and these environmental variables for North America to represent the background for our simulations. For details on how to create the reference niche, see the vignette [1. Build ellipsoid](https://castanedam.github.io/nicheR/articles/creating_ellipsoid_based_niches.html).

```{r data}
# Reference niche
data("ref_ellipse", package = "nicheR")

# Background data
data("back_data", package = "nicheR")

# Raster layers for predictions
ma_bios <- terra::rast(system.file("extdata", "ma_bios.tif",
                                   package = "nicheR"))
```

<br>

Let's take a look at the reference niche and the background data to get familiar with them before we start simulating communities. The raster layers will be used later for predictions, so we will just check them here to make sure they are loaded correctly.

```{r data_check}
# Check reference niche
print(ref_ellipse)

# Check background data
head(back_data)

# Check the raster layers
ma_bios
```

<br>

Now, let's plot the reference niche and the background data to visualize the environmental space in which we will simulate our communities. We will use the `plot_ellipsoid()` function to plot the reference niche and the background data together.

```{r data_plot}
# Pick the variables for the background data
vars <- c("bio_1", "bio_12")

# Plotting the background data to visualize the environmental space
mars <- c(4, 4, 2, 1)
par(mar = mars)  # adjust margins for better visualization

plot_ellipsoid(ref_ellipse, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797", lwd = 2,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Reference Niche and Background Environmental Space")
```

<br>

## Simulating random communities

The `random_ellipses()` function helps to create a community of species with niches that are randomly distributed in the environmental space. The size and shape of these niches are constrained by the reference niche but they vary. Below we show an general example of how to use this function.

```{r random}
# Simulating the community
rand_comm <- random_ellipses(object = ref_ellipse,
                             background = back_data[, vars],
                             n = 20)  # number of species in the community

# check the a few details from the generated community
names(rand_comm)  # elements in the community object

print(rand_comm)  # a summary of the elements in the community object
```

<br> 

Now let's plot the generated community to visualize the distribution of the random ellipses in environmental space. We will use the `plot_community()` function for this purpose.

```{r random_plot}
# Plotting the community
par(mar = mars)  # adjust margins for better visualization

plot_community(rand_comm, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Community of Random Ellipses")
```

<br>

### Effect of background density

As the `random_ellipses()` function uses the background as a reference to pick ellipse centroids, **the density of the background data can affect the distribution of these ellipses**. If the background data is dense in certain areas of the environmental space, it may lead to a higher concentration of ellipses in those areas, while sparse background data may result in fewer ellipses being generated in those regions. This can influence the overall structure and diversity of the simulated community, as well as the patterns observed in niche overlap and species interactions.

The example below shows how the density of the background can affect the distribution of random ellipses in the environmental space. We will simulate two communities: (1) using the full background as reference, and (2) using the arguments thin_background and resolution to reduce the effect of uneven point density. We will increase the number of species in the community to better visualize the effect of background density on the distribution of ellipses.

```{r random_den}
# Simulating the community with the full background
rand_comm_full <- random_ellipses(object = ref_ellipse,
                                  background = back_data[, vars], n = 25)

# Simulating the community with a thinned background
rand_comm_thin <- random_ellipses(object = ref_ellipse,
                                  background = back_data[, vars], n = 25,
                                  thin_background = TRUE, resolution = 20)
```

<br>

Let's check the distribution of ellipses in the two communities with a plot to visualize the effect of the arguments `thin_background` and `resolution`.

```{r random_den_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"}  
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community with the full background
plot_community(rand_comm_full, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community from Full Background")

## Plotting the community with the thinned background
plot_community(rand_comm_thin, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community from Thinned Background")
```

<br>

As you can see in the plots, the community generated with the full background has a higher concentration of ellipses in areas where the background data is denser, whereas the community created with thinned background has a more even distribution of ellipses. This highlights the importance of considering this factor when simulating communities using random ellipses. 

Play with the value for the argument `resolution` to see how it can affect the distribution of ellipses in the community. Keep in mind that the largest the value for this argument, the more points will be available for the function to pick centroids.

<br>

### Effect of proportion arguments

The size of the ellipses in the communities is determined by the reference niche and the arguments `smallest_proportion` and `largest_proportion` in the `random_ellipses()` function. Based on the reference niche, playing with the values for the arguments `smallest_proportion` and `largest_proportion` can help explore scenarios with . This can affect elements of the community structure such as the level of niche overlap observed among the species (ellipses).

The example below shows how these proportions can affect the coverage and overlap of random ellipses in environmental space. Considering that `smallest_proportion` must be smaller than `largest_proportion`, we will simulate four communities: (1) both arguments are small; and (2) both arguments are large; (3) `smallest_proportion` is large and `largest_proportion` is small; and (4) `smallest_proportion` is small and `largest_proportion` is large. We will thin the background in both cases to reduce the effect of uneven point density and better visualize the effect of the proportion arguments.

```{r random_prop}
# Community with both arguments small
rand_comm1 <- random_ellipses(object = ref_ellipse,
                              background = back_data[, vars], n = 25,
                              thin_background = TRUE, resolution = 20,
                              smallest_proportion = 0.1,
                              largest_proportion = 0.5)

# Community with both arguments large
rand_comm2 <- random_ellipses(object = ref_ellipse,
                              background = back_data[, vars], n = 25,
                              thin_background = TRUE, resolution = 20,
                              smallest_proportion = 0.7,
                              largest_proportion = 1.5)

# Community with smallest_proportion large and largest_proportion small
rand_comm3 <- random_ellipses(object = ref_ellipse,
                              background = back_data[, vars], n = 25,
                              thin_background = TRUE, resolution = 20,
                              smallest_proportion = 0.5,
                              largest_proportion = 0.7)

# Community with smallest_proportion small and largest_proportion large
rand_comm4 <- random_ellipses(object = ref_ellipse,
                              background = back_data[, vars], n = 25,
                              thin_background = TRUE, resolution = 20,
                              smallest_proportion = 0.1,
                              largest_proportion = 1.5)
```

<br>

Let's check the distribution of ellipses in the two communities with a plot to visualize the effect of the proportion arguments.

```{r random_prop_plot}
# Plotting the communities
par(mfrow = c(2, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Community with small proportions
plot_community(rand_comm1, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community with Small Proportions")

## Community with large proportions
plot_community(rand_comm2, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community with Large Proportions")

## Community with large smallest_proportion and small largest_proportion
plot_community(rand_comm3, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community with Large and Small")

## Community with small smallest_proportion and large largest_proportion
plot_community(rand_comm4, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community with Small and Large")
```

<br>

As shown in the plots, when both arguments are small, the ellipses are generally smaller and there is less overlap among them. When both arguments are large, the ellipses tend to be larger and there is more overlap among them. When `smallest_proportion` is large and `largest_proportion` is small, the ellipses tend to be of intermediate size. Finally, when `smallest_proportion` is small and `largest_proportion` is large, the ellipses tend to be more variable in size. This highlights how the size of the ellipses can affect the structure of the community and the "interactions" among species niches (ellipses). 

Play with the values for these arguments to explore different scenarios and pick the ones that are more convenient for your research.

<br> 

## Simulating nested communities

The `nested_ellipses()` function helps to create a community of species with niches that are nested in the environmental space. The size and shape of these niches are constrained by the reference niche but they vary. Below we show an general example of how to use this function.

```{r nested}
# Simulating the community
nest_comm <- nested_ellipses(object = ref_ellipse, n = 20)

# check the a few details from the generated community
print(nest_comm)  # a summary of the elements in the community object
```

<br>

Now let's plot the generated community to visualize the distribution of the nested ellipses in environmental space. We will use the `plot_community()` function for this purpose.

```{r nested_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"}
# Plotting the community
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community of nested ellipses
plot_community(nest_comm,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Community of Nested Ellipses")

## Plotting the community of nested ellipses with the background
plot_community(nest_comm, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested Community with Background")
```

<br>

Despite the fact that background was not used to create nested ellipses, since the reference niche is located in the same environmental space, we can see how the generated ellipses are distributed in that space. The nested structure of the community is evident in both plots, with each ellipse being contained within the previous one. 

<br>

### Effect of proportion argument

An important argument in the `nested_ellipses()` function is `smallest_proportion`, which determines the size of the smallest nested ellipse to be generated, in relation to the reference niche. A smaller value for `smallest_proportion` will result in a greater range of ellipse sizes.

```{r nested_prop}
# Simulating the community with a small smallest_proportion
nest_comm_small <- nested_ellipses(object = ref_ellipse, n = 20,
                                   smallest_proportion = 0.1)

# Simulating the community with a large smallest_proportion
nest_comm_large <- nested_ellipses(object = ref_ellipse, n = 20,
                                   smallest_proportion = 0.6)

# Lets check the volume stats for the communities for comparison
## Communities with small smallest_proportion
mean(sapply(nest_comm_small$ellipse_community, function(x) x$volume))

## Communities with large smallest_proportion
mean(sapply(nest_comm_large$ellipse_community, function(x) x$volume))
```

<br>

Below we plot the two communities to visualize the effect of the argument `smallest_proportion` on the distribution of nested ellipses in environmental space.

```{r nested_prop_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"}
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community with a small smallest_proportion
plot_community(nest_comm_small,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested with Small Proportion")

## Plotting the community with a large smallest_proportion
plot_community(nest_comm_large,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested with Large Proportion")
```

<br>

From the plot, it is clear why the mean volume of the ellipses in the community with a small `smallest_proportion` is smaller than the mean volume of the ellipses in the community with a large `smallest_proportion`.

<br>

### Effect of bias argument

Another important arguments in the `nested_ellipses()` function is `bias`, which determines the degree of bias towards generating smaller or larger ellipses in the community. A value of `bias` < 1 and close to zero clusters ellipses toward the border of the reference ellipse, whereas a value greater than 1 clusters them toward the centroid of the reference ellipse. We will keep the value for `smallest_proportion` the same in both cases to better visualize the effect of the argument `bias`.

```{r nested_bias}
# Simulating the community with a bias towards the border
nest_comm_small_bias <- nested_ellipses(object = ref_ellipse, n = 20,
                                        smallest_proportion = 0.1, bias = 0.2)

# Simulating the community with a bias towards the centroid
nest_comm_large_bias <- nested_ellipses(object = ref_ellipse, n = 20,
                                        smallest_proportion = 0.1, bias = 2)

# Lets check the volume stats for the communities for comparison
## Community with bias towards the border
mean(sapply(nest_comm_small_bias$ellipse_community, function(x) x$volume))

## Community with bias towards the centroid
mean(sapply(nest_comm_large_bias$ellipse_community, function(x) x$volume))
```

<br>

Now, let's plot the two communities to visualize the effect of the argument `bias` on the distribution of nested ellipses in environmental space.

```{r nested_bias_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"}
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community with a bias towards the border
plot_community(nest_comm_small_bias,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested with Bias Towards Border")

## Plotting the community with a bias towards the centroid
plot_community(nest_comm_large_bias,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested with Bias Towards Centroid")
```

<br>

The plots clearly show the effect of the argument `bias` on the distribution of nested ellipses. A community created with `bias` values close to zero will tend to have larger volumes than one created with values larger than one. We used these values only for reference, play with the argument to explore different scenarios and pick the one that helps you explore your questions.

<br>

## Simulating niche conservatism in communities

The function `conserved_ellipses` creates a ellipses based on a reference ellipse aiming for a set of results that are similar to the reference. To do that centroids for the new ellipses are sampled from a background with a bias towards the centroid of the reference ellipse. A general example of how to use this function is shown below.

```{r conserved}
# Simulating the community
cons_comm <- conserved_ellipses(object = ref_ellipse,
                                background = back_data[, vars],
                                n = 20)

# check the a few details from the generated community
print(cons_comm)  # a summary of the elements in the community object
```

<br>

The plot below shows the generated community.

```{r conserved_plot}
# Plotting the community withe the background
par(mar = mars)  # adjust margins for better visualization

plot_community(cons_comm, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Community of Conserved Ellipses")
```

<br>

As you can see in the plot, the ellipses generated vary in size, direction, and position, but they don't go very far from the original ellipsoid. This follows the ideas of niche conservatism, which say that closely related species tend to have similar ecological niches. 

<br>

### Effect of background density

Similar to the `random_ellipses()` function, the `conserved_ellipses()` function uses the background as a reference to pick ellipse centroids. Differences in background density can affect the bias used to generate the ellipses (bias = tend to select new centroids close to the reference centroid). If the effect of background density is not controlled, it will be like saying that niche position evolution depends on the how close you are to the ancestor niche position and how common are environments. Keeping in mind that ancestry is not considered here, the trick is that those common environments are not necessarily geographically available as geography is not considered here.

The example below shows how the density of the background can affect the distribution of random ellipses when using the function. We will simulate two communities: (1) using the full background as reference, and (2) using the arguments `thin_background` and `resolution` to reduce the effect of uneven point density.

```{r cons_den}
# Simulating the community with the full background
cons_comm_full <- conserved_ellipses(object = ref_ellipse,
                                     background = back_data[, vars], n = 20)

# Simulating the community with a thinned background
cons_comm_thin <- conserved_ellipses(object = ref_ellipse,
                                     background = back_data[, vars], n = 20,
                                     thin_background = TRUE, resolution = 10)
```

<br>

Let's check the distribution of ellipses in the two communities with a plot to visualize the effect of the arguments `thin_background` and `resolution`.

```{r cons_den_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"}  
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community with the full background
plot_community(cons_comm_full, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Conserved Community from Full Background")

## Plotting the community with the thinned background
plot_community(cons_comm_thin, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Conserved Community from Thinned Background")
```

<br>

As you can see in the plots, the community generated with the full background has a higher concentration of ellipses in areas where the background data is denser compared to the one created with thinned background. This highlights the importance of considering this factor when simulating communities using the function `conserved_ellipses()`. Play with the value for the argument `resolution` to see how it can affect the distribution of ellipses in the community (larger values, more points available).

<br>

### Effect of proportion arguments

We have explored before the effect of the arguments that control the smallest and largest proportions when generating new ellipses considering the reference. We will explore two examples with the function `conserved_ellipses()`: (1) both arguments are small; and (2) `smallest_proportion` is small and `largest_proportion` is large. We will thin the background in both examples to better visualize the effect of the proportion arguments.

```{r cons_prop}
# Community with both arguments small
cons_comm1 <- conserved_ellipses(object = ref_ellipse,
                                 background = back_data[, vars], n = 20,
                                 thin_background = TRUE, resolution = 10,
                                 smallest_proportion = 0.1,
                                 largest_proportion = 0.5)

# Community with smallest_proportion small and largest_proportion large
cons_comm2 <- conserved_ellipses(object = ref_ellipse,
                                 background = back_data[, vars], n = 20,
                                 thin_background = TRUE, resolution = 10,
                                 smallest_proportion = 0.1,
                                 largest_proportion = 1.5)
```

<br>

Let's check the distribution of ellipses in the two communities with a plot to visualize the effect of the proportion arguments.

```{r cons_prop_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"}
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Community with small proportions
plot_community(cons_comm1, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Conserved Community with Small Proportions")

## Community with small smallest_proportion and large largest_proportion
plot_community(cons_comm2, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Conserved Community with Small and Large")
```

<br>

The plot show that when both arguments are small, the ellipses are generally smaller. When `smallest_proportion` is small and `largest_proportion` is large, the ellipses are more variable in size, but larger than in the previous case because of the values we used. A larger number of ellipses generated can help visualize the effect of the arguments on the sizes of the ellipses a little better. Play with the values for these arguments to explore different scenarios and pick the ones that are more convenient for your research.

<br> 

## Predictions for communities

Once we have our communities simulated using the functions explored above, we can predict over new data to observe the patterns of Mahalanobis distance and suitability derived from the ellipses. The function `predict()` can be used with `nicheR_community` objects to obtain predictions. The difference between predict for community objects and that for `nicheR_ellipsoid` objects (see [2. Make a Prediction](https://castanedam.github.io/nicheR/articles/predict.html) is that only one type of `prediction` can be obtained at a time. The options for the `prediction` argument for `nicheR_community` objects are:

- `Mahalanobis`: Mahalanobis distance from the centroid of ellipses to every point in `newdata`.
- `suitability`: multivariate normal probability for every point in `newdata`, derived from the Mahalanobis distance, interpreted as suitability. 
- `Mahalanobis_trunc`:  the same Mahalanobis distance but truncated to the limit of the ellipses. Values outside the ellipse are returned as `NA`.
- `suitability_trunc`: the same suitability but truncated to the limit of the ellipses. Values outside the ellipse are returned as zero.

Predictions can be produced for `newdata` in the form of `data.frame` or `SpatRaster` objects.
The `newdata` can contain multiple variables, but the names of the ones used to create the reference ellipse and the communities must be included. As long as those two variable names match with the ones in `newdata`, predictions are possible (i.e., predictions to distinct areas in the world are possible). We will explore implementations for both types of objects below.

<br>

### Predict to data frames

The most intuitive way to work with ellipses is in environmental space. This is because ecological niches are defined in this space. Since environmental values can easily be organized in a data fame, predictions for those objects are easy to obtain. See a quick example below, using one of the communities generated in which we assumed niche conservatism and background data as `newdata`.

```{r predict}
# Predicting Mahalanobis distances
maha_cons_pred <- predict(cons_comm, newdata = back_data[, vars],
                          prediction = "Mahalanobis")

# Predicting suitability
suit_cons_pred <- predict(cons_comm, newdata = back_data[, vars],
                          prediction = "suitability")

# Check Mahalanobis predictions
maha_cons_pred[1:5, 1:5]

# Check suiatbility predictions
suit_cons_pred[1:5, 1:5]
```

<br>

Now let's check how the results look like in a plot. Predictions are produced for all ellipses, and in the returned results, every prediction is added as a new column to `newdata`. We will plot results for only the first and second ellipses to show clearly the differences between Mahalanobis and suitability predictions. We will use the `plot_ellipsoid()` function to plot the predictions in environmental space. The argument `col_layer` is used to specify which prediction to use for the colors in the plot. We will use different color palettes for Mahalanobis and suitability predictions to better visualize the differences between them.

```{r predict_plot}
# Plotting the results
## Plotting area parameters
par(mfrow = c(2, 2), cex = 0.6, mar = mars)  # adjust margins for visualization

## plot mahalanobis distance predictions for ellipse 1
plot_ellipsoid(object = cons_comm[[3]][[1]],  # the relevant ellipse
               prediction = maha_cons_pred,  # mahalanobis distance prediction
               col_layer = "ell_1",  # ellipse prediction to use for colors
               pal = hcl.colors(100, palette = "Oslo"),  # color palette
               rev_pal = TRUE,  # reversing the palette
               col_ell = "#e10000", lwd = 2,  # color and line for ellipse
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Mahalanobis Distance Ellipse 1")

# plot mahalanobis distance predictions for ellipse 2
plot_ellipsoid(object = cons_comm[[3]][[2]],
               prediction = maha_cons_pred,
               col_layer = "ell_2",
               pal = hcl.colors(100, palette = "Oslo"),
               rev_pal = TRUE,
               col_ell = "#e10000", lwd = 2,  # color and line for ellipse
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Mahalanobis Distance Ellipse 2")

# plot suitability predictions for ellipse 1
plot_ellipsoid(object = cons_comm[[3]][[1]],
               prediction = suit_cons_pred,  # suitability prediction
               col_layer = "ell_1",
               pal = hcl.colors(100, palette = "Viridis"),
               col_ell = "#e10000", lwd = 2,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Suitability Ellipse 1")

# plot suitability predictions for ellipse 2
plot_ellipsoid(object = cons_comm[[3]][[2]],
               prediction = suit_cons_pred,  # suitability prediction
               col_layer = "ell_2",
               pal = hcl.colors(100, palette = "Viridis"),
               col_ell = "#e10000", lwd = 2,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Suitability Ellipse 2")
```

<br>

The plots shown are simple, but show clear patterns. For reference, dark colors mean low values for both, Mahalanobis distance and suitability. We can see that distance is larger the farther from the ellipse centroid, and suitability is higher closer to the centroid. That is the back bone of the theory behind using ellipses (or ellipsoids) as models of ecological niche. 

<br>

### Predict to SpatRaster

Now, let's predict suing `SpatRaster` objects as `newdata`. The important implication of this is that it allows us to see geographic representations of our predictions. The geographic patterns of our predictions are relevant for questions about species geographic distributions. The code below shows how to produce this predictions. We will use the same community and the raster layers from which our background derives (originally form [WorldClim](https://www.worldclim.org/), included as data in the `nicheR` package).

```{r predictr}
# Predicting Mahalanobis distances
maha_cons_predr <- predict(cons_comm, newdata = ma_bios,
                           prediction = "Mahalanobis")

# Predicting suitability
suit_cons_predr <- predict(cons_comm, newdata = ma_bios,
                           prediction = "suitability")

# Check Mahalanobis predictions
maha_cons_predr

# Check suiatbility predictions
suit_cons_predr
```

<br>

Let's check how the results look like in a plot.

```{r predictr_plot, fig.width = 7, fig.height = 4.5, out.width = "95%"}
# Plotting area parameters
par(mfrow = c(2, 2), cex = 0.6)  # adjust margins for visualization
marsr <- c(0.5, 0.5, 2, 4)

# Plots
terra::plot(maha_cons_predr$ell_1,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Mahalanobis Distance Ellipse 1")

terra::plot(maha_cons_predr$ell_2,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Mahalanobis Distance Ellipse 2")

terra::plot(suit_cons_predr$ell_1,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Suitability Ellipse 1")

terra::plot(suit_cons_predr$ell_2,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Suitability Ellipse 2")
```

<br>

The raster plots show us the geographic patterns of Mahalanobis distances and suitability, which are much more complex than the ones in environmental space. Exploring this plots can help identify clusters of areas with high suitability, for instance. High suitability is expected to favor the presence of a species, which is why geographic projections of ecological niche models are used for explorations of potential distributional areas. 

<br>

### Truncating predictions

An important decision to make when studying what conditions are good for a species to maintain populations for long periods of time is when environments stop being suitable. Here is where ellipsoids shine as models of suitability, because their formulation gives us that threshold. The limit of the ellipses we created, are the theoretical limits for what is suitable and not suitable.

The options for `prediction` in our `predict()` function include results for Mahalanobis and suitability that can be "truncated" using that ellipsoid limit. For Mahalanobis distances, truncation implies that all conditions outside the ellipsoid limit become NA. Whereas, for suitability, all values outside the limit become zero. This is important when trying to translate ideas of niche into potential distributions.

Let's explore an example with suitability predictions, truncated based on the ellipse limits.

```{r trunc}
# Predicting suitability truncated using a data frame
suit_cons_predt <- predict(cons_comm, newdata = back_data[, vars],
                           prediction = "suitability_trunc")

# Predicting suitability truncated using raster data
suit_cons_predrt <- predict(cons_comm, newdata = ma_bios,
                            prediction = "suitability_trunc")

# Check predictions in data.frame
suit_cons_predt[1:5, 1:5]

# Check predictions in raster
suit_cons_predrt
```

<br>

Compared to our previous predictions of suitability, now the truncated results include zero.

Let's explore the results using plots. First, the truncated results as they are produced.

```{r trunc_plot, fig.width = 7, fig.height = 3, out.width = "95%"}
# Plotting area parameters
par(mfrow = c(1, 2), cex = 0.6, mar = mars)

## the truncated results for the background
plot_ellipsoid(object = cons_comm[[3]][[1]],
               prediction = suit_cons_predt,  # suitability prediction
               col_layer = "ell_1",
               pal = hcl.colors(100, palette = "Viridis"),
               col_ell = "NA", lwd = 0,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Suitability Trunc. E Space")

## custom color palette for raster predictions
tr_cols <- c("gray", hcl.colors(100, palette = "Viridis"))
tr_breacks <- c(0, seq(1e-6, 1, length.out = 101))

## the truncated raster
terra::plot(suit_cons_predrt$ell_1,
            col = tr_cols, breaks = tr_breacks,
            type = "continuous",
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Suitability Trunc. G Space")
```

<br>
These plots of truncated suitability look very similar to the simple suitability ones. But this is just an artifact from color assigning in plotting. 

Now, let's plot suitable vs unsuitable environments and areas. We start by transforming everything inside the ellipse into one and what is outside remains as zero.

```{r trunc_plot2, fig.width = 7, fig.height = 3, out.width = "95%"}
# Obtaining values of zero and one
## results in data.frame
bin_suit_cons_predt <- suit_cons_predt
bin_suit_cons_predt[, -(1:2)] <- (bin_suit_cons_predt[, -(1:2)] > 0) * 1

## results in raster
bin_suit_cons_predrt <- suit_cons_predrt
bin_suit_cons_predrt <- (bin_suit_cons_predrt > 0) * 1


# Plotting
## Colors for suitability
bincol <- c("#c9c9c9", "#0004d5")

## Plotting area parameters
par(mfrow = c(1, 2), cex = 0.6, mar = mars)

## the binary for the background
plot(bin_suit_cons_predt[, vars],  # plot the variables as points
     col = bincol[as.factor(bin_suit_cons_predt$ell_1)],
     pch = 16, xlab = "Bio1 (Mean Annual Temperature)",
     ylab = "Bio12 (Annual Precipitation)",
     main = "Suitability Binary E Space")

## the binary raster
terra::plot(bin_suit_cons_predrt$ell_1, col = bincol,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Suitability Binary G Space")
```

<br> 

Now it is a lot more evident what was considered inside the ellipsoid in environmental and geographic space. 

<br>

### Simple community outcomes

He have shown examples for a couple of the ellipses generated, but this can be done for all the ones in the community. Once all the elements in the community are considered, interesting community level metrics and indices can be explored. Using the binary results obtained from the truncated predictions we will derive further results.

First, we already have a **Presence-absence matrix (PAM)** with our binary results as a data.frame. These type of matrices are commonly used in ecology because they tell us which species (columns) are where (row: each row is a site). 

Let's visualize our PAM below:

```{r pam_plot}
# Exclude sites with no species to make the plot easier
pam <- bin_suit_cons_predt[, -(1:2)]
pam <- as.matrix(pam[!rowSums(pam) == 0, ])

# Plot PAM using image
par(mar = c(1, 1, 5, 1), cex = 0.6)

image(1:ncol(pam), 1:nrow(pam), t(pam[nrow(pam):1, ]),
      col = bincol, axes = FALSE)
box()
text(x = 1:ncol(pam), y = nrow(pam) * 1.01,
     labels = colnames(pam),
     srt = 90, adj = 0, xpd = TRUE)

title(main = "Species Presence-Absence Matrix", line = 3.5)
```

<br> 

Now let's plot a simple **Species Richness Map**. We need to compute richness using our binary truncated results.

```{r richness, fig.width = 7, fig.height = 4.5, out.width = "95%"}
# Compute richness
richness <- terra::app(bin_suit_cons_predrt, sum)

# Plot richness
terra::plot(richness, mar = marsr, col = rev(heat.colors(20)),
            main = "Species Richness Map")
```

<br> 

A lot of other results can be derived from the community predictions. We hope these demonstrations have helped inspire new ideas.

<br>

```{r par_reset}
# Reset plotting parameters
par(original_par)
```

<br>

## Save and import

To facilitate writing nicheR objects to local directories and importing results in later sessions, we created the `save_nicheR` and `read_nicheR` functions. These functions can be used with the results obtained from virtual community simulations as `nicheR_community` objects, as shown below. 

```{r save_import, eval=FALSE}
# file name (in a temporary directory for demonstration purposes)
temp_file <- file.path(tempdir(), "conserved_community.rds")

# Save the community object to a local directory
save_nicheR(cons_comm, file = temp_file)

# Import the community object from a local directory
read_com <- read_nicheR(temp_file)
```

<br>

Results from predictions obtained as `data.frame` or `SpatRaster` can be saved using functions that are convenient for those types of objects. For instance, `write.csv()` for the data.frame and `terra::writeRaster()` for the raster results. Importing those files can be done with `read.csv()` and `terra::rast()`.

```{r save_import2, eval=FALSE}
# file names (in a temporary directory for demonstration purposes)
temp_df_file <- file.path(tempdir(), "df_cons_com_predictions.csv")
temp_raster <- file.path(tempdir(), "raster_cons_com_predictions.tif")

# Save predictions in data.frame objects
write.csv(suit_cons_predt, file = temp_df_file, row.names = FALSE)

# Save predictions in raster objects
terra::writeRaster(suit_cons_predt, filename = temp_raster)

# Import predictions as data.frame objects
read_com_pred_df <- read.csv(temp_df_file)

# Import predictions as SpatRaster objects
read_com_pred_ras <- terra::rast(temp_raster)
```