The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Introduction to gllvm Part 2: Species correlations

Jenni Niku

2023-09-18

Introduction to gllvm

R package gllvm

# From CRAN
install.packages(gllvm)
# OR
# From GitHub using devtools package's function install_github
devtools::install_github("JenniNiku/gllvm")
Problems?

gllvm package depends on R packages TMB and mvabund, try to install these first.

Distributions

Response Distribution Method Link
Counts Poisson VA/LA log
NB VA/LA log
ZIP VA/LA log
ZINB VA/LA log
Binary Bernoulli VA/LA probit
LA logit
Ordinal Ordinal VA probit
Normal Gaussian VA/LA identity
Positive continuous Gamma VA/LA log
Non-negative continuous Exponential VA/LA log
Biomass Tweedie LA/EVA log
Percent cover beta LA/EVA probit/logit

Data input

Main function of the gllvm package is gllvm(), which can be used to fit GLLVMs for multivariate data with the most important arguments listed in the following:

gllvm(y = NULL, X = NULL, TR = NULL, family, num.lv = 2, 
 formula = NULL, method = "VA", row.eff = FALSE, n.init=1, starting.val ="res", ...)
library(gllvm)

Example: Spiders

Data fitting

Fit GLLVM with environmental variables \(g(E(y_{ij})) = \beta_{0j} + \boldsymbol{x}_i'\boldsymbol{\beta}_{j} + \boldsymbol{u}_i'\boldsymbol{\theta}_j\) using gllvm:

library(gllvm)
data("spider")
fitx <- gllvm(y = spider$abund, X=spider$x, family = "negative.binomial", num.lv = 2)
fitx
## Call: 
## gllvm(y = spider$abund, X = spider$x, family = "negative.binomial", 
##     num.lv = 2)
## family: 
## [1] "negative.binomial"
## method: 
## [1] "VA"
## 
## log-likelihood:  -588.9186 
## Residual degrees of freedom:  217 
## AIC:  1415.837 
## AICc:  1548.06 
## BIC:  1870.074

Model selection

X=spider$x
fitx1 <- gllvm(spider$abund, X, family = "negative.binomial", num.lv = 1)
fitx2 <- gllvm(spider$abund, X, family = "negative.binomial", num.lv = 2)
fitx3 <- gllvm(spider$abund, X, family = "negative.binomial", num.lv = 3)
AIC(fitx1)
## [1] 1393.837
AIC(fitx2)
## [1] 1425.35
AIC(fitx3)
## [1] 1445.35

Residual analysis

par(mfrow = c(1,2))
plot(fitx1, which = 1:2)

Exercises

E1. Load spider data from mvabund package and take a look at the dataset.

library(gllvm)
#Package **mvabund** is loaded with **gllvm** so just load with a function `data()`.
data("spider")
# more info: 
# ?spider
Show the answers.

Package mvabund is loaded with gllvm so just load with a function data().

# response matrix:
spider$abund
##       Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi Pardlugu Pardmont
##  [1,]       25       10        0        0        0        4        0       60
##  [2,]        0        2        0        0        0       30        1        1
##  [3,]       15       20        2        2        0        9        1       29
##  [4,]        2        6        0        1        0       24        1        7
##  [5,]        1       20        0        2        0        9        1        2
##  [6,]        0        6        0        6        0        6        0       11
##  [7,]        2        7        0       12        0       16        1       30
##  [8,]        0       11        0        0        0        7       55        2
##  [9,]        1        1        0        0        0        0        0       26
## [10,]        3        0        1        0        0        0        0       22
## [11,]       15        1        2        0        0        1        0       95
## [12,]       16       13        0        0        0        0        0       96
## [13,]        3       43        1        2        0       18        1       24
## [14,]        0        2        0        1        0        4        3       14
## [15,]        0        0        0        0        0        0        6        0
## [16,]        0        3        0        0        0        0        6        0
## [17,]        0        0        0        0        0        0        2        0
## [18,]        0        1        0        0        0        0        5        0
## [19,]        0        1        0        0        0        0       12        0
## [20,]        0        2        0        0        0        0       13        0
## [21,]        0        1        0        0        0        0       16        1
## [22,]        7        0       16        0        4        0        0        2
## [23,]       17        0       15        0        7        0        2        6
## [24,]       11        0       20        0        5        0        0        3
## [25,]        9        1        9        0        0        2        1       11
## [26,]        3        0        6        0       18        0        0        0
## [27,]       29        0       11        0        4        0        0        1
## [28,]       15        0       14        0        1        0        0        6
##       Pardnigr Pardpull Trocterr Zoraspin
##  [1,]       12       45       57        4
##  [2,]       15       37       65        9
##  [3,]       18       45       66        1
##  [4,]       29       94       86       25
##  [5,]      135       76       91       17
##  [6,]       27       24       63       34
##  [7,]       89      105      118       16
##  [8,]        2        1       30        3
##  [9,]        1        1        2        0
## [10,]        0        0        1        0
## [11,]        0        1        4        0
## [12,]        1        8       13        0
## [13,]       53       72       97       22
## [14,]       15       72       94       32
## [15,]        0        0       25        3
## [16,]        2        0       28        4
## [17,]        0        0       23        2
## [18,]        0        0       25        0
## [19,]        1        0       22        3
## [20,]        0        0       22        2
## [21,]        0        1       18        2
## [22,]        0        0        1        0
## [23,]        0        0        1        0
## [24,]        0        0        0        0
## [25,]        6        0       16        6
## [26,]        0        0        1        0
## [27,]        0        0        0        0
## [28,]        0        0        2        0
# Environmental variables
spider$x
##    soil.dry bare.sand fallen.leaves   moss herb.layer reflection
## 1    2.3321    0.0000        0.0000 3.0445     4.4543     3.9120
## 2    3.0493    0.0000        1.7918 1.0986     4.5643     1.6094
## 3    2.5572    0.0000        0.0000 2.3979     4.6052     3.6889
## 4    2.6741    0.0000        0.0000 2.3979     4.6151     2.9957
## 5    3.0155    0.0000        0.0000 0.0000     4.6151     2.3026
## 6    3.3810    2.3979        3.4340 2.3979     3.4340     0.6931
## 7    3.1781    0.0000        0.0000 0.6931     4.6151     2.3026
## 8    2.6247    0.0000        4.2627 1.0986     3.4340     0.6931
## 9    2.4849    0.0000        0.0000 4.3307     3.2581     3.4012
## 10   2.1972    3.9318        0.0000 3.4340     3.0445     3.6889
## 11   2.2192    0.0000        0.0000 4.1109     3.7136     3.6889
## 12   2.2925    0.0000        0.0000 3.8286     4.0254     3.6889
## 13   3.5175    1.7918        1.7918 0.6931     4.5109     3.4012
## 14   3.0865    0.0000        0.0000 1.7918     4.5643     1.0986
## 15   3.2696    0.0000        4.3944 0.6931     3.0445     0.6931
## 16   3.0301    0.0000        4.6052 0.6931     0.6931     0.0000
## 17   3.3322    0.0000        4.4543 0.6931     3.0445     1.0986
## 18   3.1224    0.0000        4.3944 0.0000     3.0445     1.0986
## 19   2.9232    0.0000        4.5109 1.6094     1.6094     0.0000
## 20   3.1091    0.0000        4.5951 0.6931     0.6931     0.0000
## 21   2.9755    0.0000        4.5643 0.6931     1.7918     0.0000
## 22   1.2528    3.2581        0.0000 4.3307     0.6931     3.9120
## 23   1.1939    3.0445        0.0000 4.0254     3.2581     4.0943
## 24   1.6487    3.2581        0.0000 4.0254     3.0445     4.0073
## 25   1.8245    3.5835        0.0000 1.0986     4.1109     2.3026
## 26   0.9933    4.5109        0.0000 1.7918     1.7918     4.3820
## 27   0.9555    2.3979        0.0000 3.8286     3.4340     3.6889
## 28   0.9555    3.4340        0.0000 3.7136     3.4340     3.6889
# Plot data using boxplot:
boxplot(spider$abund)

E2. Fit GLLVM with two latent variables to spider data with a suitable distribution. Data consists of counts of spider species.

# Take a look at the function documentation for help: 
?gllvm
Show the answers.

2. Response variables in spider data are counts, so Poisson, negative binomial and zero inflated Poisson are possible. However, ZIP is implemented only with Laplace method, so it need to be noticed, that if models are fitted with different methods they can not be compared with information criteria. Let’s try just with a Poisson and NB. NOTE THAT the results may not be exactly the same as below, as the initial values for each model fit are slightly different, so the results may

# Fit a GLLVM to data
fitp <- gllvm(y=spider$abund, family = poisson(), num.lv = 2)
fitp
## Call: 
## gllvm(y = spider$abund, family = poisson(), num.lv = 2)
## family: 
## [1] "poisson"
## method: 
## [1] "VA"
## 
## log-likelihood:  -845.8277 
## Residual degrees of freedom:  301 
## AIC:  1761.655 
## AICc:  1770.055 
## BIC:  1895.254
fitnb <- gllvm(y=spider$abund, family = "negative.binomial", num.lv = 2)
fitnb
## Call: 
## gllvm(y = spider$abund, family = "negative.binomial", num.lv = 2)
## family: 
## [1] "negative.binomial"
## method: 
## [1] "VA"
## 
## log-likelihood:  -733.6806 
## Residual degrees of freedom:  289 
## AIC:  1561.361 
## AICc:  1577.028 
## BIC:  1740.765

Based on AIC, NB distribution suits better. How about residual analysis: NOTE THAT The package uses randomized quantile residuals so each time you plot the residuals, they look a little different.

# Fit a GLLVM to data
par(mfrow = c(1,2))
plot(fitp, which = 1:2)

plot(fitnb, which = 1:2)

You could do these comparisons with Laplace method as well, using the code below, and it would give the same conclusion that NB distribution suits best:

fitLAp <- gllvm(y=spider$abund, family = poisson(), method = "LA", num.lv = 2)
fitLAnb <- gllvm(y=spider$abund, family = "negative.binomial", method = "LA", num.lv = 2)
fitLAzip <- gllvm(y=spider$abund, family = "ZIP", method = "LA", num.lv = 2)
AIC(fitLAp)
AIC(fitLAnb)
AIC(fitLAzip)
GLLVM with two latent variables can be used as a model-based approach to unconstrained ordination, as considered at the first day of the workshop.

E3. Fit GLLVM with environmental variables soil.dry and reflection to the data with suitable number of latent variables.

Show the answers.

We can extract the two columns from the environmental variable matrix or define the model using formula.

# `soil.dry` and `reflection` are in columns 1 and 6
X <- spider$x[,c(1,6)]
fitx1 <- gllvm(spider$abund, X, family = "negative.binomial", num.lv = 1)
fitx2 <- gllvm(spider$abund, X, family = "negative.binomial", num.lv = 2)
fitx3 <- gllvm(spider$abund, X, family = "negative.binomial", num.lv = 3)
AIC(fitx1)
## [1] 1449.428
AIC(fitx2)
## [1] 1471.428
AIC(fitx3)
## [1] 1506.451
# Or alternatively using formula:
fitx1 <- gllvm(spider$abund, spider$x, formula = ~soil.dry + reflection, family = "negative.binomial", num.lv = 1)
fitx1
## Call: 
## gllvm(y = spider$abund, X = spider$x, formula = ~soil.dry + reflection, 
##     family = "negative.binomial", num.lv = 1)
## family: 
## [1] "negative.binomial"
## method: 
## [1] "VA"
## 
## log-likelihood:  -664.7138 
## Residual degrees of freedom:  276 
## AIC:  1449.428 
## AICc:  1476.046 
## BIC:  1678.454
Model with one latent variable gave the lowest AIC value.

E4. Explore the model fit. Find the coefficients for environmental covariates.

Show the answers.

Estimated parameters can be obtained with coef() function. Confidence intervals for parameters are obtained with confint().

coef(fitx1)
## $Species.scores
##                   LV1
## Alopacce 1.000000e+00
## Alopcune 1.286715e+00
## Alopfabr 4.410752e-01
## Arctlute 2.261569e+00
## Arctperi 3.236902e-07
## Auloalbi 2.499151e+00
## Pardlugu 3.974207e-01
## Pardmont 9.559481e-01
## Pardnigr 2.566314e+00
## Pardpull 2.528325e+00
## Trocterr 9.659148e-01
## Zoraspin 1.751412e+00
## 
## $sigma.lv
##       LV1 
## 0.8536381 
## 
## $Intercept
##   Alopacce   Alopcune   Alopfabr   Arctlute   Arctperi   Auloalbi   Pardlugu 
##  -3.009461  -7.407015   2.078134 -21.475462  -7.919460 -12.966531   4.567866 
##   Pardmont   Pardnigr   Pardpull   Trocterr   Zoraspin 
##  -6.091019 -14.644236 -16.455489  -3.711573  -9.215034 
## 
## $Xcoef
##            soil.dry reflection
## Alopacce -0.3441060  1.5943262
## Alopcune  2.6111104  0.6246532
## Alopfabr -1.8669681  0.7467644
## Arctlute  6.3215994  0.8205912
## Arctperi -1.7567138  2.8890762
## Auloalbi  4.1813972  0.6831463
## Pardlugu -0.5453219 -1.2605150
## Pardmont  1.7924537  1.4691347
## Pardnigr  4.9041830  0.8601458
## Pardpull  5.3496432  1.2822836
## Trocterr  2.3357726  0.2234482
## Zoraspin  3.5554592  0.1198072
## 
## $inv.phi
##     Alopacce     Alopcune     Alopfabr     Arctlute     Arctperi     Auloalbi 
## 9.379320e+00 1.699533e+00 1.577633e+00 1.064240e+00 1.629747e+07 1.511543e+00 
##     Pardlugu     Pardmont     Pardnigr     Pardpull     Trocterr     Zoraspin 
## 1.543692e+00 9.784559e-01 2.517855e+00 2.339480e+00 1.471002e+01 3.625230e+00 
## 
## $phi
##     Alopacce     Alopcune     Alopfabr     Arctlute     Arctperi     Auloalbi 
## 1.066175e-01 5.883970e-01 6.338608e-01 9.396375e-01 6.135922e-08 6.615756e-01 
##     Pardlugu     Pardmont     Pardnigr     Pardpull     Trocterr     Zoraspin 
## 6.477976e-01 1.022018e+00 3.971634e-01 4.274454e-01 6.798089e-02 2.758446e-01
# Coefficients for covariates are named as `Xcoef`
# Confidence intervals for these coefficients:
confint(fitx1, parm = "Xcoef")
##               cilow       ciup
## Xcoef1  -0.89413837  0.2059263
## Xcoef2   1.60653994  3.6156809
## Xcoef3  -2.83503803 -0.8988983
## Xcoef4   5.72879548  6.9144033
## Xcoef5  -2.87611514 -0.6373125
## Xcoef6   2.77826475  5.5845296
## Xcoef7  -1.62666107  0.5360173
## Xcoef8   0.82543215  2.7594753
## Xcoef9   3.71814685  6.0902191
## Xcoef10  4.08453542  6.6147511
## Xcoef11  1.89268036  2.7788647
## Xcoef12  2.47215622  4.6387621
## Xcoef13  0.99262088  2.1960315
## Xcoef14  0.22972302  1.0195833
## Xcoef15 -0.23417207  1.7277008
## Xcoef16 -0.03218841  1.6733709
## Xcoef17  1.63801838  4.1401341
## Xcoef18  0.04128991  1.3250026
## Xcoef19 -1.72657436 -0.7944556
## Xcoef20  0.93054958  2.0077199
## Xcoef21  0.24237359  1.4779180
## Xcoef22  0.64776105  1.9168063
## Xcoef23 -0.00557523  0.4524716
## Xcoef24 -0.32437013  0.5639846
# The first 12 intervals are for soil.dry and next 12 for reflection
Problems? See hints:

I have problems in model fitting. My model converges to infinity or local maxima: GLLVMs are complex models where starting values have a big role. Choosing a different starting value method (see argument starting.val) or use multiple runs and pick up the one giving highest log-likelihood value using argument n.init. More variation to the starting points can be added with jitter.var.

My results does not look the same as in answers: The results may not be exactly the same as in the answers, as the initial values for each model fit are slightly different, so the results may also differ slightly.

Studying species correlations

Species correlations

fitnb <- gllvm(spider$abund, family = "negative.binomial", num.lv = 2)

Visualizing species correlations

fitnb <- gllvm(spider$abund, family = "negative.binomial", num.lv = 2)
ordiplot(fitnb, biplot = TRUE)
abline(h = 0, v = 0, lty=2)

fitnb <- gllvm(spider$abund, family = "negative.binomial", num.lv = 2)
cr <- getResidualCor(fitnb)
library(corrplot);
corrplot(cr, diag = FALSE, type = "lower", method = "square", tl.srt = 25)

ordiplot(fitnb, biplot = TRUE)
abline(h = 0, v = 0, lty=2)

Studying effects of covariates

Studying effects of environmental variables

rbPal <- c("#00FA9A", "#00EC9F", "#00DFA4", "#00D2A9", "#00C5AF", "#00B8B4", "#00ABB9", "#009DBF", "#0090C4", "#0083C9", "#0076CF", "#0069D4", "#005CD9", "#004EDF", "#0041E4", "#0034E9", "#0027EF", "#001AF4", "#000DF9", "#0000FF")
X <- spider$x[,c(1,6)]
par(mfrow = c(1,2), mar=c(4,4,2,2))
for(i in 1:ncol(X)){
Col <- rbPal[as.numeric(cut(X[,i], breaks = 20))]
ordiplot(fitnb, symbols = T, s.colors = Col, main = colnames(X)[i], biplot = TRUE)
abline(h=0,v=0, lty=2)

}

Coefficient plot

fitx1 <- gllvm(spider$abund, X, family = "negative.binomial", num.lv = 1)
coefplot(fitx1, mfrow = c(1,2), cex.ylab = 0.8)

Correlation matrix

crx <- getResidualCor(fitx1)
corrplot(crx, diag = FALSE, type = "lower", method = "square", tl.srt = 25)

Fourth corner models


  1. Niku, J., F.K.C. Hui, S. Taskinen, and D.I. Warton. 2019. Gllvm - Fast Analysis of Multivariate Abundance Data with Generalized Linear Latent Variable Models in R. 10. Methods in Ecology and Evolution: 2173–82↩︎

  2. Kasper Kristensen, Anders Nielsen, Casper W. Berg, Hans Skaug, Bradley M. Bell (2016). TMB: Automatic Differentiation and Laplace Approximation. Journal of Statistical Software, 70(5), 1-21↩︎

  3. Hui, F. K. C., Warton, D., Ormerod, J., Haapaniemi, V., and Taskinen, S. (2017). Variational approximations for generalized linear latent variable models. Journal of Computational and Graphical Statistics. Journal of Computational and Graphical Statistics, 26:35-43↩︎

  4. Korhonen, P., Hui, F. K. C., Niku, J., and Taskinen, S. (2021). Fast, universal estimation of latent variable models using extended variational approximations, arXiv:2107.02627 .↩︎

  5. Niku, J., Warton, D. I., Hui, F. K. C., and Taskinen, S. (2017). Generalized linear latent variable models for multivariate count and biomass data in ecology. Journal of Agricultural, Biological, and Environmental Statistics, 22:498-522.↩︎

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.