| Title: | Energy Landscape Analysis for Ecological Communities |
|---|---|
| Description: | Energy landscape analysis (ELA) is a systematic method for analyzing an energy landscape represented as a weighted network. This R package is especially designed to analyze ecological dynamics, i.e., transitions in community compositions. For details of the analysis framework, please visit Suzuki K et al. (2021) <doi:10.1002/ecm.1469>. |
| Authors: | Kenta Suzuki [aut, cre, cph], Sotaro Takano [aut], Hiroaki Fujita [aut], Zeus Sato [ctb], Fumihiko Ayabe [ctb] |
| Maintainer: | Kenta Suzuki <[email protected]> |
| License: | BSD_3_clause + file LICENSE |
| Version: | 1.0.2 |
| Built: | 2026-05-28 18:10:21 UTC |
| Source: | https://github.com/cran/ELAplus |
Energy landscape analysis (ELA) is a systematic method for analyzing an energy landscape represented as a weighted network. This R package is especially designed to analyze ecological dynamics, i.e., transitions in community compositions. For details of the analysis framework, please visit Suzuki K et al. (2021) <doi:10.1002/ecm.1469>.
This package provides tools for energy landscape analysis (ELA), a framework for analyzing dynamical systems represented as weighted networks. The main functions support construction of energy landscapes, by fitting the community composition data to a pairwise maximum entropy model (or its extension) and subsequent identification of stable states, and analysis of transitions between community compositions.
Kenta Suzuki [aut, cre, cph], Sotaro Takano [aut], Hiroaki Fujita [aut], Zeus Sato [ctb], Fumihiko Ayabe [ctb]
Maintainer: Kenta Suzuki <[email protected]>
Suzuki, K., Nakaoka, S., Fukuda, S., & Masuya, H. (2021). "Energy landscape analysis elucidates the multistability of ecological communities across environmental gradients." Ecological Monographs. doi:10.1002/ecm.1469 https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1002/ecm.1469
GitHub: https://github.com/kecosz/rELA, Portal: https://elaportal.org
# load sample dataset data("baseabtable") data("basemetadata") # Formatting dataset formatted <- Formatting(baseabtable, basemetadata, 0, c(0.05, 0.05, 0.95)) ocvecs <- formatted[[1]] envecs <- formatted[[3]] # Example (without environmental parameters) # Fitting the data to the extended maximum pairwise entropy model # without explicit environmental parameters sa <- runSA( ocvecs, enmat = NULL, rep = 1, getall = FALSE, lambda = 0.01, we = 0.001, maxlr = 0.005, totalit = 1000 ) # Using the fitted model, normal ELA is run. ela <- ELA(sa, env=NULL, threads=1,SS.itr=2000, FindingTip.itr=1000, reporting=TRUE) # Pruning stable states elap <- ELPruning(ela, th=0.1) # show disconnectivity graph showDG(elap[[1]], ocvecs, "test") # show species interaction graph showIntrGraph(elap[[1]], sa, th=0.01, annot_adj=c(0.75, 2.00)) # Example (with environmental parameters) # Find best parameter sets for the machine learning (with envecs) bpresult <- Findbp(ocvecs,enmat=envecs,rep=1,threads=1,nrepeat=1,cvmode="5fold") bp <- bpresult[[1]] bpm <- as.numeric(unlist(strsplit(names(bp)[1], split = "_"))) # Fitting the data to the extended maximum pairwise entropy model # with explicit environmental parameters sa <- runSA( ocvecs, enmat = envecs, rep = 16, getall = FALSE, lambda = bpm[[1]], we = bpm[[2]], maxlr = bpm[[3]], totalit = bpm[[4]] ) # Using the fitted model, gradient ELA is run. gela <- GradELA( sa = sa, eid = "factor.1", enmat = envecs, env = NULL, range = c(0, 1), steps = 32, th = 0.05, threads = 2 ) showSSD(gela)# load sample dataset data("baseabtable") data("basemetadata") # Formatting dataset formatted <- Formatting(baseabtable, basemetadata, 0, c(0.05, 0.05, 0.95)) ocvecs <- formatted[[1]] envecs <- formatted[[3]] # Example (without environmental parameters) # Fitting the data to the extended maximum pairwise entropy model # without explicit environmental parameters sa <- runSA( ocvecs, enmat = NULL, rep = 1, getall = FALSE, lambda = 0.01, we = 0.001, maxlr = 0.005, totalit = 1000 ) # Using the fitted model, normal ELA is run. ela <- ELA(sa, env=NULL, threads=1,SS.itr=2000, FindingTip.itr=1000, reporting=TRUE) # Pruning stable states elap <- ELPruning(ela, th=0.1) # show disconnectivity graph showDG(elap[[1]], ocvecs, "test") # show species interaction graph showIntrGraph(elap[[1]], sa, th=0.01, annot_adj=c(0.75, 2.00)) # Example (with environmental parameters) # Find best parameter sets for the machine learning (with envecs) bpresult <- Findbp(ocvecs,enmat=envecs,rep=1,threads=1,nrepeat=1,cvmode="5fold") bp <- bpresult[[1]] bpm <- as.numeric(unlist(strsplit(names(bp)[1], split = "_"))) # Fitting the data to the extended maximum pairwise entropy model # with explicit environmental parameters sa <- runSA( ocvecs, enmat = envecs, rep = 16, getall = FALSE, lambda = bpm[[1]], we = bpm[[2]], maxlr = bpm[[3]], totalit = bpm[[4]] ) # Using the fitted model, gradient ELA is run. gela <- GradELA( sa = sa, eid = "factor.1", enmat = envecs, env = NULL, range = c(0, 1), steps = 32, th = 0.05, threads = 2 ) showSSD(gela)
'baseabtable' is a data frame containing abundance for 16 taxa across 256 samples.
baseabtablebaseabtable
A data frame with 256 rows and 16 variables:
Each row corresponds to a sample, and each column represents the abundance of a specific taxon.
The values represent abundance measurements for each taxon in each sample. The exact definition of "abundance" (e.g., counts, relative abundance, or normalized values) depends on the data preprocessing pipeline used.
Generated for example purposes (replace with actual source if applicable)
data(baseabtable) head(baseabtable) formatted <- Formatting(baseabtable = baseabtable)data(baseabtable) head(baseabtable) formatted <- Formatting(baseabtable = baseabtable)
'basemetadata' is a data frame containing environmental parameters
associated with the same 256 samples as in baseabtable.
basemetadatabasemetadata
A data frame with 256 rows and 2 variables:
Each row corresponds to a sample, and each column represents an environmental variable.
The rows in this data frame correspond to the same samples and
ordering as in baseabtable. These variables can be used
for ELA with environmental parameters such as gradELA
Generated for example purposes (replace with actual source if applicable)
data(baseabtable) data(basemetadata) head(basemetadata) formatted <- Formatting(baseabtable = baseabtable, basemetadata = basemetadata)data(baseabtable) data(basemetadata) head(basemetadata) formatted <- Formatting(baseabtable = baseabtable, basemetadata = basemetadata)
Identify basin (stable state) and its energy for a given state.
Applies SteepestDescent to the input state and returns the basin (stable state) reached, together with its energy. The stable state is encoded as an integer ID using bin2id.
Bi(xi, h, j)Bi(xi, h, j)
xi |
Binary vector representing the initial state. |
h |
vector of implicit/explicit environmental effects on species. |
j |
Matrix of species interactions. |
A list with two elements: (1) integer stable state ID (basin assignment), (2) energy of the corresponding stable state.
Encode a binary vector into a base-64 id string index.
bin2id(bin)bin2id(bin)
bin |
A binary vector (0/1, community composition) to encode. |
A character scalar representing the encoded id.
Running energy landscape analysis with bootstrap sampling of estimated saresults
Performs ELA and ELPruning using SA fitting results based on the bootstrap resampled occurrence and environmental data using bootstrap_SA.
Basically, representative ELA results with original occurrence and environmental data should be given by ela=ela,
alternatively the function can select a representative parameter where commonly appearing stable states are most frequent.
The function computes bootstrap energy distributions for all stable states and tipping points
appearing in the specified pruned landscape. Energies that are not present in a given trial
are stored as NaN in the summary table.
bootstrap_ELA( bs_params, ocmat, ela = NULL, env = NULL, enmat = NULL, threads = 1, reporting = FALSE )bootstrap_ELA( bs_params, ocmat, ela = NULL, env = NULL, enmat = NULL, threads = 1, reporting = FALSE )
bs_params |
output of bootstrap_SA. list of estimated parameters using bootstrap resampled data set. |
ocmat |
community composition matrix. |
ela |
representative ELA result (output of ELA). If |
env |
target environmental factor. |
enmat |
array of the vectors of environmental factors for each sample. Normalization required. |
threads |
number of cores (threads) used for calculation. |
reporting |
enable/disable the reporting function. |
A list with two elements:
(1) representative pruned ELA result (stable states, energies, tipping points, tipping energies),
(2) a data frame of bootstrap energies (rows = state IDs, columns = trials) with NaN for missing states.
Perform runSA algorithm and estimate parameters using bootstrap resampled
data set (ocmat and enmat).
The basic configuration is same as runSA.
bootstrap_SA( ocmat, enmat = NULL, bootstrap.it = 2000, rep = 128, threads = 1, we = 0.001, totalit = 1000, lambda = 0.01, runadamW = TRUE, maxlr = 0.005, sparse = TRUE, reporting = TRUE )bootstrap_SA( ocmat, enmat = NULL, bootstrap.it = 2000, rep = 128, threads = 1, we = 0.001, totalit = 1000, lambda = 0.01, runadamW = TRUE, maxlr = 0.005, sparse = TRUE, reporting = TRUE )
ocmat |
Numeric matrix (0/1). Rows are samples and columns are species. |
enmat |
Optional numeric matrix. Environmental variables for each sample. Should be normalized/standardized beforehand. |
bootstrap.it |
number of bootstrap trials (default: 1000). |
rep |
Integer. Number of parallel runs (replicates) for optimization. |
threads |
Integer. Number of cores used when running in parallel. |
we |
Numeric. Weight decay parameter for AdamW. |
totalit |
Integer. Total number of optimization iterations. |
lambda |
Numeric. Sparsity penalty parameter (used when |
runadamW |
Logical. If |
maxlr |
Numeric. Maximum learning rate. |
sparse |
Logical. If |
reporting |
Logical. If |
a list of SA parameters estimated with bootstrap resampling (length bootstrap.it).
For each stable state, this function computes the difference between the energy values of the corresponding tipping point and the stable state across bootstrap samples, then draws a histogram of the resulting deltaE values.
deltaE_histogram(bootstrap_ene, grobj_to_plot, getData = TRUE)deltaE_histogram(bootstrap_ene, grobj_to_plot, getData = TRUE)
bootstrap_ene |
A data frame containing energies in stable states and tipping points, output of bootstrap_ELA. |
grobj_to_plot |
A data frame containing stable state - tipping point relationship, output of showDG |
getData |
Logical. If 'TRUE', return the computed deltaE values as a list. If 'FALSE', return 'NULL'. |
If 'getData = TRUE', a list of numeric vectors, one per deltaE for each stable state. If 'getData = FALSE', 'NULL'.
Draws a disconnectivity graph from the data frame returned by GraphObj.
Nodes are plotted at their assigned x-positions and energy levels, and are
connected to their parent components by paths.
If show_pvalue = TRUE, approximate p-values for the energy gaps of
stable states are displayed. The p-values are shown on a log10 scale using
the values specified in scale_labels.
DisconnectivityGraph( grobj, s, DG_sample = "DG_sample", annot_adj = c(0.75, 2), min_pval = 0.05, scale_labels = c(0.05, 0.1, 0.2, 0.5), show_pvalue = FALSE, getGraphObj = FALSE, plot = TRUE )DisconnectivityGraph( grobj, s, DG_sample = "DG_sample", annot_adj = c(0.75, 2), min_pval = 0.05, scale_labels = c(0.05, 0.1, 0.2, 0.5), show_pvalue = FALSE, getGraphObj = FALSE, plot = TRUE )
grobj |
A data frame returned by |
s |
Integer; number of species/features (used for labeling by id2bin). |
DG_sample |
Character; title of the plot. |
annot_adj |
Numeric vector of length 2; horizontal and vertical adjustment for labels. |
min_pval |
Numeric (must be > 0); minimum p-value shown in the visualization. This is usually set to the minimum detection level in bootstrap sampling (e.g., if the number of bootstrap iterations is 1000, the minimum level is 1/1000 = 0.001). |
scale_labels |
Numeric vector of values between 0 and 1; labels for the p-value scale. These values are converted to log-scaled values for plotting. |
show_pvalue |
Logical; if |
getGraphObj |
Logical; if |
plot |
Logical; if TRUE, show disconnectivity graph |
If getGraphObj = TRUE, returns a data frame containing energy values,
x-positions, and p-values for the energy gaps of the identified stable states
and tipping points. Otherwise, returns NULL.
Identifying stable states and tipping points.
This function estimates stable states by repeated steepest-descent searches
SSestimate and then checks tipping points between all pairs of stable
states using FindingTippingpoint_cpp.
The result is summarized and converted by elaconvert
into a list of stable state IDs, their energies, tipping point IDs, and
tipping energies.
ELA( sa = sa, env = NULL, SS.itr = 20000, FindingTip.itr = 10000, threads = 1, reporting = TRUE, fork = TRUE )ELA( sa = sa, env = NULL, SS.itr = 20000, FindingTip.itr = 10000, threads = 1, reporting = TRUE, fork = TRUE )
sa |
output of runSA function. |
env |
vector of environmental factors for specifying a single environmental condition (optional). |
SS.itr |
number of iterations of the process to identify the stable states. |
FindingTip.itr |
number of iterations of the process to identify the tipping points. |
threads |
number of cores (threads) used for calculation. |
reporting |
enable/disable the reporting function. |
fork |
logical; if TRUE and |
A list returned by elaconvert: (1) stable state IDs, (2) stable state energies,
(3) tipping point ID matrix/list, and (4) tipping energies matrix/list.
Pruning shallow stable states of energy landscape.
Iteratively prunes stable states (merges basins).
If type="relative", stable state whose energy depth is smaller than th times the deepest
basin depth are pruned and merged to the deeper basins.
If type="bootstrap", stable state whose bootstrap confidence interval (default: 95
cross zero is pruned.
The function updates stable states and corresponding tipping
structures, and returns both the pruned ELA object and a mapping table
from stable states before pruning to those after pruning.
ELPruning( ela, type = "relative", th = 0.05, lower_qr = 0.05, bootstrap_ene = NULL, threads = 1, reporting = TRUE, fork = TRUE )ELPruning( ela, type = "relative", th = 0.05, lower_qr = 0.05, bootstrap_ene = NULL, threads = 1, reporting = TRUE, fork = TRUE )
ela |
output of ELA function. |
type |
the type of pruning |
th |
ratio of the depth of the basin to be pruned with respect to the deepest basin. |
lower_qr |
Numeric value. Quantile used for constructing confidence intervals from bootstrap samples to check
whether the estimated energy gaps include zero. Default is |
bootstrap_ene |
a data frame of bootstrap energies generated by bootstrap_ELA. Necessary when |
threads |
number of cores (threads) used for calculation. |
reporting |
enable/disable the reporting function. |
fork |
logical; if TRUE and |
A list of two elements: (1) pruned ELA tuple (same format as ELA output),
and (2) a data frame with columns ss.before.pruning and ss.after.pruning.
Calculate energy of a given state.
Computing energy in state using given parameters,alpha and beta
based on the (extended) maximum pairwise entropy model.
Energy(state, alpha, beta)Energy(state, alpha, beta)
state |
Binary vector representing a state (community composition). |
alpha |
vector of implicit/explicit environmental effects on species. |
beta |
Matrix of species interactions. |
A numeric scalar giving the energy of the state.
Repeating cross-validation over a grid of hyperparameters for simulated annealing-based fitting, evaluates prediction/fitting errors on held-out samples, and returns candidate best parameter/timepoint combinations.
Findbp( ocmat, enmat = NULL, cvmode = c("10fold", "5fold", "loo"), nrepeat = 5, we = c(0.001), maxlrs = c(0.005), totalit = 4000, lmd = c(5e-04, 0.001, 0.0025, 0.005, 0.0075), rep = 8, intv = 100, threads = 1, tol = 0.03, runadamW = TRUE, sparse = TRUE, fastfitting = TRUE, reporting = TRUE )Findbp( ocmat, enmat = NULL, cvmode = c("10fold", "5fold", "loo"), nrepeat = 5, we = c(0.001), maxlrs = c(0.005), totalit = 4000, lmd = c(5e-04, 0.001, 0.0025, 0.005, 0.0075), rep = 8, intv = 100, threads = 1, tol = 0.03, runadamW = TRUE, sparse = TRUE, fastfitting = TRUE, reporting = TRUE )
ocmat |
Numeric matrix (0/1). occurrence matrix of the communities. |
enmat |
Optional environmental parameter matrix. |
cvmode |
Cross-validation mode. One of '"10fold"', '"5fold"', or '"loo"'. '"loo"' performs leave-one-out cross-validation. Default is '"10fold"'. |
nrepeat |
Number of repeated cross-validation runs for k-fold CV. Ignored when 'cvmode = "loo"'. Default is '10'. |
we |
Numeric vector. Candidate weight decay values. |
maxlrs |
Numeric vector. Maximum learning rate. |
totalit |
Integer. Total iterations for SA. |
lmd |
Numeric vector. Candidate |
rep |
Integer. Number of replicates for runSA. |
intv |
Integer. Recording interval. |
threads |
Integer. Number of cores. |
tol |
Numeric. Error tolerance for finding the best hyperparameters (default: 0.03) |
runadamW |
Logical. Use AdamW. |
sparse |
Logical. Use sparsity penalty. |
fastfitting |
Logical. If |
reporting |
Logical. Print progress. |
A list with:
A named list of best hyperparameter settings and scores.
All validation results for each grid point.
Generate occurrence(binary)/abundance/factor matrices and labels from input abundance/metadata tables.
Formatting( baseabtable, basemetadata = NULL, normalize = 1, parameters = c(0.05, 0.05, 0.95), grouping = 0, grouping_th = 0, reporting = TRUE )Formatting( baseabtable, basemetadata = NULL, normalize = 1, parameters = c(0.05, 0.05, 0.95), grouping = 0, grouping_th = 0, reporting = TRUE )
baseabtable |
Abundance table (samples x species) as data.frame/matrix. |
basemetadata |
Optional environmental metadata table (samples x factors); NULL/empty for none. |
normalize |
Integer 0 or 1. If 1, normalize within each sample to sum to 1. |
parameters |
Numeric vector/list of thresholds (following 3 parameters). 1. "0/1 threshold" 2. "min occurrence threshold" 3. "max occurrence threshold" (e.g., c(0.05, 0.05, 0.95) : species abundance > 0.05 is regarded as present, and species occurring in more than 0.05 and less than 0.95 in the all samples are kept.) |
grouping |
Integer 0 or 1. If 1, group similar occurrence patterns. |
grouping_th |
Numeric threshold for grouping (relative Hamming distance). |
reporting |
Logical. If TRUE, print processing summary. |
A list: (ocmatrix, abmatrix, fctmatrix, samplelabel, specieslabel, factorlabel).
Construct smoothed energy landscape surface objects from GELA output and parameters.
GELAObj(gela, sa, threads = 1)GELAObj(gela, sa, threads = 1)
gela |
GELA object/list |
sa |
Parameter list generated by runSA |
threads |
Number of parallel workers. |
A list containing surface grid data and state-specific surface points.
Estimate energy landscape across environmental gradients.
Runs ELA and ELPruning comprehensively across specified environmental parameters.
One environmental parameter can be changed by specifying eid over steps within range.
If env is not given, the mean vector of enmat is used as the reference environmental condition.
GradELA( sa = sa, pruning_type = "relative", eid = NULL, enmat = enmat, env = NULL, range = NULL, steps = 16, th = 0.05, threads = 1, SS.itr = 20000, FindingTip.itr = 10000, lower_qr = 0.05, bs_params = NULL, ocmat = NULL )GradELA( sa = sa, pruning_type = "relative", eid = NULL, enmat = enmat, env = NULL, range = NULL, steps = 16, th = 0.05, threads = 1, SS.itr = 20000, FindingTip.itr = 10000, lower_qr = 0.05, bs_params = NULL, ocmat = NULL )
sa |
output of runSA function. |
pruning_type |
the type of pruning |
eid |
target environmental factor: the specified environmental factor to be changed. |
enmat |
array/matrix of environmental factors for each sample (normalization required). |
env |
array of environmental factor values.
For the environmental parameters that are not specified for |
range |
range over which the environmental factor specified by eid is changed. |
steps |
number of steps in the range. |
th |
ratio of the depth of the basin to be pruned with respect to the deepest basin. |
threads |
number of cores (threads) used for calculation. |
SS.itr |
number of iterations to identify stable states. |
FindingTip.itr |
number of iterations to identify tipping points. |
lower_qr |
Numeric value. Quantile used for constructing confidence intervals from bootstrap samples to check
whether the estimated energy gaps include zero. Default is |
bs_params |
output of bootstrap_SA. list of estimated parameters using bootstrap resampled data set.
necessary when |
ocmat |
community composition matrix. |
A list of three elements: (1) ELA results for each step (pruned), (2) the sequence of environmental values, and (3) the reference env vector.
Estimates parameters h and J by iterative gradient updates to match
first and second moments of the observed binary data under an maximum
pairwise entropy model.
GradientDescent(data, maxit)GradientDescent(data, maxit)
data |
Numeric matrix (typically 0/1). Rows are samples and columns are species or taxa. |
maxit |
Integer. Maximum number of gradient descent iterations. |
This function calls Energy() to compute energies for candidate states.
Ensure that Energy is available in the package namespace.
A list with elements:
Numeric vector of estimated implicit environmental parameters.
Numeric matrix of estimated interaction matrix.
Calculate energy gap, energy barrier and stable state entropy.
This function evaluates stability indices for both (i) the pruned and (ii) the non-pruned
energy landscape. For each sample, it assigns the basin by Bi and computes:
- energy gap (sample energy minus basin energy),
- energy barrier (minimum tipping energy connected to the basin minus basin energy),
- stable state entropy (SSentropy).
For basin connectivity: when only one stable state exists, tipping energy is treated as Inf
and therefore energy barrier is returned as Inf.
gStability( sa, ocmat, enmat = NULL, seitr = 1000, seconvTime = 10000, pruning_type = "relative", th = 0.1, lower_qr = 0.05, threads = 1, reporting = TRUE, bs_params = NULL )gStability( sa, ocmat, enmat = NULL, seitr = 1000, seconvTime = 10000, pruning_type = "relative", th = 0.1, lower_qr = 0.05, threads = 1, reporting = TRUE, bs_params = NULL )
sa |
: output of runSA function. |
ocmat |
: matrix of presence/absence status for each sample. |
enmat |
: matrix of environmental factors for each sample. Normalization required. |
seitr |
: Number of stochastic search iterations. |
seconvTime |
: Maximum number of iterations for each search. |
pruning_type |
the type of pruning |
th |
: Numerical (0~1). Ratio of the depth of the basin to be pruned with respect to the deepest basin in each pair. |
lower_qr |
Numeric value. Quantile used for constructing confidence intervals from bootstrap samples to check
whether the estimated energy gaps include zero. Default is |
threads |
number of cores (threads) used for calculation.#' @param threads : number of cores (threads) used for calculation. |
reporting |
: enable/disable the reporting function. |
bs_params |
output of bootstrap_SA. list of estimated parameters using bootstrap resampled data set,
necessary when |
A list of four elements: (1) pruned stability data frame, (2) non-pruned stability data frame, (3) auxiliary object for pruned landscape (parameters + stable states table), and (4) auxiliary object for sample-wise evaluation (inputs and ELA objects).
Generates random implicit environmental parameters h and a interaction matrix J.
Optionally generates explicit environmental coefficients g.
The j_connectivity parameter controls the fraction of non-zero off-diagonal interactions.
hb.paramgen(nn, ne = 0, j_connectivity = 1)hb.paramgen(nn, ne = 0, j_connectivity = 1)
nn |
Integer. Number of units (species, taxa, or members). |
ne |
Integer. Number of environmental factors. If |
j_connectivity |
Numeric in [0, 1]. Fraction of possible pairwise interactions set to non-zero. |
A list containing at least:
Named numeric vector of length nn. Implicit environmental parameters
Named numeric matrix (nn x nn), symmetric. Interaction matrix
If ne > 0, also includes:
Named numeric matrix (nn x ne). Explicit environmental parameters
Generates samples of binary state vectors using heat-bath updates under
an extended maximum pairwise entropy model parameterized
by implicit environmental parameters h and interaction matrix J.
Optionally includes environmental effects g, in which case an environmental
matrix is generated internally.
HeatBath(steps, nproc, h, J, g = NULL)HeatBath(steps, nproc, h, J, g = NULL)
steps |
Integer. Number of heat-bath iterations to run. |
nproc |
Integer. Number of parallel samples (rows) to maintain. |
h |
Numeric vector of length |
J |
Numeric matrix of size |
g |
Optional numeric matrix of size |
A list with components:
Numeric matrix (nproc x nn) of sampled states (0/1).
Numeric matrix (nproc x ne) of environmental values, or NULL
if g is NULL.
set.seed(1) nn <- 5 pars <- hb.paramgen(nn) res <- HeatBath(steps = 10, nproc = 8, h = pars[[1]], J = pars[[2]]) str(res)set.seed(1) nn <- 5 pars <- hb.paramgen(nn) res <- HeatBath(steps = 10, nproc = 8, h = pars[[1]], J = pars[[2]]) str(res)
Decode a base-64 id string index to a binary vector of length s.
id2bin(id, s)id2bin(id, s)
id |
A character scalar encoded by bin2id. |
s |
Integer number of bits to return. |
An integer vector of 0/1 bits of length s.
Compute a minimum “tipping” energy path between two states.
MinTippingPath(s1, s2, hg, j, tmax)MinTippingPath(s1, s2, hg, j, tmax)
s1 |
Another integer state ID |
s2 |
Another integer state ID |
hg |
explicit/implicit environmental factors |
j |
species interaction matrix |
tmax |
Iteration limit. |
Numeric vector of energies along the selected path (or scalar if s1==s2).
Updates each row (sample) of a binary state matrix by choosing one random position per row and resampling it according to precomputed heat-bath probabilities.
OnestepHBS(y, logmo)OnestepHBS(y, logmo)
y |
Numeric matrix of size |
logmo |
Numeric matrix of probabilities compatible with |
A numeric matrix with the same dimensions as y, updated states.
Performs PCA on the community composition matrix (ocmat) and plots the first
two PCs. Points are colored by the stable state (basin) label estimated
from Bi. If pruned = TRUE, colors correspond to the
pruned/merged basin labels in ssrep (stablestates2), otherwise to the
original labels (stablestates). When export = TRUE, the plotting data
frame (stable-state labels + PCs + merged labels) is returned.
PCplot(ocmat, sa = sa, env = NULL, ssrep = NULL, pruned = TRUE, export = FALSE)PCplot(ocmat, sa = sa, env = NULL, ssrep = NULL, pruned = TRUE, export = FALSE)
ocmat |
Data frame (or matrix) of community compositions (rows = samples, columns = species/features). |
sa |
Parameter set passed to sa2params. |
env |
Optional environment object passed to sa2params. |
ssrep |
Data frame-like object providing stable state label mapping with
two columns: |
pruned |
Logical; if TRUE, use |
export |
Logical; if TRUE, return the data used for plotting. |
NULL invisibly (or a data frame when export = TRUE).
Plot mean validation scores over time for the results returned by
Findbp, including 95% confidence intervals.
plotSAtest(sa_results, ylim = c("auto", "auto"))plotSAtest(sa_results, ylim = c("auto", "auto"))
sa_results |
Output list from Findbp. |
ylim |
Numeric length-2 vector (ymin, ymax). Default: |
This function uses ggplot2 and tidyverse-style data manipulation. Ensure required packages are installed.
A ggplot2 object.
Estimate species interaction parameters and (optionally) environmental
preference parameters by stochastic approximation (SA). The function runs
multiple SA optimizations in parallel (when threads > 1) using
foreach.
runSA( ocmat, enmat = NULL, rep = 128, threads = 1, we = 0.001, totalit = 1000, lambda = 0.01, intv = 100, runadamW = TRUE, maxlr = 0.005, sparse = TRUE, getall = FALSE, reporting = TRUE )runSA( ocmat, enmat = NULL, rep = 128, threads = 1, we = 0.001, totalit = 1000, lambda = 0.01, intv = 100, runadamW = TRUE, maxlr = 0.005, sparse = TRUE, getall = FALSE, reporting = TRUE )
ocmat |
Numeric matrix (0/1). Rows are samples and columns are species. |
enmat |
Optional numeric matrix. Environmental variables for each sample. Should be normalized/standardized beforehand. |
rep |
Integer. Number of parallel runs (replicates) for optimization. |
threads |
Integer. Number of cores used when running in parallel. |
we |
Numeric. Weight decay parameter for AdamW. |
totalit |
Integer. Total number of optimization iterations. |
lambda |
Numeric. Sparsity penalty parameter (used when |
intv |
Integer. Interval of iterations at which intermediate results are recorded. |
runadamW |
Logical. If |
maxlr |
Numeric. Maximum learning rate. |
sparse |
Logical. If |
getall |
Logical. If |
reporting |
Logical. If |
If getall = TRUE, a list of SA trajectories (length rep).
Otherwise, a list with two elements, c(params, enmat) is returned.
params is the concatenated matrix of the estimated parameter set:
Implicit environmental effects on species occurrence
Species-species interaction matrix
(If enmat is not NULL) Explicit effects of the environmental factors.
Convert the output of runSA / saEndpoint into parameter
objects in (extended) maximum pairwise entropy model:
h, J, and optionally g when environmental parameter
matrix is specified.
sa2params(sa, env = NULL)sa2params(sa, env = NULL)
sa |
A list in the form |
env |
Optional numeric vector of environmental values (column names correspond to environmental factors).Default is NULL. |
A list with elements he (h), je (J), ge (g or NULL),
and hge (combined field).
Summarize the final parameter estimates across SA replicates by averaging the last recorded block of parameters.
saEndpoint(fittingMat, ocmat, enmat = NULL)saEndpoint(fittingMat, ocmat, enmat = NULL)
fittingMat |
A list of SA trajectories returned by |
ocmat |
Numeric matrix (0/1). Rows are samples and columns are species. |
enmat |
Optional environmental parameter matrix used in fitting. |
A numeric matrix of median endpoint parameters across serials. Rownames correspond to samples and Colnames correspond to species list
Computing graph objects from ela (and optionally
bootstrap energies) and draw the disconnectivity graph for a given community
composition matrix oc. More than two stable state are required.
showDG( ela, oc, label = "", bootstrap_ene = NULL, min_pval = 0.05, plot = TRUE, scale_labels = c(0.05, 0.1, 0.2, 0.5), getGraphObj = FALSE )showDG( ela, oc, label = "", bootstrap_ene = NULL, min_pval = 0.05, plot = TRUE, scale_labels = c(0.05, 0.1, 0.2, 0.5), getGraphObj = FALSE )
ela |
ELA output tuple/list used by |
oc |
Community composition matrix/data frame (used to infer |
label |
Character; plot title. |
bootstrap_ene |
Optional bootstrap energy matrix/data frame for computing statistical significance levels of energy gaps. |
min_pval |
Numeric; minimum p-value shown in the visualization. This is usually set to the minimum detection level in bootstrap sampling (e.g., if the number of bootstrap iterations is 1000, the minimum level is 1/1000 = 0.001). |
plot |
Logical; if TRUE, show disconnectivity graph |
scale_labels |
Numeric vector of values between 0 and 1; labels for the p-value scale. These values are converted to log-scaled values for plotting. |
getGraphObj |
Logical; if TRUE, returns data.frame object used for plotting DisconnectivityGraph Users can manually change the graphical configurations. |
A data frame containing energy values, xpositions, pvalues for energy gaps in identified stable states and tipping points.
Draws a 3D surface of the energy landscape and overlays line trajectories using plot3D function.
showGELA3D(gelsobj, theta = 25, phi = 50)showGELA3D(gelsobj, theta = 25, phi = 50)
gelsobj |
List with two elements: surface data frame and line data frame. |
theta |
numeric value. plot3D option: Azimuth angle (in degrees) controlling rotation around the z-axis. |
phi |
numeric value. plot3D option: Elevation angle (in degrees) controlling the vertical viewing position. |
NULL (produces a 3D plot as side effect).
Computes a 2D PCoA embedding from the interaction matrix (jest) and overlays
interaction links above a threshold th. Positive and negative links are
drawn separately, and species names are annotated at their coordinates.
showIntrGraph(ela, sa, env = NULL, th = 0, annot_adj = c(0.75, 2))showIntrGraph(ela, sa, env = NULL, th = 0, annot_adj = c(0.75, 2))
ela |
ELA output tuple/list (used to obtain stable states etc.; currently mainly used for context). |
sa |
Parameter set passed to sa2params. |
env |
target environmental factor. |
th |
Numeric threshold for drawing links (absolute value). |
annot_adj |
Numeric length-2; horizontal/vertical adjustment for labels. |
NULL.
Plots energy values of stable states across a factor (e.g., environmental or
control parameter) using the output from gela. Each unique stable state is
drawn as a separate line with points, and a legend is added.
showSSD(gela)showSSD(gela)
gela |
The output of GradELA. Object containing stable state lists, energies, the factor, and environment labels as used in the code. |
NULL (produces a plot as side effect).
Plots energy values of stable states across a factor (e.g., environmental or
control parameter) using the output from gela. Each unique stable state is
drawn as a separate line with points, and a legend is added.
showSSD_ggplot(gela, plot = TRUE, getSSDobj = TRUE)showSSD_ggplot(gela, plot = TRUE, getSSDobj = TRUE)
gela |
The output of GradELA. Object containing stable state lists, energies, the factor, and environment labels as used in the code. |
plot |
Logical; if TRUE, show SSD plot |
getSSDobj |
Logical; if TRUE, returns data.frame object used for plotting SSD diagram Users can manually change the graphical configurations. |
A data frame used for plotting SSD diagram.
Estimate stable state entropy for a given state.
Evaluating the entropy of transitions among stable states starting from a given community state.
The returned value includes the estimated entropy, the convergence time, and the number of successful convergences within the given iteration limits.
SSentropy(state, ss, alpha, beta, seitr = 1000, convTime = 10000)SSentropy(state, ss, alpha, beta, seitr = 1000, convTime = 10000)
state |
Binary vector (or matrix with one row) representing the current state. |
ss |
Matrix of stable states. |
alpha |
Vector of implicit or explicit environmental effects. |
beta |
Matrix of species interactions. |
seitr |
Number of stochastic search iterations. |
convTime |
Maximum number of iterations for each search. |
A numeric matrix with three columns: (1) stable state entropy, (2) convergence time, (3) number of converged searches.
Summary of stable states.
Builds a table of stable states and their energies from an ELA output.
The returned data frame contains stable state IDs and the corresponding binary
presence/absence vector (per species). This function is typically used for
inspecting stable states after pruning (if ela is pruned) or before pruning.
sstable(ela, ocmat)sstable(ela, ocmat)
ela |
output of ELA or the first return value of ELPruning. |
ocmat |
community composition matrix. |
A data frame with columns ID, Energy.
Calculate energy gap and stable state entropy as stability indices and returns a dataframe.
For each sample state, this function identifies the basin (stable state) by Bi,
then computes (i) the energy gap and (ii) stable state entropy
estimated by SSentropy.
If enmat is provided and sa was estimated with environmental
factors, the basin assignment is computed for each sample-specific env.
Note: the returned stable state id corresponds to the basin (stable state) the sample is assigned to.
Stability( sa, ocmat, enmat = NULL, seitr = 1000, seconvTime = 10000, threads = 1, reporting = TRUE )Stability( sa, ocmat, enmat = NULL, seitr = 1000, seconvTime = 10000, threads = 1, reporting = TRUE )
sa |
: output of runSA function. |
ocmat |
: array of the vectors of presence/absence status for each sample. |
enmat |
: array of the vectors of environmental factors for each sample. Normalization required. |
seitr |
: Number of stochastic search iterations. |
seconvTime |
: Maximum number of iterations for each search. |
threads |
: number of cores (threads) used for calculation. |
reporting |
: enable/disable the reporting function. |
A data frame with stability indices. Columns include at least
energy.gap, ss.entropy, e.realize, e.stable, state.id, and stable.state.id
(and may include conv_time and convergence depending on the branch).
Perform steepest descent on the energy landscape.
Starting from an initial binary state, the steepest descent algorithm
will find a local minimum of the energy landscape.
The computation is delegated to the C++ implementation
SteepestDescent_cpp.
SteepestDescent(state, alpha, beta)SteepestDescent(state, alpha, beta)
state |
Binary vector representing the initial state. |
alpha |
vector of implicit/explicit environmental effects. |
beta |
Matrix of species interactions. |
A numeric vector of the reached binary stable state and its energy.
Summary of tipping points.
Builds a table of tipping points (transition states) and their energies from an ELA output.
This function lists only finite tipping points
(i.e., entries not equal to Inf), and returns the corresponding binary state vectors.
Two stable states are required at least.
tptable(ela, ocmat)tptable(ela, ocmat)
ela |
|
ocmat |
community composition matrix. |
A data frame with columns TP, SS1, SS2, Energy,
and one column per species (binary state vector).
If no tipping points exist, returns NULL (invisibly) after printing a message.