Varying Demographic Rates Across Patches

[skeyser]

Apologies in advance if I have missed this in the extensive user manual documentation. I have to say the documentation for the framework is awesome!

I'm curious if it is possible to integrate variable patch-specific demographic rates (e.g., survival) in a patch-based model? For example, using a predictive model to project survival across habitat patches and then integrate this into the RangeShifter framework.

[JetteReeg]

Hi,
there is indeed a version of the R package of RangeShifter where spatially and temporally varying demographic rates can be implemented for models using habitat suitabilities.
You can install it in R using this command:

devtools::install_github("RangeShifter/RangeShiftR-pkg", ref="spatial_demography", subdir="RangeShiftR")

Please check the corresponding publication, which also includes the codes, which might be usefull for preparing the parameter master.
Please note that the input and output is different from the current release version of RangeShiftR to facilitate the use of bayesian tools. In addition, the link provided above is an updated version including two recent minor bugfixes.

I hope the description in the method section of the paper, the provided code and R help files ?ImportedLandscapes and ?StageStructure provide enough details. Let me know if you need any further technical assistance.

[skeyser]

Hi Jette,

Thanks for the prompt response and the wealth of information here.

I have looked over the publication and this looks awesome. Full transparency, I have not inspected the code thoroughly yet.

For clarification, I am curious if this approach can be translated to a discrete, patch-based model rather than a suitability-informed, cell-based model? For our work we are interested in modelling the persistence of a regional population network with ~22 discrete patches that we delineated. Currently, we are running RangeShiftr using pooled vital rates across the species range, but we know those vital rates are patch-specific and want to introduce that flexibility such that an individual in patch 1 get patch 1 vital rates, an individual in patch 2 get patch 2 vital rates, ....

My sense is that the resources you shared might hold the answer, but I wanted to clarify before I dig deeper into the code associated with that paper.

-Spencer

[JetteReeg]

Hi Spencer,

the feature is applicable for cell-based as well as patch-based models, but only using habitat suitabilities (HabPercent = T). Here are some Screenshots from the R help files:

From ?ImportedLandscape:

From ?StageStructure:

[skeyser]

Hi Jette,

I see, so assigning custom vital rates for a patch can be achieved by theoretically setting all values for the cell to the intended vital rate. For example, if patch 1 has Stage 3 Survival Rate (0.7) then the SurvLayer would have all patch 1 cells encoded as 0.7 and so forth a so on for the remaining patches? Hopefully I am interpreting this right, but in any case thank you for your help, Jette!

[JetteReeg]

Hi Spencer,

yes, you would need to set the maximal survival rate for the stages in the transition matrix and with the additional SurvLayer you can vary this value for spatially. E.g. giving all cells of patch 1 a value of 70.0. so e.g. if you set the maximal survival rate of stage 3 to 1.0 and all cells of patch 1 a value of 70 in the SurvLayer, the survival rate of stage 3 would be 0.7 in patch 1.

Best, Jette

[JetteReeg]

Hi Spencer,
could you please provide me with the code for the demogScaleLayers and the ImportedLandscape?

TLDR: the demogScaleLayers must be a list of arrays. Check if that is the case.

Explanation:
I just quickly checked where the error message is included in the code and found that it is occuring twice and is wrong/misleading once, because it doesn't assume dynamic landscapes at this point:

The screenshot shows the validation check of the ImportedLandscape function, specifically the check of the parameter demogScaleLayers:
If the parameter nrDemogScaleLayers>0 and the length(demogScaleLayers) > 0 is true, it checks whether any entry of the list demogScaleLayers is not of class array (any(sapply(object@demogScaleLayers, class) != "array")).
Following the logic of that check, your demogScaleLayers should be a list of arrays and if not all entries of that list are arrays, it prints the error message you got.
The last part of the message about DynamicLandYears is wrong at this place, because it only checks whether all entries are correct. Only in the following step it would check whether the simulation includes dynamic landscapes (and only then it is necessary that the length of the list demogScaleLayers is equal to the length of DynamicLandYears).

Hope this helps.

[skeyser]

Hi Jette,

Thanks for your reply. I found the source code per your earlier link to the Malchow et al. 2021 paper. I see there are a couple sticking issues that happen when going through the ifelse ladder. The latest one that I am stuck on is the "any(apply(ds_dims,1,var)!=0" fails because I am trying to supply only demogScaleLayers for survival. That line produces an NA since I only have one layer.

Does this mean I need to supply dummy layers for each demographic rate (Survival, Development, and Fecundity) even though I only want to vary adult survival?

My demogScaleLayer looks like this (My study region is irregular and I have NAs coded as -9999, but can change if need be):
`

class(demog)
[1]"array"

dim(demog)
[1] 388 727 4

demog[1:2, 1:2, 1:4]
, , 1

  [,1]  [,2]

[1,] -9999 -9999
[2,] -9999 -9999

, , 2

  [,1]  [,2]

[1,] -9999 -9999
[2,] -9999 -9999

, , 3

  [,1]  [,2]

[1,] -9999 -9999
[2,] -9999 -9999

, , 4

  [,1]  [,2]

[1,] -9999 -9999
[2,] -9999 -9999`

[JetteReeg]

Hi Spencer,
I had another look at the checks and found the error causing your latest issue. The model version was used/checked with a dynamic landscape model. The check causing your error message (that all arrays in DemogScaleLayers should have the same size) was only valid if users would use dynamic landscapes and thus provide multiple 3D arrays for the DemogScaleLayers. I changed this check so it will also work if users only use one 3D array. I also changed the first misleading error message and another minor bug I found while running a testscript.

Here is an example script that ran without error messages:

# # reinstall version with the bugfixes: (uncomment if you need to run the code)
# library(Rcpp)
# library(devtools)
# library(Rdpack)
# devtools::install_github("RangeShifter/RangeShiftR-pkg", ref="spatial_demography", subdir="RangeShiftR")
# 

# Example setup for demographic variability
library(RangeShiftR)

##--- create example landscape with demog scale layers
landscape <- list(matrix(c(100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100), nrow=4, byrow=TRUE))


# determine array dimensions
dimX <<- ncol(landscape[[1]])
dimY <<- nrow(landscape[[1]])
nlayers <- 1
dim_array = c(dimY,dimX,nlayers)
# create array of dynamic scalings with empty placeholder layers
tmp_scaleLayers <- rep(list(array(rep(landscape[[1]],nlayers), dim = dim_array)),1)


land <- ImportedLandscape(LandscapeFile = landscape,
                          Resolution = 100,
                          OriginCoords = c(0,0),
                          HabPercent = TRUE, K_or_DensDep = 10,
                          demogScaleLayers = tmp_scaleLayers,
                          nrDemogScaleLayers = 1)

##--- demography

nr_stages <- 3

demog <- Demography(ReproductionType = 0,
                    StageStruct = StageStructure(
                        Stages = nr_stages, 
                        TransMatrix = matrix(c(0, 0, 1,
                                              1, 0.5, 0,
                                              0,0,0.8),
                                             nrow = nr_stages,
                                             byrow = T),
                        FecLayer = c(NA,NA,1)) # the fecundity for stage 3 (adult) is spatially varying using the first layer of the provided demogScaleLayer 3D array; use NA if the default values of the TransMatrix should be used
)

##--- dispersal

disp <- Dispersal(Emigration = Emigration(EmigProb = 0.9),
                  Transfer   = DispersalKernel(Distances = 100),
                  Settlement = Settlement(Settle = 3))


##--- initialisation params

init <- Initialise(InitType = 0, 
                   FreeType = 1, InitDens = 1, PropStages = c(0,0.5,0.5))

##--- param master

s <- RSsim(batchnum = 1, land = land, demog = demog, dispersal = disp, init = init, simul = Simulation())

##--- create directory structure

dir.create("Model", showWarnings = FALSE)
dir.create("Model/Inputs", showWarnings = FALSE)
dir.create("Model/Outputs", showWarnings = FALSE)
dir.create("Model/Outputs_Maps", showWarnings = FALSE)

##--- run model

RunRS(s, dirpath="Model/")

You will see two messages while initialising the simulation:
"Dimensions of DevLayer/SurvLayer matrix do not match..."
These are just informing you that for development and survival the default values provided in the transition matrix are used because DevLayer and SurvLayer are not provided.

[skeyser]

Hi Jette,

On my end I'm getting an error when I run your example stating "Error in validityMethod(object) : object 'dims' not found". This is specific to creating the ImportedLandscape object. Perhaps the issue is a simple typo here?

if( !all(unlist(dims) == unlist(dims[[1]])) ){ msg <- c(msg, "all arrays in demogScaleLayers must have the same size.") }

needs to be

if( !all(unlist(ds_dims) == unlist(ds_dims[[1]])) ){ msg <- c(msg, "all arrays in demogScaleLayers must have the same size.") }

Also please forgive me if I am mistaken but I'm a little confused by these lines. Working with the example you provided the sapply(demogScaleLayers, dim) function returns a "matrix" "array" class type.

In the example you provides the ds_dims object looks like this to me:

      [,1]
[1,]    4
[2,]    4
[3,]    1

Seems like the if statement will fail because it would require a perfect 3-D array (ncol==nrow==slices) to match unlist(ds_dims[[1]]). Perhaps a conditional to catch if the user supplies a single element in the list would help? Additionally, if the user supplies multiple elements of the list (which I assume corresponds to the dynamic land years) this argument could check the output from the function that generates ds_dims with a function like all(apply(ds_dims, 2, function(col) identical(col, ds_dims[,1]))) to check that all of the columns are the same specs? Just a thought but I could be off-base here.

[skeyser]

I was just trying to test this a little bit more...it seems like the example you provided works with this edit:

if( !all(apply(ds_dims, 2, function(col) identical(col, ds_dims[,1]))) ){
                                        msg <- c(msg, "all arrays in demogScaleLayers must have the same size.")
                                    }

on lines 617-619 in the class_LandParams.R file. Happy to make a pull request or if you think this is a reasonable fix and are willing to implement it that works too!

The ImportedLandscape returned this on my end after the fix above:

Formal class 'ImportedLandscape' [package "RangeShiftR"] with 14 slots
  ..@ LandscapeFile     :List of 1
  .. ..$ : num [1:4, 1:4] 100 100 100 100 100 100 100 100 100 100 ...
  ..@ Resolution        : num 100
  ..@ OriginCoords      : num [1:2] 0 0
  ..@ HabPercent        : logi TRUE
  ..@ Nhabitats         : int 1
  ..@ K_or_DensDep      : num 10
  ..@ PatchFile         : list()
  ..@ CostsFile         : list()
  ..@ SpDistFile        : list()
  ..@ SpDistResolution  : num -9
  ..@ DynamicLandYears  : int 0
  ..@ nrDemogScaleLayers: num 2
  ..@ demogScaleLayers  :List of 1
  .. ..$ : num [1:4, 1:4, 1:2] 100 100 100 100 100 100 100 100 100 100 ...
  ..@ LandNum           : int 1

[JetteReeg]

Hi Spencer,
thanks for looking into it and detecting the bug. I reviewed your suggestion and committed it.
Best,
Jette

[skeyser]

Hi Jette,

Got it! Looks like that did the trick. Thanks for the quick reply.

I also noticed that when I try to retrieve individual movement paths via OutIntPaths = 1 (in our analysis we are using an SMS transfer model that uses a cost surface) in this version of RangeShiftR R doesn't play nicely and immediately aborts the session. Strangely I can produce output for the SMS heatmap but this indicates no dispersing inds despite setting density-independent emigProb = 1 for my natal dispersers. Does the spatial_demo version accommodate the SMS movement model the same way as the main branch (I was able to get dispersal to work a while back)? Any chance it could have to do with setting the demographic scaling layer values to 0 for the non-patches?

Ultimately, we want to see how sensitive our species' regional persistence is to both connectivity and spatially varying demographics so both features are important.

Best,
Spencer

[JetteReeg]

Hi Spencer,

Similar to the issue with dynamic landscapes, the functionality of this branch was only tested with dispersal kernels, but it should not make any difference because no SMS parameter is spatially varying. However, there might be issues with the input of the cost layers?!

I was only able to quickly test the SMS function and at least it runs without an error, but I was not yet able to check whether it really runs correctly.

Have you already tried running the simulation with the most recent main RangeShiftR version (and files as inputs / without adding the spatially varying demographic rates)? Could you provide me with an example code?

PS: I am not sure if I am able to have a closer look at this issue before the xmas break.

[skeyser]

Hi Jette,

Circling back to this after some time away for the holidays.

For the cost layer I used the same formatting as for the other layers.

cost <- rast("costLayer.asc")
cost <- terra::as.matrix(cost, wide = TRUE)
cost[is.nan(cost)] <- NA
dim(cost)
cost <- list(cost)

for the landscape

## Landscape params
land <- ImportedLandscape(LandscapeFile = lscp,
                          PatchFile = patch,
                          CostsFile = cost,
                          Resolution = 1000,
                          HabPercent = T,
                          K_or_DensDep = 0.004,
                          demogScaleLayers = demog, #Just adult survival
                          nrDemogScaleLayers = 1L,
                          DynamicLandYears = 0L
)

then in the dispersal module I had the SMS specified as follows

trans_sm <- SMS(PR = 1,
                PRMethod = 2,
                MemSize = 3,
                DP = 5, 
                Costs = "file",
                StepMort = 0.0)

I have tried setting all the cells to '1' for the cost layer to see if it was an issue with some of the values and I have changed the initial patch to see if it was a proximity issue but neither solved the problem. The dispersal kernel works as expected and I can get my species to colonize nearby patches with this approach. My cost layer has values between 1 and 1000, but NAs for the boundary. It's an irregular shape near a coastline.

My previous attempts with the main branch before transition to the spatial demography worked as expected with the cost layer, but I have not tested it on the latest main branch.

[JetteReeg]

Hi Spencer,

I finally had the time to look into that issue and was able to tackle down the problem causing the error in the SMS functionality.
The cost maps were not transfered correctly to the main RangeShifter model due to a redeclaration of a variable inside an if-condition.
Please reinstall the RangeShiftR-package to get the latest version.

[skeyser]

Hi Jette,

Fantastic! The latest version of RangeShiftR I installed has the SMS functionality working. Thanks for your help on all of this.

Best,
Spencer