Title:	Particle-Takens Stability
Version:	1.4
Description:	Includes a collection of functions presented in "Measuring stability in ecological systems without static equilibria" by Clark et al. (2022) <doi:10.1002/ecs2.4328> in Ecosphere. These 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). 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 estimating variability, this package does not necessarily 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).
Depends:	R (≥ 3.4)
Imports:	graphics, stats
Suggests:	BayesianTools
License:	GPL-3
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.2.3
Collate:	'bayesfun.R' 'data.R' 'fake_data.R' 'logit_funs.R' 'particlefilter.R' 'pttstability_man.R' 'smapping_functions.R'
NeedsCompilation:	no
Author:	Adam Clark [aut, cre]
Maintainer:	Adam Clark <adam.tclark@gmail.com>
Repository:	CRAN
Packaged:	2024-01-09 14:41:38 UTC; aclark
Date/Publication:	2024-01-09 15:00:02 UTC

EDM deterministic function

Description

Estimates future states of xt based on based behaviour

Usage

EDMfun0(smp_cf, yp, x, minest = 0, maxest = NULL, time)

Arguments

Value

a number or numeric vector of length xt, with predicted abundances at time t+1

Source

Adapted from Ye, Sugihara, et al. (2015), PNAS 112:E1569-E1576.

Apply S-mapping algorithm from Sugihara 1994

Description

Carries out an S-mapping analysis, following the algorithm outlined in Sugihara (1994).

Usage

S_map_Sugihara1994(Y, E, theta, X = NULL, lib = NULL, trimNA = FALSE)

Arguments

Value

a list, including the timeseries used for S-mapping (Y), the delay embedding matrix used for S-mapping (X), a vector of predictions (Y_hat), a matrix of S-mapping coefficients (C), the standard errors for the S-mapping coefficients (C_SE), and goodness of fit metrics R-squared (R2) and root mean square error (RMSE).

Source

Sugihara, G. (1994). Nonlinear forecasting for the classification of natural time-series. Philos. Trans. R. Soc. -Math. Phys. Eng. Sci., 348, 477–495.

Examples

# create an example timeseries
n = 100
set.seed(1234)
datout<-makedynamics_general(n = n+2,
                             pdet=log(c(0.8,1)),
                             proc = -2.5,
                             detfun = detfun0_sin)
plot(datout$true, type = "l") # plot timeseries
Y = datout$true # extract true values

# run s-mapping
sout = S_map_Sugihara1994(Y = Y, E = 2, theta = 0.5)
s_coef = process_scof(sout$C) # process coefficients from the S-mapping output

# find best E/theta
fitout = data.frame(E = 1:5, theta = NA, RMSE = NA)

for(i in 1:nrow(fitout)) {
   E = fitout$E[i]
   Ytmp = Y[-c(1:E)]
   optout = optimize(f = function(x) {S_map_Sugihara1994(Ytmp, E, x)$RMSE}, interval = c(0,10))
  
  fitout$theta[i] = optout$minimum # get best theta for given E
  fitout$RMSE[i] = optout$objective # get error
}
ps = which.min(fitout$RMSE)

E = fitout$E[ps] # get best E
theta = fitout$theta[ps] # get best theta
X = makeblock(Y, E) # get X for analysis
Y = Y[-c(1:E)] # trim NA values (corresponding to positions in X)
X = X[(E+1):nrow(X),] # trim NA values
sout = S_map_Sugihara1994(Y = Y, E = E,
  theta = theta, X = X) # run S-mapping for best paramter combination
sout$R2 # look at R-squared

# check fit
plot(sout$Y_hat, Y)
abline(a=0, b=1, lty=2)

default colonization function

Description

Simulates colonization events - events occur as a binomial random process with probability ilogit(p), and populations are seeded with abundance exp(A).

Usage

colfun0(co, xt)

Arguments

Value

a numeric, including number or numeric vector of length xt, with predicted abundances after colonization has occurred

Microcosm experimental data

Description

A dataset containing the abundances of Chlamydomonas terricola growing in a multi-species community. Includes 17 time steps covering 463 days in two treatments: LSA ('low' temperature, 'stable' oscillations, and 'absence' of predators) and LVA ('low' temperature, 'variable' oscillations, and 'absence' of predators).

Usage

dat

Format

A data frame with 271 rows and 4 variables:

treatment

Experimental treatment

number

Replicate number

time

Day of experiment

Chlamydomonas.terricola

Species abundance

Source

Burgmer & Hillebrand 2011, Oikos 120:922-933.

Default density function for prior

Description

Default density function, following the syntax for priors in the BayesianTools package. Uses flat priors for all paramters, within the given interval. Density function integrates to one.

Usage

density_fun0(param, minv, maxv)

Arguments

Value

returns log likelihood of parameters given priors.

default deterministic function

Description

Simulates deterministic component of Ricker model, of the form xt+1 = xt exp(exp(sdet[1])*(1-xt/exp(sdet[2])))

Usage

detfun0(sdet, xt, time = NULL, ...)

Arguments

Value

a number or numeric vector of length xt, with predicted abundances at time t+1

deterministic function with time-varying carrying capacity

Description

Simulates deterministic component of Ricker model, of the form xt+1 = xt exp(exp(sdet[1])*(1-xt/K)) where K varies with time as (sin(time/2)+exp(sdet[2])+0.5)*(2/3). Function is calibrated such that for exp(sdet[2]) = 1, mean(K) = 1.

Usage

detfun0_sin(sdet, xt, time = NULL, ...)

Arguments

Value

a number or numeric vector of length xt, with predicted abundances at time t+1

Get rates

Description

Calculates colonization rate, mortality rate, and expected mean occupancy time based on a time series

Usage

getcm(dat)

Arguments

Value

a list including colonization and mortality probability per time step (pc and pm, respectively), and pocc, the expected fraction of time that the species will be present

Inverse logit

Description

Returns the inverse logit transformation of x

Usage

ilogit(x, ...)

Arguments

Value

transformed result

Sort output of particle filter

Description

Sorts outputs of particle filter based on index - returns a sorted list of particles, based on the sampling trajectory through time. This is a somewhat more accurate estiamte of the true posterior than are the stepwise samples provided by the filter.

Usage

indexsort(fulltracemat, fulltraceindex, nsmp = NULL)

Arguments

Value

an index-sorted matrix - each column shows the trajectory of a single particle

Default inverse transormation function

Description

Takes in a matrix, where each column represents a parameter. Returns parameters in untransformed space. If length = 2, then in the order (obs1, proc1). If 3, then in the order (obs1, proc1, proc2). If 4, then in the order (obs1, obs2, proc1, proc2). If 6, then in the order (obs1, proc1, pcol1, pcol2, det1, det2) If 7, then in the order (obs1, proc1, proc2, pcol1, pcol2, det1, det2) If 8, then in the order (obs1, obs2, proc1, proc2, pcol1, pcol2, det1, det2)

Usage

inv_fun0(x)

Arguments

Value

returns back-transformed values of parameters

Default likelihood function

Description

Calculates likelihood of vector y given parameter values in param, based on the particleFilterLL function.

Usage

likelihood0(
  param,
  y = y,
  parseparam = parseparam0,
  N = 1000,
  detfun = detfun0,
  edmdat = NULL,
  obsfun = obsfun0,
  procfun = procfun0,
  neff = FALSE,
  lowerbound = (-999)
)

Arguments

Value

Log likelihood generated by particleFilterLL function

calculate likelihood for piecewise data

Description

Calculates likelihoods across several segments of data - e.g. multiple plots from a single experiment. See documentation for particleFilterLL_piecewise for examples of use.

Usage

likelihood_EDM_piecewise(
  param,
  y,
  libuse_y,
  smap_coefs,
  Euse,
  tuse,
  N,
  colpar = c(logit(1e-06), log(0.1))
)

Arguments

Value

summed log likelihood across all segments

Logit

Description

Returns the logit transformation of x

Usage

logit(x, ...)

Arguments

Value

transformed result - impossible values are replaced with NA, without warnings

Get inverse logit-normal mode

Description

Returns a mean for a logit normal such that the mode will be centered around mu

Usage

logitnormal_imode(mu, sd)

Arguments

Value

the proposed mean for the distribution

Get inverse log-normal mode

Description

Returns a mean for a lognormal such that the mode will be centered around mu

Usage

lognormal_imode(mu, sd)

Arguments

Value

the proposed mean for the distribution

Make an embedding block from timeseries data

Description

Returns a matrix X, where columns are time-delayed embeddings of Y, with number of embeddings specified by embedding dimension E. See help file for the S_map_Sugihara1994 function for examples.

Usage

makeblock(Y, E, lib = NULL)

Arguments

Value

a matrix of time-delayed embeddings

Simulate general time series

Description

Simulates a time series following a user-defined deterministic function, observation function, process noise function, and colonization function.

Usage

makedynamics_general(
  n = 1000,
  n0 = 0.1,
  pdet = c(log(3), log(1)),
  proc = c(log(1)),
  obs = c(log(1)),
  pcol = c(logit(0.2), log(1)),
  detfun = detfun0,
  procfun = procfun0,
  obsfun = obsfun0,
  colfun = colfun0,
  doplot = FALSE
)

Arguments

Value

An n-by-3 dataframe of states, including obs (observed values), truth (true values), and noproc (values without process noise)

Examples

 #run function
 datout<-makedynamics_general(n=2e4, proc = c(-2,log(1.2)))

 #show regression of variance vs. mean for binned data
 datout_ps<-datout[datout$true>0 & datout$noproc>0,]
 #bins
 sq<-seq(0, quantile(datout$true, 0.95), length=50)
 ctd<-cut(datout_ps$noproc, sq)
 #calculate mean and variance by bin
 tdat<-data.frame(mu=(sq[-1]+sq[-length(sq)])/2,
      var=tapply((datout_ps$true-datout_ps$noproc)^2, ctd, mean))
 #plot result
 plot(log(tdat$mu), log(tdat$var), xlab="mu", ylab="var")
 #show regression
 summary(mod<-lm(log(var)~log(mu), tdat)); abline(mod, col=2)

default observation noise function

Description

Two options: If inverse=FALSE, calculates the log probability density of observation yt based on true state xt and observation error. Otherwise, simulates N random observations of yt. Observation error follows a Gaussian distribution truncated at zero, using a Tobit distribution. Note that probability density is calculated based on a Tobit distribution, with lower boundary zero.

Usage

obsfun0(
  so,
  yt,
  xt = NULL,
  inverse = FALSE,
  N = NULL,
  minsd = 0.01,
  time = NULL
)

Arguments

Value

If inverse=FALSE, returns a list including LL, a number or numeric vector of length xt, with predicted log likelihoods of observation yt, and wts, a number or vector with weights corresponding to the relative likelihood of each observation (after accounting for variable continuous vs. discrete probability distributions). If inverse = FALSE, returns N random draws from the observation function.

Parse parameters

Description

Takes in a vector of 3 or 6 parameters, and puts them into a list of the format expected by the particleFilterLL function.

Usage

parseparam0(
  param,
  colparam = c(logit(0.2), log(0.1)),
  detparam = c(log(1.2), log(1))
)

Arguments

Value

a formatted list of parameters

particle filter

Description

General function for caluclating the log-likeihood of a stochastic discrete-time model, based on a noisy observation of time-series y. Returns estimates of true values of y, as well as for process noise, observation error, colonization rates, and extinction rates. Function is adapted from the R code of Knape and Valpine (2012), Ecology 93:256-263.

Usage

particleFilterLL(
  y,
  pars,
  N = 1000,
  detfun = detfun0,
  procfun = procfun0,
  obsfun = obsfun0,
  colfun = colfun0,
  edmdat = NULL,
  dotraceback = FALSE,
  fulltraceback = FALSE
)

Arguments

Value

LL (total log likelihood), LLlst (log likelihood for each time step), Nest (mean estimated state), Nsd (standard deviation of estimated state), Nest_noproc (mean estimated state at time t+1 without process error), Nsd_noproc (standard deviation of estimated state at time t+1 without process error), fulltracemat (full traceback of particle paths), fulltracemat_noproc (full traceback of particle paths at time t+1 without process noise), and fulltraceindex (index positions for the particle traces over time)

Source

Adapted from Knape and Valpine (2012), Ecology 93:256-263.

run particle filter across piecewise data

Description

Calculates likelihoods across several segments of data - e.g. multiple plots from a single experiment. Requires several implicitely defined variables to run:

Usage

particleFilterLL_piecewise(
  param,
  N,
  y,
  libuse_y,
  smap_coefs,
  Euse,
  tuse,
  colpar = c(logit(1e-06), log(0.1)),
  nsmp = 1,
  lowerbound = -999,
  maxNuse = 512000
)

Arguments

Value

results from particle filter - including mean estimates (Nest) and standard deviations (Nsd), across particles, and sample particle trajectories with (Nsmp) and without (Nsmp_noproc) process noise

Examples

# load data
data(dat)
# sort by index
dat = dat[order(dat$treatment, dat$number, dat$time),]


# make list of starting and ending positions for each replicate in the dat list
libmat<-NULL
trtmat<-data.frame(trt=as.character(sort(unique(dat$treatment))))
datnum<-1:nrow(dat)

for(i in 1:nrow(trtmat)) {
  ps1<-which(dat$treatment==trtmat$trt[i])
  replst<-sort(unique(dat$number[ps1]))

  for(j in 1:length(replst)) {
    ps2<-which(dat$number[ps1]==replst[j])
    libmat<-rbind(libmat, data.frame(trt=trtmat$trt[i], rep=replst[j],
      start=min(datnum[ps1][ps2]), end=max(datnum[ps1][ps2])))
  }
}

## run particle filter
# select treatment to analyse: enter either "LSA" or "LSP"
trtuse<-"HSP"
# extract library positions for treatment
libuse<-as.matrix(libmat[libmat$trt==trtuse,c("start", "end")])
# save abundance data to variable y
yps<-which(dat$treatment==trtuse)
y<-dat[,"Chlamydomonas.terricola"][yps]
libuse_y<-libuse-min(libuse)+1 # translate positions in dat to positions in y vector
y<-y/sd(y) # standardize to mean of one
timesteps<-dat$time[yps]

# get S-mapping parameters
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, lib = libuse_y)$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, lib = libuse_y)

# set priors (log-transformed Beta_obs, Beta_proc1, and Beta_proc2)
minvUSE_edm<-c(log(0.001), log(0.001))  # lower limits
maxvUSE_edm<-c(log(2), log(2)) # upper limits


## Not run: 
  ## Run filter
  # Commented-out code: Install BayesianTools package if needed
  #install.packages("BayesianTools")
  set.seed(2343)
  require(BayesianTools)
  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)
  niter<-5e3 # number of steps for the MCMC sampler
  N<-2e3 # number of particles
  smap_coefs<-process_scof(spred$C) # coefficients from s-mapping routine

  # likelihood and bayesian set-ups for EDM functions
  likelihood_EDM_piecewise_use<-function(x) {
    # default values for filter - see likelihood_EDM_piecewise documentation for details
    # note that colpar are set near zero because we expect no colonisation into a closed microcosm.
    likelihood_EDM_piecewise(param=x, y, libuse_y, smap_coefs, euse, tuse, N,
                             colpar = c(logit(1e-06), log(0.1)))
  }

  bayesianSetup_EDM <- createBayesianSetup(likelihood = likelihood_EDM_piecewise_use,
     prior = prior_edm)

  # run MCMC optimization (will take ~ 5 min)
  out_EDM <- runMCMC(bayesianSetup = bayesianSetup_EDM,
     settings = list(iterations=niter, consoleUpdates=20))
  burnin<-floor(niter/5) # burnin period
  plot(out_EDM, start=burnin) # plot MCMC chains
  gelmanDiagnostics(out_EDM, start=burnin) # calculate Gelman statistic
  summary(out_EDM, start=burnin) # coefficient summary

  ## extract abundance estimate from particle filter
  # use final estimate from MCMC chain
  smp_EDM<-(getSample(out_EDM, start=floor(niter/5)))
  tmp<-particleFilterLL_piecewise(param = smp_EDM[nrow(smp_EDM),], N=N, y = y, libuse_y = libuse_y,
                                  smap_coefs = smap_coefs, Euse = euse, tuse = tuse)
  # mean estimated abundance
  simout<-tmp$Nest
  # sd estimated abundance
  sdout<-tmp$Nsd
  # sample from true particle trajectory
  simout_smp<-tmp$Nsmp
  # sample from true particle trajectory pre-process noise
  simout_smp_noproc<-tmp$Nsmp_noproc

  plot(timesteps, simout, xlab="Time", ylab="Abundance")
  abline(h=0, lty=3)

## End(Not run)

Process S-mapping coefficients

Description

Processes s-mapping coefficients from S_map_Sugihara1994 into a matrix of form C1, C2, C3, ... C0, where C0 is the intercept, C1 is the current time step t, C2 is timestep t-1, C3 is timestep t-2, and so on. Rows correspond to the time step used to produce the prediction, e.g. row 4 is used to calculate predicted value for time step 5. This is the format expected by the EDMfun0 function. See help file for the S_map_Sugihara1994 function for examples.

Usage

process_scof(smap_coefs)

Arguments

Value

a matrix of s-mapping coefficients

default process noise function

Description

Simulates effects of process noise following a Gaussian perturbation. Note that process noise only influences positive abundances (i.e. process noise cannot contribute to colonization)

Usage

procfun0(sp, xt, inverse = FALSE, time = NULL)

Arguments

Value

a number or numeric vector of length xt, with predicted abundances after process noise has occurred

continuous-time process noise function

Description

Simulates effects of process noise following a Gaussian perturbation. Note that process noise only influences positive abundances (i.e. process noise cannot contribute to colonization)

Usage

procfun_ct(sp, xt, waiting_time = 1, time = NULL)

Arguments

Value

a number or numeric vector of length xt, with predicted abundances after process noise has occurred

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)

Default sampler function for prior

Description

Draws samples from a flat prior

Usage

sampler_fun0(n = 1, minv, maxv)

Arguments

Value

returns random draws from the priors

calculate estimated total variance

Description

Function for estimating stochastic variation in linar process x as a function of relative growth rate and disturbance regime standard deviation.

Usage

sdproc_abstract(sd_proc, rgr, waiting_time = 1)

Arguments

Value

standard deviation of stochastic variability in x