Package {rEDM}

Type:	Package
Title:	Empirical Dynamic Modeling ('EDM')
Version:	2.0.1
Date:	2026-07-27
Maintainer:	Joseph Park <JosephPark@IEEE.org>
Description:	An implementation of 'EDM' algorithms based on research software developed at the Sugihara Lab ('UCSD/SIO'). Primary methods include 'Simplex' projection from Sugihara & May (1990) <doi:10.1038/344734a0>, Sequential locally-weighted global linear maps 'S-map': Sugihara (1994) <doi:10.1098/rsta.1994.0106>, Convergent cross mapping described in Sugihara et al. (2012) <doi:10.1126/science.1227079>, and, 'Multiview embedding' from Ye & Sugihara (2016) <doi:10.1126/science.aag0863>.
License:	BSD_2_clause + file LICENSE
Encoding:	UTF-8
LazyData:	true
LazyLoad:	yes
Imports:	stats, utils, graphics, parallel, RANN
Suggests:	testthat (≥ 3.0.0), knitr, rmarkdown, formatR
Config/testthat/edition:	3
RoxygenNote:	8.0.0
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2026-07-30 17:38:08 UTC; jpark
Author:	Joseph Park [aut, cre], George Sugihara [aut, ccp], Ethan Deyle [aut], The Regents of the University of California [cph]
Repository:	CRAN
Date/Publication:	2026-07-31 00:00:30 UTC

Empirical dynamic modeling

Description

rEDM provides tools for Empirical Dynamic Modeling time series analyses. It is based on reconstructing multivariate state space representations from uni- or multivariate time series, then projecting state changes using various metrics applied to nearest neighbors.

Functionality includes:

• Simplex projection (Sugihara and May 1990)

• Sequential Locally Weighted Global Linear Maps (S-map) (Sugihara 1994)

• Multivariate embeddings (Dixon et. al. 1999)

• Convergent cross mapping (Sugihara et. al. 2012)

• Multiview embedding (Ye and Sugihara 2016)

Details

Main Functions:

• Simplex - simplex projection

• SMap - S-map projection

• CCM - convergent cross mapping

• Multiview - multiview forecasting

Helper Functions:

• Embed - time delay embedding

• ComputeError - forecast skill metrics

• EmbedDimension - optimal embedding dimension

• PredictInterval - optimal prediction interval

• PredictNonlinear - evaluate nonlinearity

Author(s)

Maintainer: Joseph Park

Authors: Joseph Park, Ethan Deyle, George Sugihara

References

Sugihara G. and May R. 1990. Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series. Nature, 344:734-741.

Sugihara G. 1994. Nonlinear forecasting for the classification of natural time series. Philosophical Transactions: Physical Sciences and Engineering, 348 (1688) : 477-495.

Dixon, P. A., M. Milicich, and G. Sugihara, 1999. Episodic fluctuations in larval supply. Science 283:1528-1530.

Sugihara G., May R., Ye H., Hsieh C., Deyle E., Fogarty M., Munch S., 2012. Detecting Causality in Complex Ecosystems. Science 338:496-500.

Ye H., and G. Sugihara, 2016. Information leverage in interconnected ecosystems: Overcoming the curse of dimensionality. Science 353:922-925.

Convergent cross mapping

Description

The state-space of a multivariate dynamical system (not a purely stochastic one) encodes coherent phase-space variable trajectories. If enough information is available, one can infer the presence or absence of cross-variable interactions associated with causal links between variables. CCM measures the extent to which states of variable Y can reliably estimate states of variable X. This can happen if X is causally influencing Y.

If cross-variable state predictability converges as more state-space information is provided, this indicates a causal link. CCM performs this cross-variable mapping using Simplex, with convergence assessed across a range of observational library sizes as described in Sugihara et al. 2012.

Usage

CCM(dataFrame = NULL, columns, target, E, Tp = 0, knn = 0, tau = -1,
    exclusionRadius = 0, libSizes, sample = 30, seed = NULL,
    embedded = FALSE, validLib = logical(0), includeData = FALSE,
    numProcess = 4, backend = "RANN", pathIn = "./", dataFile = "",
    pathOut = "./", predictFile = "", parameterList = FALSE, showPlot = FALSE)

Arguments

Value

A data.frame with 3 columns. The first column is LibSize specifying the subsampled library size. Columns 2 and 3 report Pearson correlation coefficients for the prediction of X from Y, and Y from X.

If includeData = TRUE and parameterList = TRUE a named list "parameters" is added.

References

Sugihara G., May R., Ye H., Hsieh C., Deyle E., Fogarty M., Munch S., 2012. Detecting Causality in Complex Ecosystems. Science 338:496-500.

Examples

data(sardine_anchovy_sst)
df = CCM( dataFrame = sardine_anchovy_sst, E = 3, Tp = 0, columns = "anchovy",
target = "np_sst", libSizes = "10 70 10", sample = 100 )

Prediction skill statistics

Description

Pearson rho, MAE, CAE and RMSE between observations and predictions.

Usage

ComputeError(obs, pred, digits = 6)

Arguments

Value

A named list with elements rho, MAE, CAE and RMSE.

Examples

ComputeError(1:10 + 0, c(1:9 + 0, 8))

Time-delay embedding

Description

Takens time-delay embedding of one or more columns.

Usage

Embed(dataFrame, E, tau, columns)

Arguments

Details

Each columns item will have E-1 time-lagged vectors created. The column name is appended with (t-n). For example, data columns X, Y, with E = 2 will have columns named X(t-0) X(t-1) Y(t-0) Y(t-1).

The returned data.frame does not have a time column. The returned data.frame is truncated by tau * (E-1) rows to remove state vectors with partial data (NaN elements).

Value

A numeric matrix of the embedding (NA in the first/last |tau|*(E-1) rows).

Examples

Embed(data.frame(x = 1:10), E = 3, tau = -1, columns = "x")

Prediction skill versus embedding dimension

Description

EmbedDimension uses Simplex to evaluate prediction accuracy as a function of embedding dimension.

Usage

EmbedDimension(dataFrame = NULL, columns, target, lib, pred, maxE = 10, Tp = 1,
               tau = -1, exclusionRadius = 0, embedded = FALSE,
               validLib = logical(0), noTime = FALSE, ignoreNan = TRUE,
               numProcess = 4, pathIn = "./", dataFile = "", pathOut = "./",
               predictFile = "", backend = "RANN",
               parameterList = FALSE, showPlot = FALSE)

Arguments

Value

A data.frame with columns E, rho.

Examples

data(TentMap)
E.rho = EmbedDimension( dataFrame = TentMap, lib = "1 100", pred = "201 500",
columns = "TentMap", target = "TentMap", showPlot = FALSE )

Water flow to NE Everglades

Description

Cumulative weekly water flow into northeast Everglades from water control structures S12C, S12D and S333 from 1980 through 2005.

Usage

EvergladesFlow

Format

A data frame with 1379 rows and 2 columns:

Date

Date.

S12CD_S333_CFS

Cumulative weekly flow (CFS).

5-D Lorenz'96

Description

5-D Lorenz'96 timeseries with F = 8.

Usage

Lorenz5D

Format

Data frame with 1000 rows and 6 columns

Time

Time.

V1

variable 1.

V2

variable 2.

V3

variable 3.

V4

variable 4.

V5

variable 5.

References

Lorenz, Edward (1996). Predictability - A problem partly solved, Seminar on Predictability, Vol. I, ECMWF.

Multiview embedding forecast

Description

Multiview embedding forecast (Ye & Sugihara 2016): rank D-column views and average the best to find optimal combinations of variables that best represent the dynamics.

Usage

Multiview(dataFrame = NULL, columns, target, lib = numeric(0), pred = numeric(0),
          D = 0, E = 1, Tp = 1, knn = 0, tau = -1, multiview = 0,
          exclusionRadius = 0, trainLib = TRUE, excludeTarget = FALSE,
          ignoreNan = TRUE, numProcess = 4, backend = "RANN",
          pathIn = "./", dataFile = "", pathOut = "./", predictFile = "",
          parameterList = FALSE, showPlot = FALSE)

Arguments

Details

Multiview embedding is a method to identify variables in a multivariate dynamical system that are most likely to contribute to the observed dynamics. It is a multistep algorithm with these general steps:

1. Compute D-dimensional variable combination forecasts.

2. Rank forecasts.

3. Compute predictions of top combinations.

4. Compute multiview averaged prediction.

If E>1, all variables are embedded to dimension E. If trainLib is TRUE initial forecasts and ranking are done in-sample (lib=pred) and predictions using the top ranked combinations use the specified lib and pred. If trainLib is FALSE initial forecasts and ranking use the specified lib and pred, the step of computing predictions of the top combinations is skipped.

Value

Named list with data.frames [[View, Predictions]].

data.frame View columns:

If parameterList = TRUE a named list "parameters" is added.

References

Ye H., and G. Sugihara, 2016. Information leverage in interconnected ecosystems: Overcoming the curse of dimensionality. Science 353:922-925.

Examples

data(block_3sp)
L = Multiview( dataFrame = block_3sp, lib = "1 100", pred = "101 190",
E = 2, columns = "x_t y_t z_t", target = "x_t" )

CCM cross-map skill versus library size, both directions.

Description

CCM cross-map skill versus library size, both directions.

Usage

PlotCCM(df, E = NULL)

Arguments

Value

Invisibly NULL.

Plot observed versus predicted values

Description

Base-graphics plot of Observations and Predictions against time.

Usage

PlotObsPred(df, dataName = "", E = NULL, Tp = NULL)

Arguments

Value

Invisibly NULL.

Plot S-map predictions and coefficients

Description

Base-graphics multi-panel plot of S-map predictions and coefficient series.

Usage

PlotSmap(SmapList, dataName = "", E = NULL, Tp = NULL)

Arguments

Value

Invisibly NULL.

rho versus a swept parameter (EmbedDimension / PredictInterval / ...).

Description

rho versus a swept parameter (EmbedDimension / PredictInterval / ...).

Usage

PlotSweep(df, xlab, title = NULL)

Prediction skill versus exclusion radius

Description

Simplex prediction skill (rho) swept over the temporal exclusion radius.

Usage

PredictExclusionRadius(dataFrame = NULL, columns, target, lib, pred,
                       exclusionRadius = NULL, E = 1, Tp = 1, tau = -1,
                       embedded = FALSE, validLib = logical(0),
                       noTime = FALSE, ignoreNan = TRUE, numProcess = 4,
                       backend = "RANN", pathIn = "./", dataFile = "",
                       pathOut = "./", predictFile = "",
                       parameterList = FALSE, showPlot = FALSE)

Arguments

Value

A data.frame with columns ExclusionRadius and rho.

Prediction skill versus forecast interval

Description

PredictInterval uses Simplex to evaluate prediction accuracy as a function of forecast interval Tp.

Usage

PredictInterval(dataFrame = NULL, columns, target, lib, pred, maxTp = 10, E = 1,
                tau = -1, exclusionRadius = 0, embedded = FALSE,
                validLib = logical(0), noTime = FALSE, ignoreNan = TRUE,
                numProcess = 4, backend = "RANN", pathIn = "./",
                dataFile = "", pathOut = "./", predictFile = "",
                parameterList = FALSE, showPlot = FALSE)

Arguments

Value

A data.frame with columns Tp, rho.

Examples

data(TentMap)
Tp.rho = PredictInterval( dataFrame = TentMap, lib = "1 100",
pred = "201 500", E = 2, columns = "TentMap", target = "TentMap",
showPlot = FALSE )

Prediction skill versus S-map localisation

Description

PredictNonlinear uses SMap to evaluate prediction accuracy as a function of the state-space localization parameter theta.

Usage

PredictNonlinear(dataFrame = NULL, columns, target, lib, pred, theta = NULL, E = 1,
                 Tp = 1, knn = 0, tau = -1, exclusionRadius = 0,
                 embedded = FALSE, validLib = logical(0), noTime = FALSE,
                 ignoreNan = TRUE, numProcess = 4, backend = "RANN",
                 pathIn = "./", dataFile = "", pathOut = "./", predictFile = "",
                 parameterList = FALSE, showPlot = FALSE, verbose = FALSE)

Arguments

Value

A data.frame with columns Theta and rho.

Examples

data(TentMapNoise)
theta.rho = PredictNonlinear( dataFrame = TentMapNoise, E = 2,
lib = "1 100", pred = "201 500", columns = "TentMap",
target = "TentMap", showPlot = FALSE )

Sequential locally weighted global linear map (S-map)

Description

SMap performs time series forecasting based on localised (or global) nearest neighbor projection in the time series phase space as described in Sugihara 1994.

Usage

SMap(dataFrame = NULL, columns, target, lib, pred, E = 0, Tp = 1, knn = 0,
     tau = -1, theta = 0, exclusionRadius = 0, embedded = FALSE,
     validLib = logical(0), noTime = FALSE, ignoreNan = TRUE,
     backend = "RANN", pathIn = "./", dataFile = "", pathOut = "./",
     predictFile = "", parameterList = FALSE, verbose = FALSE, showPlot = FALSE)

Arguments

Details

If embedded is FALSE, the data column(s) are embedded to dimension E with time lag tau. This embedding forms an n-columns * E-dimensional phase space for the SMap projection. If embedded is TRUE, the data are assumed to contain an E-dimensional embedding with E equal to the number of columns. See the Note below for proper use of multivariate data (number of columns > 1).

If ignoreNan is TRUE, the library (lib) is internally redefined to exclude nan embedding vectors. If ignoreNan is FALSE no library adjustment is made. The (lib) can be explicitly specified to exclude nan library vectors.

Predictions are made using leave-one-out cross-validation, i.e. observation rows are excluded from the prediction regression.

In contrast to Simplex, SMap uses all available neighbors and weights them with an exponential decay in phase space distance with exponent theta. theta=0 uses all neighbors corresponding to a global autoregressive model. As theta increases, neighbors closer in vicinity to the observation are considered.

Value

A named list with three data.frames [[predictions, coefficients, singularValues]]. predictions has columns Observations, Predictions. The first column contains time or index values.

coefficients data.frame has time or index values in the first column. Columns 2 through E+2 (E+1 columns) are the SMap coefficients.

singularValues data.frame has time or index values in the first column. Columns 2 through E+2 (E+1 columns) are the SVD singularValues. The first value corresponds to the SVD bias (intercept) term.

If parameterList = TRUE a named list "parameters" is added.

Note

SMap should be called with columns explicitly corresponding to dimensions E. In the univariate case (number of columns = 1) with default embedded = FALSE, the time series will be time-delay embedded to dimension E, SMap coefficients correspond to each dimension.

If a multivariate data set is used (number of columns > 1) it must use embedded = TRUE with E equal to the number of columns. This prevents the function from internally time-delay embedding the multiple columns to dimension E. If the internal time-delay embedding is performed, then state-space columns will not correspond to the intended dimensions in the matrix inversion, coefficient assignment, and prediction. In the multivariate case, the user should first prepare the embedding (using Embed for time-delay embedding), then pass this embedding to SMap with appropriately specified columns, E, and embedded = TRUE.

References

Sugihara G. 1994. Nonlinear forecasting for the classification of natural time series. Philosophical Transactions: Physical Sciences and Engineering, 348 (1688):477-495.

Simplex projection

Description

Simplex performs time series forecasting based on weighted nearest neighbors projection in the time series phase space as described in Sugihara and May.

Usage

Simplex(dataFrame = NULL, columns, target, lib, pred, E = 0, Tp = 1, knn = 0,
        tau = -1, exclusionRadius = 0, embedded = FALSE,
        validLib = logical(0), noTime = FALSE, ignoreNan = TRUE,
        backend = "RANN", pathIn = "./", dataFile = "", pathOut = "./",
        predictFile = "", parameterList = FALSE, verbose = FALSE, showPlot = FALSE)

Arguments

Details

If embedded is FALSE, the data column(s) are embedded to dimension E with time lag tau. This embedding forms an E-dimensional phase space for the Simplex projection. If embedded is TRUE, the data are assumed to contain an E-dimensional embedding with E equal to the number of columns. Predictions are made using leave-one-out cross-validation, i.e. observation vectors are excluded from the prediction simplex.

To assess an optimal embedding dimension EmbedDimension can be applied. Accuracy statistics can be estimated by ComputeError.

Value

A data.frame with columns Time, Observations, Predictions and Pred_Variance.

References

Sugihara G. and May R. 1990. Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series. Nature, 344:734-741.

Examples

data( block_3sp )
smplx = Simplex(block_3sp, "x_t", "x_t", lib = c(1, 100), pred = c(101, 195), E = 3)
ComputeError( smplx $ Predictions, smplx $ Observations )

Surrogate time series for significance testing

Description

Generate surrogate series (random_shuffle, ebisuzaki phase randomisation, or seasonal).

Usage

SurrogateData(ts, method = c("random_shuffle", "ebisuzaki", "seasonal"),
              num_surr = 100, T_period = 1, alpha = 0)

Arguments

Value

A numeric matrix with length(ts) rows and num_surr columns.

Examples

SurrogateData(sin(1:100), method = "ebisuzaki", num_surr = 5)

Time series for a tent map with mu = 2.

Description

First-differenced time series generated from the tent map recurrence relation with mu = 2.

Usage

TentMap

Format

Data frame with 999 rows and 2 columns

Time

time index.

TentMap

tent map values.

Time series of tent map plus noise.

Description

First-differenced time series generated from the tent map recurrence relation with mu = 2 and random noise.

Usage

TentMapNoise

Format

Data frame with 999 rows and 2 columns

Time

time index.

TentMap

tent map values.

Apple-blossom Thrips time series

Description

Seasonal outbreaks of Thrips imaginis.

References

Davidson and Andrewartha, Annual trends in a natural population of Thrips imaginis Thysanoptera, Journal of Animal Ecology, 17, 193-199, 1948.

Time series for a three-species coupled model.

Description

Time series generated from a discrete-time coupled Lotka-Volterra model exhibiting chaotic dynamics.

Usage

block_3sp

Format

A data frame with 198 rows and 10 columns:

time

time index (# of generations)

x_t

abundance of simulated species x at time t

x_t-1

abundance of simulated species x at time t-1

x_t-2

abundance of simulated species x at time t-2

y_t

abundance of simulated species y at time t

y_t-1

abundance of simulated species y at time t-1

y_t-2

abundance of simulated species y at time t-2

z_t

abundance of simulated species z at time t

z_t-1

abundance of simulated species z at time t-1

z_t-2

abundance of simulated species z at time t-2

2-D timeseries of a circle.

Description

Time series of of circle in 2-D (sin and cos).

Usage

circle

Format

A data frame with 200 rows and 3 columns:

Time

time index.

x

sin component.

y

cos component.

Time series for the Paramecium-Didinium laboratory experiment

Description

Time series of Paramecium and Didinium abundances (#/mL) from an experiment by Veilleux (1979)

Usage

paramecium_didinium

Time series for the California Current Anchovy-Sardine-SST system

Description

Time series of Pacific sardine landings (CA), Northern anchovy landings (CA), and sea-surface temperature (3-year average) at the SIO pier and Newport pier

Usage

sardine_anchovy_sst

Format

year

year of measurement

anchovy

anchovy landings, scaled to mean = 0, sd = 1

sardine

sardine landings, scaled to mean = 0, sd = 1

sio_sst

3-year running average of sea surface temperature at SIO pier, scaled to mean = 0, sd = 1

np_sst

3-year running average of sea surface temperature at Newport pier, scaled to mean = 0, sd = 1