In KevCaz/seedlingsRecruitmentAnalysis: Research compendium of Solarik et al. (2019) DOI: 10.1111/1365-2745.13311

library(seedlingsRecruitment)

This vignette aims at guiding the reader through the different steps of the analysis carried out in Priority effects will prevent range shifts of temperate tree species into the boreal forest.[^not1].

Data

Data are currently archived on Dryad at the following URL https://doi.org/10.5061/dryad.q573n5tdx and can be download using retrieveData() as follows

retrieveData()

This create a folder (by default it is input) with all datasets as .csv or .xlsx files. Once done, trees locations and regeneration files can readily be formatted with formatData(). Note that, as favuorability data are strored as .docx files,in order to use the code below it requires to convert the .docx into .csv files.

dat_abitibi <- formatData("input/abitibitrees.csv", "input/abitibiregen.csv",
  "input/abitibisubstrate.csv")
dat_lebic <- formatData("input/lebictrees.csv", "input/lebicregen.csv",
  "input/lebicsubstrate.csv",  abbr_site = "bic")
dat_sutton <- formatData("input/suttontrees.csv", "input/suttonregen.csv", "input/suttonsubstrate.csv", abbr_site = "sut")

By default all the files created will be exported as .rds file in output. Note that we stored all these data as .rda files in data. For every site data are list of five objects:

names(dat_abitibi)

• trees that includes identity and position of trees;
• regen2015 and regen2016 regeneration and favourability data for 2015 and 2016 respectively (favourability values are the same for both years);
• lsdist a list of distances: for every plot it gives the distance of all surrounding tree including in the site;
• iddist a list of logicals of the same format as lsdist that states whether the tree is included in the buffer.

See, for instance, ?dat_abitibi

Sites of the study

As explained in the original study, recruitment data have been measured on three different sites:

• 'abi': Abitibi
• 'bic': Le Bic
• 'sut': Sutton

The map below shows the location of those sites. For each site, extra information are available on click.

datS <- data.frame(
  id = c('abi', 'bic', 'sut'),
  name = c('Abitibi-Temiscaminque', 'Le Bic', 'Sutton'),
  lon = c(-79.401219, -68.817719, -72.541297),
  lat = c(48.162539, 48.333619, 45.112803),
  plots= c(400, 1024, 2000)
)
##
rownames(datS) <- datS$id
spSite <- sp::SpatialPointsDataFrame(
  datS[,c('lon', 'lat')],
  data = datS,
  proj4string = sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
 )
 mapview::mapview(spSite)@map

Tree species of the study

| abbr1 | abbr2 | Scientific Names | TSN |abi |bic |sut | |:------|:------|:----------------------|:-------|:-----|:-----|:-----| | BF | ABBA | Abies balsamea | 18032 |TRUE |TRUE |TRUE | | RM | ACRU | Acer rubrum | 28728 |TRUE |TRUE |FALSE | | SM | ACSA | Acer saccharum | 28731 |TRUE |TRUE |TRUE | | YB | BEAL | Betula alleghaniensis | 19481 |FALSE |FALSE |TRUE | | WB | BEPA | Betula papyrifera | 19489 |TRUE |TRUE |TRUE | | AB | FAGR | Fagus grandifolia | 19462 |FALSE |FALSE |TRUE | | AS | POTR | Populus tremuloides | 195773 |TRUE |TRUE |TRUE |

: Abbreviations, names, TSN (Taxonomic Serial Number, visit the Integrated Taxonomic Information System website for more details about TSN) of the tree species considered and their respective presence in the three sites.

Recruitment data

A link to the dataset would be added once data are released.

Statistical models

Below are detailed the steps to reproduce our analysis.

Models

General equation

For a given species and a given year, the seeds recruitment depends on:

1. the fecundity of surrounding individuals of the same species;
2. the substrate favourability;
3. the neighborhood effect.

$$R_i = \text{STR}\sum_{j=1}^Sc_{i,j}f_{i,j}\sum_{k=1}^T\left(\frac{\text{DBH}k}{30}\right)^2 \frac{1}{K}e^{-\text{DM}{i,k}}e^{-\frac{h_i}{1-h_i}\sum_{l=1}^Ba_{i,l}}$$

Number of parameters per process

| Ecological process |abbr |number of parameters| |:--------------------|:-----|:-------------------| | dispersion |disp |2 | | neighborhood effect |neigh |5 | | favorability |favo |1 |

Kernels

| numb | disp | which | clip | |:-----|:-----|:-----------------|:------| |0 |-- |non-significant | - | |1 |dc |lognormal |clip | |2 |d- |lognormal |no-clip| |3 |Dc |power-exponential |clip | |4 |D- |power-exponential |no-clip|

Simulations

Overall strategy

For each model, we estimate the parameters through the maximization of the likelihood (actually the minimization of the likelihood log-transformed) for the corresponding data. Every model is a combination :

• the site and the tree species considered;
• the age of seedling taken into account (1 or 2);
• whether we used a zero-inflated poisson or a regular poisson distribution;
• the kernel dispersal used;
• the distribution used to model seedling recruitment (zero-inflated poisson or a regular poisson).

As an example, below are the 15 first models out of the r nrow(simuDesign)/7 tested for Abies balsamea:

data(simuDesign)
# NB: each row is a model
knitr::kable(head(simuDesign, 15))

As the models we used are not linear and highly customized we optimized the likelihood with the simulated annealing algorithm implemented in GenSA[^gensa].

Numerical implementation

For each row of simuDesign, columns are arguments passed to recruitmentAnalysis() in which:

1. parameters are computed for the simulated annealing (getParameters());

2. simulated annealing is performed (GenSA() from GenSA package) to minimize the likelihood (getLikelihood()).

knitr::kable(getParameters(TRUE, TRUE, TRUE, lognormal = FALSE))

Example

The line below (not run) allows to run for the 16^th model.

res <- launchIt(iter = 16, simu = simuDesign, mxt = 100,
  record = "output/test.txt", quiet = TRUE, path = "output/")
# value (likelihood)
res$value
# parameters value
res$par

Note that test.txt records all the values tested during the simulated algorithm and the corresponding likelihood based on which we obtain a confident interval.

getConfInt("inst/res/test.txt", simuDesign, 16, basename_out = "inst/resf_")

Figure 5

Once all models run and the best models identified (which took days on a server and we thank Calcul Quebec for making such computations possible), we draw figure 5 that is the dispersal kernels inferred for one year seedlings

figKernels(res_bestModels)

Similarly, we used this function to draw the same figure for 2 years seedlings (see Supplementary Information of the paper):

figKernels(res_bestModels, age = 2)

References

[^not1]: Solarik et al. DOI: 10.1111/1365-2745.13311 [^gensa]: Y. Xiang, S. Gubian. B. Suomela, J. Hoeng (2013). Generalized Simulated Annealing for Efficient Global Optimization: the GenSA Package for R. The R Journal, Volume 5/1, June 2013. https://journal.r-project.org/archive/2013/RJ-2013-002/index.html

KevCaz/seedlingsRecruitmentAnalysis documentation built on Oct. 23, 2019, 12:58 p.m.