---
title: 'Analysing sparse ecological percent cover data using gllvm'
author: "Pekka Korhonen"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Analysing sparse ecological percent cover data using gllvm}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, size="footnotesize", fig.width=7, fig.height=7, 
fig.align="center",dev="png", code.frame = TRUE, warning = FALSE, fig.pos='H')
```

This short document showcases how to use the **gllvm** package to analyse multivariate percent cover data. Namely, we show how to apply the hurdle beta GLLVM (with logistic link), as detailed in @Korhonen2024, to analyse the kelp forest dataset from the Santa Barbara Coastal Long-Term Ecological Research project, available from [https://doi.org/10.6073/pasta/0af1a5b0d9dde5b4e5915c0012ccf99c](https://doi.org/10.6073/pasta/0af1a5b0d9dde5b4e5915c0012ccf99c).

A multivariate percent cover dataset often comes in the form of a $n \times m$ matrix $\bf{Y}$, with $n$ being the number of observational units/sites, and $m$ being the number of species of plants, macroalgae, sessile invertebrates, et cetera. The response $y_{ij}$ is then the recorded percentage of the relative area covered by species $j$ on unit/site $i$. Typically such datasets contain considerable proportion of zero observations, as a given obs. unit is often populated by only a small subset of the $m$ species in total. More rarely, it might even be the case that one of the species covers some obs. unit completely.

## Hurdle beta GLLVM

Traditionally, beta regression has been used to model responses that take the form of percentages. If $y^*_{ij}$ is in the (open) interval $(0,1)$, and is distributed according to beta distribution, $y^*_{ij} \sim \text{Beta}(\mu_{ij}, \phi_j)$, then it has the following density:
$$ f_{\text{beta}}(y_{ij}^*; \mu_{ij}, \phi_j) = \frac{\Gamma(\phi_j)}{\Gamma(\mu_{ij}\phi_j)\Gamma(\phi_j-\phi_j\mu_{ij})} (y_{ij}^*)^{\mu_{ij}\phi_j-1}(1-y_{ij}^*)^{\phi_j(1-\mu_{ij}-1)}.$$

The mean parameter $\mu_{ij}$ can be connected to a set of covariates and latent variables through some link function, by the equation $g(\mu_{ij})=\eta_{ij} = \beta_{0j} + \boldsymbol{\beta}_j^\top\bf{x_i} + \boldsymbol{\lambda}_j^\top\bf{u_i}$.

However, such a model is ill-suited for most percent cover datasets, due to the fact that the beta distribution isn't capable of handling zero (or $100\%$) responses. To make use of the beta regression model in such scenario, one needs to use some transformation in order to 'push' the responses away from the boundaries. This procedure might provide reasonable results when the numbers of zeros or ones in the data are relatively low. On the other hand, by transforming the zeros and ones, we might lose some important information.

A more sophisticated way of tackling such issue is by considering a zero-accommodating model, couple of which have been proposed recently. One such model is the hurdle beta model, which models the two (can be extended to three if $100\%$ covers are recorded, see @Korhonen2024) classes the response $y_{ij}$ can take, i.e., $\{0\}$ and $(0,1)$, by separate processes. Namely, the zeros are assumed to be generated by a Bernoulli process. Conditional on $y_{ij}\in(0,1)$, the response is modeled using the standard beta distribution as presented above. The likelihood function for a response $Y_{ij}$ following the hurdle beta distribution is of the form:
$$
P(Y_{ij};\mu_{ij}, \mu_{ij}^0, \phi_j) = \begin{cases}
        1-\mu_{ij}^0, & Y_{ij} = 0,\\
         \mu_{ij}^0 \cdot f_{\text{beta}}(Y_{ij};\mu_{ij},\phi_j), & Y_{ij} \in (0,1). \\
        \end{cases}
$$
where $g(\mu_{ij}^0) = \eta_{ij}^0 = \beta_{0j}^0 + \bf{x_i}^\top\boldsymbol{\beta}_j^0 + \bf{u_i}^\top\boldsymbol{\lambda}_j^0$ for and $g(\mu_{ij}) = \eta_{ij} = \beta_{0j} + \bf{x_i}^\top\boldsymbol{\beta}_j + \bf{u_i}^\top\boldsymbol{\lambda}_j$ and $g(\cdot)$ can be either probit or logistic link function. Note, that here, the separate linear predictors share the same environmental covariates $\bf{x}_i$ and latent variable scores $\bf{u}_i$, while the coefficients and loadings are allowed to differ.

The **gllvm** package implements the hurdle beta GLLVM with two different estimation methods available. First, accessed by the argument ``method="VA"`` when calling \texttt{gllvm()} uses a hybrid approach, where the method of variational approximations, or VA, is applied to the Bernoulli-process part of the data (only probit link allowed), while the method of **extended** variational approximations, see @Korhonen2023, is applied to the beta distributed part. By instead specifying ``method="EVA"``, the **EVA** method is applied to both parts of the likelihood.

## Example

In the following, we show how to fit the logistic hurdle beta GLLVM using **EVA** on the SBC LTER marine macroalgae (i.e., seaweed) percent cover dataset. The data has been collected in 2000-2020 along 44 permanent transect lines along coastal southern California. We will specify a model with two latent variables (for ordination) and will include the rockiness of the seabed and the average number of stripes of giant kelp as environmental covariates in the model.

```{r}
library(gllvm)
data("kelpforest")
Yabund <- kelpforest$Y
Xenv <- kelpforest$X
SPinfo <- kelpforest$SPinfo

# Data contains both algae and sessile invertebrates
table(SPinfo$GROUP)

# Select only the macroalgae:
Yalg <- Yabund[,SPinfo$GROUP=="ALGAE"]

# Remove empty rows and rows with missing values from the data
rem0 <- which((rowSums(Yalg)==0) | is.na(rowSums(Yalg)))
Yalg <- Yalg[-rem0,]
Xenv <- Xenv[-rem0,]

# Use only the data from the year 2016:

Yalg <- Yalg[Xenv$YEAR==2016,]
Xenv <- Xenv[Xenv$YEAR==2016,]

# Remove species which have no observations
Yalg <- Yalg[,-which(colSums(Yalg>0)==0)]
# Number of obs. and species:
dim(Yalg)

# Specify the covariates in the linear predictor
Xformulai = ~KELP_FRONDS + PERCENT_ROCKY
```

```{r, echo=FALSE}
rownames(Yalg) <- interaction(Xenv$SITE, Xenv$TRANSECT)
```

After setting up the data, LV design and the covariates, the model is estimated by
```{r}
fit <- gllvm(Yalg, X=Xenv, formula = Xformulai, family = "betaH", method="EVA", 
             num.lv = 2, link="logit", control=list(reltol=1e-12))
```

To inspect e.g., the covariate effects, use
```{r}
fit$params$Xcoef
```
In the above, the prefix \texttt{H01} indicates that the coefficient relates to the Bernoulli part of the hurdle model.

Ordination plot can then be generated as per usual:
```{r}
ordiplot(fit, jitter = TRUE, s.cex = .8)
```


## References

---
references:
- id: Korhonen2023
  title: Fast and universal estimation of latent variable models using extended variational approximations.
  author: 
  - family: Korhonen
    given: P.
  - family: Hui
    given: F.K.C.
  - family: Niku
    given: J.
  - family: Taskinen
    given: S.
  publisher: Statistics and Computing
  volume: 33
  page: 26
  type: article-journal
  issued: 
    year: 2023
  doi: 10.1007/s11222-022-10189-w
  
- id: Korhonen2024
  title: A Comparison of Joint Species Distribution Models for Percent Cover Data.
  author: 
  - family: Korhonen
    given: P.
  - family: Hui
    given: F.K.C.
  - family: Niku
    given: J.
  - family: Taskinen
    given: S.
  - family: van der Veen
    given: B.
  publisher:
  volume:
  page:
  type: article-journal
  issued: 
    year: 2024
  doi: 10.48550/arXiv.2403.11562

---






















