In cascade-wp6/socialecological: The economics of catastrophic shifts in arid rangelands.

This package vignette contains examples for the use of the functions provided in the package. It has the following main features:

provide a spatially explicit livestock resilience model for the cellular automata framework
initialise a list of landscape objects
update the list of landscape objects over a given period with a given set of parameters
summarize and analyse the distribution of the landscape states found in a list

the livestock model

The livestock model has been developed as a spatially explicit version of a graphical model deviced by Noy-Meir 1975^[Noy-Meir, I. 1975. Stability of grazing systems: an application of predator-prey graphs. Journal of Ecology 63, 459-481.] which explained alternative stable states in rangelands from the logistic growth of plants and the non-linear feeding response of herbivore grazers. The intersections of plant growth and mortality define the potential steady states of the system.

The livestock model adds spatially explicit positive and negative interactions to the terms of growth and mortality which are based on assumptions of local facilitation and local competition, as well as plant defensive traits against grazing and attractant decoy effects.

While these assumptions define the shape of growth and mortality along the dual gradient of vegetation cover and spatial structure, the livestock rates and basic productivity of the system define the absolute extend of mortality and growth respectively.

There is a range of parameter combinations that allow for alternative stable states of the system, that is, the vegetation trends towards a different attractor if starting from high cover and clustering than it would at low cover and little clustering. The history of the ecosystem defines about the direction of its further development (Hysteresis).

The basins of attraction can be approached in a pair-approximation simplification of the model. However, only the stochastic noise of the cellular automata allows for an investigation of likelihood of future outcomes. Thus, by adding further stochastic elements to the model, we can assess uncertainties in terms of future revenue or success of a management method.

More details on the model can be found in the vignette example.

generate landscapes and run simulations using caspr

The caspr package provides a framework to run cellular automata simulations in R. Many of the functions provided in this example refer to functions of caspr. So first, both packages have to be loaded into your current R session.

library(socialecological)
library(caspr)

Now we can generate landscapes using the function init_landscape():

l <- init_landscape(
          states= c("1","0"), # potential cell states: 1 = vegetated, 0 = empty ground
          cover = c(0.7,0.3), # initial random cover 
          width = 50         # width and height of the landscape object
          )
plot(l)

The function ca() is the core function of caspr and applies a cellular automata simulation to the given landscape with a set of parameters:

p <- list(  
  r = 1.0,  # max. regeneration rate of plants
  b = 0.3,  # environmental quality
  sigma = 0, # environmental noise
  f = 0.6,  # local facilitation
  alpha = 0, # water runoff
  K = 0.9, # carrying capacity of the system
  c = 0.2, # local competition
  m = 0.05, # intrinsic mortality of plants (inverse of av. lifespan)
  v = 0.8, # attractant-decoy
  p = 0, # associational resistance
  L = 1, # Livestock density
  q = 0, # hill exponent of functional response
  h = 10, # handling time 
  a = 5 # attack rate of livestock
) 

r <- ca(l, model = livestock, parms = p, t_max = 100)

The plot function knows how to handle output of this model run. It automatically plots a timeseries of vegetation cover. Alternatively, single snapshots of landscapes can be assessed at any timestep.

plot(r) # plot timeseries of vegetation cover

par(mfrow = c(1,4), mar = c(1,1,1,1)) # plot snapshots of landscapes at certain timesteps
plot(r$landscapes[[1]])
plot(r$landscapes[[25]])
plot(r$landscapes[[50]])
plot(r$landscapes[[100]])

add random environmental noise

The parameters b and sigma define the environmental quality and the random variability of it across years. Thus, if specifying sigma, the climate effects vary randomly over time. In the following example code, the longtime outcome of the simulation is only little affected by high variation in environmental quality (red case). However, in situations close to the tipping point, high stochasticity may be inducing the shift earlier.

p$sigma = 0.01

r <- ca(l, model = livestock, parms = p, t_max = 100)
par(mar=c(4,4,1,4)); plot(r); axis(4); mtext("environmental quality", side = 4, line = 2)
temp <- with(p, b*rnorm(100, 1, sigma))
lines(temp)     

p$sigma = 0.5

run <- ca(l, livestock, p, t_max = 100, saveeach = 100)
lines(run$time, run$cover[[1]], col = "red")

temp <- with(p, b*rnorm(100, 1, sigma))
lines(temp, col = "red")

If facilitation effects are strong, the climate variation does not have much of an effect since local conditions are defined by parameter $f$ rather than environmental quality, $b$. However, if facilitation is of minor importance, the variation over time can be strongly affected by climatic noise. Generally, climatic stochasticity is expected to be of relevance particularly in systems close to the tipping point. The likelihood of degradation might be enhanced by an increased amplitude of the climate extremes, i.e. a large sigma.

set parameter time series

Manipulations to the parameters over time will be required for the intended research questions. Either because the questions address the impacts of ongoing climate change on arid landscape stability without intervention, or because the consequences of intervention by management methods are to be assessed. In both cases, a predefined time sequence of a parameter needs to be fed into the updating procedure of the cellular automata. This will be achieved by providing the model with vectors instead of single value parameters. Each vector has to have the length of the requested timeseries or it will be recycled, e.g. for rotational management measures.

Thus, as required by the caspr package, the parameter set is given as a list of parameters, but each list entry can be either a single value or a vector representing the sequence of parameter values over time (in years). This is inflated by the requested time to evaluate using the function parms_timeseries() of caspr.

p <- list(  
  r = 1.0,  # max. regeneration rate of plants
  b = seq(0.2,0.05, length = 100)*rnorm(100, 1, 0.2),  # environmental quality
  sigma = 0, # random annual variation of environmental quality
  f = 0.9,  # local facilitation
  alpha = 0, # water runoff
  K = 0.7, # carrying capacity of the system
  c = 0.2, # local competition
  m = 0.05, # intrinsic mortality of plants (inverse of av. lifespan)
  v = 0.2, # attractant-decoy
  p = 0.99, # associational resistance
  L = rep(c(0,3.5), each = 2), # Livestock density
  q = 0, # hill exponent of functional response
  h = 10, # handling time 
  a = 5 # attack rate of livestock
) 


plist <- parms_timeseries(p, 100)

l <- init_landscape(c("1","0"), cover = c(0.8,0.2), width = 100)

r <- ca(l, livestock, plist, t_max = 100, saveeach = 5)

par(mar = c(4,4,1,4), las = 1)
plot(r)
axis(4); mtext("environmental quality (blue)", side = 4, line = 2, las = 0)
lines(1:100, plist[,"b"], col = "blue")

lists of landscape objects

In a next step, we want to explore the space of possible outcomes under a given set of parameters. This is achieved by running the model iteratively, starting from steady state situations, and by mapping the trajectories and analysing the final state of the model after a couple of years.

timeseries simulations and steady state

I wrote a set of functions which take a list of landscape objects (created by function init_list()) and updates it with the given parameter set for the requested time period (by applying function update_list()). If a parallel backend is provided, the single runs are processesed in parallel. The output is an updated list of landscape objects which can be summarized into a report of the distribution of outcomes.

All simulations have to start out from a pre-formed landscape with realistic patch structure, since the spatial structure is what determines the future development of the landscape. That is particularly important if looking at only short time spans such as 5 or 50 years. For now, we assume that steady state is achieved after a period of 100 years.

L0 <- init_list(100, runif_range = c(0.6,0.99), width = 25)

p <- list(  
  r = 1.0,  # max. regeneration rate of plants
  b = 0.3,  # environmental quality
  sigma = 0, #
  f = 0.7,  # local facilitation
  alpha = 0, # water runoff
  K = 0.9, # carrying capacity of the system
  c = 0.0, # local competition
  m = 0.05, # intrinsic mortality of plants (inverse of av. lifespan)
  v = 0.0, # attractant-decoy
  p = 0.0, # associational resistance
  L = 2, # Livestock density
  q = 0, # hill exponent of functional response
  h = 10, # handling time 
  a = 5 # attack rate of livestock
) 

L100 <- update_list(L0, 100, p)

summary(L100)

Now, we can think of many different scenarios, e.g. the application of a management method for restoration. Let's say after a long period of intensive grazing, the landscape is to be grazed only by one tenth of the livestock rates.

p$b <- 0.1
p$L <- 2

L110 <- update_list(L100, 10, p)
L120 <- update_list(L110, 10, p)
L130 <- update_list(L120, 10, p)
L140 <- update_list(L130, 10, p)
L150 <- update_list(L140, 10, p)

summary(L150)

Note that the summary function returns a set of statistics for a list of landscape objects. A call like summary(L100) returns the most obvious overview statistics. The single values for each landscape can be called using the $ selector, e.g. summary(L100)$cover will return a vector with the vegetation cover of each landscape object in the list. In the background, more information is computed that will be explained in the following sections.

landscape productivity

To translate vegetation cover into a measure of productivity in terms of input for the livestock production, we need to calculate the amount of plants that is likely to die of grazing in the upcoming timestep, as a proxy for the quantity of plant biomass that is eaten.

The function that produces this output is feed() which can be applied to a single lanscape object and returns the number of plants, normalised per hectar, that is about to fall to livestock feeding. Remember to enter the desired set of parameters, otherwise the function will calculate the feeding for the default parameterset!

feed(L100[[1]], parms = p)

The summary function internally calculates this for a given list of landscape objects. It can be accessed as a vector by calling

summary(L100, parms = p)$feed

visualisation and analysis

Violin plots represent the approximate distribution of the iterations current cover and can be interpreted as probabilities. The distribution pattern varies with the number of replicates, the amount of environmental noise and size of the landscape.

library(vioplot) 

plot(NA,NA, 
     xlim = c(0,50), ylim  = c(0,1), 
     xlab = "time (years)", 
     ylab = "vegetation cover" )

vioplot(summary(L100)$cover, 
        at = 0, add = T, col = 'white', 
        wex = 3, pchMed = "-", colMed = "white")

# continue running in steps of 5 years but at increased 
# environmental quality

L <- L100
p$b <- 0.15
p$sigma = 0.2
p$L <- 1.7

for(i in seq(5,50, 5)) {
L <- update_list(L, 5, p)

vioplot(summary(L)$cover,
        at = i, add = T, col = 'white',
        wex = 3, pchMed = "-", colMed = "white")
}

translation into biomass production for livestock

More interesting than analysing cover would be analysing the output of the landscape for livestock production, i.e. the amount of plants that will be dying due to grazing as estimated by the function feed() (see above).

plot(NA,NA, 
     xlim = c(0,50), ylim  = c(0,0.3), 
     xlab = "time (years)", 
     ylab = "plant output" )

# continue running in steps of 5 years but at increased 
# environmental quality

L <- L100
p$b <- 0.15
p$f <- 0.6
p$sigma <- 0.2
p$L <- 1.1

vioplot(summary(L, parms = p)$feed, 
        at = 0, add = T, col = 'white', 
        wex = 3, pchMed = "-", colMed = "white")

for(i in seq(5,50, 5)) {
L <- update_list(L, 5, p)

vioplot(summary(L, parms = p)$feed,
        at = i, add = T, col = 'white',
        wex = 3, pchMed = "-", colMed = "white")
}

This distorts the meaning of plant productivity in the light of plant biomass that is accessible to livestock feeding. It impressively illustrates the importance of the associational resistance which has its maximal effect at high cover.

density Kernel estimate

The violin plot is a combination of a boxplot and a distribution kernel. The vioplot function applies a gaussian data structure to approximate the density estimate. It cuts the violin plot at 0 and 1. This is not mathematically reflecting the structure of the proportional landscape cover data.

Instead, the summary function will calculate the kernel based on a beta distribution, which assumes proportional data limited between 0 and 1.

The kernel itself is stored in the output created by the summary function and can be assessed for further analysis, e.g. to match it with economic revenue of a particular landscape state. The boundary corrected kernel density estimation is derived using the function dbckden() of package evmix and is wrapped into the function betakernel() provided with this package. It defaults to parameters: bcmethod = "beta1", lambda = 0.2, xmax = 1.

plot(summary(L100)$kernel, type = "l", xlim = c(0,1), ylim = c(0,13), col = "grey60")

lines(summary(L110)$kernel, type = "l", col = "grey30")
lines(summary(L120)$kernel, type = "l", col = "grey10")

The package also provides a vioplot function that plots the beta distributed density kernel estimate: vioplot_beta(). Note that different from the original vioplot() function, this does not assume a cutoff at a minimum and maximum of the data range.

vioplot_beta(summary(L100)$feed, summary(L110)$feed,
             summary(L120)$feed, summary(L130)$feed,
             summary(L140)$feed, summary(L150)$feed,
        names = c("0","10","20","30","40","50"),
         add = F, col = 'white', 
        wex = 0.6, pchMed = "-", colMed = "white")