Title:	Hierarchical Model of Species Communities
Type:	Package
Version:	3.0-13
Description:	Hierarchical Modelling of Species Communities (HMSC) is a model-based approach for analyzing community ecological data. This package implements it in the Bayesian framework with Gibbs Markov chain Monte Carlo (MCMC) sampling (Tikhonov et al. (2020) <doi:10.1111/2041-210X.13345>).
License:	GPL-3 | file LICENSE
URL:	https://www.helsinki.fi/en/researchgroups/statistical-ecology/software/hmsc
BugReports:	https://github.com/hmsc-r/HMSC/issues/
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.2.1
Suggests:	R.rsp, testthat, corrplot
VignetteBuilder:	R.rsp
Depends:	R (≥ 3.0.2), coda
Imports:	abind, ape, BayesLogit, fields, FNN, ggplot2, MASS, Matrix, matrixStats, MCMCpack, methods, nnet, rlang, parallel, pracma, pROC, sp, statmod, truncnorm
NeedsCompilation:	no
Packaged:	2022-08-11 13:37:09 UTC; jarioksa
Author:	Gleb Tikhonov [aut], Otso Ovaskainen [aut, cre], Jari Oksanen [aut], Melinda de Jonge [aut], Oystein Opedal [aut], Tad Dallas [aut]
Maintainer:	Otso Ovaskainen <otso.ovaskainen@helsinki.fi>
Repository:	CRAN
Date/Publication:	2022-08-11 14:10:14 UTC

Hmsc: Hierarchical Model of Species Communities

Description

Hierarchical Modelling of Species Communities (Hmsc) is a flexible framework for Joint Species Distribution Modelling (JSDMs). The framework can be used to relate species occurrences or abundances to environmental covariates, species traits and phylogenetic relationships. JSDMs are a special case of species distribution models (SDMs) that take into account the multivariate nature of communities which allows us to estimate community level responses as well capture biotic interactions and the influence of missing covariates in residual species associations. The Hmsc package contains functions to fit JSDMs, analyze the output and to generate predictions with these JSDMs.

Workflow

A typical workflow for a Hmsc analysis constists of 5 steps:

Step 1: Setting the model structure and fitting the model

The obligatory data for a Hmsc analysis includes a matrix of species occurrences or abundances (Y) and a data frame of environmental covariates (XData). The species matrix Y consists of ns columns representing the species and np rows representing the sampling units. Species data for different species does not have to be on the same scale, i.e. species 1 may be recorded as presence/absence while species 2 is recorded as abundance. XData consists of nc columns representing the environmental variables and np rows representing the sampling units. Y and XData need to have the same amount of rows (sampling units).

Optionally, the user can include information species traits, phylogenetic relationships and information on the spatiotemporal or hierarchical context of the sampling design to account for dependencies among the sampling units.

The model structure is created using the Hmsc function. As input, this function needs at least the matrix with species data (Y) and the dataframe with environmental variables. This is also the place where the data model is specified, as a default Hmsc assumes normally distributed species data. Other options are 'probit' for binary data, 'Poisson' and 'overdispersedPoisson' for count data. Additionally, you can specify the study design, the random effects structure, the species traits to be used, species phylogeny, how the covariate should be scales and how viable selection should be applied.

The random levels supplied to Hmsc are generated using HmscRandomLevel. Here, the user can specify the spatial or temporal information for the units as well as the covariates for covariate-depedent associations. If there is no spatial, temporal or coviate information for this level, the user should provide here the names of the units for this level.

After setting the model structure, the model is fitted by running sampleMcmc to sample from the posterior distributions of the model parameters.

Step 2: Examining MCMC convergence

After fitting the model, the MCMC convergence needs to be evaluated. The easiest way to do this is generate a coda object from the fitted model using convertToCodaObject. Convergence of the chains should be assessed qualitatively using trace plots of the model parameters and quantitatively using gelman diagnostics gelman.diag and by calculating the effective sample size with effectiveSize.

Step 3: Evaluating model fit

If MCMC convergence is satisfacotry, the model fit can be evaluated. Using the computePredictedValues function. It is also recommended to also compute predictive power using cross validation. Cross validation is done by supplying the partitioning created using using createPartition to computePredictedValues.

Step 4: Exploring parameter estimates

When the fitted model has satisfactory convergence and fit, the next step is generally to investigate the parameter estimates of the model. The posterior distribution of the parameter of choice is extracted from the model object using getPostEstimate. plotBeta and plotGamma can be used to visuzalize the beta and gamma parameter estimates. Additionally, the function biPlot can be used to construct an ordination plot from the eta and lamda parameters

Additionally, at this stage of the analysis, the user may want to look at how the variance explained by the model is partitioned using the computeVariancePartitioning and the plotVariancePartitioning function.

Step 5: Making predictions

Hmsc comes with a generic predict function that can be used to predict species occurences or counts in new units. These predictions can be unconditional or conditional on the occurence of other species. For conditional predictions, the user needs to supply species data for at least some species in the new units. Additionally, the package specific tools to make predictions along environmental gradients. These predictions can be vizualized at both the species and the community level (constructGradient, plotGradient).

Citing the Package

Tikhonov, G., Opedal, O.H., Abrego, N., Lehikoinen, A., de Jonge, M.M.J, Oksanen, J. and Ovaskainen, O. (2020) Joint species distribution modelling with the R-package Hmsc. Methods in Ecology and Evolution 11, 442–447. doi:10.1111/2041-210X.13345

See Also

Useful links:

• https://www.helsinki.fi/en/researchgroups/statistical-ecology/software/hmsc

• Report bugs at https://github.com/hmsc-r/HMSC/issues/

Hmsc

Description

Creates an Hmsc-class object

Usage

Hmsc(
  Y,
  XFormula = ~.,
  XData = NULL,
  X = NULL,
  XScale = TRUE,
  XSelect = NULL,
  XRRRData = NULL,
  XRRRFormula = ~. - 1,
  XRRR = NULL,
  ncRRR = 2,
  XRRRScale = TRUE,
  YScale = FALSE,
  studyDesign = NULL,
  ranLevels = NULL,
  ranLevelsUsed = names(ranLevels),
  TrFormula = NULL,
  TrData = NULL,
  Tr = NULL,
  TrScale = TRUE,
  phyloTree = NULL,
  C = NULL,
  distr = "normal",
  truncateNumberOfFactors = TRUE
)

Arguments

Details

Matrix Y may contain missing values, but it is not recommended to add a species/sampling unit with fully missing data, since those do not bring any new additional information.

Only one of XFormula-XData and X arguments can be specified. Similar requirement applies to TrFormula-TrData and Tr. It is recommended to use the specification with formula, since that information enables additional features for postprocessing of the fitted model.

As default, scaling is applied for X and Tr matrices, but not for Y matrix. If the X and/or Tr matrices are scaled, the estimated parameters are back-transformed so that the estimated parameters correspond to the original X and Tr matrices, not the scaled ones. In contrast, if Y is scaled, the estimated parameters are not back-transformed because doing so is not possible for all model parameters. Thus, the estimated parameters correspond to the scaled Y matrix, not the original one. If the Y matrix is scaled, the predictions generated by predict are back-transformed, so that the predicted Y matrices are directly comparable to the original Y matrix. If default priors are assumed, it is recommended that all matrices (X, Tr and Y) are scaled.

The object XSelect is a list. Each object of the list Xsel = XSelect[[i]] is a named list with objects Xsel$covGroup, Xsel$spGroup and Xsel$q. The parameter covGroup is a vector containing the columns of the matrix X for which variable selection is applied. The parameter spGroup is a vector of length equal to the number of species ns, with values 1,...,ng, where ng is the number of groups of species for which variable selection is applied simultanously. The parameter q is a vector of length ng, containing the prior probabilities by which the variables are to be included. For example, choosing covGroup = c(2,3), spGroup = rep(1,ns) and q=0.1 either includes or excludes both of the covariates 2 and 3 simultaneously for all species. For another example, choosing covGroup = c(2,3), spGroup = 1:ns and q=rep(0.1,ns) either includes or excludes both of the covariates 2 and 3 separately for each species.

The included random levels are specified by the ranLevels and ranLevelsUsed arguments. The correspondence between units of each random level and rows of Y must be specified by a column of studyDesign, which corresponds to the name of a list item in ranLevels. It is possible to provide an arbitrary number of columns in studyDesign that are not listed in ranLevels. These do not affect the model formulation or fitting scheme, but can be utilized during certain functions postprocessing the results of statistical model fit.

The distr argument may be either a matrix, a string literal, or a vector of string literals. In the case of a matrix, the dimension must be n_s \times 2, where the first column defines the family of the observation model and the second argument defines the dispersion property. The elements of the first column must take values 1-normal, 2-probit and 3-Poisson with log link function. The second argument stands for the dispersion parameter being fixed (0) or estimated (1). The default fixed values of the dispersion parameters are 1 for normal and probit, and 0.01 for Poisson (implemented as a limiting case of lognormally-overdispersed Poisson). Alternatively, a string literal shortcut can be given as a value to the distr argument, simultaniously specifying similar class of observation models for all species. The available shortcuts are "normal", "probit", "poisson", "lognormal poisson". If distr is a vector of string literals, each element corresponds to one species, should be either "normal", "probit", "poisson", "lognormal poisson", and these can be abbreviated as long as they are unique strings. The matrix argument and the vector of string literals allows specifying different observation models for different species.

By default this constructor assigns default priors to the latent factors. Those priors are designed to be reasonably flat assuming that the covariates, species traits and normally distributed responses are scaled. In case when other priors needed to be specified, a call of setPriors.Hmsc methods should be made, where the particular priors may be specified.

Value

An object of Hmsc class without any posterior samples.

See Also

HmscRandomLevel, sampleMcmc, setPriors.Hmsc

Examples

# Creating a Hmsc object without phylogeny, trait data or random levels
m = Hmsc(Y=TD$Y, XData=TD$X, XFormula=~x1+x2)

# Creating a Hmsc object with phylogeny and traits
m = Hmsc(Y=TD$Y, XData=TD$X, XFormula=~x1+x2,
TrData=TD$Tr, TrFormula=~T1+T2, phyloTree=TD$phylo)

# Creating a Hmsc object with 2 nested random levels (50 sampling units in 20 plots)
studyDesign = data.frame(sample = as.factor(1:50), plot = as.factor(sample(1:20,50,replace=TRUE)))
rL1 = HmscRandomLevel(units=levels(TD$studyDesign$plot))
rL2 = HmscRandomLevel(units=levels(TD$studyDesign$sample))
m = Hmsc(Y=TD$Y, XData=TD$X, XFormula=~x1+x2,
studyDesign=studyDesign,ranLevels=list("sample"=rL1,"plot"=rL2))

Create an Hmsc random level

Description

Specifies the structure of a random factor, including whether the random factor is assumed to be spatially explicit or not, the spatial coordinates and the potential structure of covariate-dependent random factors.

Usage

HmscRandomLevel(
  sData = NULL,
  sMethod = "Full",
  distMat = NULL,
  xData = NULL,
  units = NULL,
  N = NULL,
  nNeighbours = 10,
  sKnot = NULL,
  longlat = FALSE
)

Arguments

Details

Only one of sData, distMat, xData, units and N arguments can be provided.

As a good practice, we recommend to specify all available units for a random level, even if some of those are not used for training the model.

Value

a HmscRandomLevel-class object that can be used for Hmsc-class object construction

See Also

setPriors.Hmsc to change the default priors of an existing HmscRandomLevel object.

Examples

# Setting a random level with 50 units
rL = HmscRandomLevel(units=TD$studyDesign$sample)

# Setting a spatial random level
rL = HmscRandomLevel(sData=TD$xycoords)

# Setting a covariate-dependent random level.
rL = HmscRandomLevel(xData=data.frame(x1=rep(1,length(TD$X$x1)),x2=TD$X$x2))

Simulated data and a fitted Hmsc model for a small species community.

Description

This dataset contains simulated occurence data for 4 species in 50 sampling units. The data is based on a hierarchical study design consisting of 50 sampling units in 10 georeferenced plots. Occurences of 4 species were simulated using one continuous environmental variable (x1) and spatial autocorrelation between the plots. Response of species to the environment are related to one species trait which is fully phylogenetically structured. This dataset is used for the examples and package testing. The variables are as follows:

Usage

TD

Format

A list of 12 objects

ns

Number of species in the dataset

units

Number of sampling units

plots

Number of plots

X

A 3 by 50 environmental matrix consisting of one continuous and one categorical variable. Also includes intercept column

phy

A list containing the simulated phylogenetic tree for 4 species

C

A 4 by 4 phylogenetic variance covariance matrix

Tr

A 4 by 3 trait matrix with one phylogenetically phylogenetically structured continuous trait, one categroical trait and an intercept

xycoords

simulated 2 dimensional coordinates

studyDesign

Sampling unit and plot IDs

Y

Simulated species occurences

m

A fitted Hmsc object with 100 posterior samples

alignPosterior

Description

Aligns posterior in terms of variables susceptible to label switching

Usage

alignPosterior(hM)

Arguments

Value

an Hmsc model object that is identical to the input except the posterior being aligned

Examples

# Align the posterior for a previously fitted HMSC model
m = alignPosterior(TD$m)

biPlot

Description

Constructs an ordination biplot based on the fitted model

Usage

biPlot(
  hM,
  etaPost,
  lambdaPost,
  factors = c(1, 2),
  colVar = NULL,
  colors = NULL,
  spNames = hM$spNames,
  ...
)

Arguments

Examples

# Construct an ordination biplot using two chosen latent factors from a previously fitted HMSC model
etaPost = getPostEstimate(TD$m, "Eta")
lambdaPost=getPostEstimate(TD$m, "Lambda")
biPlot(TD$m, etaPost = etaPost, lambdaPost=lambdaPost, factors=c(1,2))

Combine Posterior Samples of Several Hmsc Models

Description

Function combines posterior samples of several sampled Hmsc models (see sampleMcmc) as new chains in the first fitted model. The combined models must be comparable, and there are some tests for detecting non-equal models. These tests will only give warning, and it is at user deliberation to decide which models and which posterior samples can be combined. You should be careful not start two models from the same random number seed, because these will only duplicate your data instead of providing new independent samples.

Usage

## S3 method for class 'Hmsc'
c(...)

Arguments

Value

An Hmsc model with chains of posterior samples.

Examples

## Fit a toy model with two chains
m1 <- sampleMcmc(TD$m, samples=10, transient=5, nChains=2, verbose=0)
## Need more data? Add chains: check carefully that these are
## sampled exactly like the previous model
m2 <- sampleMcmc(TD$m, nChains=2, samples=10, transient=5, verbose=0)
## Now four chains
m4 <- c(m1, m2)
m4

computeAssociations

Description

Computes the species association matrices associated with each random level

Usage

computeAssociations(hM, start = 1, thin = 1)

Arguments

Value

list of association matrices (\omega) corresponding to each random level in the model

Examples

# Compute the associations (residual correlations) between species from a HMSC model
assoc = computeAssociations(TD$m)

computeDataParameters

Description

Computes initial values before the sampling starts

Usage

computeDataParameters(hM)

Arguments

Value

a list including pre-computed matrix inverses and determinants (for phylogenetic and spatial random effects) needed in MCMC sampling

computeInitialParameters

Description

Computes initial parameter values before the sampling starts

Usage

computeInitialParameters(hM, initPar)

Arguments

Value

a list of Hmsc model parameters

computePredictedValues

Description

Computes predicted values from the fitted Hmsc model

Usage

computePredictedValues(
  hM,
  partition = NULL,
  partition.sp = NULL,
  start = 1,
  thin = 1,
  Yc = NULL,
  mcmcStep = 1,
  expected = TRUE,
  initPar = NULL,
  nParallel = 1,
  nChains = length(hM$postList),
  updater = list(),
  verbose = hM$verbose,
  alignPost = TRUE
)

pcomputePredictedValues(
  hM,
  partition = NULL,
  partition.sp = NULL,
  start = 1,
  thin = 1,
  Yc = NULL,
  mcmcStep = 1,
  expected = TRUE,
  initPar = NULL,
  nParallel = 1,
  useSocket = .Platform$OS.type == "windows",
  nChains = length(hM$postList),
  updater = list(),
  verbose = nParallel == 1,
  alignPost = TRUE
)

Arguments

Details

There are two alternative functions computePredictedValues and pcomputePredictedValues. Function pcomputePredictedValues uses more aggressive parallelization and can be much faster when partition is used. Function computePredictedValues can run chains of each sampleMcmc partition in parallel, but pcomputePredictedValues can run each partition fold times chain in parallel (if hardware and operating systems permit). Function pcomputePredictedValues is still experimental, and therefore we provide both the old and new functions, but the old functions is scheduled to be removed in the future. Species partitions are not yet parallelized, and they can be very slow, especially with many mcmcSteps.

If the option partition is not used, the posterior predictive distribution is based on the model fitted to the full data. If the option partition is used but partition.sp is not used, the posterior predictive distribution is based on cross-validation over the sampling units. If partition.sp is additionally used, then, when predictions are made for each fold of the sampling units, the predictions are done separately for each fold of species. When making the predictions for one fold of species, the predictions are conditional on known occurrences (those observed in the data) of the species belonging to the other folds. If partition.sp is used, the parameter mcmcStep should be set high enough to obtain appropriate conditional predictions. The option Yc can be used alternatively to partition.sp if the conditioning is to be done based on a fixed set of data (independently of which sampling unit and species the predictions are made for).

Value

an array of model predictions, made for each posterior sample

See Also

predict.Hmsc

Examples

# Compute predicted values using a previously fitted HMSC model
preds = computePredictedValues(TD$m)

## Not run: 
# Compute predicted values for a previously fitted HMSC model using 2 folds
partition = createPartition(TD$m, nfolds = 2)
predsCV1 = computePredictedValues(TD$m,partition=partition)

# Compute conditional predictions for a previously fitted HMSC model using 2 folds
partition = createPartition(TD$m, nfolds = 2)
predsCV2 = computePredictedValues(TD$m, partition = partition,
partition.sp = 1:TD$m$ns, mcmcStep = 100)

## End(Not run)

computeVariancePartitioning

Description

Computes variance components with respect to given grouping of fixed effects and levels of random effects

Usage

computeVariancePartitioning(
  hM,
  group = NULL,
  groupnames = NULL,
  start = 1,
  na.ignore = FALSE
)

Arguments

Details

The vector group has one value for each column of the matrix hM$X, describing the index of the group in which this column is to be included. The names of the group are given by groupnames. The output object VP$vals gives the variance proportion for each group and species. The output object VP$R2T gives the variance among species explained by traits, measured for species' responses to covariates (VP$R2T$Beta) and species occurrences (VP$R2T$Y)

Value

returns an object VP with components VP$vals, VP$R2T, VP$group and VP$groupnames.

See Also

Use plotVariancePartitioning to display the result object.

Examples

# Partition the explained variance for a previously fitted model
# without grouping environmental covariates
VP = computeVariancePartitioning(TD$m)

# Partition the explained variance for a previously fitted model
# while grouping the two environmental variables together
VP = computeVariancePartitioning(TD$m, group=c(1,1), groupnames = c("Habitat"))

computeWAIC

Description

Computes the value of WAIC (Widely Applicable Information Criterion) for the Hmsc model

Usage

computeWAIC(hM, ghN = 11, byColumn = FALSE)

Arguments

Details

The result is exact for normal and probit observational models. For Poisson-type observational model the result is obtained through numerical integration using Gauss-Hermite quadrature.

Value

the scalar WAIC

Examples

# Compute WAIC of previously sampled Hmsc object
WAIC = computeWAIC(TD$m)

constructGradient

Description

Constructs an environmental gradient over one of the variables included in XData

Usage

constructGradient(
  hM,
  focalVariable,
  non.focalVariables = list(),
  ngrid = 20,
  coordinates = list()
)

Arguments

Details

In basic form, non.focalVariables is a list, where each element is on the form variable=list(type,value), where variable is one of the non-focal variables, and the value is needed only if type = 3. Alternatives type = 1 sets the values of the non-focal variable to the most likely value (defined as expected value for covariates, mode for factors), type = 2 sets the values of the non-focal variable to most likely value, given the value of focal variable, based on a linear relationship, and type = 3 fixes to the value given. As a shortcut, a single number 1 or 2 can be given as a type used for all non-focal variables. If a non.focalVariable is not listed, type=2 is used as default. Note that if the focal variable is continuous, selecting type 2 for a non-focal categorical variable can cause abrupt changes in response.

The function needs access to the original XData data frame, and cannot be used if you defined your model with X model matrix. In that case you must construct your gradient manually.

Value

a named list with slots XDataNew, studyDesignNew and rLNew

See Also

plotGradient, predict.

Examples

# Construct gradient for environmental covariate called 'x1'.
Gradient = constructGradient(TD$m, focalVariable="x1")

# Construct gradient for environmental covariate called 'x1'
# while setting the other covariate to its most likely values
Gradient = constructGradient(TD$m, focalVariable="x1",non.focalVariables=list(x2=list(1)))

constructKnots

Description

Construct a Regular Grid of Knot Locations for Spatial GPP Model

Usage

constructKnots(sData, nKnots = NULL, knotDist = NULL, minKnotDist = NULL)

Arguments

Details

This is a helper function for spatial Hmsc models with the spatial method set to GPP where user must provide knot locations. Knot locations with a distance greater than minKnotDist to the nearest data point are dropped from the grid. If the input locations are SpatialPoints data, these are treated like Euclidean coordinates, and if the points are not projected, a warning is issued.

Only one of nKnots and minKnotDist arguments can be provided.

Value

a data frame with knot locations

Examples

#Creating knots for some 2 dimensional spatial data
n = 100
xycoords = matrix(runif(2*n),ncol=2)
xyKnots = constructKnots(xycoords,knotDist = 0.2, minKnotDist = 0.5)

convertToCodaObject

Description

Converts the Hmsc posterior into a named list of mcmc.list objects

Usage

convertToCodaObject(
  hM,
  start = 1,
  spNamesNumbers = c(TRUE, TRUE),
  covNamesNumbers = c(TRUE, TRUE),
  trNamesNumbers = c(TRUE, TRUE),
  Beta = TRUE,
  Gamma = TRUE,
  V = TRUE,
  Sigma = TRUE,
  Rho = TRUE,
  Eta = TRUE,
  Lambda = TRUE,
  Alpha = TRUE,
  Omega = TRUE,
  Psi = TRUE,
  Delta = TRUE
)

Arguments

Value

A named list that can be analysed with coda functions.

Examples

# Convert recorded posterior samples in \code{Hmsc} object to coda object
codaObject = convertToCodaObject(TD$m)

# Convert recorded posterior samples, starting from sample 100, in m object to coda object
codaObject = convertToCodaObject(TD$m, start=100)

createPartition

Description

Constructs a partition vector given the number of folds and column of study design

Usage

createPartition(hM, nfolds = 10, column = NULL)

Arguments

Value

a vector describing the fold of each sampling unit

Examples

# Create 3 folds for a HMSC object
partition = createPartition(TD$m, nfolds = 3)

evaluateModelFit

Description

Computes measures of model fit for a Hmsc model

Usage

evaluateModelFit(hM, predY)

Arguments

Details

All measures of model fit are based on comparing the posterior predictive distribution (predY)) to the observed values (hM$Y). The predicted distribution is first summarized to a single matrix of predicted values by taking the posterior mean (for normal and probit models) or posterior median (for Poisson models). All measures of model fit are given as vectors with one value for each species.

The kinds of measures of model fit depend on the type of response variable. For all types of response variables, root-mean-square error (RMSE) between predicted and observed values is computed. For normal models, R2 is computed as squared pearson correlation between observed and predicted values, times the sign of the correlation. For probit models, Tjur R2 and AUC are computed. For Poisson models, a pseudo-R2 is computed as squared spearman correlation between observed and predicted values, times the sign of the correlation (SR2). For Poisson models, the observed and predicted data are also truncated to occurrences (presence-absences), for which the same measures are given as for the probit models (O.RMSE, O.AUC and O.TjurR2). For Poisson models, the observed and predicted data are also subsetted to conditional on presence, for which the root-mean-square error and pseudo-R2 based on squared spearman correlation are computed (C.RMSE, C.SR2).

The measures O.RMSE, O.AUC, O.TjurR2, C.RMSE and C.SR2 can be computed only if the option expected=FALSE has been used when making the predictions

If the model includes a mixture of response variable types, the resulting measures of model fit contain NA's for those response variables for which they cannot be computed.

Value

a list of measures of model fit

Examples

# Evaluate model fit
preds = computePredictedValues(TD$m)
MF = evaluateModelFit(hM=TD$m, predY=preds)

# Evaluate model performance based on cross validation: this will be slow
## Not run: 
partition = createPartition(TD$m, nfolds = 2)
predsCV1 = computePredictedValues(TD$m, partition=partition)
MF = evaluateModelFit(hM=TD$m, predY=predsCV1)

## End(Not run)

getPostEstimate

Description

Calculates mean, support and other posterior quantities for a specified model parameter

Usage

getPostEstimate(
  hM,
  parName,
  r = 1,
  x = NULL,
  q = c(),
  chainIndex = 1:length(hM$postList),
  start = 1,
  thin = 1
)

Arguments

Value

A named list of posterior quantities.

Examples

# Get posterior mean and support for species' responses to environmental covariates
postBeta = getPostEstimate(TD$m, parName='Beta')

# Get posterior mean and support for species' responses to latent factors for the first random level
postLambda = getPostEstimate(TD$m, parName='Lambda', r=1)

plotBeta

Description

Plots heatmaps of parameter estimates or posterior support values of species' environmental responses, i.e. how species in Y responds to covariates in X

Usage

plotBeta(
  hM,
  post,
  param = "Support",
  plotTree = FALSE,
  SpeciesOrder = "Original",
  SpVector = NULL,
  covOrder = "Original",
  covVector = NULL,
  spNamesNumbers = c(TRUE, TRUE),
  covNamesNumbers = c(TRUE, TRUE),
  supportLevel = 0.9,
  split = 0.3,
  cex = c(0.7, 0.7, 0.8),
  colors = colorRampPalette(c("blue", "white", "red")),
  colorLevels = NULL,
  mar = NULL,
  marTree = c(6, 0, 2, 0),
  mgp = c(3, 2, 0),
  main = NULL,
  smallplot = NULL,
  bigplot = NULL,
  newplot = TRUE
)

Arguments

Examples

# Plot posterior support values of species' environmental responses
betaPost=getPostEstimate(TD$m, "Beta")
plotBeta(TD$m, post=betaPost, param="Support")

# Plot parameter estimates of species' environmental responses together with the phylogenetic tree
betaPost=getPostEstimate(TD$m, "Beta")
plotBeta(TD$m, post=betaPost, param="Mean", plotTree=TRUE)

plotGamma

Description

Plots heatmaps of parameter estimates or posterior support values of trait effects on species' environmental responses, i.e. how environmental responses in Beta responds to covariates in X

Usage

plotGamma(
  hM,
  post,
  param = "Support",
  trOrder = "Original",
  trVector = NULL,
  covOrder = "Original",
  covVector = NULL,
  trNamesNumbers = c(TRUE, TRUE),
  covNamesNumbers = c(TRUE, TRUE),
  supportLevel = 0.9,
  main = NULL,
  cex = c(0.8, 0.8, 0.8),
  colors = colorRampPalette(c("blue", "white", "red")),
  colorLevels = NULL,
  mar = c(6, 9, 2, 0),
  smallplot = NULL,
  bigplot = NULL,
  newplot = TRUE
)

Arguments

Examples

# Plot posterior support values of trait effects on environmental responses
gammaPost=getPostEstimate(TD$m, "Gamma")
plotGamma(TD$m, post=gammaPost, param="Support")

# Plot parameter estimates of trait effects on environmental responses
gammaPost=getPostEstimate(TD$m, "Gamma")
plotGamma(TD$m, post=gammaPost, param="Mean")

plotGradient

Description

Plots an environmental gradient over one of the variables included in XData

Usage

plotGradient(
  hM,
  Gradient,
  predY,
  measure,
  xlabel = NULL,
  ylabel = NULL,
  index = 1,
  q = c(0.025, 0.5, 0.975),
  cicol = rgb(0, 0, 1, alpha = 0.5),
  pointcol = "lightgrey",
  pointsize = 1,
  showData = FALSE,
  jigger = 0,
  yshow = NA,
  showPosteriorSupport = TRUE,
  main,
  ...
)

Arguments

Details

For measure="Y", index selects which species to plot from hM$spNames. For measure="T", index selects which trait to plot from hM$trNames. With measure="S" the row sum of pred is plotted, and thus the interpretation of "species richness" holds only for probit models. For Poisson models "S" shows the total count, whereas for normal models it shows the summed response. For measure="T", in probit model the weighting is over species occurrences, whereas in count models it is over individuals. In normal models, the weights are exp-transformed predictions to avoid negative weights

Value

For the case of a continuous covariate, returns the posterior probability that the plotted variable is greater for the last sampling unit of the gradient than for the first sampling unit of the gradient. For the case of a factor, returns the plot object.

See Also

constructGradient, predict

Examples

# Plot response of species 2 over the gradient of environmental variable x1
Gradient = constructGradient(TD$m, focalVariable="x1")
predY = predict(TD$m, Gradient=Gradient)
plotGradient(TD$m, Gradient, pred=predY, measure="Y", index = 2, showData = TRUE, jigger = 0.05)
# Plot modelled species richness over the gradient of environmental variable x1
Gradient = constructGradient(TD$m, focalVariable="x1")
predY = predict(TD$m, Gradient=Gradient)
plotGradient(TD$m, Gradient, pred=predY, measure="S")

plotVariancePartitioning

Description

Plots the results of variance partitioning of a Hmsc model produced by computeVariancePartitioning as a barplot

Usage

plotVariancePartitioning(
  hM,
  VP,
  cols = NULL,
  main = "Variance Partitioning",
  ...
)

Arguments

Examples

# Plot how the explained variance of a previously fitted model is partitioned
VP = computeVariancePartitioning(TD$m)
plotVariancePartitioning(TD$m, VP)

poolMcmcChains

Description

Combines a list of single or several MCMC chains into a single chain

Usage

poolMcmcChains(postList, chainIndex = 1:length(postList), start = 1, thin = 1)

Arguments

Value

a list with combined MCMC samples

Examples

# Combine the posteriors from all chains in a Hmsc object
postList = TD$m$postList
pooledPost = poolMcmcChains(postList)

predict

Description

Calculates predicted values from a fitted Hmsc model.

Usage

## S3 method for class 'Hmsc'
predict(
  object,
  post = poolMcmcChains(object$postList),
  XData = NULL,
  X = NULL,
  XRRRData = NULL,
  XRRR = NULL,
  studyDesign = object$studyDesign,
  ranLevels = object$ranLevels,
  Gradient = NULL,
  Yc = NULL,
  mcmcStep = 1,
  expected = FALSE,
  predictEtaMean = FALSE,
  predictEtaMeanField = FALSE,
  nParallel = 1,
  useSocket = TRUE,
  ...
)

Arguments

Details

In mcmcStep,the number of extra mcmc steps used for updating the random effects for the Eta parameters, starting from the samples of the fitted Hmsc model in order to account for the conditional infromation provided in the Yc argument. The higher this number is, the more the obtained updated samples are unaffected by the posterior estimates of latent factors in the model fitted to the training data and more resembles the true conditional posterior. However, the elapsed time for conditional prediction grows approximately linearly as this parameter increases. The exact number for sufficient is problem-dependent and should be assessed by e.g. gradually increasing this parameter till the stationarity of the produced predictions.

Value

A list of length length(post), each element of which contains a sample from the posterior predictive distribution (given the sample of the Hmsc model parameters in the corresponding element of the post argument)

See Also

predictLatentFactor

predictLatentFactor

Description

Draws samples from the conditional predictive distribution of latent factors

Usage

predictLatentFactor(
  unitsPred,
  units,
  postEta,
  postAlpha,
  rL,
  predictMean = FALSE,
  predictMeanField = FALSE
)

Arguments

Details

Length of units vector and number of rows in postEta matrix shall be equal. The method assumes that the i-th row of postEta correspond to i-th element of units.

This method uses only the coordinates rL$s field of the rL$s argument. This field shall be a matrix with rownames covering the union of unitsPred and units factors. Alternatively, it can use distance matrix rL$distMat which is a symmetric square matrix with similar row names as the coordinate data (except for the GPP models that only can use coordinates).

In case of spatial random level, the computational complexity of the generic method scales cubically as the number of unobserved units to be predicted. Both predictMean=TRUE and predictMeanField=TRUE options decrease the asymptotic complexity to linear. The predictMeanField=TRUE option also preserves the uncertainty in marginal distribution of predicted latent factors, but neglects the inter-dependece between them.

Value

a list of length length(postEta) containing samples of random factors at unitsPred from their predictive distribution conditional on the values at units

prepareGradient

Description

prepares a user-made environmental and/or spatial gradient to be used for prediction

Usage

prepareGradient(hM, XDataNew, sDataNew)

Arguments

Details

The dataframe XDataNew is the for output as for input. The main purpose of this function is to prepare the study design and random levels so that predictions can be made with the predict function. Note that the difference between constructGradient and prepareGradient is that while prepareGradient takes as input the new environmental and spatial data, constructGradient generates those data to represent a new environmental gradient.

Value

a named list with members XDataNew, studyDesignNew and rLNew

See Also

constructGradient, predict

sampleMCMC

Description

Samples the posterior with block-conditional Gibbs MCMC sampler

Usage

sampleMcmc(
  hM,
  samples,
  transient = 0,
  thin = 1,
  initPar = NULL,
  verbose,
  adaptNf = rep(transient, hM$nr),
  nChains = 1,
  nParallel = 1,
  useSocket = TRUE,
  dataParList = NULL,
  updater = list(),
  fromPrior = FALSE,
  alignPost = TRUE
)

Arguments

Details

The exact number of samples to be recorded in order to get a proper estimate of the full posterior with Gibbs MCMC algorithms, as well as the required thinning and cut-off of transient is very problem-specific and depends both on the model structure and the data itself. Therefore, in general it is very challenging to a priori provide an informed recommendation on what values should be used for a particular problem. A common recommended strategy involves executing the posterior sampling with MCMC with some guess of the values for these arguments, checking the properties of the obtained samples (primarily potential scale reduction factor and effective sample size), and adjusting the guess accordingly.

The value of 1 for thin argument means that at each MCMC step after the transient a sample is recorded.

Typically, the value of nParallel equal to nChains leads to most efficient usage of available parallelization capacities. However, this may be not the case if R is configured with multi-tread linear algebra libraries. For debug and test purposes, the nParallel should be set to 1, since only in this case a details of the potentially encountered errors would be available.

The dataParList argument may be handy for large problems that needs to be refitted multiple times, e.g. with different prior values. In that case, the data parameters that are precomputed for the Hmsc sampling scheme may require an undesirably lot of storage space if they are saved for each of the model. Instead, they could be computed only once and then directly reused, therefore reducing the storing redundancy.

Some of the available conditional updaters partially duplicate each other. In certain cases, the usage of all of them may lead to suboptimal performance, compared to some subset of those. Then, it is possible to manually disable some of them, by adding a $UPDATER_NAME=FALSE pair to the updater argument. Another usage of this argument involves cases when some of the model parameters are known and have to be fixed. However, such tweaks of the sampling scheme should be done with caution, as if compromized they would lead to erroneuos results.

Value

An Hmsc-class object with chains of posterior samples added to the postList field

See Also

Hmsc

Examples

## you need 1000 or more samples, but that will take too long
## in an example
m = sampleMcmc(TD$m, samples=10)

## Not run: 
## Record 1000 posterior samples while skipping 1 MCMC step between samples
## from 2 chains after discarding the first 500 MCMC steps
m = sampleMcmc(TD$m, samples=1000, transient=500, thin=2, nChains=2, nParallel=1)

## End(Not run)

samplePrior

Description

Samples the parameter vector from prior

Usage

samplePrior(hM, dataParList = NULL)

Arguments

Value

A named list containing the Hmsc model parameters

setPriors

Description

Sets or resets priors to objects

Usage

setPriors(...)

Arguments

Value

Object of same type as first input

See Also

setPriors.Hmsc, setPriors.HmscRandomLevel

Examples

# Set priors for random level so that there is minimum of 2 latent factors and maximum of 3
rL1 = HmscRandomLevel(units=TD$studyDesign$plot)
rL1 = setPriors(rL1, nfMax=3, nfMin=2)

# Set shrinkage parameters for priors of random level
rL1 = HmscRandomLevel(units=TD$studyDesign$plot)
rL1 = setPriors(rL1, a1=10, a2=10, b1=1, b2=1)

setPriors.Hmsc

Description

Sets or resets priors to the Hmsc object

Usage

## S3 method for class 'Hmsc'
setPriors(
  hM,
  V0 = NULL,
  f0 = NULL,
  mGamma = NULL,
  UGamma = NULL,
  aSigma = NULL,
  bSigma = NULL,
  nuRRR = NULL,
  a1RRR = NULL,
  b1RRR = NULL,
  a2RRR = NULL,
  b2RRR = NULL,
  rhopw = NULL,
  setDefault = FALSE,
  ...
)

Arguments

Value

Modified Hmsc object

setPriors.HmscRandomLevel

Description

Sets or resets priors to the Hmsc object

Usage

## S3 method for class 'HmscRandomLevel'
setPriors(
  rL,
  nu = NULL,
  a1 = NULL,
  a2 = NULL,
  b1 = NULL,
  b2 = NULL,
  alphapw = NULL,
  nfMax = NULL,
  nfMin = NULL,
  setDefault = FALSE,
  ...
)

Arguments

Value

Modified HmscRandomLevel object