| Title: | Generalized Linear Latent Variable Models |
|---|---|
| Description: | Analysis of multivariate data using generalized linear latent variable models (gllvm). Estimation is performed using either the Laplace method, variational approximations, or extended variational approximations, implemented via TMB (Kristensen et al. (2016), <doi:10.18637/jss.v070.i05>). |
| Authors: | Jenni Niku [aut, cre], Wesley Brooks [aut], Riki Herliansyah [aut], Francis K.C. Hui [aut], Pekka Korhonen [aut], Sara Taskinen [aut], Bert van der Veen [aut], David I. Warton [aut] |
| Maintainer: | Jenni Niku <[email protected]> |
| License: | GPL-2 |
| Version: | 2.0 |
| Built: | 2024-11-27 01:50:34 UTC |
| Source: | CRAN |
Calculates corrected Akaike information criterion for small sample sizes, and extracts number of observations.
## S3 method for class 'gllvm' AICc(object, ...) ## S3 method for class 'gllvm' nobs(object, ...)## S3 method for class 'gllvm' AICc(object, ...) ## S3 method for class 'gllvm' nobs(object, ...)
object |
an object of class 'gllvm'. |
... |
Not used. |
Jenni Niku, Bert van der Veen
Computes an analysis of deviance table for two generalized linear latent variable model fits.
## S3 method for class 'gllvm' anova(object, ...)## S3 method for class 'gllvm' anova(object, ...)
object |
an object of class 'gllvm'. |
... |
one or more objects of class 'gllvm' |
Computes likelihood-ratio test for two or more gllvm models. Test results makes sense only for nested models. Notice also that this test is not designed for testing models which have degrees of freedom difference larger than 20. For such models the P-value should be treated as very approximate.
Jenni Niku
## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- antTraits$abund X <- antTraits$env TR <- antTraits$traits # Fit gllvm model fit1 <- gllvm(y, X, TR, formula = ~ Bare.ground + Shrub.cover, family = poisson()) fit2 <- gllvm(y, X, TR, formula = ~ Bare.ground + Shrub.cover + (Bare.ground + Shrub.cover) : Webers.length, family = poisson()) # Test if the model with fourth corner interaction terms is significantly # better using likelihood-ratio test: anova(fit1, fit2)## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- antTraits$abund X <- antTraits$env TR <- antTraits$traits # Fit gllvm model fit1 <- gllvm(y, X, TR, formula = ~ Bare.ground + Shrub.cover, family = poisson()) fit2 <- gllvm(y, X, TR, formula = ~ Bare.ground + Shrub.cover + (Bare.ground + Shrub.cover) : Webers.length, family = poisson()) # Test if the model with fourth corner interaction terms is significantly # better using likelihood-ratio test: anova(fit1, fit2)
These data describe the abundance of ground beetle assemblages from the Scottish landscape alongside with the environmental data and species traits. The data includes abundances of 68 species of ground beetle species.
data(beetle)data(beetle)
A data frame of species composition of the ground-beetle assemblages.
A data frame of study design variables and environmental data.
Unique id for sample.
Land use type.
Grid where sample were collected.
sampling area in Scotland.
Sampling year
1, peat; 2, peaty loam; 3, loamy sand; 4, sandy loam; 5, sandy clay loam; 6, sandy silt loam; 7, silty clay
organic content (% loss of organic content on ignition), log10 transformed
soil pH
available P (mg/L), log10 transformed
available K (mg/L)
percentage moisture content
percentage cover estimate of bare ground in 11 1-m2 quadrats, arcsine transformed
percentage cover estimate of litter cover in 11 1-m2 quadrats, log10 transformed
percentage cover estimate of bryophytes in 11 1-m2 quadrats, arcsine transformed
number of reproductive stems (flowering or fruiting) in 11 1-m2 quadrats
height canopy height (cm) in 11 1-m2 quadrats
number of stems (ramets) in 100 cm2
dry mass (g) of biomass 0–5 cm from soil surface in 400 cm2
dry mass (g) of biomass .5 cm from soil surface in 400 cm2, log10 transformed
biomass of reproductive parts (flowers and fruits) in 100 cm2, log10 transformed
elevation in m a.s.l.
management intensity index (see Materials and Methods: Environmental data)
A data frame of the species names and species morphological and life trait characteristics.
Species names.
Species codes corresponding species names in abundance matrix
diameter of the eye, measured from above
length of the antenna
maximum width of the pronotum
maximum depth (“vaulting”) of the pronotum
maximum width of the elytra
length of the metafemur (with the articulation segments), from the coxa to the apex
length of the metatrochanter
length of the metatarsi
maximum width of the metafemur
total length (length of the pronotum in the medial line plus length of the elytra, from the medial ridge of the scutellum to the apex)
color of the legs (1, pale; 2, black; 3, metallic)
color of the body (1, pale; 2, black; 3, metallic)
wing development (1, apterous or brachypterous; 2, dimorphic; 3, macropterous)
shape of the pronotum (1, oval; 2, cordiform; 3, trapezoidal)
overwintering (1, only adults; 2, adults and larvae or only larvae)
food of the adult (1, mostly Collembola; 2, generalist predator; 3, mostly plant material)
daily activity (1, only diurnal; 2, nocturnal)
breeding season (1, spring; 2, summer; 3, autumn or winter)
main period of emergence of the adults (1, spring; 2, summer; 3, autumn)
main period of adult activity (1, autumn; 2, summer only)
Beetles were sampled with two parallel rows of nine pitfall traps (diameter 7.5 cm, 2 m apart) at each site, starting in early May.
Detailed description of the data available in the reference and in the Ecological Archives E082-012.
Ignacio Ribera, Sylvain Dolédec, Iain S. Downie, and Garth N. Foster. 2001. Effect of land disturbance and stress on species traits of ground beetle assemblages. Ecology 82:1112-1129.
## Not run: data(beetle) # Abundance matrix Y <- beetle$Y # Environmental data X <- beetle$X # Species traits TR <- beetle$TR ## End(Not run)## Not run: data(beetle) # Abundance matrix Y <- beetle$Y # Environmental data X <- beetle$X # Species traits TR <- beetle$TR ## End(Not run)
Plots covariate coefficients and their confidence intervals.
## S3 method for class 'gllvm' coefplot( object, y.label = TRUE, which.Xcoef = NULL, order = TRUE, cex.ylab = 0.5, cex.xlab = 1.3, mfrow = NULL, mar = c(4, 6, 2, 1), xlim.list = NULL, ... )## S3 method for class 'gllvm' coefplot( object, y.label = TRUE, which.Xcoef = NULL, order = TRUE, cex.ylab = 0.5, cex.xlab = 1.3, mfrow = NULL, mar = c(4, 6, 2, 1), xlim.list = NULL, ... )
object |
an object of class 'gllvm'. |
y.label |
logical, if |
which.Xcoef |
vector indicating which covariate coefficients will be plotted. Can be vector of covariate names or numbers. Default is |
order |
logical, whether or not coefficients are ordered, defaults to |
cex.ylab |
the magnification to be used for axis annotation relative to the current setting of cex. |
cex.xlab |
the magnification to be used for axis annotation. |
mfrow |
same as |
mar |
vector of length 4, which defines the margin sizes: |
xlim.list |
list of vectors with length of two to define the intervals for an x axis in each covariate plot. Defaults to NULL when the interval is defined by the range of point estimates and confidence intervals |
... |
additional graphical arguments. |
Jenni Niku <[email protected]>, Francis K.C. Hui, Sara Taskinen, Bert van der Veen
# Extract subset of the microbial data to be used as an example data(microbialdata) X <- microbialdata$Xenv y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, X, formula = ~ pH + Phosp, family = poisson()) coefplot(fit) ## Not run: # Fit gllvm model with environmental covariances and reduced rank fitRR <- gllvm(y = y, X = X, num.RR = 2, family = "negative.binomial") coefplot(fitRR) ## End(Not run)# Extract subset of the microbial data to be used as an example data(microbialdata) X <- microbialdata$Xenv y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, X, formula = ~ pH + Phosp, family = poisson()) coefplot(fit) ## Not run: # Fit gllvm model with environmental covariances and reduced rank fitRR <- gllvm(y = y, X = X, num.RR = 2, family = "negative.binomial") coefplot(fitRR) ## End(Not run)
Computes confidence intervals for parameters in a fitted gllvm model.
## S3 method for class 'gllvm' confint(object, parm = NULL, level = 0.95, ...)## S3 method for class 'gllvm' confint(object, parm = NULL, level = 0.95, ...)
object |
an object of class 'gllvm'. |
parm |
a specification of which parameters are to be given confidence intervals, a vector of names. Examples of options are "beta0", "Xcoef",theta", "phi". If missing, all parameters are considered. |
level |
the confidence level. Scalar between 0 and 1. |
... |
not used. |
Jenni Niku <[email protected]>
## Not run: ## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- as.matrix(antTraits$env[,1:2]) # Fit gllvm model fit <- gllvm(y = y, X = X, family = poisson()) # 95 % confidence intervals for coefficients of X variables confint(fit, level = 0.95, parm = "Xcoef") ## End(Not run)## Not run: ## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- as.matrix(antTraits$env[,1:2]) # Fit gllvm model fit <- gllvm(y = y, X = X, family = poisson()) # 95 % confidence intervals for coefficients of X variables confint(fit, level = 0.95, parm = "Xcoef") ## End(Not run)
Extracts species optima and tolerances, potentially with standard errors (derived with the Delta method).
## S3 method for class 'gllvm' optima(object, sd.errors = TRUE, ...) ## S3 method for class 'gllvm' tolerances(object, sd.errors = TRUE, ...)## S3 method for class 'gllvm' optima(object, sd.errors = TRUE, ...) ## S3 method for class 'gllvm' tolerances(object, sd.errors = TRUE, ...)
object |
an object of class 'gllvm'. |
sd.errors |
logical. If |
... |
Not used. |
Currently no separate method for calculating species maxima or gradient length are implemented.
Gradient length can be inferred from the standard deviation of the latent variables, which is reported by summary.gllvm.
Bert van der Veen
Extended dataset of counts of hunting spiders in a dune area in the Netherlands, measured at 100 pitfall traps.
data(eSpider)data(eSpider)
A data frame with abundances of 12 hunting spider species measured at 100 sites.
A matrix of 26 predictor variables at 28 of the 100 sites.
An vector of indices indicating at which sites the predictor variables were measured.
Counts of hunting spiders in a dune area in the Netherlands, measured with 100 different pitfall traps. This dataset was published with permission from the CANOCO FORTRAN package (version 4 or higher) example datasets.
Species names have been abbreviated, corresponding to: Alopacce = Alopecosa accentuata, Alopcune = Alopecosa cuneata, Alopfabr = Alopecosa fabrilis, Arctlute = Arctosa lutetiana, Arctperi = Arctosa perita, Auloalbi = Aulonia albimana, Pardlugu = Pardosa lugubris, Pardmont = Pardosa monticola, Pardnigr = Pardosa nigriceps, Pardpull = Pardosa pullata, Trocterr = Trochosa terricola, Zoraspin = Zora spinimana.
Environmental measurements were taken at 28 of the 100 pitfall traps measuring soil properties(Water content: "conWate", Humus content: "conHumu", Acidity (pH-KCl)), vegetation ("BareSand": percentage bare sand, "FallTwig": cover on the ground by leaves and twigs, "CovMoss": cover by mosses and lichens, "CovHerb": cover by the herb and grass layer (including maximum height, minimum height, "CovCala": cover by Calamagrostis epigejos, cover by Carex arenaria, "CovFest": cover by Festuca ovina, "CovCory": cover by Corynephorus canescens, "CovUrti": cover by Urtica dioica, "CovMoeh": cover by Moehringia trinervia), "CovShru": cover by the shrub layer (minimum and maximum height, and "CovLigu": cover by Ligustrum vulgare), "CovTree": cover by the tree layer (including maximum height, cover by Populus tremula or Crataegus monogyna)), and light properties ("LuxGrey": lux at equal grey sky, "LuxSun": lux at cloudless sky, "LuxRef": lux by reflection of the soil surface)
The original publication of Canonical Correspondence Analysis used standardized versions of the log and log1p transformed predictors "ConWate", "BareSand", "FallTwig", "CovMoss", "CovHerb", "LuxRef".
ter Braak, C.J.F. and Smilauer, P. (1998). CANOCO reference manual and user's guide to CANOCO for Windows: software for canonical community ordination (version 4). Microcomputer Power, New York, New York, USA. ter Braak, C.J.F. (1986). Canonical correspondence analysis: a new eigenvector technique for multivariate direct gradient analysis. Ecology, 67(5), 1167-1179. 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), 1-45.
data(eSpider) Y <- eSpider$abund[eSpider$nonNA, ] X <- eSpider$X[eSpider$nonNA, ] model <- gllvm(y = Y, X = X, lv.formula = ~log(ConWate) + log1p(BareSand) + log1p(FallTwig) + log1p(CovMoss) + log1p(CovHerb) + log(LuxRef), num.RR = 2, family = "poisson")data(eSpider) Y <- eSpider$abund[eSpider$nonNA, ] X <- eSpider$X[eSpider$nonNA, ] model <- gllvm(y = Y, X = X, lv.formula = ~log(ConWate) + log1p(BareSand) + log1p(FallTwig) + log1p(CovMoss) + log1p(CovHerb) + log(LuxRef), num.RR = 2, family = "poisson")
Dataset of 1666 binary observations for 215 fungal species, in different 53 European Beech forests spread across 8 regions.
data(fungi)data(fungi)
A data frame with the presence-absences of 215 fungal species measured at 1666 logs.
A data frame of 8 predictor variables.
Diameter at breast height (cm)
Decay stage of logs on a 1-5 scale
Connectivity of the surrounding forest at 10km scale
Annual temperature range in degrees Celsius
Annual precipitation in milimeters
ln(area in hectares) of reserves
Site groups identified based on spatial clusters
Site name
A data frame of the traits used in Abrego et al. (2022).
The phylogenetic tree.
The phylogenetic covariance matrix.
The phylogenetic distance matrix.
Observations of fungi species inhabiting European beech logs, in different European countries. The countries have been grouped in eight different regions. Logs were surveyed in 53 different reserves (or sites). Included environment and trait covariates are limited to those analyzed in the original article, though more are available in the published dataset on datadryad.org.
Abrego, N., Bässler, C., Christensen, M., and Heilmann‐Clausen, J. (2022). Traits and phylogenies modulate the environmental responses of wood‐inhabiting fungal communities across spatial scales. Journal of Ecology, 110(4), 784-798.
Abrego, N., Bässler, C., Christensen, M., and Heilmann-Clausen, J. (2022). Data and code from: Traits and phylogenies modulate the environmental responses of wood-inhabiting fungal communities across spatial scales [Dataset]. Dryad. https://doi.org/10.5061/dryad.t76hdr82r
## Not run: data(fungi) Y <- fungi$Y X <- fungi$X TR <- fungi$TR C <- fungi$C dist <- fungi$dist #model <- gllvm(y = Y, X = cbind(int = 1, X), TR = TR, # formula = ~DBH.CM + AVERDP + I(AVERDP^2) + CONNECT10 + TEMPR + PRECIP + # log.AREA + (DBH.CM + AVERDP + I(AVERDP^2) + CONNECT10 + TEMPR + PRECIP + # log.AREA):(FB.type + Sp.log.vol.µ3 + Lifestyle), # family = "binomial", num.lv = 0, studyDesign = X[,c("REGION", "RESERVE")], # colMat = list(C, dist = dist), colMat.rho.struct = "term", # row.eff = ~(1 | REGION/RESERVE), sd.errors = FALSE, # randomX = ~int + DBH.CM + AVERDP + I(AVERDP^2) + # CONNECT10 + TEMPR + PRECIP + log.AREA, # beta0com = TRUE, nn.colMat = 10, maxit = 20000) ## End(Not run)## Not run: data(fungi) Y <- fungi$Y X <- fungi$X TR <- fungi$TR C <- fungi$C dist <- fungi$dist #model <- gllvm(y = Y, X = cbind(int = 1, X), TR = TR, # formula = ~DBH.CM + AVERDP + I(AVERDP^2) + CONNECT10 + TEMPR + PRECIP + # log.AREA + (DBH.CM + AVERDP + I(AVERDP^2) + CONNECT10 + TEMPR + PRECIP + # log.AREA):(FB.type + Sp.log.vol.µ3 + Lifestyle), # family = "binomial", num.lv = 0, studyDesign = X[,c("REGION", "RESERVE")], # colMat = list(C, dist = dist), colMat.rho.struct = "term", # row.eff = ~(1 | REGION/RESERVE), sd.errors = FALSE, # randomX = ~int + DBH.CM + AVERDP + I(AVERDP^2) + # CONNECT10 + TEMPR + PRECIP + log.AREA, # beta0com = TRUE, nn.colMat = 10, maxit = 20000) ## End(Not run)
Calculates the species environment covariance matrix for a gllvm object.
## S3 method for class 'gllvm' getEnvironCov(object, x = NULL, ...)## S3 method for class 'gllvm' getEnvironCov(object, x = NULL, ...)
object |
an object of class 'gllvm'. |
x |
(optional) vector of covariate values to calculate the covariance for. Defaults to a vector of 1s. If both 'randomX' and random species effects are present in the model this should be a list of length two. |
... |
not used |
Species covariance matrix due to the environment is calculated.
Covariances due to the covariates can only be calculated when random effects are included in the model, and are thus limited to reduced rank models (those including constrained and concurrent ordinations) fitted with random slopes, models fitted with random effects via the formula interface, or the fourth corner model fitted with random slopes. For full rank models with random slopes, i.e., with the formula interface or the fourth corner model, the covariances of species are formulated as:
where is a correlation matrix for the columns in the response (e.g., a Phylogenetic matrix), the signal parameter, and R the covariance matrix for the random effects. Here,
, with x a vector of covariate values for each of the random effects, which defaults to a vector of 1s. when there are covariate-specific phylogenetic signal parameters in the model, this is instead:
where is a block-diagonal lower triangular matrix, and each the lower triangular matrix of the covariance matrix for each covariate.
For reduced rank models, the covariance is separately defined for the different variance structures of the canonical coefficients in the package. With LV-specific variances, we have:
where is the matrix of loadings, and S the (diagonal) covariance matrix for the canonical coefficients. With predictor-specific variances, we instead have:
with I_d an identity matrix for the number of constrained and informed latent variables, and the variance per predictor for the canonical coefficients. When correlations are included, we have:
Expressions for the quadratic models in the package are determined similarly but not documented here for brevity.
Function returns the following components:
cov |
species covariances due to covariates |
trace.randomB |
trace of the covariance matrix due to random canonical coefficients |
trace.randomB.quad |
trace of the covariance matrix components due to squared model terms |
trace.col.eff |
trace of the covariance matrix due to random column (species) effects |
Bert van der Veen
getEnvironCor ,getResidualCov.gllvm, getResidualCor.gllvm,.
## Not run: # Example with the spider dataset data(eSpider) y = eSpider$abund[eSpider$nonNA,] X = eSpider$X[eSpider$nonNA,] fit <- gllvm(y, X = scale(X), num.RR = 2, randomB = "P", family = "negative.binomial") envcov <- getEnvironCov(fit) envcov$trace.randomB # As proportion of variance in the model envcov$trace.randomB/sum(envcov$trace.randomB) ## End(Not run)## Not run: # Example with the spider dataset data(eSpider) y = eSpider$abund[eSpider$nonNA,] X = eSpider$X[eSpider$nonNA,] fit <- gllvm(y, X = scale(X), num.RR = 2, randomB = "P", family = "negative.binomial") envcov <- getEnvironCov(fit) envcov$trace.randomB # As proportion of variance in the model envcov$trace.randomB/sum(envcov$trace.randomB) ## End(Not run)
Extract loadings (species scores) from a gllvm object.
## S3 method for class 'gllvm' getLoadings(object, ...)## S3 method for class 'gllvm' getLoadings(object, ...)
object |
an object of class 'gllvm'. |
... |
not used |
Function retrieves the loadings a.k.a. species scores for a GLLVM. For the optima of a quadratic response model, see optima.gllvm
Extract latent variables from gllvm object.
## S3 method for class 'gllvm' getLV(object, type = NULL, ...)## S3 method for class 'gllvm' getLV(object, type = NULL, ...)
object |
an object of class 'gllvm'. |
type |
type of latent variable scores to retrieve from a gllvm object. For models with unconstrained latent variables, defaults to "residual". For models with constrained latent variables, defaults to conditional. Alternatively, "marginal" returns linear combination scores without residual error. |
... |
not used |
Function retrieves the site scores for a GLLVM. Each type corresponds to a separate term of the model. For a GLLVM with unconstrained latent variables the default is "residual". "Residual" scores represent the error term in concurrent ordination, and are not available for constrained ordination.
For GLLVMs with informed latent variables, "conditional" returns the complete site scores, due to both fixed- and latent effects, where the latent effect is always scaled by the diagonal of the species loadings so that it can be small relative to the fixed-effects. "Conditional" here means conditional on the random-effect i.e. the residual.
Type "marginal" returns linear combination scores, i.e. the site scores only due to fixed-effects. These are available for constrained and concurrent ordination.
If both unconstrained and constrained latent variables are included in the model, type "marginal" returns linear combination scores for constrained latent variables but "residual" scores for unconstrained latent variables.
Calculates the prediction errors for latent variables and random effects for gllvm model.
## S3 method for class 'gllvm' getPredictErr(object, CMSEP = TRUE, cov = FALSE, ...)## S3 method for class 'gllvm' getPredictErr(object, CMSEP = TRUE, cov = FALSE, ...)
object |
an object of class 'gllvm'. |
CMSEP |
logical, if |
cov |
if |
... |
not used |
Calculates conditional mean squared errors for predictions. If variational approximation is used, prediction errors can be based on covariances of the variational distributions, and therefore they do not take into account the uncertainty in the estimation of (fixed) parameters.
Function returns following components:
lvs |
prediction errors for latent variables |
row.effects |
prediction errors for random row effects if included |
Francis K.C. Hui, Jenni Niku, David I. Warton
## Not run: # Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # prediction errors for latent variables: getPredictErr(fit) ## End(Not run)## Not run: # Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # prediction errors for latent variables: getPredictErr(fit) ## End(Not run)
Calculates the residual correlation matrix for gllvm model.
## S3 method for class 'gllvm' getResidualCor(object, adjust = 1, x = NULL, ...)## S3 method for class 'gllvm' getResidualCor(object, adjust = 1, x = NULL, ...)
object |
an object of class 'gllvm'. |
adjust |
The type of adjustment used for negative binomial and binomial distribution when computing residual correlation matrix. Options are 0 (no adjustment), 1 (the default adjustment) and 2 (alternative adjustment for NB distribution). See details. |
x |
(optional) vector of covariate values to calculate the covariance for, when applicable. |
... |
not used |
Residual correlation matrix is calculated based on the residual covariance matrix, see details from getResidualCov.gllvm.
Francis K.C. Hui, Jenni Niku, David I. Warton
#'# Extract subset of the microbial data to be used as an example data(microbialdata) y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, family = poisson()) fit$logL cr <- getResidualCor(fit) cr[1:5,1:5] ## Not run: # Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # residual correlations: cr <- getResidualCor(fit) # Plot residual correlations: install.packages("corrplot", "gclus") library(corrplot) library(gclus) corrplot(cr[order.single(cr), order.single(cr)], diag = F, type = "lower", method = "square", tl.cex = 0.8, tl.srt = 45, tl.col = "red") ## End(Not run)#'# Extract subset of the microbial data to be used as an example data(microbialdata) y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, family = poisson()) fit$logL cr <- getResidualCor(fit) cr[1:5,1:5] ## Not run: # Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # residual correlations: cr <- getResidualCor(fit) # Plot residual correlations: install.packages("corrplot", "gclus") library(corrplot) library(gclus) corrplot(cr[order.single(cr), order.single(cr)], diag = F, type = "lower", method = "square", tl.cex = 0.8, tl.srt = 45, tl.col = "red") ## End(Not run)
Calculates the residual covariance matrix for gllvm model.
## S3 method for class 'gllvm' getResidualCov(object, adjust = 1, x = NULL, ...)## S3 method for class 'gllvm' getResidualCov(object, adjust = 1, x = NULL, ...)
object |
an object of class 'gllvm'. |
adjust |
The type of adjustment used for negative binomial, binomial and normal distribution when computing residual correlation matrix. Options are 0 (no adjustment), 1 (the default adjustment) and 2 (alternative adjustment for NB distribution), see details. |
x |
(optional) vector of covariate values to calculate the covariance for, when applicable. |
... |
not used. |
Residual covariance matrix, storing information on species co-occurrence that is not explained by the environmental variables (if included), is calculated using the matrix of latent variables loadings, that is, , and the dispersion parameter related to the distribution of choice, is applicable (e.g. in the case of negative-binomial distributed responses).
When the responses are modelled using the negative binomial distribution, the residual variances for each species must be adjusted for overdispersion. The two possible adjustment terms are (adjust = 1) and (adjust = 2), where is the trigamma function.
The negative binomial model can be written using different parameterizations.
The residual covariance with adjust = 1 can be obtained using the lognormal-Poisson parametrization, that is,
where and and . Now and variance , which are the same as for the NB distribution.
Therefore, on linear predictor scale, we have the variance
which leads to the residual covariance matrix , where is the diagonal matrix with as diagonal elements (adjust = 1).
Or, for a GLLVM where species are a quadratic function of the latent variables, we instead have
which leads to the residual covariance matrix , where holds the quadratic coefficients.
Since the quadratic coefficients are constrained to be positive, the residual covariance in the latter case is, given the same coefficients on the linear term, equal or more positive than in the linear case.
The residual covariance matrix with adjust = 2 can be obtained by using Poisson-Gamma parametri-zation
where and is as above. The mean and the variance are of similar form as above and we have that
where is the trigamma function.
In the case of binomial distribution, the adjustment terms (adjust = 1) are 1 for probit link and for logit link.
These are obtained by treating binomial model as latent variable model. Assume
where for probit model, and for logit model.
Then binary response is defined as , if and 0 otherwise.
Now we have that which leads to probit and logit models.
On linear predictor scale we then have that
For the probit model, the residual covariance matrix is then , and for the logit model .
Similarly as above, for a GLLVM where species are a quadratic function of the latent variables, the term is added to the residual covariance matrix.
For normal distribution, we can write
where and thus we have that
For the gaussian model, the residual covariance matrix is then .
Function returns following components:
cov |
residual covariance matrix |
trace |
trace of the residual covariance matrix, the total variance explained |
var.q |
trace of the residual covariance matrix per latent variable, variance explained per latent variable |
var.q2 |
trace of the squared term of the residual covariance matrix per latent variable, for quadratic responses. Variance explained per latent variable by the quadratic term |
Francis K.C. Hui, Jenni Niku, David I. Warton, Bert van der Veen
## Not run: # Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # residual covariance: rescov <- getResidualCov(fit) rescov$cov # Trace of the covariance matrix rescov$trace # Variance explained per latent variable rescov$var.q ## End(Not run)## Not run: # Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # residual covariance: rescov <- getResidualCov(fit) rescov$cov # Trace of the covariance matrix rescov$trace # Variance explained per latent variable rescov$var.q ## End(Not run)
Fits generalized linear latent variable model for multivariate data. The model can be fitted using Laplace approximation method or variational approximation method.
gllvm( y = NULL, X = NULL, TR = NULL, data = NULL, formula = NULL, family, num.lv = NULL, num.lv.c = 0, num.RR = 0, lv.formula = NULL, lvCor = NULL, studyDesign = NULL, dist = list(matrix(0)), distLV = matrix(0), colMat = NULL, colMat.rho.struct = "single", corWithin = FALSE, corWithinLV = FALSE, quadratic = FALSE, row.eff = FALSE, sd.errors = TRUE, offset = NULL, method = "VA", randomB = FALSE, randomX = NULL, beta0com = FALSE, zeta.struc = "species", plot = FALSE, link = "probit", Ntrials = 1, Power = 1.1, seed = NULL, scale.X = TRUE, return.terms = TRUE, gradient.check = FALSE, disp.formula = NULL, control = list(reltol = 1e-08, reltol.c = 1e-08, TMB = TRUE, optimizer = ifelse((num.RR + num.lv.c) == 0 | randomB != FALSE, "optim", "alabama"), max.iter = 6000, maxit = 6000, trace = FALSE, optim.method = NULL, nn.colMat = 10), control.va = list(Lambda.struc = "unstructured", Ab.struct = ifelse(is.null(colMat), "blockdiagonal", "MNunstructured"), Ab.struct.rank = 1, Ar.struc = "diagonal", diag.iter = 1, Ab.diag.iter = 0, Lambda.start = c(0.3, 0.3, 0.3), NN = 3), control.start = list(starting.val = "res", n.init = 1, n.init.max = 10, jitter.var = 0, jitter.var.br = 0, start.fit = NULL, start.lvs = NULL, randomX.start = "res", quad.start = 0.01, start.struc = "LV", scalmax = 10, MaternKappa = 1.5, rangeP = NULL, zetacutoff = NULL), setMap = NULL, ... )gllvm( y = NULL, X = NULL, TR = NULL, data = NULL, formula = NULL, family, num.lv = NULL, num.lv.c = 0, num.RR = 0, lv.formula = NULL, lvCor = NULL, studyDesign = NULL, dist = list(matrix(0)), distLV = matrix(0), colMat = NULL, colMat.rho.struct = "single", corWithin = FALSE, corWithinLV = FALSE, quadratic = FALSE, row.eff = FALSE, sd.errors = TRUE, offset = NULL, method = "VA", randomB = FALSE, randomX = NULL, beta0com = FALSE, zeta.struc = "species", plot = FALSE, link = "probit", Ntrials = 1, Power = 1.1, seed = NULL, scale.X = TRUE, return.terms = TRUE, gradient.check = FALSE, disp.formula = NULL, control = list(reltol = 1e-08, reltol.c = 1e-08, TMB = TRUE, optimizer = ifelse((num.RR + num.lv.c) == 0 | randomB != FALSE, "optim", "alabama"), max.iter = 6000, maxit = 6000, trace = FALSE, optim.method = NULL, nn.colMat = 10), control.va = list(Lambda.struc = "unstructured", Ab.struct = ifelse(is.null(colMat), "blockdiagonal", "MNunstructured"), Ab.struct.rank = 1, Ar.struc = "diagonal", diag.iter = 1, Ab.diag.iter = 0, Lambda.start = c(0.3, 0.3, 0.3), NN = 3), control.start = list(starting.val = "res", n.init = 1, n.init.max = 10, jitter.var = 0, jitter.var.br = 0, start.fit = NULL, start.lvs = NULL, randomX.start = "res", quad.start = 0.01, start.struc = "LV", scalmax = 10, MaternKappa = 1.5, rangeP = NULL, zetacutoff = NULL), setMap = NULL, ... )
y |
(n x m) matrix of responses. |
X |
matrix or data.frame of environmental covariates. |
TR |
matrix or data.frame of trait covariates. |
data |
data in long format, that is, matrix of responses, environmental and trait covariates and row index named as "id". When used, model needs to be defined using formula. This is alternative data input for y, X and TR. |
formula |
an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted (for fixed-effects predictors). |
family |
distribution function for responses. Options are |
num.lv |
number of latent variables, d, in gllvm model. Non-negative integer, less than number of response variables (m). Defaults to 2, if |
num.lv.c |
number of latent variables, d, in gllvm model to inform, i.e., with residual term. Non-negative integer, less than number of response (m) and equal to, or less than, the number of predictor variables (k). Defaults to 0. Requires specification of "lv.formula" in combination with "X" or "datayx". Can be used in combination with num.lv and fixed-effects, but not with traits. |
num.RR |
number of latent variables, d, in gllvm model to constrain, without residual term (reduced rank regression). Cannot yet be combined with traits. |
lv.formula |
an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted (for latent variables). |
lvCor |
correlation structure for latent variables, defaults to |
studyDesign |
variables related to eg. sampling/study design, used for defining correlation structure of the latent variables and row effects. |
dist |
list of length equal to the number of row effects with correlation structure |
distLV |
matrix of coordinates or time points used for LV correlation structure |
colMat |
either a list of length 2 with matrix of similarity for the column effects and named matrix "dist" of pairwise distances (of columns, to use in selecting nearest neighbours) for a sparse approximation of the matrix inverse in the likelihood, or only a (p.d.) matrix of similarity for the column effects for a normal inverse calculation. |
colMat.rho.struct |
either |
corWithin |
logical. Vector of length equal to the number of row effects. For structured row effects with correlation, If |
corWithinLV |
logical. For LVs with correlation, If |
quadratic |
either |
row.eff |
|
sd.errors |
logical. If |
offset |
vector or matrix of offset terms. |
method |
model can be fitted using Laplace approximation method ( |
randomB |
either |
randomX |
formula for species specific random effects of environmental variables in fourth corner model. Defaults to |
beta0com |
logical. If |
zeta.struc |
structure for cut-offs in the ordinal model. Either "common", for the same cut-offs for all species, or "species" for species-specific cut-offs. For the latter, classes are arbitrary per species, each category per species needs to have at least one observations. Defaults to "species". |
plot |
logical. If |
link |
link function for binomial family if |
Ntrials |
number of trials for binomial family. |
Power |
fixed power parameter in Tweedie model. Scalar from interval (1,2). Defaults to 1.1. If set to NULL it is estimated (note: experimental). |
seed |
a single seed value if |
scale.X |
logical. If |
return.terms |
logical. If |
gradient.check |
logical. If |
disp.formula |
a vector of indices, or alternatively formula, for the grouping of dispersion parameters (e.g. in a negative-binomial distribution, ZINB, tweedie), shape parameters (gamma, Beta, ordered Beta, hurdle Beta models) or variance parameters (gaussian distribution). Defaults to NULL so that all species have their own dispersion parameter. Is only allowed to include categorical variables. If a formula, data should be included as named rows in y. |
control |
A list with the following arguments controlling the optimization:
|
control.va |
A list with the following arguments controlling the variational approximation method:
|
control.start |
A list with the following arguments controlling the starting values:
|
setMap |
under development, not properly tested, except for ordinal beta cutoffs (zeta) and for rho_lvc. a list of a set of parameters to be fixed. Parameters to be fixed need to be defined with factors. Other arguments may overwrite these definitions. |
... |
Not used. |
Fits generalized linear latent variable models as in Hui et al. (2015 and 2017) and Niku et al. (2017).
Method can be used with two types of latent variable models depending on covariates. If only
site related environmental covariates are used, the expectation of response is determined by
where is a known link function, are -variate latent variables (<<), is an optional community level row effect
at site , and it can be fixed or random effect (also other structures are possible, see below), is an intercept term for species , and are column
specific coefficients related to covariates and the latent variables, respectively.
Alternatively, a more complex version of the model can be fitted with quadratic = TRUE, where species are modeled as a quadratic function of the latent variables:
. Here, D_j is a diagonal matrix of positive only quadratic coefficients, so that the model generates concave shapes only. This implementation follows the ecological theoretical model where species are generally recognized to exhibit non-linear response curves. For a model with quadratic responses, quadratic coefficients are assumed to be the same for all species:
. This model requires less information
per species and can be expected to be more applicable to most datasets. The quadratic coefficients D can be used to calculate the length of
ecological gradients.
For quadratic responses, it can be useful to provide the latent variables estimated with a GLLVM with linear responses, or estimated with (Detrended) Correspondence Analysis.
The latent variables can then be passed to the start.lvs argument inside the control.start list, which in many cases gives good results.
For GLLVMs with both linear and quadratic response model, a series of predictors can be included to explain the latent variables:
where are latent variables informed by the predictors, but not constrained compared to unconstrained ordination as in methods such as CCA or RDA.
Omitting the predictors results in an unconstrained ordination, and omitting in the usual constrained ordination, which can also be fitted.
An alternative model is the fourth corner model (Brown et al., 2014, Warton et al., 2015) which will be fitted if also trait covariates
are included. The expectation of response is
where g(.), , and are defined as above. Vectors and are the main effects
or coefficients related to environmental and trait covariates, respectively, matrix includes interaction terms. Vectors are
optional species-specific random slopes for environmental covariates.
The interaction/fourth corner terms are optional as well as are the main effects of trait covariates.
In addition to the sample-specific community level random effects, , it is also possible to set arbitrary structure/design for the row effects.
That is, assume that observations / rows in the data matrix are from groups , so that each row belongs to one of the groups, denote . Each group has a number of observations , so that .
Now we can set random intercept for each group , (see argument 'row.eff'):
There is also a possibility to set correlation structure for the random intercepts between groups, so that . That might be the case, for example, when the groups are spatially or temporally dependent.
Another option is to set row specific random intercepts , but to set the correlation structure for the observations within groups, (see argument 'corWithin'). That is, we can set according to some correlation function , when .
This model is restricted to the case, where each group has equal number of observations (rows), that is for all .
The correlation structures available in the package are
corAR1 autoregressive process of order 1.
corExp exponentially decaying, see argument 'dist'.
corCS compound symmetry.
The method is sensitive for the choices of initial values of the latent variables. Therefore it is
recommendable to use multiple runs and pick up the one giving the highest log-likelihood value (see argument 'n.init').
However, sometimes this is computationally too demanding, and default option
starting.val = "res" is recommended. For more details on different starting value methods, see Niku et al., (2018).
Models are implemented using TMB (Kristensen et al., 2015) applied to variational approximation (Hui et al., 2017), extended variational approximation (Korhonen et al., 2021) and Laplace approximation (Niku et al., 2017).
With ordinal family response classes must start from 0 or 1.
Mean and variance for distributions are defined as follows.
family = poisson():Expectation , variance , or
family = "negative.binomial": Expectation , variance , or
family = "ZIP": Expectation , variance .
family = "ZINB": Expectation , variance .
family = binomial(): Expectation , variance .
family = "beta": Expectation , variance .
family = "gamma":Expectation , variance , where is species specific shape parameter.
family = "exponential":Expectation , variance .
family = "tweedie"
Expectation , variance , where is a power parameter of Tweedie distribution. See details Dunn and Smyth (2005).
family = "ordinal":Cumulative probit model, see Hui et.al. (2016).
family = gaussian(): Expectation , variance
An object of class "gllvm" includes the following components:
call |
function call. |
y |
(n x m) matrix of responses. |
X |
matrix or data.frame of environmental covariates. |
X.design |
design matrix of environmental covariates. |
lv.X |
design matrix or data.frame of environmental covariates for latent variables. |
lv.X.design |
design matrix or data.frame of environmental covariates for latent variables. |
TR |
Trait matrix. |
formula |
Formula for predictors. |
lv.formula |
Formula of latent variables in constrained and concurrent ordination. |
randomX |
Formula for species specific random effects in fourth corner model. |
Xd |
design matrix for species specific random effects in fourth corner model. |
randomB |
Boolean flag for random slopes in constrained and concurrent ordination. |
num.lv |
Number of unconstrained latent variables. |
num.lv.c |
Number of latent variables in concurrent ordination. |
num.RR |
Number of latent variables in constrained ordination. |
Ntrials |
Number of trials in a binomial model. |
method |
Method used for integration. |
family |
Response distribution. |
row.eff |
Type of row effect used. |
n.init |
Number of model runs for best fit. |
disp.group |
Groups for dispersion parameters. |
sd |
List of standard errors. |
lvs |
Latent variables. |
params |
List of parameters
|
Power |
power parameter |
sd |
list of standard errors of parameters |
prediction.errors |
list of prediction covariances for latent variables and variances for random row effects when method |
A, Ar, Ab_lv, spArs
|
covariance matrices for variational densities of latent variables, random row effects, random slopes, and colum effects respectively |
seed |
Seed used for calculating starting values |
TMBfn |
TMB objective and derivative functions |
logL |
log likelihood |
convergence |
convergence code of optimizer |
quadratic |
flag for quadratic model |
Hess |
List holding matrices of second derivatives |
beta0com |
Flag for common intercept |
cstruc |
Correlation structure for row effects |
cstruclv |
Correlation structure for LVs |
dist |
Matrix of coordinates or time points used for row effects |
distLV |
Matrix of coordinates or time points used for LVs |
col.eff |
list of components for column random effects |
variational covariance structure of fitted model
fitted rank of variational covariance matrix
flag indicating if column random effects are included
design matrix
character vector for signal parameter
terms |
Terms object for main predictors |
start |
starting values for model |
optim.method |
Optimization method when using 'optim', 'alabama', or 'nloptr' |
If function gives warning: 'In f(x, order = 0) : value out of range in 'lgamma”, optimizer have visited an area where gradients become too big. It is automatically fixed by trying another step in the optimization process, and can be ignored if errors do not occur.
Jenni Niku <[email protected]>, Wesley Brooks, Riki Herliansyah, Francis K.C. Hui, Pekka Korhonen, Sara Taskinen, Bert van der Veen, David I. Warton
Brown, A. M., Warton, D. I., Andrew, N. R., Binns, M., Cassis, G., and Gibb, H. (2014). The fourth-corner solution - using predictive models to understand how species traits interact with the environment. Methods in Ecology and Evolution, 5:344-352.
Dunn, P. K. and Smyth, G. K. (2005). Series evaluation of tweedie exponential dispersion model densities. Statistics and Computing, 15:267-280.
Hui, F. K. C., Taskinen, S., Pledger, S., Foster, S. D., and Warton, D. I. (2015). Model-based approaches to unconstrained ordination. Methods in Ecology and Evolution, 6:399-411.
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.
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.
Korhonen, P., Hui, F. K. C., Niku, J., and Taskinen, S. (2021). Fast, universal estimation of latent variable models using extended variational approximations. Stat Comput 33, 26 (2023).
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.
Niku, J., Brooks, W., Herliansyah, R., Hui, F. K. C., Taskinen, S., and Warton, D. I. (2018). Efficient estimation of generalized linear latent variable models. PLoS One, 14(5):1-20.
Warton, D. I., Guillaume Blanchet, F., O'Hara, R. B., Ovaskainen, O., Taskinen, S., Walker, S. C. and Hui, F. K. C. (2015). So many variables: Joint modeling in community ecology. Trends in Ecology & Evolution, 30:766-779.
coefplot.gllvm, confint.gllvm, ordiplot.gllvm, plot.gllvm, summary.gllvm.
# Extract subset of the microbial data to be used as an example data(microbialdata) X <- microbialdata$Xenv y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, X, formula = ~ pH + Phosp, family = poisson()) fit$logL ordiplot(fit) coefplot(fit) # Inclusion of structured random row effect sDesign<-data.frame(Site = microbialdata$Xenv$Site) fit <- gllvm(y, X, formula = ~ pH + Phosp, family = poisson(), studyDesign=sDesign, row.eff=~(1|Site)) ## Load a dataset from the mvabund package library(mvabund) data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- as.matrix(antTraits$env) TR <- antTraits$traits # Fit model with environmental covariates Bare.ground and Shrub.cover fit <- gllvm(y, X, formula = ~ Bare.ground + Shrub.cover, family = poisson()) ordiplot(fit) coefplot.gllvm(fit) ## Example 1: Fit model with two unconstrained latent variables # Using variational approximation: fitv0 <- gllvm(y, family = "negative.binomial", method = "VA") ordiplot(fitv0) plot(fitv0, mfrow = c(2,2)) summary(fitv0) confint(fitv0) ## Example 1a: Fit concurrent ordination model with two latent variables and with # quadratic response model # We scale and centre the predictors to improve convergence fity1 <- gllvm(y, X = scale(X), family = "negative.binomial", num.lv.c=2, method="VA") ordiplot(fity1, biplot = TRUE) #'## Example 1b: Fit constrained ordination model with two latent variables and with # random canonical coefficients fity2 <- gllvm(y, X = scale(X), family = "negative.binomial", num.RR=2, randomB="LV", method="VA") # Using Laplace approximation: (this line may take about 30 sec to run) fitl0 <- gllvm(y, family = "negative.binomial", method = "LA") ordiplot(fitl0) # Poisson family: fit.p <- gllvm(y, family = poisson(), method = "LA") ordiplot(fit.p) # Use poisson model as a starting parameters for ZIP-model, this line # may take few minutes to run fit.z <- gllvm(y, family = "ZIP", method = "LA", control.start = list(start.fit = fit.p)) ordiplot(fit.z) ## Example 2: gllvm with environmental variables # Fit model with two latent variables and all environmental covariates, fitvX <- gllvm(formula = y ~ X, family = "negative.binomial") ordiplot(fitvX, biplot = TRUE) coefplot.gllvm(fitvX) # Fit model with environmental covariates Bare.ground and Shrub.cover fitvX2 <- gllvm(y, X, formula = ~ Bare.ground + Shrub.cover, family = "negative.binomial") ordiplot(fitvX2) coefplot.gllvm(fitvX2) # Use 5 initial runs and pick the best one fitvX_5 <- gllvm(y, X, formula = ~ Bare.ground + Shrub.cover, family = "negative.binomial", control.start=list(n.init = 5, jitter.var = 0.1)) ordiplot(fitvX_5) coefplot.gllvm(fitvX_5) ## Example 3: Data in long format # Reshape data to long format: datalong <- reshape(data.frame(cbind(y,X)), direction = "long", varying = colnames(y), v.names = "y") head(datalong) fitvLong <- gllvm(data = datalong, formula = y ~ Bare.ground + Shrub.cover, family = "negative.binomial") ## Example 4: Fourth corner model # Fit fourth corner model with two latent variables fitF1 <- gllvm(y = y, X = X, TR = TR, family = "negative.binomial") coefplot.gllvm(fitF1) # Fourth corner can be plotted also with next lines #fourth = fitF1$fourth.corner #library(lattice) #a = max( abs(fourth) ) #colort = colorRampPalette(c("blue","white","red")) #plot.4th = levelplot(t(as.matrix(fourth)), xlab = "Environmental Variables", # ylab = "Species traits", col.regions = colort(100), # at = seq( -a, a, length = 100), scales = list( x = list(rot = 45))) #print(plot.4th) # Specify model using formula fitF2 <- gllvm(y = y, X = X, TR = TR, formula = ~ Bare.ground + Canopy.cover * (Pilosity + Webers.length), family = "negative.binomial") ordiplot(fitF2) coefplot.gllvm(fitF2) ## Include species specific random slopes to the fourth corner model fitF3 <- gllvm(y = y, X = X, TR = TR, formula = ~ Bare.ground + Canopy.cover * (Pilosity + Webers.length), family = "negative.binomial", randomX = ~ Bare.ground + Canopy.cover, control.start = list(n.init = 3)) ordiplot(fitF3) coefplot.gllvm(fitF3) ## Example 5: Fit Tweedie model # Load coral data data(tikus) ycoral <- tikus$abund # Let's consider only years 1981 and 1983 ycoral <- ycoral[((tikus$x$time == 81) + (tikus$x$time == 83)) > 0, ] # Exclude species which have observed at less than 4 sites ycoral <- ycoral[-17, (colSums(ycoral > 0) > 4)] # Fit Tweedie model for coral data (this line may take few minutes to run) fit.twe <- gllvm(y = ycoral, family = "tweedie", method = "EVA", seed=111) fit.twe ## Example 6: Random row effects fitRand <- gllvm(y, family = "negative.binomial", row.eff = "random") ordiplot(fitRand, biplot = TRUE)# Extract subset of the microbial data to be used as an example data(microbialdata) X <- microbialdata$Xenv y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, X, formula = ~ pH + Phosp, family = poisson()) fit$logL ordiplot(fit) coefplot(fit) # Inclusion of structured random row effect sDesign<-data.frame(Site = microbialdata$Xenv$Site) fit <- gllvm(y, X, formula = ~ pH + Phosp, family = poisson(), studyDesign=sDesign, row.eff=~(1|Site)) ## Load a dataset from the mvabund package library(mvabund) data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- as.matrix(antTraits$env) TR <- antTraits$traits # Fit model with environmental covariates Bare.ground and Shrub.cover fit <- gllvm(y, X, formula = ~ Bare.ground + Shrub.cover, family = poisson()) ordiplot(fit) coefplot.gllvm(fit) ## Example 1: Fit model with two unconstrained latent variables # Using variational approximation: fitv0 <- gllvm(y, family = "negative.binomial", method = "VA") ordiplot(fitv0) plot(fitv0, mfrow = c(2,2)) summary(fitv0) confint(fitv0) ## Example 1a: Fit concurrent ordination model with two latent variables and with # quadratic response model # We scale and centre the predictors to improve convergence fity1 <- gllvm(y, X = scale(X), family = "negative.binomial", num.lv.c=2, method="VA") ordiplot(fity1, biplot = TRUE) #'## Example 1b: Fit constrained ordination model with two latent variables and with # random canonical coefficients fity2 <- gllvm(y, X = scale(X), family = "negative.binomial", num.RR=2, randomB="LV", method="VA") # Using Laplace approximation: (this line may take about 30 sec to run) fitl0 <- gllvm(y, family = "negative.binomial", method = "LA") ordiplot(fitl0) # Poisson family: fit.p <- gllvm(y, family = poisson(), method = "LA") ordiplot(fit.p) # Use poisson model as a starting parameters for ZIP-model, this line # may take few minutes to run fit.z <- gllvm(y, family = "ZIP", method = "LA", control.start = list(start.fit = fit.p)) ordiplot(fit.z) ## Example 2: gllvm with environmental variables # Fit model with two latent variables and all environmental covariates, fitvX <- gllvm(formula = y ~ X, family = "negative.binomial") ordiplot(fitvX, biplot = TRUE) coefplot.gllvm(fitvX) # Fit model with environmental covariates Bare.ground and Shrub.cover fitvX2 <- gllvm(y, X, formula = ~ Bare.ground + Shrub.cover, family = "negative.binomial") ordiplot(fitvX2) coefplot.gllvm(fitvX2) # Use 5 initial runs and pick the best one fitvX_5 <- gllvm(y, X, formula = ~ Bare.ground + Shrub.cover, family = "negative.binomial", control.start=list(n.init = 5, jitter.var = 0.1)) ordiplot(fitvX_5) coefplot.gllvm(fitvX_5) ## Example 3: Data in long format # Reshape data to long format: datalong <- reshape(data.frame(cbind(y,X)), direction = "long", varying = colnames(y), v.names = "y") head(datalong) fitvLong <- gllvm(data = datalong, formula = y ~ Bare.ground + Shrub.cover, family = "negative.binomial") ## Example 4: Fourth corner model # Fit fourth corner model with two latent variables fitF1 <- gllvm(y = y, X = X, TR = TR, family = "negative.binomial") coefplot.gllvm(fitF1) # Fourth corner can be plotted also with next lines #fourth = fitF1$fourth.corner #library(lattice) #a = max( abs(fourth) ) #colort = colorRampPalette(c("blue","white","red")) #plot.4th = levelplot(t(as.matrix(fourth)), xlab = "Environmental Variables", # ylab = "Species traits", col.regions = colort(100), # at = seq( -a, a, length = 100), scales = list( x = list(rot = 45))) #print(plot.4th) # Specify model using formula fitF2 <- gllvm(y = y, X = X, TR = TR, formula = ~ Bare.ground + Canopy.cover * (Pilosity + Webers.length), family = "negative.binomial") ordiplot(fitF2) coefplot.gllvm(fitF2) ## Include species specific random slopes to the fourth corner model fitF3 <- gllvm(y = y, X = X, TR = TR, formula = ~ Bare.ground + Canopy.cover * (Pilosity + Webers.length), family = "negative.binomial", randomX = ~ Bare.ground + Canopy.cover, control.start = list(n.init = 3)) ordiplot(fitF3) coefplot.gllvm(fitF3) ## Example 5: Fit Tweedie model # Load coral data data(tikus) ycoral <- tikus$abund # Let's consider only years 1981 and 1983 ycoral <- ycoral[((tikus$x$time == 81) + (tikus$x$time == 83)) > 0, ] # Exclude species which have observed at less than 4 sites ycoral <- ycoral[-17, (colSums(ycoral > 0) > 4)] # Fit Tweedie model for coral data (this line may take few minutes to run) fit.twe <- gllvm(y = ycoral, family = "tweedie", method = "EVA", seed=111) fit.twe ## Example 6: Random row effects fitRand <- gllvm(y, family = "negative.binomial", row.eff = "random") ordiplot(fitRand, biplot = TRUE)
These data describe the cover of sessile invertebrates, understory macroalgae, and bottom substrate types as determined by a uniform point contact method. The presence of over 150 taxa of sessile invertebrates and macroalgae are recorded at 80 uniformly spaced points along permanent 40m x 2m transects. Multiple species can be recorded at any given point. Percent cover of a given species on a transect can be estimated from UPC observations as the fraction of total points at which that species was present x 100. The total percent cover of all species combined using this method can exceed 100%; however, the percent cover of any single species cannot exceed 100%. This specific version of the data includes 61 species of macroalgae, 69 species of sessile invertebrates, and two species of plants.
data(kelpforest)data(kelpforest)
A data frame with the percent cover records of 132 sessile invertebrates and understory macroalgae measured at 836 permanent transects.
A data frame of study design variables and predictors.
Kelp forest site
Sampling year
Permanent transect identifying number (unique within site), nested within sites
A number of stipes of giant kelp
percent rock coverage
A data frame of the species information including species names, group and taxonomy.
Species codes corresponding species names in abundance matrix
Species group; algae, invertebrate or plant
Species' common and scientific names
Species taxonomic information
These data are part of SBC LTERs kelp forest monitoring program, which began in 2000 and was designed to track long-term patterns in species abundance and diversity of reef-associated organisms in the Santa Barbara Channel, California, USA. The sampling locations in this dataset include nine reef sites along the mainland coast of the Santa Barbara Channel and at two sites on the north side of Santa Cruz Island. These sites reflect several oceanographic regimes in the channel and vary in distance from sources of terrestrial runoff. Data collection began in 2000 and this dataset is updated annually.
The time period of data collection varied among the 11 kelp forest sites. Sampling at BULL, CARP, and NAPL began in 2000, sampling at the other 6 mainland sites (AHND, AQUE, IVEE, GOLB, ABUR, MOHK) began in 2001 (transects 3, 5, 6, 7, 8 at IVEE were added in 2011). Data collection at the two Santa Cruz Island sites (SCTW and SCDI) began in 2004.
Detailed description of the data available in the reference and the website https://sbclter.msi.ucsb.edu/data/catalog/package/?package=knb-lter-sbc.15
Reed, D, R. Miller. 2023. SBC LTER: Reef: Kelp Forest Community Dynamics: Cover of sessile organisms, Uniform Point Contact ver 33. Environmental Data Initiative. https://doi.org/10.6073/pasta/0af1a5b0d9dde5b4e5915c0012ccf99c. (Accessed: 2023-12-01).
## Not run: data(kelpforest) Y <- kelpforest$Y X <- kelpforest$X SPinfo <- kelpforest$SPinfo ## End(Not run)## Not run: data(kelpforest) Y <- kelpforest$Y X <- kelpforest$X SPinfo <- kelpforest$SPinfo ## End(Not run)
Extracts Log-likelihood from 'gllvm' objects.
## S3 method for class 'gllvm' logLik(object, ...)## S3 method for class 'gllvm' logLik(object, ...)
object |
an object of class 'gllvm'. |
... |
not used. |
David I. Warton, Jenni Niku
## Not run: ## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # log-Likelihood: logLik(fit) ## End(Not run)## Not run: ## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # log-Likelihood: logLik(fit) ## End(Not run)
Microbial community data consist of abundances of 985 bacteria species measured at 56 soil sample sites from three regions, Kilpisjarvi (Finland), Ny-Alesund (Norway), and Mayrhofen (Austria). In addition to bacteria counts, three continuous environmental variables (pH, available phosphorous and soil organic matter) were measured from each soil sample.
data(microbialdata)data(microbialdata)
A data frame with abundances of 985 bacteria species measured at 56 soil sample sites
Environmental variables SOM: soil organic matter, pH: soil pH value and Phosp: available phosphorus and information from the samples, including Region: sampling region (Kilpisjarvi (Finland), Ny-Alesund (Norway), and Mayrhofen (Austria).), Site: sampling site and Soiltype: soil sample type (top soil (T) or bottom soil (B))
Kumar, M., Brader, G., Sessitsch, A., Mäki, A., van Elsas, J.D., and Nissinen, R. (2017). Plants Assemble Species Specific Bacterial Communities from Common Core Taxa in Three Arcto-Alpine Climate Zones. Frontiers in Microbiology, 8:12.
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.
Plots latent variables and their corresponding coefficients (biplot).
## S3 method for class 'gllvm' ordiplot( object, biplot = FALSE, ind.spp = NULL, alpha = 0.5, main = NULL, which.lvs = c(1, 2), predict.region = FALSE, level = 0.95, jitter = FALSE, jitter.amount = 0.2, s.colors = 1, s.cex = 1.2, symbols = FALSE, cex.spp = 0.7, spp.colors = "blue", arrow.scale = 0.8, arrow.spp.scale = 0.8, arrow.ci = TRUE, arrow.lty = "solid", fac.center = FALSE, spp.arrows = NULL, spp.arrows.lty = "dashed", cex.env = 0.7, lab.dist = 0.1, lwd.ellips = 0.5, col.ellips = 4, lty.ellips = 1, type = NULL, rotate = TRUE, ... )## S3 method for class 'gllvm' ordiplot( object, biplot = FALSE, ind.spp = NULL, alpha = 0.5, main = NULL, which.lvs = c(1, 2), predict.region = FALSE, level = 0.95, jitter = FALSE, jitter.amount = 0.2, s.colors = 1, s.cex = 1.2, symbols = FALSE, cex.spp = 0.7, spp.colors = "blue", arrow.scale = 0.8, arrow.spp.scale = 0.8, arrow.ci = TRUE, arrow.lty = "solid", fac.center = FALSE, spp.arrows = NULL, spp.arrows.lty = "dashed", cex.env = 0.7, lab.dist = 0.1, lwd.ellips = 0.5, col.ellips = 4, lty.ellips = 1, type = NULL, rotate = TRUE, ... )
object |
an object of class 'gllvm'. |
biplot |
|
ind.spp |
the number of response variables (usually, species) to include on the biplot. The default is none, or all if |
alpha |
a numeric scalar between 0 and 1 that is used to control the relative scaling of the latent variables and their coefficients, when constructing a biplot. |
main |
main title. |
which.lvs |
indices of two latent variables to be plotted if number of the latent variables is more than 2. A vector with length of two. Defaults to |
predict.region |
if |
level |
level for prediction regions. |
jitter |
if |
jitter.amount |
numeric, positive value indicating an amount of jittering for each point, defaults to 0.2 (jitter range). |
s.colors |
colors for sites |
s.cex |
size of site labels |
symbols |
logical, if |
cex.spp |
size of species labels in biplot |
spp.colors |
colors for sites, defaults to |
arrow.scale |
positive value, to scale arrows |
arrow.spp.scale |
positive value, to scale arrows of species |
arrow.ci |
represent statistical uncertainty for arrows in constrained or concurrent ordination using confidence or prediction interval? Defaults to |
arrow.lty |
linetype for arrows in constrained |
fac.center |
logical. If |
spp.arrows |
plot species scores as arrows if outside of the range of the plot? Defaults to |
spp.arrows.lty |
linetype for species arrows |
cex.env |
size of labels for arrows in constrained ordination |
lab.dist |
distance between label and arrow heads. Value between 0 and 1 |
lwd.ellips |
line width for prediction ellipses. See graphical parameter lwd. |
col.ellips |
colors for prediction ellipses. |
lty.ellips |
line type for prediction ellipses. See graphical parameter lty. |
type |
which type of ordination plot to construct. Options are "residual", "conditional", and "marginal". Defaults to "residual" for GLLVMs with unconstrained latent variables and "conditional" otherwise. |
rotate |
logical, if |
... |
additional graphical arguments. |
Function constructs a scatter plot of two latent variables, i.e. an ordination plot. Latent variables are re-rotated to their principal direction using singular value decomposition, so that the first plotted latent variable does not have to be the first latent variable in the model. If only one latent variable is in the fitted model, latent variables are plotted against their corresponding row indices. The latent variables are labeled using the row index of the response matrix y.
Coefficients related to latent variables are plotted in the same figure with the latent
variables if biplot = TRUE. They are labeled using the column names of y. The number
of latent variable coefficients to be plotted can be controlled by ind.spp. An argument alpha
is used to control the relative scaling of the latent variables and their coefficients.
If alpha = 0.5, the latent variables and their coefficients are on the same scale.
For details for constructing a biplot, see Gabriel (1971).
For a quadratic response model, species optima are plotted. Any species scores that are outside the range of the predicted site scores are not directly plotted, but their main direction is indicated with arrows instead. This ensures that the plot remains on a reasonable scale.
Effects of environmental variables in constrained ordination are indicated with arrows. If any of the arrows exceeds the range of the plot, arrows are scaled to 80 but so that the relative contribution of predictors is maintained. If standard errors are available in the provided model, the slopes of environmental variables for which the 95 are slightly less intensely coloured.
For constrained ordination, a conditional plot includes both fixed- and random-effects to
optimally represent species co-occurrence patterns, corresponding to "conditional" site scores in getLV.gllvm.
Marginal corresponds to an ordination plot that excludes residual patterns (i.e. excluding the random-effect),
so that it is only available with num.lv.c>0 or num.RR>0. A conditional plot requires num.lv.c>0.
The "residual" type corresponds to an ordination diagram of only residual patterns.
See getLV.gllvm for details.
- If error is occurred when using ordiplot(), try full name of the function ordiplot.gllvm() as functions named 'ordiplot' might be found in other packages as well.
Jenni Niku <[email protected]>, Francis K.C. Hui, Bert van der Veen
Gabriel, K. R. (1971). The biplot graphic display of matrices with application to principal component analysis. Biometrika, 58, 453-467.
#'# Extract subset of the microbial data to be used as an example data(microbialdata) y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, family = poisson()) fit$logL ordiplot(fit, predict.region = TRUE) ## Not run: #'## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) fit <- gllvm(y, family = poisson()) # Ordination plot: ordiplot(fit) # Biplot with 10 species ordiplot(fit, biplot = TRUE, ind.spp = 10) ## End(Not run)#'# Extract subset of the microbial data to be used as an example data(microbialdata) y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, family = poisson()) fit$logL ordiplot(fit, predict.region = TRUE) ## Not run: #'## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) fit <- gllvm(y, family = poisson()) # Ordination plot: ordiplot(fit) # Biplot with 10 species ordiplot(fit, biplot = TRUE, ind.spp = 10) ## End(Not run)
Plots phylogenetic random effects with the phylogeny, and community effects
## S3 method for class 'gllvm' phyloplot( object, tree, comm.eff = TRUE, row.eff = FALSE, which.Xcoef = NULL, xlim = NULL, level = 0.95, col = c("#E69F00", "white", "#009E73"), col.sym = TRUE, mar.spec = c(3, 2, 0, 0), mar.phy = c(0, 2, 2, 0), mar.comm = c(3, 0.5, 2, 1.5), cex = 0.6, lwd = 1, col.edge = "black", pch = "x", heights = c(0.55, 0.35), widths = c(0.64, 0.1), phy.place = "top", ... )## S3 method for class 'gllvm' phyloplot( object, tree, comm.eff = TRUE, row.eff = FALSE, which.Xcoef = NULL, xlim = NULL, level = 0.95, col = c("#E69F00", "white", "#009E73"), col.sym = TRUE, mar.spec = c(3, 2, 0, 0), mar.phy = c(0, 2, 2, 0), mar.comm = c(3, 0.5, 2, 1.5), cex = 0.6, lwd = 1, col.edge = "black", pch = "x", heights = c(0.55, 0.35), widths = c(0.64, 0.1), phy.place = "top", ... )
object |
an object of class 'gllvm'. |
tree |
an object of class ” |
comm.eff |
logical, defaults to |
row.eff |
logical, defaults to |
which.Xcoef |
List of length 2 with names to subset the effects to plot. The first vector is for the species plot, the second for community effects. |
xlim |
vector of length two. Limits for the x-axis of the caterpillar plot. Defaults to NULL, in which case the limits are chosen based on the confidence intervals. |
level |
the confidence level. Scalar between 0 and 1. |
col |
vector of three colors (defaults to |
col.sym |
logical, defaults to |
mar.spec |
vector of length 4, which defines the margins sizes for the species random effects plot. Defaults to |
mar.phy |
vector of length 4, which defines the margins sizes for plotting the phylogeny. Defaults to |
mar.comm |
vector of length 4, which defines the margins sizes for the caterpillar plot. Defaults to |
cex |
the magnification to be used for text in the plot. Defaults to 0.6. |
lwd |
line thickness for the branches in the phylogeny and the confidence intervals in the caterpillar plot. Defaults to 1. |
col.edge |
character. Color of branches in the phylogeny. |
pch |
symbol used in the catter pillar plot. Defaults to "x". |
heights |
vector of length two. Relative row heights, defaults to |
widths |
vector of length two. Relative column widths, defaults to |
phy.place |
not (yet) in use. |
... |
additional not in use. |
Plots phylogenetically structured random effects together with the phylogeny, and with community-level effects (i.e., effects that are the same across species). If standard errors have been calculated for the model, the prediction intervals for species random effects are checked, and crossed out (i.e., displayed as white) if they cross zero.
Bert van der Veen
van der Veen, B., O'Hara, R.B. (2024). Fast fitting of Fast fitting of phylogenetic mixed effects models. arXiv.
## Not run: # Load dataset data(fungi) Y <- fungi$Y # Scale the predictor X <- fungi$X X[,"DBH.CM"] <- scale(X[, "DBH.CM"]) tree <- fungi$tree # the tree colMat <- fungi$C # e.g., from ape::vcv(tree) dist <- fungi$dist # e.g., from ape::cophenetic.phylo(tree) order <- gllvm:::findOrder(covMat = colMat, distMat = dist, nn = 15, order = order(dist[1:length(tree$tip.label), nrow(dist)], decreasing = TRUE))$order order <- tree$tip.label[order] model <- gllvm(y = Y[,order], X = X, formula = ~(DBH.CM|1), beta0com = TRUE, family = "binomial", num.lv = 0, nn.colMat = 15, colMat = list(colMat[order,order], dist = dist[order,order]), colMat.rho.struct = "term") phyloplot(model, tree) ## End(Not run)## Not run: # Load dataset data(fungi) Y <- fungi$Y # Scale the predictor X <- fungi$X X[,"DBH.CM"] <- scale(X[, "DBH.CM"]) tree <- fungi$tree # the tree colMat <- fungi$C # e.g., from ape::vcv(tree) dist <- fungi$dist # e.g., from ape::cophenetic.phylo(tree) order <- gllvm:::findOrder(covMat = colMat, distMat = dist, nn = 15, order = order(dist[1:length(tree$tip.label), nrow(dist)], decreasing = TRUE))$order order <- tree$tip.label[order] model <- gllvm(y = Y[,order], X = X, formula = ~(DBH.CM|1), beta0com = TRUE, family = "binomial", num.lv = 0, nn.colMat = 15, colMat = list(colMat[order,order], dist = dist[order,order]), colMat.rho.struct = "term") phyloplot(model, tree) ## End(Not run)
Five plots (selectable by which) are currently available: a plot of residuals against linear predictors of fitted values, a Normal Q-Q plot of residuals with a simulated point-wise 95% confidence interval envelope, residuals against row index and column index and scale location plot.
## S3 method for class 'gllvm' plot( x, which = 1:5, caption = c("Residuals vs linear predictors", "Normal Q-Q", "Residuals vs row", "Residuals vs column", "Scale-Location"), var.colors = NULL, add.smooth = TRUE, envelopes = TRUE, reps = 150, envelope.col = c("blue", "lightblue"), n.plot = NULL, ... )## S3 method for class 'gllvm' plot( x, which = 1:5, caption = c("Residuals vs linear predictors", "Normal Q-Q", "Residuals vs row", "Residuals vs column", "Scale-Location"), var.colors = NULL, add.smooth = TRUE, envelopes = TRUE, reps = 150, envelope.col = c("blue", "lightblue"), n.plot = NULL, ... )
x |
an object of class 'gllvm'. |
which |
if a subset of the plots is required, specify a subset of the numbers 1:5, see caption below. |
caption |
captions to appear above the plots. |
var.colors |
colors for responses, vector with length of number of response variables or 1. Defaults to NULL, when different responses have different colors. |
add.smooth |
logical with default |
envelopes |
logical, indicating if simulated point-wise confidence interval envelope will be added to Q-Q plot, defaults to |
reps |
number of replications when simulating confidence envelopes for normal Q-Q plot |
envelope.col |
colors for envelopes, vector with length of two |
n.plot |
number of species (response variables) to be plotted. Defaults to |
... |
additional graphical arguments. |
plot.gllvm is used for model diagnostics. Dunn-Smyth residuals (randomized quantile residuals) (Dunn and Smyth, 1996) are used in plots. Colors indicate different species.
Jenni Niku <[email protected]>
Dunn, P. K., and Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5, 236-244.
Hui, F. K. C., Taskinen, S., Pledger, S., Foster, S. D., and Warton, D. I. (2015). Model-based approaches to unconstrained ordination. Methods in Ecology and Evolution, 6:399-411.
## Not run: # Fit gllvm model with Poisson family data(microbialdata) X <- microbialdata$Xenv y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, X, formula = ~ pH + Phosp, family = poisson()) # Plot residuals plot(fit, mfrow = c(3,2)) ## End(Not run)## Not run: # Fit gllvm model with Poisson family data(microbialdata) X <- microbialdata$Xenv y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, X, formula = ~ pH + Phosp, family = poisson()) # Plot residuals plot(fit, mfrow = c(3,2)) ## End(Not run)
Obtains predictions from a fitted generalized linear latent variable model object.
## S3 method for class 'gllvm' predict( object, newX = NULL, newTR = NULL, newLV = NULL, type = "link", level = 1, offset = TRUE, ... )## S3 method for class 'gllvm' predict( object, newX = NULL, newTR = NULL, newLV = NULL, type = "link", level = 1, offset = TRUE, ... )
object |
an object of class 'gllvm'. |
newX |
A new data frame of environmental variables. If omitted, the original matrix of environmental variables is used. |
newTR |
A new data frame of traits for each response taxon. If omitted, the original matrix of traits is used. |
newLV |
A new matrix of latent variables. If omitted, the original matrix of latent variables is used. Note that number of rows/sites must be the same for |
type |
the type of prediction required. The default ( |
level |
specification for how to predict. Level one ( |
offset |
specification whether of not offset values are included to the predictions in case they are in the model, defaults to |
... |
not used. |
If newX, newTR and newLV are omitted the predictions are based on the data used for fitting the model. Notice that newTR need to match with the number of species in the original data.
Instead, new sites can be specified in newX. If predictors newX (and newTR) are given, and newLV is not, latent variables are not used in the predictions.
A matrix containing requested predictor types.
Jenni Niku <[email protected]>, David Warton
# Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- scale(antTraits$env[, 1:3]) # Fit gllvm model fit <- gllvm(y = y, X, family = poisson()) # fitted values predfit <- predict(fit, type = "response") # linear predictors predlin <- predict(fit) # Predict new sites: # Generate matrix of environmental variables for 10 new sites xnew <- cbind(rnorm(10), rnorm(10), rnorm(10)) colnames(xnew) <- colnames(X) predfit <- predict(fit, newX = xnew, type = "response", level = 0) TR <- (antTraits$tr[, 1:3]) fitt <- gllvm(y = y, X, TR, family = poisson()) # linear predictors predlin <- predict(fitt) # Predict new sites: # Generate matrix of environmental variables for 10 new sites xnew <- cbind(rnorm(10), rnorm(10), rnorm(10)) colnames(xnew) <- colnames(X) # Generate matrix of traits for species trnew <- data.frame(Femur.length = rnorm(41), No.spines = rnorm(41), Pilosity = factor(sample(0:3, 41, replace = TRUE))) predfit <- predict(fitt, newX = xnew, newTR = trnew, type = "response", level = 0)# Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- scale(antTraits$env[, 1:3]) # Fit gllvm model fit <- gllvm(y = y, X, family = poisson()) # fitted values predfit <- predict(fit, type = "response") # linear predictors predlin <- predict(fit) # Predict new sites: # Generate matrix of environmental variables for 10 new sites xnew <- cbind(rnorm(10), rnorm(10), rnorm(10)) colnames(xnew) <- colnames(X) predfit <- predict(fit, newX = xnew, type = "response", level = 0) TR <- (antTraits$tr[, 1:3]) fitt <- gllvm(y = y, X, TR, family = poisson()) # linear predictors predlin <- predict(fitt) # Predict new sites: # Generate matrix of environmental variables for 10 new sites xnew <- cbind(rnorm(10), rnorm(10), rnorm(10)) colnames(xnew) <- colnames(X) # Generate matrix of traits for species trnew <- data.frame(Femur.length = rnorm(41), No.spines = rnorm(41), Pilosity = factor(sample(0:3, 41, replace = TRUE))) predfit <- predict(fitt, newX = xnew, newTR = trnew, type = "response", level = 0)
Obtains predictions for latent variables from a fitted generalized linear latent variable model object. Currently works only for the variational approximation method.
## S3 method for class 'gllvm' predictLVs(object, newX = NULL, newY = object$y, ...)## S3 method for class 'gllvm' predictLVs(object, newX = NULL, newY = object$y, ...)
object |
an object of class 'gllvm'. |
newX |
A new data frame of environmental variables. If omitted, the original matrix of environmental variables is used. |
newY |
A new response data. Defaults to the dataset used for original model fit. |
... |
not used. |
Obtains predictions for latent variables from a fitted generalized linear latent variable model object.
A matrix containing requested predictor types.
David Warton, Jenni Niku <[email protected]>
# Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- scale(antTraits$env[, 1:3]) # Fit gllvm model fit <- gllvm(y = y, X, family = poisson()) # fitted values predLVs <- predictLVs.gllvm(fit)# Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- scale(antTraits$env[, 1:3]) # Fit gllvm model fit <- gllvm(y = y, X, family = poisson()) # fitted values predLVs <- predictLVs.gllvm(fit)
Plots random slopes and their prediction intervals.
## S3 method for class 'gllvm' randomCoefplot( object, y.label = TRUE, which.Xcoef = NULL, cex.ylab = 0.5, mfrow = NULL, mar = c(4, 6, 2, 1), xlim.list = NULL, order = FALSE, ... )## S3 method for class 'gllvm' randomCoefplot( object, y.label = TRUE, which.Xcoef = NULL, cex.ylab = 0.5, mfrow = NULL, mar = c(4, 6, 2, 1), xlim.list = NULL, order = FALSE, ... )
object |
an object of class 'gllvm'. |
y.label |
logical, if |
which.Xcoef |
fector indicating which covariate coefficients will be plotted. Can be vector of covariate names or numbers. Default is NULL when all covariate coefficients are plotted. |
cex.ylab |
the magnification to be used for axis annotation relative to the current setting of cex. |
mfrow |
same as |
mar |
vector of length 4, which defines the margin sizes: |
xlim.list |
list of vectors with length of two to define the intervals for x axis in each covariate plot. Defaults to NULL when the interval is defined by the range of point estimates and confidence intervals |
order |
logical, if |
... |
additional graphical arguments. |
Jenni Niku <[email protected]>, Francis K.C. Hui, Bert van der Veen, Sara Taskinen,
## Not run: ## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- as.matrix(antTraits$env) TR <- antTraits$traits # Fit model with random slopes fitF <- gllvm(y = y, X = X, TR = TR, formula = ~ Bare.ground + Bare.ground : Webers.length, family = poisson(), randomX = ~ Bare.ground) randomCoefplot(fitF) ## End(Not run)## Not run: ## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- as.matrix(antTraits$env) TR <- antTraits$traits # Fit model with random slopes fitF <- gllvm(y = y, X = X, TR = TR, formula = ~ Bare.ground + Bare.ground : Webers.length, family = poisson(), randomX = ~ Bare.ground) randomCoefplot(fitF) ## End(Not run)
Calculates Dunn-Smyth residuals for gllvm model.
## S3 method for class 'gllvm' residuals(object, ...)## S3 method for class 'gllvm' residuals(object, ...)
object |
an object of class 'gllvm'. |
... |
not used. |
Computes Dunn-Smyth residuals (randomized quantile residuals, Dunn and Smyth, 1996) for gllvm model.
For the observation Dunn-Smyth residuals are defined as
where and are the cumulative probability functions of the standard normal
distribution, is the limit as is approached from the negative side, and has been
generated at random from the standard uniform distribution.
residuals |
matrix of residuals |
linpred |
matrix of linear predictors |
Jenni Niku <[email protected]>
Dunn, P. K., and Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5, 236-244.
Hui, F. K. C., Taskinen, S., Pledger, S., Foster, S. D., and Warton, D. I. (2015). Model-based approaches to unconstrained ordination. Methods in Ecology and Evolution, 6:399-411.
## Not run: # Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # residuals res <- residuals(fit) ## End(Not run)## Not run: # Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) # residuals res <- residuals(fit) ## End(Not run)
Calculates Hessian and standard errors for gllvm model.
## S3 method for class 'gllvm' se(object, ...)## S3 method for class 'gllvm' se(object, ...)
object |
an object of class 'gllvm'. |
... |
not used. |
Computes Hessian and standard errors for gllvm model.
sd |
list of standard errors of parameters |
Hess |
list including Hessian matrix and approximative covariance matrix of parameters |
Jenni Niku <[email protected]>
Dunn, P. K., and Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5, 236-244.
Hui, F. K. C., Taskinen, S., Pledger, S., Foster, S. D., and Warton, D. I. (2015). Model-based approaches to unconstrained ordination. Methods in Ecology and Evolution, 6:399-411.
data(eSpider) mod <- gllvm(eSpider$abund, num.lv = 2, family = "poisson", sd.errors = FALSE) # Calculate standard errors after fitting sdErr <- se(mod) # Store the standard errors in the right place mod$sd <-sdErr$sd # Store the Hessian in the right place mod$Hess <- sdErr$Hessdata(eSpider) mod <- gllvm(eSpider$abund, num.lv = 2, family = "poisson", sd.errors = FALSE) # Calculate standard errors after fitting sdErr <- se(mod) # Store the standard errors in the right place mod$sd <-sdErr$sd # Store the Hessian in the right place mod$Hess <- sdErr$Hess
Generate new data using the fitted values of the parameters
## S3 method for class 'gllvm' simulate(object, nsim = 1, seed = NULL, conditional = FALSE, ...)## S3 method for class 'gllvm' simulate(object, nsim = 1, seed = NULL, conditional = FALSE, ...)
object |
an object of class 'gllvm'. |
nsim |
an optional positive integer specifying the number of simulated datasets. Defaults to 1. |
seed |
an optional integer to set seed number, passed to set.seed. Defaults to a random seed number. |
conditional |
if |
... |
not used. |
simulate function for gllvm objects.
A matrix containing generated data.
David Warton, Jenni Niku <[email protected]>
# Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- scale(antTraits$env[, 1:3]) # Fit gllvm model fit <- gllvm(y = y, X, family = poisson()) # Simulate data newdata <- simulate(fit)# Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) X <- scale(antTraits$env[, 1:3]) # Fit gllvm model fit <- gllvm(y = y, X, family = poisson()) # Simulate data newdata <- simulate(fit)
Dataset of ordinal observations of plants, on the island Skabbholmen in the Stocholm archipelago. Includes 65 unique sites and 70 species, surveyed in two different years.
data(Skabbholmen)data(Skabbholmen)
A data frame with ordinal of 70 plant species measured at 126 plots.
A matrix of 2 predictor variables at 126 plots.
A matrix of full species names and abbreviations used in the community data (Y).
Observations of vascular plant cover in 126 one-square-meter plots divided over four transects. The ordinal responses are on a five-degree Hult-Sernander-Du Rietz scale, and were originally recorded by Wolfgang and Cramer (1987) and additionally analyzed by ter Braak (1987). There is a total of 64 unique sites, that were surveyed in two different years (1978 and 1984), but two plots were only surveyed in one year (thus bringing the total number of rows in the data to 126). The plots were located on an elevation gradient, running from the shoreline to the edge of old-growth forest. Elevation to the shoreline was recorded in centimeters during the sampling in 1978.
This dataset was published with permission from the CANOCO FORTRAN package example datasets.
ter Braak, C.J.F. and Smilauer, P. (1998). CANOCO reference manual and user's guide to CANOCO for Windows: software for canonical community ordination (version 4). Microcomputer Power, New York, New York, USA.
Jongman, E., & Jongman, S. R. R. (1995). Data analysis in community and landscape ecology. Cambridge university press.
ter Braak, C.J.F. (1987). The analysis of vegetation-environment relationships by canonical correspondence analysis. Vegetatio, 69(1), 69-77.
Cramer, W. & Hytteborn, H. (1987). The separation of fluctuation and long-term change in vegetation dynamics of a rising seashore. Vegetatio, 69, 157–167.
# Uncomment the example #data(Skabbholmen) #Y <- Skabbholmen$Y #X <- Skabbholmen$X #model <- gllvm(y = Y, X = X, # num.RR = 2, # family = "ordinal", # zeta.struc="common", # row.eff=~(1|transectID))# Uncomment the example #data(Skabbholmen) #Y <- Skabbholmen$Y #X <- Skabbholmen$X #model <- gllvm(y = Y, X = X, # num.RR = 2, # family = "ordinal", # zeta.struc="common", # row.eff=~(1|transectID))
A summary of the fitted 'gllvm' object, including function call, distribution family and model parameters.
## S3 method for class 'gllvm' summary( object, by = "all", digits = max(3L, getOption("digits") - 3L), signif.stars = getOption("show.signif.stars"), dispersion = FALSE, spp.intercepts = FALSE, row.intercepts = FALSE, Lvcoefs = FALSE, rotate = TRUE, type = NULL, ... ) ## S3 method for class 'summary.gllvm' print(x, ...) ## S3 method for class 'summary.gllvm' plot(x, component = NULL, ...)## S3 method for class 'gllvm' summary( object, by = "all", digits = max(3L, getOption("digits") - 3L), signif.stars = getOption("show.signif.stars"), dispersion = FALSE, spp.intercepts = FALSE, row.intercepts = FALSE, Lvcoefs = FALSE, rotate = TRUE, type = NULL, ... ) ## S3 method for class 'summary.gllvm' print(x, ...) ## S3 method for class 'summary.gllvm' plot(x, component = NULL, ...)
object |
an object of class 'gllvm' |
by |
By = "all" (default) will return a Wald statistics per predictor and LV if the ordination includes predictors, by = "terms" will return a multivariate Wald statistic per predictor (displayed at first LV), and by = "LV" will do the same but per dimension (displayed at first predictors). |
digits |
the number of significant digits to use when printing |
signif.stars |
If |
dispersion |
option to return dispersion parameters, defaults to |
spp.intercepts |
option to return species intercepts, defaults to |
row.intercepts |
option to return row intercepts, defaults to |
Lvcoefs |
option to return species scores in the ordination, defaults to |
rotate |
defaults to |
type |
to match "type" in |
... |
not used. |
x |
a summary object |
component |
component to be plotted |
Various options are available to include extra parameter estimates in the summary, which have been excluded by default, for readability.
Jenni Niku <[email protected]>, Bert van der Veen
## Not run: ## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) summary(fit) ## End(Not run)## Not run: ## Load a dataset from the mvabund package data(antTraits, package = "mvabund") y <- as.matrix(antTraits$abund) # Fit gllvm model fit <- gllvm(y = y, family = poisson()) summary(fit) ## End(Not run)
Calculates variance partitioning for gllvm object with function varPartitioning().
Function plotVarPartitioning() (alias plotVP()) plots the results of variance partitioning of a fitted gllvm.
## S3 method for class 'gllvm' varPartitioning( object, group = NULL, groupnames = NULL, adj.cov = TRUE, grouplvs = FALSE, ... ) plotVarPartitioning( VP, main = "Variance Partitioning", xlab = "Response", ylab = "Variance proportion", legend.text = NULL, args.legend = list(cex = 0.7, x = "topright", bty = "n", inset = c(0, -0.15)), mar = c(4, 4, 6, 2), ... ) plotVP(VP, ...)## S3 method for class 'gllvm' varPartitioning( object, group = NULL, groupnames = NULL, adj.cov = TRUE, grouplvs = FALSE, ... ) plotVarPartitioning( VP, main = "Variance Partitioning", xlab = "Response", ylab = "Variance proportion", legend.text = NULL, args.legend = list(cex = 0.7, x = "topright", bty = "n", inset = c(0, -0.15)), mar = c(4, 4, 6, 2), ... ) plotVP(VP, ...)
object |
an object of class 'gllvm'. |
group |
a vector of integers identifying grouping of X covariates, the default is to use model terms formula and lv.formula. |
groupnames |
a vector of strings given as names for the groups defined in group |
adj.cov |
logical, whether or not to adjust co-variation within the group |
grouplvs |
logical, whether or not to group latent variables to one group |
... |
additional graphical arguments passed to the barplot function |
VP |
a variance partitioning object for a gllvm produced by function varPartitioning. |
main |
main title |
xlab |
a label for the x axis. |
ylab |
a label for the y axis. |
legend.text |
a vector of names for the groups, as a default 'groupnames' from varPartitioning. If FALSE, legend not printed. |
args.legend |
a list of additional arguments to pass to |
mar |
Margins of the plot. Default |
Variance for the linear predictor for response j can be calculated as
where is a vector consisting of predictor/latent variable/row effect etc values for all sampling units i.
If s are not correlated, covariance term is 0 and thus the variance explained of a response j for predictor is given as .
In case of correlated predictors, it is advised to group them into a same group. The variance explained is calculated for the correlated group of predictors together and adjusted with the covariance term.
Jenni Niku <[email protected]>
# Extract subset of the microbial data to be used as an example data(microbialdata) X <- microbialdata$Xenv y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, X[,1:3], formula = ~ pH + Phosp, family = poisson(), studyDesign = X[,4:5], row.eff = ~(1|Site)) VP <- varPartitioning(fit) plotVarPartitioning(VP) ## Not run: # Plot the result of variance partitioning plotVP(VP, col = palette(hcl.colors(5, "Roma"))) ## End(Not run)# Extract subset of the microbial data to be used as an example data(microbialdata) X <- microbialdata$Xenv y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]] fit <- gllvm(y, X[,1:3], formula = ~ pH + Phosp, family = poisson(), studyDesign = X[,4:5], row.eff = ~(1|Site)) VP <- varPartitioning(fit) plotVarPartitioning(VP) ## Not run: # Plot the result of variance partitioning plotVP(VP, col = palette(hcl.colors(5, "Roma"))) ## End(Not run)
Returns the variance-covariance matrix of the parameters from a GLLVM. If the variance-covariance matrix was not calculated after model fitting, this function will have to calculate the variance-covariance matrix, which may be computational intensive for a large number of species and/or sites.
## S3 method for class 'gllvm' vcov(object, ...)## S3 method for class 'gllvm' vcov(object, ...)
object |
an object of class 'gllvm'. |
... |
not used. |
Calculates the variance-covariance matrix of a GLLVM object using se.gllvm, which may be computational intensive with many parameters.The parameters might have unintuitive names. Fixed-effects coefficients are labeled "b", and are ordered per species as: 1) intercepts 2) fixed-effects slopes. Coefficients of the latent variables are labled "lambda" (linear coefficients) or "lambda2".
Bert van der Veen