| Title: | Statistical Scaling Functions for Ecological Systems |
|---|---|
| Description: | Implementation of the scaling functions presented in "General statistical scaling laws for stability in ecological systems" by Clark et al in Ecology Letters <DOI:10.1111/ele.13760>. Includes functions for extrapolating variability, resistance, and resilience across spatial and ecological scales, as well as a basic simulation function for producing time series, and a regression routine for generating unbiased parameter estimates. See the main text of the paper for more details. |
| Authors: | Adam Clark [aut, cre] |
| Maintainer: | Adam Clark <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1 |
| Built: | 2025-02-17 05:00:08 UTC |
| Source: | https://github.com/adamtclark/ecostatscale |
Helper function for symdynN, used to simulate dynamics between disturbance events. Note that K is set to 1 for all species.
df_col(time, state, pars)df_col(time, state, pars)
time |
time step - ignored, but required for consistency with ode function |
state |
vector of current states |
pars |
parameter list, including matrix A with full interaction matrix (and growth rates along the diagonal), and a value Ifrac, which is the dispersal rate (D in the main text), and Ksim, which is the carrying capacity |
rate of change for each species
Helper function for symdynN, used to simulate dynamics between disturbance events. Note that K is set to 1 for all species.
df_col_loss(time, state, pars)df_col_loss(time, state, pars)
time |
time step - ignored, but required for consistency with ode function |
state |
vector of current states |
pars |
parameter list, including matrix A with full interaction matrix (and growth rates along the diagonal), and a value Ifrac, which is the dispersal rate (D in the main text), and Ksim, which is the carrying capacity, and Iloss which is the loss rate. |
rate of change for each species
Helper function for symdynN, used to simulate dynamics between disturbance events.
df0(time, state, pars)df0(time, state, pars)
time |
time step - ignored, but required for consistency with ode function |
state |
vector of current states |
pars |
parameter list, including matrix A with full interaction matrix (and growth rates along the diagonal) |
rate of change for each species
Extrapolate resilience observed at the scale of a single spatial or ecological scale b (e.g. a patch or species) to a larger scale, B (e.g. functional group or landscape).
res_scale( mvar_b, murho_b_abundance, mucov_b_abundance = NULL, msd_b, murho_b_disturbance, mucov_b_disturbance = NULL, b = 1, B, lambda )res_scale( mvar_b, murho_b_abundance, mucov_b_abundance = NULL, msd_b, murho_b_disturbance, mucov_b_disturbance = NULL, b = 1, B, lambda )
mvar_b |
Mean abundance variance observed at scale b |
murho_b_abundance |
Mean Pearson correlation coefficient of abundance values observed at scale b |
mucov_b_abundance |
Mean covariance of abundance values observed at scale b. Ignored unless murho_b_abundance is NULL Defaults to NULL. |
msd_b |
Mean disturbance standard deviation observed at scale b |
murho_b_disturbance |
Mean Pearson correlation coefficient of disturbances observed at scale b |
mucov_b_disturbance |
Mean covariance of disturbances observed at scale b. Ignored unless murho_b_abundance is NULL Defaults to NULL. |
b |
Size of observed scale. Defaults to 1. |
B |
Larger scale being extrapolated to (e.g. total number of species, or size of patch B relative to b) |
lambda |
Mean disturbance frequency. |
Extrapolated median resilience at scale of M species.
# extrapolate from scale of 1 species to 10 species res_scale(mvar_b = 0.25, murho_b_abundance = -0.034, msd_b = sqrt(0.1), murho_b_disturbance = 0, B = 30, lambda=1) # plot relationship for groups of 1 to 30 species plot(1:30, res_scale(mvar_b = 0.25, murho_b_abundance = -0.034, msd_b = sqrt(0.1), murho_b_disturbance = 0, B = 1:30, lambda=1), xlab="ecological scale", ylab="resilience, r", type="b")# extrapolate from scale of 1 species to 10 species res_scale(mvar_b = 0.25, murho_b_abundance = -0.034, msd_b = sqrt(0.1), murho_b_disturbance = 0, B = 30, lambda=1) # plot relationship for groups of 1 to 30 species plot(1:30, res_scale(mvar_b = 0.25, murho_b_abundance = -0.034, msd_b = sqrt(0.1), murho_b_disturbance = 0, B = 1:30, lambda=1), xlab="ecological scale", ylab="resilience, r", type="b")
Extrapolate disturbance standard deviation observed at spatial or ecological scale b to a different scale, B (inversely related to resistance). Equivalent to Eq.7b in the main text.
sd_scale(msd_b, murho_b, mucov_b = NULL, b = 1, B)sd_scale(msd_b, murho_b, mucov_b = NULL, b = 1, B)
msd_b |
Mean disturbance standard deviation observed at scale b. |
murho_b |
Mean Pearson correlation coefficient of disturbances observed at scale b, calculated as mucov_b/mvar_b. If NULL, mucov_b is used instead. |
mucov_b |
Mean covariane of disturbances observed at scale b. Ignored if mrho_b is not NULL. Defaults to NULL. |
b |
Size of observed scale. Defaults to 1. |
B |
Size of desired scale for extrapolation. |
Extrapolated disturbance standard deviation at scale B.
#extrapolate from scale of 1 to 10 - e.g. from a 1m2 patch to a 10m2 patch sd_scale(msd_b = 1, murho_b = 0.5, b = 1, B = 10)#extrapolate from scale of 1 to 10 - e.g. from a 1m2 patch to a 10m2 patch sd_scale(msd_b = 1, murho_b = 0.5, b = 1, B = 10)
Function for simulating dynamics from Eq.1 in the main text.
symdyn( r, f, d, d_sd, sf, tmax, stochd = TRUE, stocht = TRUE, as.matrix = FALSE, oscillate_dist = FALSE )symdyn( r, f, d, d_sd, sf, tmax, stochd = TRUE, stocht = TRUE, as.matrix = FALSE, oscillate_dist = FALSE )
r |
per-capita growth rate (r in Eq.1) |
f |
the waiting time (or average waiting time) between disturbance events (equal to 1/lambda in Eq.1) |
d |
mean size of disturbance function (mu in Eq.1) |
d_sd |
standard deviation of disturbance function (sigma in Eq.1) |
sf |
waiting time between sampling events |
tmax |
the time series length to be simulated |
stochd |
a logical variable, indicating whether disturbance size should be stochastic - otherwise, all disturbances are of magnitude d - defaults to TRUE |
stocht |
a logical variable, indicating whether waiting time between disturbance events should be stochastic - otherwise, waiting time is always f - defaults to TRUE |
as.matrix |
indicates whether results should be returned as matrix (potentially faster for some applications) - defaults to FALSE |
oscillate_dist |
a logical variable indicating whether the sign of the disturbance should oscillate between positive and negative - ignored if stochd==TRUE - defaults to FALSE |
a matrix or data.frame with columns for sampling times, abundances, and number of disturbances for each time interval
# see xt2fun# see xt2fun
Function for simulating dynamics from Eq.2-3 in the main text.
symdynN( r, amu, asd, f, d, d_sd, d_cov, N, sf, tmax, stochd = TRUE, stocht = TRUE, as.matrix = FALSE, amax = 0, amin = -Inf, Ifrac = NULL, Iloss = NULL, dffun = df0, fullout = FALSE, xstart = NULL, Ksim = 1 )symdynN( r, amu, asd, f, d, d_sd, d_cov, N, sf, tmax, stochd = TRUE, stocht = TRUE, as.matrix = FALSE, amax = 0, amin = -Inf, Ifrac = NULL, Iloss = NULL, dffun = df0, fullout = FALSE, xstart = NULL, Ksim = 1 )
r |
per-capita growth rate (r in Eq.2-3) |
amu |
the mean interaction strength |
asd |
standard deviation used to generate interaction strengths |
f |
the waiting time (or average waiting time) between disturbance events (equal to 1/lambda in Eq.2-3) |
d |
mean size of disturbance function (mu in Eq.2-3) |
d_sd |
standard deviation of disturbance function (sigma in Eq.5) |
d_cov |
the covariance for generating disturbances |
N |
number of species or patches |
sf |
waiting time between sampling events |
tmax |
the time series length to be simulated |
stochd |
a logical variable, indicating whether disturbance size should be stochastic - otherwise, all disturbances are of magnitude d - defaults to TRUE |
stocht |
a logical variable, indicating whether waiting time between disturbance events should be stochastic - otherwise, waiting time is always f - defaults to TRUE |
as.matrix |
indicates whether results should be returned as matrix (potentially faster for some applications) - defaults to FALSE |
amax |
the maximum value allowed for interaction coefficients - defaults to zero |
amin |
the minimum value allowed for interaction coefficients - defaults to -Inf |
Ifrac |
dispersal rate (D in Eq. 2) - defaults to NULL (i.e. no dispersal) |
Iloss |
loss rate from dispersal that falls outside of the patch - defaults to NULL |
dffun |
the function handed to the ODE solver - should be df_col for spatial simulations, and df0 for multi-species simulations - defaults to df0 |
fullout |
a logical, determining whether the full output or just a summary is returned - defaults to fullout |
xstart |
optional vector of starting abundances - defaults to NULL (i.e. no values) |
Ksim |
carrying capacities - defaults to 1 |
a matrix or data.frame with columns for sampling times, abundances, and number of disturbances for each time interval
### Example 1: 10 patches r<-1 #rate of recovery d<-(0) #mean size of disturbance (mu in text) d_sd<-sqrt(0.1) #SD of disturbances (sigma in text) f<-1 #average time between disturbances (1/lambda in text) sf<-0.1 #sampling interval tmax<-120 #maximum time for simulation d_cov<-d_cov0<-(d_sd)^2/2 #covariance in disturbance size among patches xtNpatches<-symdynN(r = r, amu=0, asd=0, f=f, d=d, d_sd=d_sd, d_cov=d_cov, N=10, sf=sf, tmax=tmax, Ifrac=0, dffun = df_col) ### Example 2: 30 species r<-1 #rate of recovery d<-(0) #mean size of disturbance (mu in text) d_sd<-sqrt(0.1) #SD of disturbances (sigma in text) f<-1 #average time between disturbances (1/lambda in text) sf<-0.1 #sampling interval tmax<-120 #maximum time for simulation d_cov<-0 #covariance in disturbances among species amu<-(-r/2) #average interaction coefficient asd<-0.1 #standard deviation of interaction coefficient xtNsp<-symdynN(r = r, amu=amu, asd=asd, f=f, d=d, d_sd=d_sd, d_cov=d_cov, N=30, sf=sf, tmax=tmax)### Example 1: 10 patches r<-1 #rate of recovery d<-(0) #mean size of disturbance (mu in text) d_sd<-sqrt(0.1) #SD of disturbances (sigma in text) f<-1 #average time between disturbances (1/lambda in text) sf<-0.1 #sampling interval tmax<-120 #maximum time for simulation d_cov<-d_cov0<-(d_sd)^2/2 #covariance in disturbance size among patches xtNpatches<-symdynN(r = r, amu=0, asd=0, f=f, d=d, d_sd=d_sd, d_cov=d_cov, N=10, sf=sf, tmax=tmax, Ifrac=0, dffun = df_col) ### Example 2: 30 species r<-1 #rate of recovery d<-(0) #mean size of disturbance (mu in text) d_sd<-sqrt(0.1) #SD of disturbances (sigma in text) f<-1 #average time between disturbances (1/lambda in text) sf<-0.1 #sampling interval tmax<-120 #maximum time for simulation d_cov<-0 #covariance in disturbances among species amu<-(-r/2) #average interaction coefficient asd<-0.1 #standard deviation of interaction coefficient xtNsp<-symdynN(r = r, amu=amu, asd=asd, f=f, d=d, d_sd=d_sd, d_cov=d_cov, N=30, sf=sf, tmax=tmax)
Extrapolate variance observed at spatial or ecological scale b to a different scale, B. Equivalent to Eq.7a in the main text.
var_scale(mvar_b, murho_b, mucov_b = NULL, b = 1, B)var_scale(mvar_b, murho_b, mucov_b = NULL, b = 1, B)
mvar_b |
Mean variance of abundance values observed at scale b. |
murho_b |
Mean Pearson correlation coefficient of abundance values observed at scale b, calculated as mucov_b/mvar_b. If NULL, mucov_b is used instead. |
mucov_b |
Mean covariance of abundances observed at scale b. Ignored if mrho_b is not NULL. Defaults to NULL. |
b |
Size of observed scale. Defaults to 1. |
B |
Size of desired scale for extrapolation. |
Extrapolated variance at scale B.
# extrapolate from scale of 1 to 10 - e.g. from a 1m2 patch to a 10m2 patch var_scale(mvar_b = 1, murho_b = 0.5, b = 1, B = 10) # example with 100 simulated species nsp<-100 # number of species var_b<-1 # species-level abundance variance cov_b<-(-0.01) # between-specie abundance covariance # note - if nsp is large, cov_b must be near zero # this is because, e.g. many variables cannot all be # simultaneously negatively correlated # make a covariance matrix based on var_b and cov_b sigmamat<-diag(nsp)*var_b+(1-diag(nsp))*cov_b # simulate 1000 observations of 100 species sim_x<-mvtnorm::rmvnorm(n=1e3, mean = rep(0,100), sigma = sigmamat) # calculate mean variance, covariance, and correlation from sim_x cvmat<-cov(sim_x) mvar_b<-mean(diag(cvmat)) mucov_b<-mean(cvmat[row(cvmat)!=col(cvmat)]) murho_b<-mucov_b/mvar_b # test function vs. observation # note - answers match exactly var(rowSums(sim_x)) var_scale(mvar_b, murho_b = murho_b, b=1, B=100)# extrapolate from scale of 1 to 10 - e.g. from a 1m2 patch to a 10m2 patch var_scale(mvar_b = 1, murho_b = 0.5, b = 1, B = 10) # example with 100 simulated species nsp<-100 # number of species var_b<-1 # species-level abundance variance cov_b<-(-0.01) # between-specie abundance covariance # note - if nsp is large, cov_b must be near zero # this is because, e.g. many variables cannot all be # simultaneously negatively correlated # make a covariance matrix based on var_b and cov_b sigmamat<-diag(nsp)*var_b+(1-diag(nsp))*cov_b # simulate 1000 observations of 100 species sim_x<-mvtnorm::rmvnorm(n=1e3, mean = rep(0,100), sigma = sigmamat) # calculate mean variance, covariance, and correlation from sim_x cvmat<-cov(sim_x) mvar_b<-mean(diag(cvmat)) mucov_b<-mean(cvmat[row(cvmat)!=col(cvmat)]) murho_b<-mucov_b/mvar_b # test function vs. observation # note - answers match exactly var(rowSums(sim_x)) var_scale(mvar_b, murho_b = murho_b, b=1, B=100)
Helper function for symdyn, used to simulate dynamics between disturbance events.
xt(t0, t1, B0, r)xt(t0, t1, B0, r)
t0 |
initial time step |
t1 |
desired time step |
B0 |
initial abundance |
r |
relative growth rate |
predicted value of x at time t1
# see xt2fun# see xt2fun
Function for solving for stability paramter values from observed time series. Equivalent to Eq.5 in the main text.
xt2fun(x0, r, d, d_sd, dt, ndist)xt2fun(x0, r, d, d_sd, dt, ndist)
x0 |
value of x^2 at time t (x(t) in Eq.5) |
r |
per-capita growth rate (r in Eq.5) |
d |
mean size of disturbance function (mu in Eq.5) |
d_sd |
standard deviation of disturbance function (sigma in Eq.5) |
dt |
time step (i.e. time between x0 and x1) - can be a vector of the same length as x0, or a number if all time steps are of equal length |
ndist |
number of disturbances in each time step (equivalent to p(t+tau) in Eq.5) - must be same length as x0 |
predicted value of x^2 at time t+dt
# simulate dynamics, with r=1, d=0, and d_sd=0.1 xtout<-symdyn(r=1, f=1, d=0, d_sd=0.1, sf=0.1, tmax=100) # abundance in current time step x0<-xtout$state[1:(nrow(xtout)-1)] # abundance at t+1 x1<-xtout$state[2:nrow(xtout)] dt<-diff(xtout$time) ndist<-xtout$disturbed[-1] # fit model - note square root transform of response variable, # and log transform of parameter values mod<-nls(sqrt(x1^2)~sqrt(xt2fun(x0, r=exp(lr), d=0, d_sd=exp(ld_sd), dt, ndist)), start=c(lr=log(1), ld_sd=log(0.1))) exp(coef(mod)) # model estimates# simulate dynamics, with r=1, d=0, and d_sd=0.1 xtout<-symdyn(r=1, f=1, d=0, d_sd=0.1, sf=0.1, tmax=100) # abundance in current time step x0<-xtout$state[1:(nrow(xtout)-1)] # abundance at t+1 x1<-xtout$state[2:nrow(xtout)] dt<-diff(xtout$time) ndist<-xtout$disturbed[-1] # fit model - note square root transform of response variable, # and log transform of parameter values mod<-nls(sqrt(x1^2)~sqrt(xt2fun(x0, r=exp(lr), d=0, d_sd=exp(ld_sd), dt, ndist)), start=c(lr=log(1), ld_sd=log(0.1))) exp(coef(mod)) # model estimates
Helper function for symdynN. Simulate an ODE given parameters, starting value, and times.
xtN(t0, t1, B0, odepars, dffun, nsteps = 2)xtN(t0, t1, B0, odepars, dffun, nsteps = 2)
t0 |
initial time step |
t1 |
desired time step |
B0 |
initial abundance |
odepars |
parameter list |
dffun |
function for calculating derivatives |
nsteps |
number of time steps to return - defaults to 2 |
a matrix of species abundances