R/setup.R
In gabackus/rcomms: Simulate Randomized Metacommunities under Stochastic Climate Change
Defines functions
commSetup
manageSetup
nSetup
Documented in
commSetup
manageSetup
nSetup
#' Setup Model Structure
#'
#' This sets up the basic structure of the model.
#' It creates a list of biological parameters in \code{P} (\code{S} randomized species and their parameters) and it creates a list of environmental parameters in \code{X}.
#' All parameters have default values, so input is not necessary.
#'
#' @param S Total number of species created with hetSetup.
#' @param L Total number of patches in the metacommunity.
#' @param W Number of microhabitats in each patch.
#' @param zo A vector of pre-defined optimal temperature values. Only works if length(zo) is S.
#' @param gam A vector of pre-defined mean dispersal distances. Only works if length(gam) is S. If no vector is specified, gamMean and gamSD are used to randomly generate the vector gam.
#' @param sig A vector of pre-defined thermal tolerance breadths. Only works if length(sig) is S. If no vector is specified, sigMean and sigSD are used to randomly generate the vector sig.
#' @param A Matrix of the relative competition coefficients between species.
#' @param m A vector of pre-defined mortality probabilities (or a single value that will be shared for all species). These are probabilies, so the value must be between 0 and 1 (inclusive). Only works if length(m) is 1 or S.
#' @param gamMean Mean dispersal distance for randomized species. Default is based on Urban et al. 2012.
#' @param gamSD Standard deviation of dispersal distance for randomized species. Default is based on Urban et al. 2012.
#' @param sigMean Mean thermal tolerance breadth for randomized species. Default is based on Urban et al. 2012.
#' @param sigSD Standard deviation of thermal tolerance breadth for randomized species. Default is based on Urban et al. 2012.
#' @param lam Skewness in thermal tolerance. Default is based on Urban et al. 2012. (to have a mild decrease moving toward colder temperatures and a sharp decrease moving toward warmer temperatures).
#' @param B Area of integrated birth rate over all T for each species.
#' @param compType The type of competition in a string. Must be either "lottery" or "temp"
#' @param XW Window of analysis (to remove edge effects)
#' @param temp2d A matrix of pre-defined temperatures over x. Only works if nrow(temp2d) is L and ncol(temp2d) is W.
#' @param tempLow Lowest mean temperature on linear temperature gradient. temp1d(L)=tempLow
#' @param tempHigh Highest mean temperature on linear temperature gradient. temp1d(1)=tempHigh
#' @param tempRev If tempRev=T, then temp1d(1)=tempLow and temp1d(L)=tempHigh.
#' @param tempGH Hurst exponent for global temperature heterogeneity. H=0.5 is Brownian motion; 0.5<H<=1 is long-term positive autocorrelation and 0<=H<0.5 is long-term negative autocorrelation.
#' @param tempGSD Controls the magnitude of the global temperature heterogeneity.
#' @param tempLSD Standard deviation in temperature between microhabitats in each patch.
#' @param years Maximum number of years for initialization + climate change.
#' @param y Current year.
#' @param tempY A vector of pre-defined temperatures over time. Only works if length(tempY) is years.
#' @param tau Average temperature change per year.
#' @param tempYAC Temperature autocorrelation over time. Default is global temperature AC from 1880-1979.
#' @param tempYSD Temperature standard deviation over time. Default is global temperature SD from 1880-1979.
#' @param Tau Average temperature change per year (in a matrix). If
#' @param tempYXGH Hurst exponent for global heterogeneity in temperature change over time.
#' @param tempYXGSD Magnitude of global heterogeneity in temperature change over time.
#' @param tempYXLSD Standard deviation in temperature change between microhabitats in each patch.
#' @param tempYXGLH Hurst exponent for global heterogeneity in temperature change over time and space.
#' @param tempYXGLSD Standard deviation in temperature change between microhabitats in each patch over space.
#' @param Q A matrix of pre-defined habitat quality over x. Only works if nrow(Q) is L and ncol(Q) is W.
#' @param QMean Average habitat quality for any microhabitat.
#' @param QGH Hurst exponent for global habitat quality heterogeneity. H=0.5 is Brownian motion; 0.5<H<=1 is long-term positive autocorrelation and 0<=H<0.5 is long-term negative autocorrelation.
#' @param QGSD Controls the magnitude of the global habitat quality heterogeneity.
#' @param QLType How is the habitat quality lost? Must be "none", "subtract", or "multiplicative".
#' @param QLSD Standard deviation in habitat quality between microhabitats in each patch.
#'
#' @return A list that contains P (a list of biotic paraters) and X (a list of abiotic parameters)
#' @export
#'
#' @examples
#' cSetup<-commSetup(S=4,L=512,W=4,years=202)
#' P<-cSetup$P
#' X<-cSetup$X
#' n0 <- nSetup(8,4,512,4)
#' initSim <- commSimulate(n0,P,X,y=1,years=100,init=T)
#' n1 <- initSim$n
#' ccSim <- commSimulate(n1,P,X,y=101,years=100,init=F)
#' vComTime(ccSim$N)
#'
commSetup <- function(S=32, L=512, W=8,
zo=NULL, gam=NULL, sig=NULL, A=NULL, m=1,
gamMean=2.5, gamSD=2.5,
sigMean=5, sigSD=5,
lam=-2.7, B=10,
compType="lottery",
XW=seq(129,384),
temp2d=NULL, tempLow=9.78, tempHigh=30.22, tempRev=F, tempGH=1, tempGSD=0, tempLSD=1, tempGLH=0.81,tempGLSD=2, tempLSDRed=0,
years=1000, y=1,
tempY=NULL, tau=0.04, tempYAC=0.767, tempYSD=0.1639,
Tau=NULL, tempYXGH=1, tempYXGSD=0, tempYXLSD=0, tempYXGLH=1, tempYXGLSD=0,
Q=NULL, QMean=8, QGH=1, QGSD=0, QLSD=0, QLType="none", QLoss=0){
##########################################################
# First, specify biological parameter values for all of the S species.
# Optimal temperature for each species. These can be randomly picked from a uniform distribution or pre-defined.
if(is.null(zo)){
zo <- runif(S,9.9,30.1)
} else {
if(length(zo)!=S){
stop("zo does not match the number of species!")
}
}
# Dispersal distance for each species. These can be randomly picked from a lognormal distribution or pre-defined.
if(is.null(gam)){
# To use the lognormal distribution, we need to convert mean and SD values
gamMu <- log(gamMean/sqrt(1+gamSD^2/gamMean^2))
gamSig <- sqrt(log(1+gamSD^2/gamMean^2))
gam <- rlnorm(S,gamMu,gamSig)
} else {
if(length(gam)!=S){
stop("gam does not match the number of species!")
}
}
# Thermal tolerance breadth for each species. These can be randomly picked from a lognormal distribution or pre-defined.
if(is.null(sig)){
# To use the lognormal distribution, we need to convert mean and SD values
sigMu <- log(sigMean/sqrt(1+sigSD^2/sigMean^2))
sigSig <- sqrt(log(1+sigSD^2/sigMean^2))
sig <- rlnorm(S,sigMu,sigSig)
} else {
if(length(sig)!=S){
stop("sig does not match the number of species!")
}
}
# Competition coefficients between each pair of species species. By default, all coefficients are 1 (lottery competition), but A can be pre-defined.
if(is.null(A)){
A <- matrix(1,S,S)
}else{
if(!(nrow(A)==S & ncol(A)==S)){
stop("A does not match the number of species!")
}
}
# Yearly mortality probability for each species.
if(any(m<0) | any(m>1)){
stop("m should be between 0 and 1 (inclusive)")
}
if(length(m)==1){
# If only one value is provided, all species will have the same mortality
m <- rep(m,S)
} else if(length(m)!=S){
stop("m does not match the number of species!")
}
# Mortality is more convenient in this form
M <- rep(m,W*L)
# Make sure that compType is used correctly
if(!(compType=="lottery" | compType=="temp")){
stop("compType must be either 'lottery' or 'temp'!")
}
##########################################################
# Using the randomized species parameters, we derive other variables needed for computation.
# zo helps define where the species' thermal optimum is, but mathematically this is not completely correct.
# If we let zo be the true optimum, z is the value we plug into the reproduction function so that argmax(b)==zo
# We apply the function zAdjust to all values of zo to calculate z
z <- mapply(zAdjust, sig, zo, lam, 2^13)
# To speed up computation time, we define full dispersal kernels now.
# The dispersal kernel uses q, a transformation of gam
q <- sapply(1:S, function(i) 1+1/gam[i]-sqrt(1+1/gam[i]^2))
l <- (-L):(L)
k <- t(sapply(1:S, function(i) doubGeom(l,q[i])))
k[k<10^(-15)] <- 0.0
# K is a list of separate L by 2L dispersal kernel matrices
# K[[s]] is the dispersal kernel matrix of species s
# K[[s]][i,] is a vector of probabilities for a propagule in patch i to spread to patch j-L/2 (this is extra long to account for a propagule spreading beyond the limits of the ecosystem)
K <- rep(list(matrix(0,L,2*L)),S)
for(i in 1:S){
Ki<-matrix(0,2*L,4*L)
for(j in 1:(2*L)){
Ki[j,j:(j+2*L)]<-k[i,]
Ki[j,3/2*L]<-sum(Ki[j,1:(3/2*L)])
Ki[j,5/2*L+1]<-sum(Ki[j,(5/2*L+1):(4*L)])
}
K[[i]]<-Ki[(L/2+1):(3*L/2),(3/2*L):(5/2*L)]
}
# Tolerance and reproductive strength have a tradeoff
# Birth rate is adjusted so each species has roughly equal birth rate when integrated over all T
# ro is a constant that adjusts to this reproductive output
ro <- sapply(sig,function(x) rAdjust(x,B,lam,1e-06,2^13))
##########################################################
# Put the biological parameters together into a single list, P
P=list(S=S,
z=z,
gam=gam,
sig=sig,
lam=lam,
A=A,
M=M,
ro=ro,
zo=zo,
K=K,
compType=compType)
##########################################################
# Next, we define environmental parameters.
# Discrete spatial domain: from 1 to L (integers)
x <- seq(1,L)
# temp2d is the current temperature over all x and microhabitats
if(is.null(temp2d)){
temp1dr <- seq(tempHigh,tempLow,length=L)
if(tempRev){
temp1dr<-rev(temp1dr)
}
tempG <- tempVarH(L,tempGH)
temp1d <- temp1dr+tempG*tempGSD/sd(tempG)
if(W>1){
tempLSDx <- tempVarH(L,tempGLH)
tempLSDx2 <- tempLSD*( (1-tempLSDRed) + tempLSDRed/(1+exp(-tempLSDx*tempGLSD/sd(tempLSDx))))
temp2dr <- matrix(temp1d,L,W)
tempLH <- matrix(rnorm(L*W,0,1),L,W)
tempLH2 <- t(sapply(1:L, function(i) mean(tempLH[i,])+(tempLH[i,]-mean(tempLH[i,]))*tempLSDx2[i]/sd(tempLH[i,])))
tempLH3 <- t(matrix(sapply(1:L, function(x) sort(tempLH2[x,]-mean(tempLH2[x,]))),W,L))
temp2d <- temp2dr+tempLH3
} else{
temp2d<-temp1d
}
} else {
if(!(nrow(temp2d)==L & ncol(tempsd)==W)){
stop("temp2d does not match environment size!")
}
temp1d<-rowMeans(temp2d)
}
# tempY is a vector of the temperature over time
if(is.null(tempY)){
tempY<-0:years
for (i in 1:years){
# A new epsi is calculated for each time step
tempY[i+1] <- tempYAC*tempY[i]+rnorm(1,0,tempYSD)*sqrt(1-tempYAC^2)
}
} else{
if(length(tempY)!=years+1){
stop("tempY does not match years!")
}
}
# Tau is the temperature change over time in each patch and subpatch
if(is.null(Tau)){
tauLW <- matrix(tau*100,L,W)
tauG <- tempVarH(L,tempYXGH)
tauLW2 <- tauLW + tauG*tempYXGSD/sd(tauG)
if(W>1){
tauLH <- matrix(rnorm(L*W,0,tempYXLSD),L,W)
tauLH <- t(matrix(sapply(1:L, function(x) tauLH[x,]-mean(tauLH[x,])),W,L))
} else{
tauLH <- 0
}
Tau <- tauLW2+tauLH
Tau <- Tau/100
} else{
if(!(nrow(Tau)==L & ncol(Tau)==W)){
stop("Tau does not match environment size!")
}
}
# Habitat quality could differ in space or with species, but we will keep it constant for now
if(is.null(Q)){
if(QLType=="none"){
Q<-matrix(QMean,L,W)
} else{
QG <- tempVarH(L,QGH)
QLH <- matrix(rnorm(L*W,0,QLSD),L,W)
if(QLType=="subtract"){
Q1<-QMean-QMean/(QLoss^-1*(1+exp(-QLSD*QG/sd(QG))))
} else if(QLType=="multiplicative"){
if(QLoss<=0){
QLoss<-0.01
}
Q1<-(QMean*QLoss^-1)/(1+exp(-QG*QGSD/sd(QG)))
}
Q <- array(Q1,c(L,W))+QLH
}
} else {
if(!(nrow(Q)==L & ncol(A)==W)){
stop("Q does not match environment size!")
}
}
##########################################################
# Put the abiotic parameters together into a single list, X
X=list(L=L,
x=x,
XW=XW,
temp1d=temp1d,
tau=tau,
Tau=Tau,
Q=Q,
tempY=tempY,
tempYAC=tempYAC,
tempYSD=tempYSD,
temp2d=temp2d,
tempLSDx=tempLSDx,
W=W)
##########################################################
# Export it all as a list
return(list(P=P,X=X))
}
#' Setup Management Structure
#'
#' This sets up the management options for a simulation.
#' It creates a list of several lists that include paramaters for each management option.
#' Many of these don't bother to check for errors, so make sure you're careful.
#' @param P List of biotic parameters
#' @param X List of abiotic parameters
#' @param years The number of years that the model is going to be run.
#' @param AM Assisted migration. TRUE or FALSE determine whether the model runs through AM during each time step. If left NULL, this function will not define any parameters for AM.
#' @param corridor Corridor establishment. TRUE or FALSE determine whether the model runs through corridor management during each time step. If left NULL, this function will not define any parameters for corridor.
#' @param habQual Change habitat quality. TRUE or FALSE determine whether the model runs through changing habitat quality during each time step. If left NULL, this function will not define any parameters for habQual.
#' @param locHet Change local heterogeneity. TRUE or FALSE determine whether the model runs through changing local heterogeneity during each time step. If left NULL, this function will not define any parameters for locHet.
#' @param velocHet Change heterogeneity climate velocity. TRUE or FALSE determine whether the model runs through changing heterogeneity in climate velocity during each time step. If left NULL, this function will not define any parameters for velocHet.
#' @param targs Which species are specifically targeted for assisted migration? A list of 0s (for not targeted) and 1s (for targeted). Must be a vector of length S. Could also be a vector of species number (index from P parameters). Lastly, can just be "all" to include all species.
#' @param eta Relocation threshold. Any population of a target species that falls below this is relocated. Should be a nonnegative single integer or a vector of length S.
#' @param tCD Cooldown time between relocations. AM is not repeated if tCD has not passed since the last AM event. Should be a single nonnegative integer or a vector of length S.
#' @param rho Proportion of the total population moved when AM is triggered. Should be a single number from 0 to 1 or a vector of length S.
#' @param mu Probability of surviving AM. Should be a single number from 0 to 1 or a vector of length S.
#' @param zEst Estimate of species' thermal optimum. Should be a vector of length S.
#' @param xLoc How far ahead (in the direction of the leading edge) of the location closest to the thermal optimum estimate are we centering the relocation? Should be a single integer or a vector of length S.
#' @param recRad How wide is the recipient location? The recipient location consists of the center location (xLoc in front of the thermal optimum) and each patch +/- recRad from this center. Should be a single positive integer or a vector of length S.
#' @param donor Which individuals do we take from the donor community? 1: randomly (from the "top"), 2: from the trailing edge, 3: from the leading edge. Possible future extensions are 4: from the middle and 5: randomly from the old donor population. Should be a single integer.
#' @param recipient What shape should the recipient population look like? 1: an evenly distributed box. Possible future extensions are 2: rounded, 3: triangular, 4: thermal tolerance shaped. Should be a single integer.
#' @param tLR When was the last time the species was relocated? Should be an integer vector of length S.
#' @param randPick Do you want to pick exactly rho of the population (rounded down) or would you rather randomly pick from a binomial distribution with rho as the probability? If you want to do this randomly, set randPick to T. By default, it is set to F.
#'
#' @return A list of parameters that tell the model when and how to enact various managment techniques.
#' @export
#' @examples
#' cSetup<-commSetup(S=32,L=512,W=4,years=202)
#' P<-cSetup$P
#' X<-cSetup$X
#' mSetup<-manageSetup(P,X,AM=T,targs = "all",eta=40,rho=0.6)
#' n0 <- nSetup(8,32,512,4)
#' initSim <- commSimulate(n0,P,X,y=1,years=100,init=T)
#' n1 <- initSim$n
#' ccSim <- commSimulate(n1,P,X,y=101,years=100,init=F)
#' vComTime(log(ccSim$N[rowSums(flatComm(n1)[,X$XW])>0,]))
#' lines(c(0,100),c(log(mSetup$AM$eta[1]),log(mSetup$AM$eta[1])))
manageSetup <- function(P,X,years=100,
AM=NULL,corridor=NULL,habQual=NULL,locHet=NULL,velocHet=NULL,
targs=NULL,eta=50,tCD=5,rho=0.8,mu=0.8,zEst=P$zo,xLoc=10,recRad=2,donor=1,recipient=1,tLR=NULL,randPick=F){
# Now we set all of this up
S <- P$S
if(is.null(AM)){
# If you have no plans for AM at all, there is no need to define anything but AM$AM.
AM<-list(AM=F)
} else{
if(targs=="all"){targs <- 1:S} else{
if(is.null(targs)){targs<-rep(0,S)}
if(length(targs)>S){stop("targ too long")} else if(length(targs)<1){stop("targ too small")}
if(all(targs %in% c(0,1))){if(!(length(targs)==S | length(targs)==1)){stop("If using a 0/1 vector, targs needs to be length S")} else if(length(targs)==1){tInds<-targs; targs<-rep(0,S); targs[tInds]<-1}}else
if(all(targs %in% 1:S)){tInds<-targs; targs<-rep(0,S); targs[tInds]<-1} else{stop("If using targs as index vector, at least one of your indices were not valid.")}
}
if(!(length(eta)==S | length(eta)==1)){stop("Wrong length")} else if(length(eta)==1){eta<-rep(eta,S)}
if(!(length(tCD)==S | length(tCD)==1)){stop("Wrong length")} else if(length(tCD)==1){tCD<-rep(tCD,S)}
if(!(length(rho)==S | length(rho)==1)){stop("Wrong length")} else if(length(rho)==1){rho<-rep(rho,S)}
if(!(length(mu)==S | length(mu)==1)){stop("Wrong length")} else if(length(mu)==1){mu<-rep(mu,S)}
if(!(length(xLoc)==S | length(xLoc)==1)){stop("Wrong length")} else if(length(xLoc)==1){xLoc<-rep(xLoc,S)}
if(!(length(recRad)==S | length(recRad)==1)){stop("Wrong length")} else if(length(recRad)==1){recRad<-rep(recRad,S)}
relTimes<-matrix(0,S,years)
if(is.null(tLR)){tLR<-rep(0,S)}
AM <- list(AM=AM,targs=targs,eta=eta,tCD=tCD,rho=rho,mu=mu,zEst=zEst,xLoc=xLoc,recRad=recRad,donor=donor,recipient=recipient,tLR=tLR,randPick=randPick,relTimes=relTimes)
}
# CORRIDOR PARAMETERS
if(is.null(corridor)){
# If you have no plans for corridor establishment at all, there is no need to define anything but corridor$corridor.
corridor<-list(corridor=F)
}
# HABITAT QUALITY PARAMETERS
if(is.null(habQual)){
# If you have no plans for increasing habitat quality at all, there is no need to define anything but habQual$habQual.
habQual<-list(habQual=F)
}
# LOCAL HETEROGENEITY PARAMETERS
if(is.null(locHet)){
# If you have no plans for local heterogeneity change at all, there is no need to define anything but locHet$locHet.
locHet<-list(locHet=F)
}
# CLIMATE VELOCITY HETEROGENEITY PARAMETERS
if(is.null(velocHet)){
# If you have no plans for changing climate velocity heterogeneity at all, there is no need to define anything but velocHet$velocHet.
velocHet<-list(velocHet=F)
}
manage <- list(AM=AM,
corridor=corridor,
locHet=locHet,
velocHet=velocHet)
return(manage)
}
#' Setup Population Array
#'
#' Easy way to create initial population array that is not initialized to environment.
#' @param inds Initial population size in all species by location, by subpatch combo. Must be an integer!
#' @param S Number of species
#' @param L Length (number of patches)
#' @param W Width (number of subpatches)
#'
#' @return \code{S}x\code{L}x\code{W} array of population sizes over species, location, and subpatch
#' @export
#'
#' @examples
#' n <- nSetup(4,4,32,2)
nSetup <- function(inds=20,S,L,W){
n<-array(inds,c(S,L,W))
return(n)
}
gabackus/rcomms documentation built on Nov. 4, 2019, 1 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions commSetup manageSetup nSetup
Documented in commSetup manageSetup nSetup
#' Setup Model Structure
#'
#' This sets up the basic structure of the model.
#' It creates a list of biological parameters in \code{P} (\code{S} randomized species and their parameters) and it creates a list of environmental parameters in \code{X}.
#' All parameters have default values, so input is not necessary.
#'
#' @param S Total number of species created with hetSetup.
#' @param L Total number of patches in the metacommunity.
#' @param W Number of microhabitats in each patch.
#' @param zo A vector of pre-defined optimal temperature values. Only works if length(zo) is S.
#' @param gam A vector of pre-defined mean dispersal distances. Only works if length(gam) is S. If no vector is specified, gamMean and gamSD are used to randomly generate the vector gam.
#' @param sig A vector of pre-defined thermal tolerance breadths. Only works if length(sig) is S. If no vector is specified, sigMean and sigSD are used to randomly generate the vector sig.
#' @param A Matrix of the relative competition coefficients between species.
#' @param m A vector of pre-defined mortality probabilities (or a single value that will be shared for all species). These are probabilies, so the value must be between 0 and 1 (inclusive). Only works if length(m) is 1 or S.
#' @param gamMean Mean dispersal distance for randomized species. Default is based on Urban et al. 2012.
#' @param gamSD Standard deviation of dispersal distance for randomized species. Default is based on Urban et al. 2012.
#' @param sigMean Mean thermal tolerance breadth for randomized species. Default is based on Urban et al. 2012.
#' @param sigSD Standard deviation of thermal tolerance breadth for randomized species. Default is based on Urban et al. 2012.
#' @param lam Skewness in thermal tolerance. Default is based on Urban et al. 2012. (to have a mild decrease moving toward colder temperatures and a sharp decrease moving toward warmer temperatures).
#' @param B Area of integrated birth rate over all T for each species.
#' @param compType The type of competition in a string. Must be either "lottery" or "temp"
#' @param XW Window of analysis (to remove edge effects)
#' @param temp2d A matrix of pre-defined temperatures over x. Only works if nrow(temp2d) is L and ncol(temp2d) is W.
#' @param tempLow Lowest mean temperature on linear temperature gradient. temp1d(L)=tempLow
#' @param tempHigh Highest mean temperature on linear temperature gradient. temp1d(1)=tempHigh
#' @param tempRev If tempRev=T, then temp1d(1)=tempLow and temp1d(L)=tempHigh.
#' @param tempGH Hurst exponent for global temperature heterogeneity. H=0.5 is Brownian motion; 0.5<H<=1 is long-term positive autocorrelation and 0<=H<0.5 is long-term negative autocorrelation.
#' @param tempGSD Controls the magnitude of the global temperature heterogeneity.
#' @param tempLSD Standard deviation in temperature between microhabitats in each patch.
#' @param years Maximum number of years for initialization + climate change.
#' @param y Current year.
#' @param tempY A vector of pre-defined temperatures over time. Only works if length(tempY) is years.
#' @param tau Average temperature change per year.
#' @param tempYAC Temperature autocorrelation over time. Default is global temperature AC from 1880-1979.
#' @param tempYSD Temperature standard deviation over time. Default is global temperature SD from 1880-1979.
#' @param Tau Average temperature change per year (in a matrix). If
#' @param tempYXGH Hurst exponent for global heterogeneity in temperature change over time.
#' @param tempYXGSD Magnitude of global heterogeneity in temperature change over time.
#' @param tempYXLSD Standard deviation in temperature change between microhabitats in each patch.
#' @param tempYXGLH Hurst exponent for global heterogeneity in temperature change over time and space.
#' @param tempYXGLSD Standard deviation in temperature change between microhabitats in each patch over space.
#' @param Q A matrix of pre-defined habitat quality over x. Only works if nrow(Q) is L and ncol(Q) is W.
#' @param QMean Average habitat quality for any microhabitat.
#' @param QGH Hurst exponent for global habitat quality heterogeneity. H=0.5 is Brownian motion; 0.5<H<=1 is long-term positive autocorrelation and 0<=H<0.5 is long-term negative autocorrelation.
#' @param QGSD Controls the magnitude of the global habitat quality heterogeneity.
#' @param QLType How is the habitat quality lost? Must be "none", "subtract", or "multiplicative".
#' @param QLSD Standard deviation in habitat quality between microhabitats in each patch.
#'
#' @return A list that contains P (a list of biotic paraters) and X (a list of abiotic parameters)
#' @export
#'
#' @examples
#' cSetup<-commSetup(S=4,L=512,W=4,years=202)
#' P<-cSetup$P
#' X<-cSetup$X
#' n0 <- nSetup(8,4,512,4)
#' initSim <- commSimulate(n0,P,X,y=1,years=100,init=T)
#' n1 <- initSim$n
#' ccSim <- commSimulate(n1,P,X,y=101,years=100,init=F)
#' vComTime(ccSim$N)
#'
commSetup <- function(S=32, L=512, W=8,
zo=NULL, gam=NULL, sig=NULL, A=NULL, m=1,
gamMean=2.5, gamSD=2.5,
sigMean=5, sigSD=5,
lam=-2.7, B=10,
compType="lottery",
XW=seq(129,384),
temp2d=NULL, tempLow=9.78, tempHigh=30.22, tempRev=F, tempGH=1, tempGSD=0, tempLSD=1, tempGLH=0.81,tempGLSD=2, tempLSDRed=0,
years=1000, y=1,
tempY=NULL, tau=0.04, tempYAC=0.767, tempYSD=0.1639,
Tau=NULL, tempYXGH=1, tempYXGSD=0, tempYXLSD=0, tempYXGLH=1, tempYXGLSD=0,
Q=NULL, QMean=8, QGH=1, QGSD=0, QLSD=0, QLType="none", QLoss=0){
##########################################################
# First, specify biological parameter values for all of the S species.
# Optimal temperature for each species. These can be randomly picked from a uniform distribution or pre-defined.
if(is.null(zo)){
zo <- runif(S,9.9,30.1)
} else {
if(length(zo)!=S){
stop("zo does not match the number of species!")
}
}
# Dispersal distance for each species. These can be randomly picked from a lognormal distribution or pre-defined.
if(is.null(gam)){
# To use the lognormal distribution, we need to convert mean and SD values
gamMu <- log(gamMean/sqrt(1+gamSD^2/gamMean^2))
gamSig <- sqrt(log(1+gamSD^2/gamMean^2))
gam <- rlnorm(S,gamMu,gamSig)
} else {
if(length(gam)!=S){
stop("gam does not match the number of species!")
}
}
# Thermal tolerance breadth for each species. These can be randomly picked from a lognormal distribution or pre-defined.
if(is.null(sig)){
# To use the lognormal distribution, we need to convert mean and SD values
sigMu <- log(sigMean/sqrt(1+sigSD^2/sigMean^2))
sigSig <- sqrt(log(1+sigSD^2/sigMean^2))
sig <- rlnorm(S,sigMu,sigSig)
} else {
if(length(sig)!=S){
stop("sig does not match the number of species!")
}
}
# Competition coefficients between each pair of species species. By default, all coefficients are 1 (lottery competition), but A can be pre-defined.
if(is.null(A)){
A <- matrix(1,S,S)
}else{
if(!(nrow(A)==S & ncol(A)==S)){
stop("A does not match the number of species!")
}
}
# Yearly mortality probability for each species.
if(any(m<0) | any(m>1)){
stop("m should be between 0 and 1 (inclusive)")
}
if(length(m)==1){
# If only one value is provided, all species will have the same mortality
m <- rep(m,S)
} else if(length(m)!=S){
stop("m does not match the number of species!")
}
# Mortality is more convenient in this form
M <- rep(m,W*L)
# Make sure that compType is used correctly
if(!(compType=="lottery" | compType=="temp")){
stop("compType must be either 'lottery' or 'temp'!")
}
##########################################################
# Using the randomized species parameters, we derive other variables needed for computation.
# zo helps define where the species' thermal optimum is, but mathematically this is not completely correct.
# If we let zo be the true optimum, z is the value we plug into the reproduction function so that argmax(b)==zo
# We apply the function zAdjust to all values of zo to calculate z
z <- mapply(zAdjust, sig, zo, lam, 2^13)
# To speed up computation time, we define full dispersal kernels now.
# The dispersal kernel uses q, a transformation of gam
q <- sapply(1:S, function(i) 1+1/gam[i]-sqrt(1+1/gam[i]^2))
l <- (-L):(L)
k <- t(sapply(1:S, function(i) doubGeom(l,q[i])))
k[k<10^(-15)] <- 0.0
# K is a list of separate L by 2L dispersal kernel matrices
# K[[s]] is the dispersal kernel matrix of species s
# K[[s]][i,] is a vector of probabilities for a propagule in patch i to spread to patch j-L/2 (this is extra long to account for a propagule spreading beyond the limits of the ecosystem)
K <- rep(list(matrix(0,L,2*L)),S)
for(i in 1:S){
Ki<-matrix(0,2*L,4*L)
for(j in 1:(2*L)){
Ki[j,j:(j+2*L)]<-k[i,]
Ki[j,3/2*L]<-sum(Ki[j,1:(3/2*L)])
Ki[j,5/2*L+1]<-sum(Ki[j,(5/2*L+1):(4*L)])
}
K[[i]]<-Ki[(L/2+1):(3*L/2),(3/2*L):(5/2*L)]
}
# Tolerance and reproductive strength have a tradeoff
# Birth rate is adjusted so each species has roughly equal birth rate when integrated over all T
# ro is a constant that adjusts to this reproductive output
ro <- sapply(sig,function(x) rAdjust(x,B,lam,1e-06,2^13))
##########################################################
# Put the biological parameters together into a single list, P
P=list(S=S,
z=z,
gam=gam,
sig=sig,
lam=lam,
A=A,
M=M,
ro=ro,
zo=zo,
K=K,
compType=compType)
##########################################################
# Next, we define environmental parameters.
# Discrete spatial domain: from 1 to L (integers)
x <- seq(1,L)
# temp2d is the current temperature over all x and microhabitats
if(is.null(temp2d)){
temp1dr <- seq(tempHigh,tempLow,length=L)
if(tempRev){
temp1dr<-rev(temp1dr)
}
tempG <- tempVarH(L,tempGH)
temp1d <- temp1dr+tempG*tempGSD/sd(tempG)
if(W>1){
tempLSDx <- tempVarH(L,tempGLH)
tempLSDx2 <- tempLSD*( (1-tempLSDRed) + tempLSDRed/(1+exp(-tempLSDx*tempGLSD/sd(tempLSDx))))
temp2dr <- matrix(temp1d,L,W)
tempLH <- matrix(rnorm(L*W,0,1),L,W)
tempLH2 <- t(sapply(1:L, function(i) mean(tempLH[i,])+(tempLH[i,]-mean(tempLH[i,]))*tempLSDx2[i]/sd(tempLH[i,])))
tempLH3 <- t(matrix(sapply(1:L, function(x) sort(tempLH2[x,]-mean(tempLH2[x,]))),W,L))
temp2d <- temp2dr+tempLH3
} else{
temp2d<-temp1d
}
} else {
if(!(nrow(temp2d)==L & ncol(tempsd)==W)){
stop("temp2d does not match environment size!")
}
temp1d<-rowMeans(temp2d)
}
# tempY is a vector of the temperature over time
if(is.null(tempY)){
tempY<-0:years
for (i in 1:years){
# A new epsi is calculated for each time step
tempY[i+1] <- tempYAC*tempY[i]+rnorm(1,0,tempYSD)*sqrt(1-tempYAC^2)
}
} else{
if(length(tempY)!=years+1){
stop("tempY does not match years!")
}
}
# Tau is the temperature change over time in each patch and subpatch
if(is.null(Tau)){
tauLW <- matrix(tau*100,L,W)
tauG <- tempVarH(L,tempYXGH)
tauLW2 <- tauLW + tauG*tempYXGSD/sd(tauG)
if(W>1){
tauLH <- matrix(rnorm(L*W,0,tempYXLSD),L,W)
tauLH <- t(matrix(sapply(1:L, function(x) tauLH[x,]-mean(tauLH[x,])),W,L))
} else{
tauLH <- 0
}
Tau <- tauLW2+tauLH
Tau <- Tau/100
} else{
if(!(nrow(Tau)==L & ncol(Tau)==W)){
stop("Tau does not match environment size!")
}
}
# Habitat quality could differ in space or with species, but we will keep it constant for now
if(is.null(Q)){
if(QLType=="none"){
Q<-matrix(QMean,L,W)
} else{
QG <- tempVarH(L,QGH)
QLH <- matrix(rnorm(L*W,0,QLSD),L,W)
if(QLType=="subtract"){
Q1<-QMean-QMean/(QLoss^-1*(1+exp(-QLSD*QG/sd(QG))))
} else if(QLType=="multiplicative"){
if(QLoss<=0){
QLoss<-0.01
}
Q1<-(QMean*QLoss^-1)/(1+exp(-QG*QGSD/sd(QG)))
}
Q <- array(Q1,c(L,W))+QLH
}
} else {
if(!(nrow(Q)==L & ncol(A)==W)){
stop("Q does not match environment size!")
}
}
##########################################################
# Put the abiotic parameters together into a single list, X
X=list(L=L,
x=x,
XW=XW,
temp1d=temp1d,
tau=tau,
Tau=Tau,
Q=Q,
tempY=tempY,
tempYAC=tempYAC,
tempYSD=tempYSD,
temp2d=temp2d,
tempLSDx=tempLSDx,
W=W)
##########################################################
# Export it all as a list
return(list(P=P,X=X))
}
#' Setup Management Structure
#'
#' This sets up the management options for a simulation.
#' It creates a list of several lists that include paramaters for each management option.
#' Many of these don't bother to check for errors, so make sure you're careful.
#' @param P List of biotic parameters
#' @param X List of abiotic parameters
#' @param years The number of years that the model is going to be run.
#' @param AM Assisted migration. TRUE or FALSE determine whether the model runs through AM during each time step. If left NULL, this function will not define any parameters for AM.
#' @param corridor Corridor establishment. TRUE or FALSE determine whether the model runs through corridor management during each time step. If left NULL, this function will not define any parameters for corridor.
#' @param habQual Change habitat quality. TRUE or FALSE determine whether the model runs through changing habitat quality during each time step. If left NULL, this function will not define any parameters for habQual.
#' @param locHet Change local heterogeneity. TRUE or FALSE determine whether the model runs through changing local heterogeneity during each time step. If left NULL, this function will not define any parameters for locHet.
#' @param velocHet Change heterogeneity climate velocity. TRUE or FALSE determine whether the model runs through changing heterogeneity in climate velocity during each time step. If left NULL, this function will not define any parameters for velocHet.
#' @param targs Which species are specifically targeted for assisted migration? A list of 0s (for not targeted) and 1s (for targeted). Must be a vector of length S. Could also be a vector of species number (index from P parameters). Lastly, can just be "all" to include all species.
#' @param eta Relocation threshold. Any population of a target species that falls below this is relocated. Should be a nonnegative single integer or a vector of length S.
#' @param tCD Cooldown time between relocations. AM is not repeated if tCD has not passed since the last AM event. Should be a single nonnegative integer or a vector of length S.
#' @param rho Proportion of the total population moved when AM is triggered. Should be a single number from 0 to 1 or a vector of length S.
#' @param mu Probability of surviving AM. Should be a single number from 0 to 1 or a vector of length S.
#' @param zEst Estimate of species' thermal optimum. Should be a vector of length S.
#' @param xLoc How far ahead (in the direction of the leading edge) of the location closest to the thermal optimum estimate are we centering the relocation? Should be a single integer or a vector of length S.
#' @param recRad How wide is the recipient location? The recipient location consists of the center location (xLoc in front of the thermal optimum) and each patch +/- recRad from this center. Should be a single positive integer or a vector of length S.
#' @param donor Which individuals do we take from the donor community? 1: randomly (from the "top"), 2: from the trailing edge, 3: from the leading edge. Possible future extensions are 4: from the middle and 5: randomly from the old donor population. Should be a single integer.
#' @param recipient What shape should the recipient population look like? 1: an evenly distributed box. Possible future extensions are 2: rounded, 3: triangular, 4: thermal tolerance shaped. Should be a single integer.
#' @param tLR When was the last time the species was relocated? Should be an integer vector of length S.
#' @param randPick Do you want to pick exactly rho of the population (rounded down) or would you rather randomly pick from a binomial distribution with rho as the probability? If you want to do this randomly, set randPick to T. By default, it is set to F.
#'
#' @return A list of parameters that tell the model when and how to enact various managment techniques.
#' @export
#' @examples
#' cSetup<-commSetup(S=32,L=512,W=4,years=202)
#' P<-cSetup$P
#' X<-cSetup$X
#' mSetup<-manageSetup(P,X,AM=T,targs = "all",eta=40,rho=0.6)
#' n0 <- nSetup(8,32,512,4)
#' initSim <- commSimulate(n0,P,X,y=1,years=100,init=T)
#' n1 <- initSim$n
#' ccSim <- commSimulate(n1,P,X,y=101,years=100,init=F)
#' vComTime(log(ccSim$N[rowSums(flatComm(n1)[,X$XW])>0,]))
#' lines(c(0,100),c(log(mSetup$AM$eta[1]),log(mSetup$AM$eta[1])))
manageSetup <- function(P,X,years=100,
AM=NULL,corridor=NULL,habQual=NULL,locHet=NULL,velocHet=NULL,
targs=NULL,eta=50,tCD=5,rho=0.8,mu=0.8,zEst=P$zo,xLoc=10,recRad=2,donor=1,recipient=1,tLR=NULL,randPick=F){
# Now we set all of this up
S <- P$S
if(is.null(AM)){
# If you have no plans for AM at all, there is no need to define anything but AM$AM.
AM<-list(AM=F)
} else{
if(targs=="all"){targs <- 1:S} else{
if(is.null(targs)){targs<-rep(0,S)}
if(length(targs)>S){stop("targ too long")} else if(length(targs)<1){stop("targ too small")}
if(all(targs %in% c(0,1))){if(!(length(targs)==S | length(targs)==1)){stop("If using a 0/1 vector, targs needs to be length S")} else if(length(targs)==1){tInds<-targs; targs<-rep(0,S); targs[tInds]<-1}}else
if(all(targs %in% 1:S)){tInds<-targs; targs<-rep(0,S); targs[tInds]<-1} else{stop("If using targs as index vector, at least one of your indices were not valid.")}
}
if(!(length(eta)==S | length(eta)==1)){stop("Wrong length")} else if(length(eta)==1){eta<-rep(eta,S)}
if(!(length(tCD)==S | length(tCD)==1)){stop("Wrong length")} else if(length(tCD)==1){tCD<-rep(tCD,S)}
if(!(length(rho)==S | length(rho)==1)){stop("Wrong length")} else if(length(rho)==1){rho<-rep(rho,S)}
if(!(length(mu)==S | length(mu)==1)){stop("Wrong length")} else if(length(mu)==1){mu<-rep(mu,S)}
if(!(length(xLoc)==S | length(xLoc)==1)){stop("Wrong length")} else if(length(xLoc)==1){xLoc<-rep(xLoc,S)}
if(!(length(recRad)==S | length(recRad)==1)){stop("Wrong length")} else if(length(recRad)==1){recRad<-rep(recRad,S)}
relTimes<-matrix(0,S,years)
if(is.null(tLR)){tLR<-rep(0,S)}
AM <- list(AM=AM,targs=targs,eta=eta,tCD=tCD,rho=rho,mu=mu,zEst=zEst,xLoc=xLoc,recRad=recRad,donor=donor,recipient=recipient,tLR=tLR,randPick=randPick,relTimes=relTimes)
}
# CORRIDOR PARAMETERS
if(is.null(corridor)){
# If you have no plans for corridor establishment at all, there is no need to define anything but corridor$corridor.
corridor<-list(corridor=F)
}
# HABITAT QUALITY PARAMETERS
if(is.null(habQual)){
# If you have no plans for increasing habitat quality at all, there is no need to define anything but habQual$habQual.
habQual<-list(habQual=F)
}
# LOCAL HETEROGENEITY PARAMETERS
if(is.null(locHet)){
# If you have no plans for local heterogeneity change at all, there is no need to define anything but locHet$locHet.
locHet<-list(locHet=F)
}
# CLIMATE VELOCITY HETEROGENEITY PARAMETERS
if(is.null(velocHet)){
# If you have no plans for changing climate velocity heterogeneity at all, there is no need to define anything but velocHet$velocHet.
velocHet<-list(velocHet=F)
}
manage <- list(AM=AM,
corridor=corridor,
locHet=locHet,
velocHet=velocHet)
return(manage)
}
#' Setup Population Array
#'
#' Easy way to create initial population array that is not initialized to environment.
#' @param inds Initial population size in all species by location, by subpatch combo. Must be an integer!
#' @param S Number of species
#' @param L Length (number of patches)
#' @param W Width (number of subpatches)
#'
#' @return \code{S}x\code{L}x\code{W} array of population sizes over species, location, and subpatch
#' @export
#'
#' @examples
#' n <- nSetup(4,4,32,2)
nSetup <- function(inds=20,S,L,W){
n<-array(inds,c(S,L,W))
return(n)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.