Introduction to gllvm Part 1: Ordination
Jenni Niku
2023-09-19
Source:vignettes/vignette3.rmd
vignette3.rmdIntroduction to gllvm
R package gllvm
R package gllvm fits Generalized linear latent variable models (GLLVM) for multivariate data1.
Developed by J. Niku, W.Brooks, R. Herliansyah, F.K.C. Hui, S. Taskinen, D.I. Warton, B. van der Veen.
-
The package available in
Package installation:
# 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.
GLLVMs are computationally intensive to fit due the integral in log-likelihood.
gllvm package overcomes computational problems by applying closed form approximations to log-likelihood and using automatic differentiation in C++ to accelerate computation times (TMB2).
Estimation is performed using either variational approximation (VA3), extended variational approximation method (EVA4) or Laplace approximation (LA5) method implemented via R package TMB.
VA method is faster and more accurate than LA, but not applicable for all distributions and link functions.
-
Using gllvm we can fit
- GLLVM without covariates gives model-based ordination and biplots
- GLLVM with environmental covariates for studying factors explaining species abundance
- Fourth corner models with latent variables for studying environmental-trait interactions
- GLLVM without latent variables fits basic multivariate GLMs
Additional tools: model checking, model selection, inference, visualization.
Distributions
| Response | Distribution | Method | Link |
|---|---|---|---|
| Counts | Poisson | VA/LA | log |
| NB | VA/LA | log | |
| ZIP | LA | log | |
| Binary | Bernoulli | VA/LA | probit |
| EVA/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", ...)- y: matrix of abundances
- X: matrix or data.frame of environmental variables
- TR: matrix or data.frame of trait variables
- family: distribution for responses
- num.lv: number of latent variables
- method: approximation used “VA” or “LA”
- row.eff: type of row effects
- n.init: number of random starting points for latent variables
- starting.val: starting value method
Example: Spiders
- Abundances of 12 hunting spider species measured as a count at 28 sites6.
- Six environmental variables measured at each site.
-
soil.dry: Soil dry mass -
bare.sand: cover of bare sand -
fallen.leaves: cover of fallen leaves/twigs -
moss: cover of moss -
herb.layer: cover of herb layer -
reflection: reflection of the soil surface with a cloudless sky
-
Data fitting
Fitting basic GLLVM \(g(E(y_{ij})) = \beta_{0j} + \boldsymbol{u}_i'\boldsymbol{\theta}_j\) with gllvm:
library(mvabund)
data("spider")
library(gllvm)
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.765Residual analysis
Residual analysis can be used to assess the appropriateness of the fitted model (eg. in terms of mean-variance relationship).
Randomized quantile/Dunn-Smyth residuals7 are used in the package, as they provide standard normal distributed residuals, even for discrete responses, in the case of a proper model.
Exercises
Try to do these exercises for the next 10 minutes, as many as time is enough for.
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:
# ?spiderShow the answers.
1. Print the data and covariates and draw a boxplot of the 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 to spider data with a suitable distribution. Data consists of counts of spider species.
# Take a look at the function documentation for help:
?gllvmShow 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 also differ slightly.
# 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.765Based 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
plot(fitp)
plot(fitnb)
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:
E3. Explore the fitted model. Where are the estimates for parameters? What about predicted latent variables? Standard errors?
Show the answers.
3. Lets explore the fitted model:
# Parameters:
coef(fitnb)
## $Species.scores
## LV1 LV2
## Alopacce 1.0000000 0.00000000
## Alopcune -0.2631399 1.00000000
## Alopfabr 0.9015777 -0.48406268
## Arctlute -0.7036398 2.03035434
## Arctperi 0.5939498 -1.58944410
## Auloalbi -0.6526969 1.53311872
## Pardlugu -0.8886192 -0.02803854
## Pardmont 0.7498069 0.58843699
## Pardnigr -0.5748737 1.78025729
## Pardpull -0.4146911 1.94808390
## Trocterr -0.4694572 0.78799274
## Zoraspin -0.6964251 1.15697052
##
## $sigma.lv
## LV1 LV2
## 1.803174 1.829265
##
## $Intercept
## Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi Pardlugu
## 0.7740337 0.5514772 -0.1420830 -3.5927662 -3.3659862 -0.7311799 0.3187352
## Pardmont Pardnigr Pardpull Trocterr Zoraspin
## 1.6247762 -0.2756119 -0.1351194 2.6556796 0.2459100
##
## $inv.phi
## Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi Pardlugu
## 1.3061971 1.4686553 0.7051402 0.8975848 2.3318838 1.5391773 1.1528673
## Pardmont Pardnigr Pardpull Trocterr Zoraspin
## 1.6178940 2.5410879 2.1962396 12.9778835 2.5728375
##
## $phi
## Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi Pardlugu
## 0.76558124 0.68089495 1.41815778 1.11410081 0.42883784 0.64969775 0.86740248
## Pardmont Pardnigr Pardpull Trocterr Zoraspin
## 0.61808745 0.39353224 0.45532374 0.07705417 0.38867593
# Where are the predicted latent variable values? just fitp$lvs or
getLV(fitnb)
## LV1 LV2
## Row1 0.92852820 1.270518399
## Row2 -0.68267483 0.793812287
## Row3 0.67620039 1.274447740
## Row4 -0.24715948 1.149854311
## Row5 -0.36104903 1.217237043
## Row6 -0.35232477 1.008692032
## Row7 0.02923371 1.426060354
## Row8 -1.52132300 -0.070700584
## Row9 0.91074458 0.001240765
## Row10 1.03735488 -0.493620603
## Row11 1.48873297 0.303011200
## Row12 1.32413991 0.807370007
## Row13 0.11539506 1.356631020
## Row14 -0.43682462 1.014590269
## Row15 -1.31924491 -0.544659051
## Row16 -1.17467393 -0.242685488
## Row17 -1.11253204 -0.512783657
## Row18 -1.15894416 -0.542942463
## Row19 -1.34067365 -0.502769859
## Row20 -1.36734900 -0.592073248
## Row21 -1.13969788 -0.475614140
## Row22 0.76439926 -1.345943726
## Row23 0.92078686 -1.380618543
## Row24 0.97904291 -1.395238182
## Row25 0.63127498 0.576843413
## Row26 0.27406930 -1.964097643
## Row27 1.11557505 -1.309365584
## Row28 1.01882165 -0.825794682
# Standard errors for parameters:
fitnb$sd
## $theta
## LV1 LV2
## Alopacce 0.0000000 0.0000000
## Alopcune 0.2736279 0.0000000
## Alopfabr 0.3301706 0.2043807
## Arctlute 0.6749238 0.7849031
## Arctperi 0.4757364 0.5592237
## Auloalbi 0.4407884 0.3892597
## Pardlugu 0.2275838 0.2005392
## Pardmont 0.2196723 0.1816359
## Pardnigr 0.4734857 0.3799201
## Pardpull 0.4966159 0.4106411
## Trocterr 0.2098677 0.1568714
## Zoraspin 0.3480341 0.2796558
##
## $sigma.lv
## LV1 LV2
## 0.4163199 0.4257772
##
## $beta0
## Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi Pardlugu Pardmont
## 0.4518413 0.4630228 0.5478096 1.7104771 1.4645880 0.7893086 0.4391421 0.3996804
## Pardnigr Pardpull Trocterr Zoraspin
## 0.7870448 0.8340612 0.3339019 0.5830956
##
## $inv.phi
## Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi Pardlugu Pardmont
## 0.6025083 0.6026042 0.3957454 0.5787042 1.8042603 0.8383739 0.5810509 0.7756862
## Pardnigr Pardpull Trocterr Zoraspin
## 1.3150293 1.2472530 6.8923550 1.2850988
##
## $phi
## Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi Pardlugu
## 0.35313891 0.27937811 0.79591192 0.71829958 0.33180688 0.35388362 0.43717515
## Pardmont Pardnigr Pardpull Trocterr Zoraspin
## 0.29633701 0.20365546 0.25858012 0.04092229 0.19413856E4. Fit model with different numbers of latent variables.
Show the answers.
4. Default number of latent variables is 2. Let’s try 1 and 3 latent variables as well:
# In exercise 2, we fitted GLLVM with two latent variables
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
# How about 1 or 3 LVs
fitnb1 <- gllvm(y=spider$abund, family = "negative.binomial", num.lv = 1)
fitnb1
## Call:
## gllvm(y = spider$abund, family = "negative.binomial", num.lv = 1)
## family:
## [1] "negative.binomial"
## method:
## [1] "VA"
##
## log-likelihood: -759.4012
## Residual degrees of freedom: 300
## AIC: 1590.802
## AICc: 1599.712
## BIC: 1728.218
getLV(fitnb1)
## LV1
## Row1 -0.81680896
## Row2 -0.94960668
## Row3 -0.91687412
## Row4 -1.10569740
## Row5 -1.22148426
## Row6 -1.03402453
## Row7 -1.26770805
## Row8 -0.53474588
## Row9 0.36604841
## Row10 1.00893873
## Row11 0.31290080
## Row12 -0.31672766
## Row13 -1.17621724
## Row14 -1.03498983
## Row15 -0.09007961
## Row16 -0.27457284
## Row17 -0.02445335
## Row18 -0.02078993
## Row19 -0.12968110
## Row20 -0.06887949
## Row21 -0.04322797
## Row22 1.54084587
## Row23 1.65026381
## Row24 1.68381227
## Row25 -0.28900450
## Row26 1.92918533
## Row27 1.62340700
## Row28 1.20021556
fitnb3 <- gllvm(y=spider$abund, family = "negative.binomial", num.lv = 3)
fitnb3
## Call:
## gllvm(y = spider$abund, family = "negative.binomial", num.lv = 3)
## family:
## [1] "negative.binomial"
## method:
## [1] "VA"
##
## log-likelihood: -733.6806
## Residual degrees of freedom: 279
## AIC: 1581.361
## AICc: 1605.145
## BIC: 1798.937
getLV(fitnb3)
## LV1 LV2 LV3
## Row1 0.92853206 1.270518792 1.442656e-05
## Row2 -0.68270629 0.793766171 -8.316391e-06
## Row3 0.67620584 1.274441623 1.274558e-05
## Row4 -0.24717658 1.149822689 -1.650495e-06
## Row5 -0.36104853 1.217208844 -2.658857e-06
## Row6 -0.35231766 1.008662441 -8.923418e-06
## Row7 0.02922554 1.426035357 1.303585e-06
## Row8 -1.52132542 -0.070757571 -1.614037e-05
## Row9 0.91074989 0.001246445 -1.715171e-05
## Row10 1.03737383 -0.493619564 -1.267990e-05
## Row11 1.48872754 0.303020073 -1.235111e-05
## Row12 1.32415144 0.807380178 -1.107074e-05
## Row13 0.11538557 1.356608431 5.898143e-06
## Row14 -0.43680735 1.014565392 1.794549e-07
## Row15 -1.31923869 -0.544698780 9.511407e-06
## Row16 -1.17465539 -0.242723147 4.184919e-06
## Row17 -1.11252716 -0.512820326 1.006938e-05
## Row18 -1.15895359 -0.542988543 1.150770e-05
## Row19 -1.34064255 -0.502799998 3.281065e-06
## Row20 -1.36731532 -0.592102183 4.863251e-06
## Row21 -1.13968268 -0.475648971 -2.590011e-06
## Row22 0.76441614 -1.345945278 1.463777e-06
## Row23 0.92082974 -1.380606945 1.898252e-06
## Row24 0.97908360 -1.395215043 -2.739204e-06
## Row25 0.63126723 0.576832440 1.211589e-05
## Row26 0.27409901 -1.964094875 3.665146e-06
## Row27 1.11558557 -1.309359352 -2.385965e-06
## Row28 1.01883916 -0.825784128 1.569796e-06E5. Include environmental variables to the GLLVM and explore the model fit.
Show the answers.
5. Environmental variables can be included
with an argument X:
fitnbx <- gllvm(y = spider$abund, X = spider$x, family = "negative.binomial", seed = 123, num.lv = 2)
fitnbx
## Call:
## gllvm(y = spider$abund, X = spider$x, family = "negative.binomial",
## num.lv = 2, seed = 123)
## family:
## [1] "negative.binomial"
## method:
## [1] "VA"
##
## log-likelihood: -588.9187
## Residual degrees of freedom: 217
## AIC: 1415.837
## AICc: 1548.06
## BIC: 1870.074
coef(fitnbx)
## $Species.scores
## LV1 LV2
## Alopacce 1.000000e+00 0.00000000
## Alopcune 2.201376e+00 1.00000000
## Alopfabr 1.600175e+00 -0.02790578
## Arctlute 2.864153e+00 0.85491601
## Arctperi -2.611968e-07 -2.41258495
## Auloalbi 1.491118e+00 3.19928366
## Pardlugu 3.944288e-01 -0.02065962
## Pardmont 8.907405e-01 0.15136216
## Pardnigr 2.417516e+00 -3.48347325
## Pardpull 1.680954e+00 2.24819374
## Trocterr 1.030635e+00 -0.67087884
## Zoraspin 1.570957e+00 -2.45793904
##
## $sigma.lv
## LV1 LV2
## 4.769211e-01 1.750845e-08
##
## $Intercept
## Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi
## -1.8064046 -4.6628013 -4.2140438 -15.5681238 -14.9088011 -16.3558605
## Pardlugu Pardmont Pardnigr Pardpull Trocterr Zoraspin
## 7.9303148 -6.3995949 -5.6547891 -13.6114971 -0.3336265 -4.1372776
##
## $Xcoef
## soil.dry bare.sand fallen.leaves moss herb.layer
## Alopacce -0.8878738 -0.009315833 -0.291126554 0.274014193 0.65118653
## Alopcune 1.3901366 -0.516772190 0.007261267 -0.136642394 0.48662594
## Alopfabr -0.7974089 0.823100137 -0.145041818 0.531349894 0.42565781
## Arctlute 6.2604058 -0.006996697 -1.824637605 0.502257514 -0.09695928
## Arctperi -1.7801727 0.238572359 -2.157239869 0.495821007 0.02126666
## Auloalbi 0.3204160 -0.134778257 0.608841892 0.129651631 4.18697844
## Pardlugu -2.5394813 -0.732449351 0.364909249 -0.358249299 0.60351399
## Pardmont 1.6381786 0.137509982 -0.546863657 0.917166651 0.70949096
## Pardnigr 2.6361823 -0.161228714 -0.857110295 -0.350929172 0.84806179
## Pardpull 3.0903528 -0.454461209 -0.478489258 0.447283079 2.07099025
## Trocterr 1.2578039 -0.263436888 -0.228903368 -0.158047787 0.47711404
## Zoraspin 2.1516927 0.046123722 -0.577840336 -0.005416612 0.67469059
## reflection
## Alopacce 0.69602739
## Alopcune 0.31252195
## Alopfabr 0.52347182
## Arctlute -0.85619380
## Arctperi 4.00643750
## Auloalbi -0.66909955
## Pardlugu -1.19124727
## Pardmont 0.04132679
## Pardnigr -0.68171846
## Pardpull -0.48741030
## Trocterr -0.33629434
## Zoraspin -1.00338851
##
## $inv.phi
## Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi
## 2.063410e+01 3.532395e+00 9.535779e+00 2.266454e+00 4.640004e+07 1.096939e+01
## Pardlugu Pardmont Pardnigr Pardpull Trocterr Zoraspin
## 3.466149e+01 2.199583e+00 9.423731e+00 1.904354e+01 2.691070e+01 6.434174e+00
##
## $phi
## Alopacce Alopcune Alopfabr Arctlute Arctperi Auloalbi
## 4.846347e-02 2.830941e-01 1.048682e-01 4.412180e-01 2.155171e-08 9.116278e-02
## Pardlugu Pardmont Pardnigr Pardpull Trocterr Zoraspin
## 2.885047e-02 4.546317e-01 1.061151e-01 5.251125e-02 3.715994e-02 1.554201e-01
# confidence intervals for parameters:
confint(fitnbx)
## 2.5 % 97.5 %
## theta.LV1.1 3.412564e-01 6.125858e-01
## theta.LV1.2 -2.252097e-05 2.255599e-05
## theta.LV1.3 1.000000e+00 1.000000e+00
## theta.LV1.4 2.197535e+00 2.205217e+00
## theta.LV1.5 1.594203e+00 1.606147e+00
## theta.LV1.6 2.861504e+00 2.866801e+00
## theta.LV1.7 -8.936370e-04 8.931146e-04
## theta.LV1.8 1.485669e+00 1.496567e+00
## theta.LV1.9 3.907375e-01 3.981202e-01
## theta.LV1.10 8.840093e-01 8.974717e-01
## theta.LV1.11 2.400409e+00 2.434624e+00
## theta.LV1.12 1.641161e+00 1.720746e+00
## theta.LV2.1 9.720619e-01 1.089209e+00
## theta.LV2.2 1.563848e+00 1.578067e+00
## theta.LV2.3 0.000000e+00 0.000000e+00
## theta.LV2.4 1.000000e+00 1.000000e+00
## theta.LV2.5 -2.790578e-02 -2.790578e-02
## theta.LV2.6 8.549160e-01 8.549160e-01
## theta.LV2.7 -2.412585e+00 -2.412585e+00
## theta.LV2.8 3.199284e+00 3.199284e+00
## theta.LV2.9 -2.065962e-02 -2.065962e-02
## theta.LV2.10 1.513622e-01 1.513622e-01
## theta.LV2.11 -3.483473e+00 -3.483473e+00
## theta.LV2.12 2.248194e+00 2.248194e+00
## sigma.LV1 -6.708788e-01 -6.708788e-01
## sigma.LV2 -2.457939e+00 -2.457939e+00
## Intercept.Alopacce -1.812280e+00 -1.800529e+00
## Intercept.Alopcune -4.674639e+00 -4.650964e+00
## Intercept.Alopfabr -4.224739e+00 -4.203349e+00
## Intercept.Arctlute -1.558161e+01 -1.555464e+01
## Intercept.Arctperi -1.491537e+01 -1.490223e+01
## Intercept.Auloalbi -1.636441e+01 -1.634732e+01
## Intercept.Pardlugu 7.922342e+00 7.938288e+00
## Intercept.Pardmont -6.409584e+00 -6.389605e+00
## Intercept.Pardnigr -5.664810e+00 -5.644768e+00
## Intercept.Pardpull -1.361989e+01 -1.360310e+01
## Intercept.Trocterr -3.410847e-01 -3.261682e-01
## Intercept.Zoraspin -4.147216e+00 -4.127339e+00
## Xcoef.soil.dry:Alopacce -9.237258e-01 -8.520218e-01
## Xcoef.soil.dry:Alopcune 1.349472e+00 1.430801e+00
## Xcoef.soil.dry:Alopfabr -8.197262e-01 -7.750915e-01
## Xcoef.soil.dry:Arctlute 6.218875e+00 6.301936e+00
## Xcoef.soil.dry:Arctperi -1.787646e+00 -1.772700e+00
## Xcoef.soil.dry:Auloalbi 2.865838e-01 3.542481e-01
## Xcoef.soil.dry:Pardlugu -2.558985e+00 -2.519978e+00
## Xcoef.soil.dry:Pardmont 1.570186e+00 1.706171e+00
## Xcoef.soil.dry:Pardnigr 2.595635e+00 2.676729e+00
## Xcoef.soil.dry:Pardpull 3.043656e+00 3.137049e+00
## Xcoef.soil.dry:Trocterr 1.224275e+00 1.291333e+00
## Xcoef.soil.dry:Zoraspin 2.116890e+00 2.186495e+00
## Xcoef.bare.sand:Alopacce -1.074288e-01 8.879716e-02
## Xcoef.bare.sand:Alopcune -5.297806e-01 -5.037638e-01
## Xcoef.bare.sand:Alopfabr 7.961501e-01 8.500502e-01
## Xcoef.bare.sand:Arctlute -1.648518e-02 2.491786e-03
## Xcoef.bare.sand:Arctperi 2.142026e-01 2.629421e-01
## Xcoef.bare.sand:Auloalbi -1.700558e-01 -9.950076e-02
## Xcoef.bare.sand:Pardlugu -7.513930e-01 -7.135057e-01
## Xcoef.bare.sand:Pardmont 6.414102e-02 2.108789e-01
## Xcoef.bare.sand:Pardnigr -2.007248e-01 -1.217326e-01
## Xcoef.bare.sand:Pardpull -5.101206e-01 -3.988018e-01
## Xcoef.bare.sand:Trocterr -3.811233e-01 -1.457505e-01
## Xcoef.bare.sand:Zoraspin 2.892192e-02 6.332553e-02
## Xcoef.fallen.leaves:Alopacce -2.949788e-01 -2.872743e-01
## Xcoef.fallen.leaves:Alopcune -9.378432e-02 1.083068e-01
## Xcoef.fallen.leaves:Alopfabr -1.467799e-01 -1.433038e-01
## Xcoef.fallen.leaves:Arctlute -1.837483e+00 -1.811792e+00
## Xcoef.fallen.leaves:Arctperi -2.157240e+00 -2.157240e+00
## Xcoef.fallen.leaves:Auloalbi 4.899677e-01 7.277161e-01
## Xcoef.fallen.leaves:Pardlugu 2.962014e-01 4.336171e-01
## Xcoef.fallen.leaves:Pardmont -6.051007e-01 -4.886266e-01
## Xcoef.fallen.leaves:Pardnigr -9.480451e-01 -7.661755e-01
## Xcoef.fallen.leaves:Pardpull -5.683111e-01 -3.886675e-01
## Xcoef.fallen.leaves:Trocterr -2.939830e-01 -1.638237e-01
## Xcoef.fallen.leaves:Zoraspin -6.752642e-01 -4.804165e-01
## Xcoef.moss:Alopacce 2.429067e-01 3.051217e-01
## Xcoef.moss:Alopcune -1.658938e-01 -1.073910e-01
## Xcoef.moss:Alopfabr 5.010307e-01 5.616691e-01
## Xcoef.moss:Arctlute 4.825618e-01 5.219532e-01
## Xcoef.moss:Arctperi 4.759735e-01 5.156686e-01
## Xcoef.moss:Auloalbi 1.078977e-01 1.514056e-01
## Xcoef.moss:Pardlugu -3.851409e-01 -3.313577e-01
## Xcoef.moss:Pardmont 8.346161e-01 9.997172e-01
## Xcoef.moss:Pardnigr -3.981595e-01 -3.036989e-01
## Xcoef.moss:Pardpull 3.261368e-01 5.684293e-01
## Xcoef.moss:Trocterr -2.748672e-01 -4.122841e-02
## Xcoef.moss:Zoraspin -2.378061e-02 1.294739e-02
## Xcoef.herb.layer:Alopacce 6.000648e-01 7.023083e-01
## Xcoef.herb.layer:Alopcune 4.402765e-01 5.329754e-01
## Xcoef.herb.layer:Alopfabr 3.818734e-01 4.694422e-01
## Xcoef.herb.layer:Arctlute -1.558862e-01 -3.803234e-02
## Xcoef.herb.layer:Arctperi 5.889576e-03 3.664374e-02
## Xcoef.herb.layer:Auloalbi 4.151275e+00 4.222682e+00
## Xcoef.herb.layer:Pardlugu 5.295818e-01 6.774462e-01
## Xcoef.herb.layer:Pardmont 6.420815e-01 7.769004e-01
## Xcoef.herb.layer:Pardnigr 8.043680e-01 8.917556e-01
## Xcoef.herb.layer:Pardpull 2.017884e+00 2.124096e+00
## Xcoef.herb.layer:Trocterr 4.176353e-01 5.365928e-01
## Xcoef.herb.layer:Zoraspin 6.309830e-01 7.183982e-01
## Xcoef.reflection:Alopacce 6.718149e-01 7.202399e-01
## Xcoef.reflection:Alopcune 2.611065e-01 3.639374e-01
## Xcoef.reflection:Alopfabr 4.880824e-01 5.588612e-01
## Xcoef.reflection:Arctlute -8.874129e-01 -8.249747e-01
## Xcoef.reflection:Arctperi 3.979206e+00 4.033669e+00
## Xcoef.reflection:Auloalbi -7.274869e-01 -6.107122e-01
## Xcoef.reflection:Pardlugu -1.243111e+00 -1.139383e+00
## Xcoef.reflection:Pardmont -1.403512e-02 9.668869e-02
## Xcoef.reflection:Pardnigr -7.292709e-01 -6.341661e-01
## Xcoef.reflection:Pardpull -5.525356e-01 -4.222850e-01
## Xcoef.reflection:Trocterr -4.643246e-01 -2.082641e-01
## Xcoef.reflection:Zoraspin -1.046829e+00 -9.599477e-01
## inv.phi.Alopacce 2.062497e+01 2.064323e+01
## inv.phi.Alopcune 3.519971e+00 3.544819e+00
## inv.phi.Alopfabr 9.528301e+00 9.543258e+00
## inv.phi.Arctlute 2.265933e+00 2.266974e+00
## inv.phi.Arctperi 4.640004e+07 4.640004e+07
## inv.phi.Auloalbi 1.096598e+01 1.097280e+01
## inv.phi.Pardlugu 3.464960e+01 3.467337e+01
## inv.phi.Pardmont 2.192376e+00 2.206789e+00
## inv.phi.Pardnigr 9.408806e+00 9.438655e+00
## inv.phi.Pardpull 1.903077e+01 1.905630e+01
## inv.phi.Trocterr 2.674175e+01 2.707965e+01
## inv.phi.Zoraspin 6.420768e+00 6.447580e+00Problems? 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.
Ordination
GLLVM as a model based ordination method
GLLVMs can be used as a model-based approach to unconstrained ordination by including two latent variables in the model: \(g(E(y_{ij})) = \beta_{0j} + \boldsymbol{u}_i'\boldsymbol{\theta}_j\)
The latent variable term try to capture the underlying factors driving species abundances at sites.
Predictions for the two latent variables, \(\boldsymbol{\hat u}_i=(\hat u_{i1}, \hat u_{i2})\), then provide coordinates for sites in the ordination plot and then provides a graphical representation of which sites are similar in terms of their species composition.
Ordination plot
-
ordiplot()produces ordination plots based on fitted GLLVMs. - Uncertainty of the ordination points in model based ordination can be assessed with prediction regions based on the prediction errors of latent variables.
- Prediction regions may also help interpreting which differences between ordination points are really a sign of the difference between species composition at those sites.
Biplot
- Between species correlations can be visualized with biplot8 by adding
latent variable loadings \(\boldsymbol{\theta}_j\) to the ordination
of sites, by producing a biplot, (argument
biplot = TRUEinordiplot()). - In a biplot latent variables and their loadings are rotated so that the LV loadings of the species are in the same direction with the sites where they are most abundant.
- Biplot can be used for finding groups of correlated species or finding indicator species common at specific sites.
- For example, species named Pardlugu is common in the group of sites located on the bottom (eg sites 8, 19, 20 and 21).
- x and y axes may help the interpretation.
Environmental gradients
- The potential impact of environmental variables on species communities can be viewed by coloring ordination points according to the variables.
- For example, species named Pardlugu seems to prefer sites with lot of dry soil mass and low reflection of the soil surface with a cloudless sky.
# Arbitrary color palette, a vector length of 20. Can use, for example, colorRampPalette from package grDevices
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
par(mfrow = c(3,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)
}
Here environmental gradients stand out quite clearly, indicating that, at least, some of the differences in species compositions at sites can be explained by the differences in environmental conditions.
The next step would be to include covariates to the model to study more precisely the effects of environmental variables: \(g(E(y_{ij})) = \beta_{0j} + \boldsymbol{x}_i'\boldsymbol{\beta}_{j} + \boldsymbol{u}_i'\boldsymbol{\theta}_j\)
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↩︎
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↩︎
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↩︎
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 .↩︎
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.↩︎
van der Aart, P. J. M., and Smeenk-Enserink, N. (1975) Correlations between distributions of hunting spiders (Lycosidae, Ctenidae) and environmental characteristics in a dune area. Netherlands Journal of Zoology 25, 1-45.↩︎
Dunn, P. K., and Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5, 236-244.↩︎
Gabriel, K. R. (1971). The biplot graphic display of matrices with application to principal component analysis. Biometrika, 58, 453-467↩︎