pttstability: pttstability: Methods for Measuring Stability in Systems...

In pttstability: Particle-Takens Stability

pttstability: Methods for Measuring Stability in Systems Without Static Equilibria

Description

The pttstability ("ParTicle-Takens Stability") package is a collection of functions that can be used to estimate the parameters of a stochastic state space model (i.e. a model where a time series is observed with error).

Applications

The goal of this package is to estimate the variability around a deterministic process, both in terms of observation error - i.e. variability due to imperfect observations that does not influence system state - and in terms of process noise - i.e. stochastic variation in the actual state of the process. Unlike classical methods for estmiating variability, this package does not necesarilly assume that the deterministic state is fixed (i.e. a fixed-point equilibrium), meaning that variability around a dynamic trajectory can be estimated (e.g. stochastic fluctuations during predator-prey dynamics).

By combining information about both the estimated deterministic state of the system and the estimated effects of process noise, this package can be used to compute a dynamic analog of various stability metrics - e.g. coefficient of variation (CV) or invariability. Estimated extinction rates and colonization rates can also be estimated.

Contents

This package builds on three existing toolkits. First, it applies an updated version of the "particleFilterLL" particle filter function of Knape and Valpine (2012) to calculate likelihoods of observed time series given modeled dynamics. Second, it applies empirical dynamic modeling (EDM) methods to estimate deterministic dynamics even in cases where the underlying equations governing system behavior are not known. These models are based on Takens delay-embedding theorem, from which this package takes its name. Finally, it (optionally) uses the MCMC fitting methods from the BayesianTools package to estimate paramter values for the observation error, process noise, and (optionally) deterministic functions underlying observed dynamics.

The default observation error and process noise functions in this package (obsfun0 and procfun0) take advantage of the Taylor Power law to separate noise components for relatively short time series. Observation error is assumed to scale with the observed state as sd_obs(x) = x*exp(obs), Process noise is either a constant (i.e. sd_proc(x) = exp(proc)), or, if two variables are given, process noise scales as a power function of the observed value as sd_proc(x) = sqrt(exp(proc1)*x^exp(proc2))

Note that although we include default functions in this package, users are able (and encouraged!) to write their own (including for observation error, process noise, deterministic dynamics, priors, and likelihoods).

Source

Adam T. Clark, Lina K. Mühlbauer, Helmut Hillebrand, and Canan Karakoç. (2022). Measuring Stability in Ecological Systems Without Static Equilibria. Ecosphere 13(12):e4328.

Knape, J., and Valpine, P. (2012). Fitting complex population models by combining particle filters with Markov chain Monte Carlo. Ecology 93:256-263.

Ye, H., Sugihara, G., et al. (2015). Equation-free ecosystem forecasting. PNAS 112:E1569-E1576.

Ye, H., et al. (2019). rEDM: Applications of Empirical Dynamic Modeling from Time Series. R package version 0.7.2.

Hartig, F., et al. (2019). BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics. R package version 0.1.6.

Examples

#Set seed
set.seed(5826)

## Simulate data
pars_true<-list(obs=log(0.15),
            proc=c(log(0.1)),
            pcol=c(logit(0.2), log(0.1)),
            det=c(log(1.2),log(1)))
#parameters for the filter
pars_filter<-pars_true

#generate random parameter values
datout<-makedynamics_general(n = 100, n0 = exp(pars_true$det[2]),
                             pdet=pars_true$det, proc = pars_true$proc,
                             obs = pars_true$obs, pcol = pars_true$pcol,
                             detfun = detfun0_sin, procfun = procfun0,
                             obsfun=obsfun0, colfun=colfun0)
y<-datout$obs
plot(y, type = "l", xlab="time", ylab="observed abundance")

#get parameters for S-mapping
sout<-data.frame(E = 2:4, theta = NA, RMSE = NA)
for(i in 1:nrow(sout)) {
  optout = optimize(f = function(x) {S_map_Sugihara1994(Y = y, E = sout$E[i],
    theta = x)$RMSE}, interval = c(0,10))
  sout$theta[i] = optout$minimum
  sout$RMSE[i] = optout$objective
}
tuse<-sout$theta[which.min(sout$RMSE)] # find theta (nonlinerity) parameter
Euse<-sout$E[which.min(sout$RMSE)] # find embedding dimension
spred<-S_map_Sugihara1994(Y = y, E = Euse, theta = tuse) # fit S-mapping for best paramter set
plot(spred$Y, spred$Y_hat); abline(a=0, b=1, lty=2) # observed vs. predicted

## Run filter with "correct" parameter values
N = 1e3 # number of particles
#based on detful0
filterout_det<-particleFilterLL(y, pars=pars_filter, N,
                                detfun = detfun0_sin, procfun = procfun0,
                                dotraceback = TRUE, fulltraceback = TRUE)
#based on EDM
filterout_edm<-particleFilterLL(y, pars=pars_filter, N, detfun = EDMfun0,
                                edmdat = list(E=Euse, theta=tuse),
                                procfun = procfun0, dotraceback = TRUE,
                                fulltraceback = TRUE)

#plot filter output
op = par(mar=c(4,4,2,2), mfrow=c(3,1))
#plot 30 of the 1000 particles to show general trend
# correct deterministic function
matplot(1:length(y), filterout_det$fulltracemat[,1:30],
        col=adjustcolor(1,alpha.f = 0.5), lty=3,
        type="l", xlab="time", ylab="abund", main="detfun0")
lines(1:length(y), y, col=2, lwd=1.5) #observations
lines(1:length(y), datout$true, col="blue", lwd=1.5, lty=2) #true values
lines(1:length(y), filterout_det$Nest, col=3, lwd=1.5) #mean filter estimate

# EDM function
matplot(1:length(y), filterout_edm$fulltracemat[,1:30],
        col=adjustcolor(1,alpha.f = 0.5), lty=3,
        type="l", xlab="time", ylab="abund", main="EDM")
lines(1:length(y), y, col=2, lwd=1.5)
lines(1:length(y), datout$true, col="blue", lwd=1.5, lty=2)
lines(1:length(y), filterout_edm$Nest, col=3, lwd=1.5)

plot(filterout_det$Nest, datout$true, xlim=range(c(filterout_det$Nest,
     filterout_edm$Nest)),
     xlab="predicted", ylab="true value", col=4)
points(filterout_edm$Nest, datout$true, col=2)
points(y, datout$true, col=3)
abline(a=0, b=1, lty=2)
legend("topleft", c("detfun0", "EDM", "obs"), pch=1, col=c(4,2,3), bty="n")

#note improvement in fit, for both filters
cor(datout$true, datout$obs)^2 #observations
cor(datout$true, filterout_det$Nest)^2 #deterministic filter
cor(datout$true, filterout_edm$Nest)^2 #EDM filter
par(op) # reset plotting parameters

## Not run: 
# Commented-out code: Install BayesianTools if needed
#install.packages("BayesianTools")
require(BayesianTools)
## Run optimizers
#create priors
minvUSE<-c(-4, -4) #minimum interval for obs and proc
maxvUSE<-c(0, 0) #maximum interval for obs and proc

minvUSE_edm<-c(-4, -4) #minimum interval for obs and proc
maxvUSE_edm<-c(0, 0) #maximum interval for obs and proc

#density, sampler, and prior functions for deterministic function
density_fun_USE<-function(param) density_fun0(param = param, minv = minvUSE,
                         maxv=maxvUSE)
sampler_fun_USE<-function(x) sampler_fun0(n = 1, minv = minvUSE,
                         maxv=maxvUSE)
prior_USE <- createPrior(density = density_fun_USE,
                         sampler = sampler_fun_USE,
                         lower = minvUSE, upper = maxvUSE)

#density, sampler, and prior functions for EDM function
density_fun_USE_edm<-function(param) density_fun0(param = param,
                         minv = minvUSE_edm, maxv=maxvUSE_edm)
sampler_fun_USE_edm<-function(x) sampler_fun0(n = 1, minv = minvUSE_edm,
                         maxv=maxvUSE_edm)
prior_edm <- createPrior(density = density_fun_USE_edm,
                         sampler = sampler_fun_USE_edm,
                         lower = minvUSE_edm, upper = maxvUSE_edm)
## Run filter
niter<-5000 #number of steps for the MCMC sampler
N<-1e3 #number of particles
Euse<-Euse #number of embedding dimensions

#likelihood and bayesian set-ups for deterministic functions
likelihood_detfun0<-function(x) likelihood0(param=x, y=y,
                         parseparam = parseparam0, detfun = detfun0_sin,
                         procfun = procfun0, N = N)
bayesianSetup_detfun0 <- createBayesianSetup(likelihood = likelihood_detfun0,
                         prior = prior_USE)

#likelihood and bayesian set-ups for EDM functions
likelihood_EDM<-function(x) {
  likelihood0(param = x, y=y, parseparam = parseparam0, procfun = procfun0,
              detfun = EDMfun0, edmdat = list(E=Euse, theta=tuse), N = N)
}

bayesianSetup_EDM <- createBayesianSetup(likelihood = likelihood_EDM,
                         prior = prior_edm)
#run MCMC optimization
out_detfun0 <- runMCMC(bayesianSetup = bayesianSetup_detfun0,
                   settings = list(iterations=niter, consoleUpdates=20))
out_EDM <- runMCMC(bayesianSetup = bayesianSetup_EDM,
                   settings = list(iterations=niter, consoleUpdates=20))

#plot results, with a 1000-step burn-in
plot(out_detfun0, start = 1000, thin = 2)
plot(out_EDM, start = 1000, thin = 2)

## extract and plot parameter distributions
smp_detfun0<-getSample(out_detfun0, start = 1000, thin = 2)
smp_EDM<-getSample(out_EDM, start=1000, thin = 2)

op = par(mfrow=c(2,2))
hist(exp(smp_detfun0[,1]), xlim=c(exp(minvUSE[1]), exp(maxvUSE[1])),
                    main="det. function", xlab="obs", breaks = 20)
abline(v=exp(pars_true$obs), col=2) # true value
abline(v=c(exp(minvUSE[1]), exp(maxvUSE[1])), col=1, lty=2)# Priors

hist(exp(smp_detfun0[,2]), xlim=c(exp(minvUSE[2]), exp(maxvUSE[2])),
                    main="det. function", xlab="proc", breaks = 20)
abline(v=exp(pars_true$proc), col=2) # true value
abline(v=c(exp(minvUSE[2]), exp(maxvUSE[2])), col=1, lty=2)# Priors

hist(exp(smp_EDM[,1]), xlim=c(exp(minvUSE_edm[1]), exp(maxvUSE_edm[1])),
                    main="EDM function", xlab="obs", breaks = 20)
abline(v=exp(pars_true$obs), col=2) # true value
abline(v=c(exp(minvUSE_edm[1]), exp(maxvUSE_edm[1])), col=1, lty=2)# Priors

hist(exp(smp_EDM[,2]), xlim=c(exp(minvUSE_edm[2]), exp(maxvUSE_edm[2])),
                    main="EDM function", xlab="proc", breaks = 20)
abline(v=exp(pars_true$proc), col=2) # true value
abline(v=c(exp(minvUSE_edm[2]), exp(maxvUSE_edm[2])), col=1, lty=2)# Priors
# note slight over-estimate of process noise due to EDM error
par(op) # reset plotting parameters

## End(Not run)

pttstability documentation built on May 29, 2024, 1:28 a.m.