Skip to content

minor improvements of vignette #40

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
BertvanderVeen wants to merge 1 commit into
base: master
Choose a base branch
from
Open

minor improvements of vignette #40

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 12 additions & 11 deletions vignettes/SAM-example.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,21 +20,21 @@ knitr::opts_chunk$set(
```

### Species Archetype Models
The Species Archetype Models (SAMs) are variant of 'Mixture-of-Regressions' and quantify how groups of **species** vary with the environment. These groups are referred to as Species Archetypes [@dunstan_model_2011]. Each Archetype represents a group of species which jointly response to environmental data in the model (covariates). For example, one group of species might like warm water, while another might like cold. Typically, we model observations of species at sites. The covariates are used to describe the variation (and response) of each Archetype to the physical and environmental gradients. We can describe the data as $i = 1...n$ sampling sites, $j = 1...S$ species (or desired taxonomic unit) and $k = 1...K$ Archetypes. The model conditional mean observations (occurrence/count/biomass) of species is, $\mathbb{E}(y_{ij}|\textrm{z}_{k})$, on $g_{k}(X_i)$ Archetype covariates. In `ecomix`, the intercepts are species-specific, this is an update from the original 'SpeciesMix' R package [@dunstan_speciesmix_2013], where the intercepts were only species-specific for the Negative Binomial and Tweedie distributions, as described in @dunstan_finite_2013. For error distributions with dispersion or variance parameters (Negative Binomial, Tweedie and Gaussian) these parameters are also species-specific. The model can be described as follows:
The Species Archetype Models (SAMs) are variant of 'Mixture-of-Regressions' and quantify how groups of **species** vary with the environment. These groups are referred to as Species Archetypes [@dunstan_model_2011]. Each Archetype represents a group of species which jointly response to environmental data in the model (covariates). For example, one group of species might like warm water, while another might like cold. Typically, we model observations of species at sites. The covariates are used to describe the variation (and response) of each Archetype to the physical and environmental gradients. We can describe the data as $i = 1...n$ sampling sites, $j = 1...S$ species (or desired taxonomic unit) and $k = 1...K$ Archetypes. The model for the conditional mean observations (occurrence/count/biomass) of species is, $\mathbb{E}(y_{ij}|\textrm{z}_{k})$, on $g_{k}(X_i)$ Archetype covariates. In `ecomix`, the intercepts are species-specific, this is an update from the original 'SpeciesMix' R package [@dunstan_speciesmix_2013], where the intercepts were only species-specific for the Negative Binomial and Tweedie distributions, as described in @dunstan_finite_2013. For error distributions with dispersion or variance parameters (Negative Binomial, Tweedie and Gaussian) these parameters are also species-specific. The model can be described as follows:

$$
\begin{equation}
h[\mathbb{E}(y_{ij}|z_{k})] = \alpha_j + g_k(X_{i}^\top\beta_{k}) + \nu_i \tag{1}\label{eq:one}
h[\mathbb{E}(y_{ij}|z_{k})] = \alpha_j + g_k(X_{i})^\top\beta_{k} + \nu_i \tag{1}\label{eq:one}
\end{equation}
$$

where $Pr(z_{k}) = \pi_k,$ and $\sum^K_{k=1}{\pi_k=1}$. The functional form of $g_k(.)$ can be specified to be any function commonly used within a Generalized Linear Model framework. Including linear, quadratic, spline and interaction terms. Additionally an offset term $\nu_i$ can be included to account for sampling artefacts and are included into the model on a log-scale (e.g. log(area sampled)). We refer to the model as the `species_mix` in the `ecomix` package.
where $Pr(z_{k}) = \pi_k,$ and $\sum^K_{k=1}{\pi_k=1}$. The functional form of $g_k(.)$ can be specified to be any function commonly used within a Generalized Linear Model framework. Including linear, quadratic, spline and interaction terms. Additionally an offset term $\nu_i$ can be included to account for sampling artefacts which is included into the model on a log-scale (e.g. log(area sampled)). We refer to the model as the `species_mix` in the `ecomix` package.

```{r setup}
library(ecomix)
```

Here we will demonstrate how to fit and interpret a Species Archetype Models (SAMs) from our `ecomix` package. We will present a simulation study to demonstrate the functionality of `species_mix`.
Here we will demonstrate how to fit and interpret Species Archetype Models (SAMs) from our `ecomix` package. We will present a simulation study to demonstrate the functionality of `species_mix`.

#### Simulate Environmental Data
We generate a set of simulated environment predictors using Gaussian random fields with nugget effects on some of the variables.
Expand Down Expand Up @@ -129,7 +129,7 @@ plot_covariates(env.df)
```

#### Simulate Biological Data
We simulated a set of synthetic species to be fitted using the `species_mix` function. We generated the expected species intercepts $\alpha_j$ from a beta distribution, and assign known group level covariates $\beta_k$. $\beta_k$ represents the archetype (group) response to each covariate in the model. We simulated species archetypical responses using the `species_mix.simulate` function. If no known parameters are provided random parameters will be generated for the formula and data provided. Here we provided parameters for the species intercepts (alphas) and the archetype mean responses (betas).
We simulated a set of synthetic species to be fitted using the `species_mix` function. We generate the species intercepts $\alpha_j$ from a beta distribution, and assign pre-defined group level covariate effects $\beta_k$. $\beta_k$ represents the archetype (group) response to each covariate in the model. We then simulate species archetypical responses using the `species_mix.simulate` function. If no known parameters are provided, random parameters will be generated for the formula and data provided. Here we provide parameters for the species intercepts (alphas) and the archetype mean responses (betas).

```{r, eval=TRUE, echo=TRUE}
set.seed(42)
Expand Down Expand Up @@ -171,9 +171,9 @@ simulated_data200 <- species_mix.simulate(archetype_formula=sam_form,


#### Model fitting and evaluation
We have generated occurrence records for 100 synthetic species across the 200 randomly surveyed sites. Because we have presence and absence data we can fit a SAM with a Bernoulli family. In this example, we may chose to remove the rare species (< 10 occurrences across all sites). These species could be potentially included in the model, but would likely create noise and unexplained variance, making the model harder to fit and estimate [@hui_mix_2013].
We have generated occurrence records for 100 synthetic species across the 200 randomly surveyed sites. Because we have presence and absence data we fit a SAM with a Bernoulli family. In this example, we may chose to remove the rare species (< 10 occurrences across all sites). These species could be potentially included in the model, but would likely create noise and unexplained variance, making the model harder to fit and estimate [@hui_mix_2013].

```{r, echo=FALSE, fig.cap='Figure 2. The simulated occurrences for the 100 species. We will remove all species with less than 10 presences across all 200 sites. '}
```{r, echo=FALSE, fig.cap='Figure 2. The simulated occurrences for the 100 species. We will remove all species with less than 10 presences across all 200 sites. ', dpi = 150}
count <- table(colSums(simulated_data200[,1:100]))
occur <- as.numeric(names(count))
mat <- as.data.frame(cbind(occur,count))
Expand All @@ -185,7 +185,7 @@ bpdf <- cbind(bp,df1)
axis(1, at=bp[c(0,25,50,75,100,nrow(bpdf)-1)+1,1], labels = c(0,25,50,75,100,""),las=1)
```

In this, example we fit independent polynomials with two degrees of freedom for the four simulated covariates Temperature, Oxygen, Depth & Productivity. We do this to demonstrate a simple simulated example where can show the known response of covariates to simulated data. We include time as a factor, which represents the time in which the simulated samples were recorded, for these data the two factors are 'Day' or 'Night'. But additional discrete times could be included such as years. Below is a small code block with a basic example on how to fit a single species_mix model.
In this example we fit a model with independent polynomials with two degrees of freedom for the four simulated covariates Temperature, Oxygen, Depth & Productivity. We do this to demonstrate a simple simulated example where can show the known response of covariates to simulated data. We include time ("Day" or "Night") as a factor, which represents when the simulated samples were recorded. Additional temporal covariates could be included, such as year. Below is a small code block with a basic example on how to fit a single `species_mix` model.

```{r, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE}
## load the ecomix package
Expand Down Expand Up @@ -213,9 +213,9 @@ sam_fit <- species_mix(archetype_formula = archetype_formula, # Archetype formul
```

#### Group selection
One challenge when developing SAMs (or any finite mixture model) is selecting $k$; the number of archetypes (groups) in the model. The number of archetypes is latent, so must be estimated from the data and the functional form of the covariates. In this example, we know that the optimal number of groups for these data is three, because we simulated the data with these characteristics. However, $k$ is generally not known in real world applications. So one import part of the fitting process is finding $k$, having said that it is also totally fine to define $k$ for say a management objective. We can estimate the "best" number of groups based on the most parsimonious fit to the data. We can do this based on the model log-likelihood, and information criterion such as BIC. We provide a function `species_mix.multifit` which can assist in group selection if a vector of archetypes is provided. Below is an example of how one might do this. Note that with only a single random start per archetype count (`nstart=1`), some fits -- especially for larger $k$ -- can land on a degenerate solution where an archetype collapses to (near) zero species; `species_mix.multifit`/`plot.species_mix.multifit` detect and drop these, which can leave gaps in the BIC curve. Using a handful of random starts per archetype count (`nstart=3` here) makes it far more likely that at least one converges to a sensible solution; for real applications, increase `nstart` further still.
One challenge when fitting SAMs (or any finite mixture model) is selecting $k$; the number of archetypes (groups) in the model. The number of archetypes is latent, so must be estimated from the data and the functional form of the covariates. In this example, we know that the optimal number of groups for these data is three, because we simulated the data with these characteristics. However, $k$ is generally not known in real world applications. So one import part of the fitting process is finding $k$. Having said that, $k$ may also be defined a-priori for, say, a management objective. We can estimate the "best" number of groups based on the most parsimonious fit to the data. We can do this based on the model log-likelihood, and information criterion such as BIC. We provide a function `species_mix.multifit` which can assist in group selection if a vector of archetypes is provided. Below is an example of how one might do this. Note that with only a single random start per archetype count (`nstart=1`), some fits -- especially for larger $k$ -- can land on a degenerate solution where an archetype collapses to (near) zero species; `species_mix.multifit`/`plot.species_mix.multifit` detect and drop these, which can leave gaps in the BIC curve. Using a handful of random starts per archetype count (`nstart=3` here) makes it far more likely that at least one converges to a sensible solution; for real applications, increase `nstart` further still.

```{r, eval = TRUE, fig.width=5,fig.height=4,fig.cap="Figure 3. Group selection from multiple fit function, we can see that three archetypes is the best fit to these data based on Bayesian Information Criterion (BIC).", cache=TRUE}
```{r, eval = TRUE, fig.width=5,fig.height=4,fig.cap="Figure 3. Group selection from multiple fit function, we can see that three archetypes is the best fit to these data based on the Bayesian Information Criterion (BIC).", cache=TRUE}
nArchetypes <- 1:6
sam_multifit <- species_mix.multifit(archetype_formula = archetype_formula, # Archetype formula
species_formula = species_formula, # Species formula
Expand Down Expand Up @@ -252,7 +252,8 @@ plot(x = eff.df, object = sam_fit, boot.object=sam_boot,ylim=c(0,1))
```

#### Model prediction
We can generate predictions for the each archetype in ecomix using `predict()` function. The default is used to generate a point mean estimate for each archetype. With the inclusion of a `bootstrap` object we can also provide estimates of uncertainty or standard error for each prediction. We can see spatial predict the point mean and standard error in the predicted distributions of each archetype (Fig. 7). We can also see the spatial response of each archetype is strongly correlated to it's environmental response to each covariate (Fig. 6). For example, archetype two has a strong response to depth (Fig. 7), which is evident in the archetype responses (Fig. 3d).
We can generate predictions for the each archetype in ecomix using the `predict()` function. The default is used to generate a point mean estimate for each archetype. With the inclusion of a `bootstrap` object we can also provide estimates of uncertainty or standard error for each prediction. The response of each archetype is strongly related to its environmental response to each covariate (Fig. 6). For example, archetype two has a strong positive response to depth (Fig. 6). We can also examine the spatial predictions of point mean estimates standard error in the predicted distributions of each archetype visually (Fig. 7).

```{r predict}
env.df$Time <- factor("Day",levels=c("Day","Night"))
sam3_pred <- predict(sam_fit, sam_boot, newdata=env.df)
Expand Down