#' DispFunMP R code for dispersal function #' #' # parameters: #' @param N <-para$N Integer = no. of Individuals #' @param G <-para$G vector of Integers from 1: x = number of generations #' @param SLbD <-para$SLbD vector of Integers from 1: x = Season Length = number of days before dispersal #' @param SLaD <-para$SLaD vector of Integers from 1: x = Season Length = number of days after dispersal #' @param PL <-para$PL Predation Lethality between 0 and 1 #' @param PP <-para$PP Predation Probability between 0 and 1; #' @param RLF <-para$RLF # RLF <- 1 # Resource Learning Factor (the higher its value the faster indies learn or more precisely the less items they need to collect to reduce the HT; if==1 the pop$L[k] is not modified; if >1 they learn extra fast; < 1 the get slower ) #' @param PLF <-para$PLF # PLF <- 0.5 # Predation Learning Factor the higher the factor the more L helps to reduce PL #' @param CF <-para$CF # CF <- 4 # if the average E of the population is E=0.5 (as in the beginning) the expected amount of Rmax = Rmax*0.5 or Rmax/2. CF==2 means no competition in a population with meanpop$E=0.5 #' @param a <- specifies the cost of learning - the higher the value the lower the costs. standard value = 1.4 #' @param patchID <- # number of patches e.g.1:20 #' @param R <-para$R vector Resource Types c("R0"= 0,"R1"=1,"R2"=2,"R3"=3,"R4"=4,"R5"=5,"R6"=6,"R7"=7,"R8"=8,"R9"=9,"R10"=10), #' @param Rval <-para$Rval vector of Resource values c("R0"= 0,"R1"=1,"R2"=5,"R3"=10,"R4"=15,"R5"=20,"R6"=25,"R7"=30,"R8"=35,"R9"=40,"R10"=45) #' @param RHT <-para$RHT vector of Resource HandlingTimes RHT= c("R0"= 0,"R1"=1,"R2"=5,"R3"=10,"R4"=15,"R5"=20,"R6"=25,"R7"=30,"R8"=35,"R9"=40,"R10"=45) #' @param Env <- list of different Environments from which the metapop will be created #' @param OF <-para$OF vector of ResourceType OverlookingFactor c("R0"= 0.0,"R1"=0.0,"R2"=0.9,"R3"=0.9,"R4"=0.9,"R5"=0.9,"R6"=0.4,"R7"=0.3,"R8"=0.2,"R9"=0.1,"R10"=0.05) #' @param dC <-dispersal cost (0-1) the higher the more likely individuals die during dispersal, any values > 1 means certain death #' @param qL <-para$qL mutation probability for trait "L" (0-1); standard value = 0.1 #' @param qE <-para$qE mutation probability for trait "E" (0-1); standard value = 0.1 #' @param qD <-para$qD mutation probability for trait "dispersal" (0-1); standard value = 0.1 #' @param Ext <- para$Ext specifies whether extinction occurs or not ("y/n") #' @param EX <- EXtinction counter (integer) gives how often a extinction takes place: if(g%%EX==0) Thus number gives the number of generations when an event occurs #' @param EXn <- gives how many patches will get extinct when this generation extinction takes place. #' @param saveplot <-para$saveplot # used to indicate whether a plot should be saved #' @param printPlot <-para$printPlot # used to indicate whether a plot should be made #' @param writeDF <-para$writeDF indicate whether (1) or not (0) the dataframe should be saved #' @param Cyc <-para$Cyc # provide Index for write.csv in season and for .png in function= "plotIndi" # not implemented yet #' @param path_to <-para$path_to directory were the df files are saved #' #' @examples #' para <- list(N = 500, G = 1:150, SLbD=1:40, SLaD=1:40, PL = 0.5, PP = 0.1, #' RLF=1, PLF=1, CF=8, a=1.4, qL=0.1, qE=0.1, qD=0.1,Ext="y",EX=1,EXn=2, #' #' R= c("R0"= 0,"R1"=1,"R2"=2,"R3"=3,"R4"=4,"R5"=5,"R6"=6,"R7"=7,"R8"=8,"R9"=9,"R10"=10), #' Rval= c("R0"= 0,"R1"=1,"R2"=5,"R3"=10,"R4"=15,"R5"=20,"R6"=25,"R7"=30,"R8"=35,"R9"=40,"R10"=45), #' RHT= c("R0"= 0,"R1"=1,"R2"=5,"R3"=10,"R4"=15,"R5"=20,"R6"=25,"R7"=30,"R8"=35,"R9"=40,"R10"=45), #' OF= c("R0"= 0.0,"R1"=0.0,"R2"=0.9,"R3"=0.9,"R4"=0.9,"R5"=0.9,"R6"=0.4,"R7"=0.3,"R8"=0.2,"R9"=0.1,"R10"=0.05), #' Cyc=1, #' printPlot=0, #' writeDF=0, #' path_to= c("~user/Desktop/Dispersal/DispersalR/Results/") #' ) #' system.time(dispFunMP(para,Env)) #' #' @export dispFunMP <- function(para,Env){ print(Sys.time()) # unpack parameter from list "para" N<-para$N G<-para$G SLbD<-para$SLbD SLaD<-para$SLaD PL<-para$PL PP<-para$PP RLF<-para$RLF PLF<-para$PLF CF<-para$CF a<-para$a dC<-para$dC patchID <- para$patchID patchEnv <- para$patchEnv R<-para$R Rval<-para$Rval RHT<-para$RHT OF<-para$OF qL<-para$qL qE<-para$qE qD<-para$qD Ext<-para$Ext EX<-para$EX EXn<-para$EXn printPlot<-para$printPlot saveplot<- para$saveplot writeDF<-para$writeDF Cyc<-para$Cyc path_to <- para$path_to # list to store population information for all patches mpop <- list() # metapopulation list ### Create "dfLE" list with data.frames for individual results of L, E and D in generation[g] per patch #### dfLE <- list() for(p in patchID){ dfLE[[p]]<-as.data.frame(matrix(0, ncol = length(G)+1, nrow = N*3)) dfLE[[p]][,1]<-c(rep("L",N),rep("E",N),rep("D",N)) } names(dfLE)<-patchID # length indicators for putting results into dfLE len1 <- 1 len2 <- N len3 <- N+1 len4 <- N*2 len5 <- N*2+1 len6 <- N*3 #### start-population for each patch #### for(e in 1:length(patchID)){ mpop[[e]] <- data.frame(matrix(0, nrow = N*length(patchID)+N, ncol = 78)) # nrow = N*length(patchID)+N because during dispersal # in theory every individual of the metapop might settle in this patch and thus the number of rows should be that long +N so that individuals of all patches can be located # to a specific row. And thus indis from patch 1 will not be transported to row N*length(patchID) which is = N and thus is occupied by a resident already but to N*length(patchID)+N, which is empty mpop[[e]][1:N,1] = 1:N # Individuals ID mpop[[e]][1:N,2] = runif(N) # Individuals L mpop[[e]][1:N,3] = runif(N) # Individuals E mpop[[e]][1:N,4] = runif(N) # Individuals D mpop[[e]][1:N,5] = 1 # alive yes/no mpop[[e]][1:N,76] = 0 # in the beginning it is yet not determined whether an indi will disperse or not names(mpop[[e]])<- c( "ID", "L" , # Individual trait values for learning in the first generation comes from a normal distribution "E", # Individual trait values for exploration in the first generation comes from a normal distribution "D", # Dispersal; i.e. how likely individuals will disperse (0-1) "alive", # indicates whether a individual is alive or not "dead", # indicate the day on which an indi died "sumS", # represents the gained resourced values per individual over time (whole season) # Number of collected resources "RCol0","RCol1","RCol2","RCol3","RCol4","RCol5","RCol6","RCol7","RCol8","RCol9","RCol10" , # Number of collected resources (corrected for efficency) "RCol0Cor","RCol1Cor","RCol2Cor","RCol3Cor","RCol4Cor","RCol5Cor", "RCol6Cor","RCol7Cor","RCol8Cor","RCol9Cor","RCol10Cor", # Number of collected resources values (corrected for efficiency) "RCol0V","RCol1V","RCol2V","RCol3V","RCol4V","RCol5V", "RCol6V","RCol7V","RCol8V","RCol9V","RCol10V", # Number of collected resources after dispersal "RCol0aD","RCol1aD","RCol2aD","RCol3aD","RCol4aD","RCol5aD","RCol6aD","RCol7aD","RCol8aD","RCol9aD","RCol10aD" , # Number of collected resources (corrected for efficiency) "RCol0CoraD","RCol1CoraD","RCol2CoraD","RCol3CoraD","RCol4CoraD","RCol5CoraD", "RCol6CoraD","RCol7CoraD","RCol8CoraD","RCol9CoraD","RCol10CoraD", # Number of collected resources Values (corrected for efficiency) "RCol0VaD","RCol1VaD","RCol2VaD","RCol3VaD","RCol4VaD","RCol5VaD", "RCol6VaD","RCol7VaD","RCol8VaD","RCol9VaD","RCol10VaD", "orignPatch","Resident","Disperser","patchType", "Daysalive") } #### creating vectors for plots and co #### meanL <- list() # output list for mean population L meanE <- list() # output list for mean population E meanD <- list() # output list for mean population D in the beginning of season (i.e. before dispersal) RColTab <- list() # output list for number of collected resources per patch and RType mOut <- list() #output list for tab_Coeff library(geepack) # for running GEEs patchDensBefore<-list() #output list for number of individuals in patch[[e]] before dispersal patchDensAfter<-list() # output list for number of individuals in patch[[e]] after dispersal suc_residents<-list() # output list for number of successful residents in patch[[e]] for(p in patchID){ meanL[[p]]<-vector("double",length = length(G)) meanE[[p]]<-vector("double",length = length(G)) meanD[[p]]<-vector("double",length = length(G)) patchDensBefore[[p]]<-vector("double",length = length(G)) patchDensAfter[[p]]<-vector("double",length = length(G)) suc_residents[[p]]<-vector("double",length = length(G)) RColTab[[p]] <- data.frame(matrix(0, nrow = 5, ncol = 11)) # for control rownames(RColTab[[p]]) <- c("Cor_OF_E_dP", "Cor_Eff","Cor_Competi","CCFbD","CCFaD") colnames(RColTab[[p]]) <- c("R0","R1","R2","R3","R4","R5","R6","R7","R8","R9","R10") } names(meanL)<-names(meanE)<-names(meanD)<-names(patchDensAfter)<-names(patchDensBefore)<-names(suc_residents)<-patchID #### Loop Generation #### for(g in 1:length(G)) { # for each generation go through all patches and all individuals (rows) for(e in 1:length(patchID)) { # go through all patches in metapop("mpop") # mean L in each generation meanL[[e]][g] <- sum(mpop[[e]][,2])/ sum(mpop[[e]][,1] > 0) # as there are many "placeholder" rows with zeros, thus I cannot take simply mean() # mean E in each generation meanE[[e]][g] <- sum(mpop[[e]][,3])/ sum(mpop[[e]][,1] > 0) # mean D in each generation meanD[[e]][g] <- sum(mpop[[e]][,4])/ sum(mpop[[e]][,1] > 0) ### put individual results for generation [g] into "dfLE" dfLE[[e]][len1:len2, g+1] <- mpop[[e]][1:N,2] dfLE[[e]][len3:len4, g+1] <- mpop[[e]][1:N,3] dfLE[[e]][len5:len6, g+1] <- mpop[[e]][1:N,4] ### erase RcolTable form previous generation # RcolTable is a summarize Table for control, which counts how many resources have been collected # and tables the different number of itmes corrected for efficiency, competition etc. RColTab[[e]][,] <- 0 #### start loop for each Indi in patch mpop[[e]] before dispersal #### for(k in 1:nrow(mpop[[e]])){ # start loop for each Indi in patch mpop[[e]] if(mpop[[e]][k,1] > 0){ # only if there is a ID for an individual in this row, the individual exists and thus can collect resources mpop[[e]][k,74] <- e # specify patch origin mpop[[e]][k,77] <- patchEnv[e] # specify patchType of patch origin ##### RCollect ##### # RCollect gives the number of Items per RType an individual comes across based on the abundancies of RTypes givin in Envi[e] # and the individuals exploration tendency - the more it explores the more sites it visits and thus the more resources it comes across # however, if it finds the resources and whether it can exploit them will be calculated below. RCollect <- 5*max(SLbD) * Env[[patchEnv[e]]] + 5*max(SLbD)*mpop[[e]][k,3] * Env[[patchEnv[e]]] # basal amount of resources calculated by: # 5 per day before dispersal (SLbD) multiplied with the abundance factor of each RType specified as in Env[patchEnv[e]] multiplied with the individuals exploration tendency (pop$E[k]) # (The multipilcation factor "5" is not included in the publication due to readability.The same amount of resoruces could also be obtained by simply increasing the values in the "Envi" vector.) # Amount of resources individuals gain depending on their exploration tendency. Thus, half of the total amount of resources every individual gains (it explores at least every second timestep). # An individual gains the more of other half of the total amount the more it explores. Thus, with an E of 0 you gain max 50% with an E of=1 you gain max 100% depending on the OF of resources. # The assumption that an individual will move every second time step is in line with findings from our previous models in which E was around 40-50%. #### calculating death before dispersal #### Pevent<-rbinom(length(SLbD), size = 1, prob = PP*mpop[[e]][k,3]) # gives for each day of SL if a predator encounter happens; the higher E the higher the probability becomes cSumP<-cumsum(Pevent) # gives cumulative predator encounters lethal<-(1/(1+(mpop[[e]][k,2]*PLF)*(cSumP-1))) * PL # gives the probabilities of lethality per encounter corrected for anti-predation learning lethality<-lethal*Pevent death<-runif(length(SLbD)) # random death-threshold dead<-lethality>death # specifies if lethality is higher than threshold if(is.finite(match(1,dead)/length(SLbD))){ mpop[[e]][k,5] <- 0 mpop[[e]][k,6] <- match(1,dead)/length(SLbD) # dp= deathPenalty:dead: gives the index of the first time that lethality > "random" death. # This divided by SL gives the length of life as a fraction of 1. Thus this parameter can be used to correct the Number of collected resources mpop[[e]][k,78] <- match(1,dead) } else { mpop[[e]][k,6] <- 1 mpop[[e]][k,78] <- length(SLbD) } #### correcting resources #### # corrected for individuals E and OF per RType and by dP (time of death) for(i in 1:length(R)){ mpop[[e]][k,7+i] <- RCollect[i]*(1-mpop[[e]][k,3]*OF[i]) * mpop[[e]][k,6] # sum of collected items per RType corrected for individuals E and OF per RType and by dP (time of death) RColTab[[e]][1,i] <- RCollect[i]*(1-mpop[[e]][k,3]*OF[i]) * mpop[[e]][k,6] + RColTab[[e]][1,i] # for control only } ## correcting for efficiency depending on sum of collected itmes per RType, HT, and L[k] nR <- vector("double", length = length(R)) for(i in 1:length(R)){ # start loop correcting values and account for individual handling efficiency nR[i]<- mpop[[e]][k,7+i] if(nR[i] > 0 && R[i] > 1){ # (R[i] > 1) for "R0" and "R1" no efficiency correction needed, since R1 has a HT=1 if(mpop[[e]][k,2] > 0){ # only if individuals have an L > 0 it will be calculated if they gain something from RTypes>R1 eff_rf <- 1-sqrt(RHT[i]-1)/(sqrt(1:nR[i])*(mpop[[e]][k,2]*RLF)) # eff_rf = efficiency_resource_factor eff_rf <- ifelse(eff_rf >= 0 , eff_rf, 0) mpop[[e]][k,18+i]<-sum(eff_rf) # gives number of R itmes exploited corrected for efficiency mpop[[e]][k,29+i]<-sum(eff_rf*Rval[i]) # gives sum of value for specific RType exploited corrected for efficiency RColTab[[e]][2,i] <- sum(eff_rf) + RColTab[[e]][2,i] # for control only } else { # if indi L==0, O will be assignd # This step should be unnecessary as there should be zeros assigned already mpop[[e]][k,18+i] <- 0 # } # end of correcting Resoruce sum(resource values) for all > 0 count and not "R0" or "R1" } else { # start "correction" for "R1" = fill in Column "pop$RCol1V" and "pop$RCol1Cor" if(R[i]==1){ # take out this "if" then R0 would be also counted mpop[[e]][k,18+i]<-nR[i] mpop[[e]][k,29+i]<-mpop[[e]][k,18+i]*Rval[i] RColTab[[e]][2,i] <- nR[i] + RColTab[[e]][2,i] # for control only }}# end of counting "R1" } } }# end of collecting R itmes for each individual #### Competition #### # in the environment of this population there should be around that many resources in total per RType Rmax <- 10 * max(SLbD) * N * Env[[patchEnv[e]]] # 10 per day = 5 before dispersal and 5 per day after dispersal) fields visited per day and individual (if all individuals have E=1) times the abundancy of items per RTypes. Thus, this is the maximum number of resource items per Rytpe in the population # I use this number as the max amount of resources to be found per type in this environment. # To induce competition I divide this number by a factor (Competition Factor = CF) and take this number as the maximum of items per RType available RmaxRtype <- Rmax / CF # if then the total amount of collected items per RType in the population is higher, the individuals count of items will be corrected. # This correction is done in relation of RmaxRtype to how many items have been collected by a specific individual. e.g. by: ### Competition Correction with the proxi number of itmes an indi exploited CCF <- ifelse(colSums(mpop[[e]][,19:29]) <= 0, 0, RmaxRtype/colSums(mpop[[e]][,19:29])) # for all RTypes CCF=CompetitionCorrectioFactor = vector with CCF for each RType RColTab[[e]][4,] <- CCF # for control only # pop[k,19:29] # If CCF < 1 it means that more resource items have been exploited than the maximum number of items avaiable in this environment. # Thus, for all individuals the number/values of collected resources for RTypes in which the maximum has been exceeded need to be corrected # by the CCF. To do so multiply the RCol4V ect. vector with the CCF ### now for all RTypes for which CCF is < 1 (since when less items has been exploited than available = CCF > 1 and does not need to be corrected) for(i in 1:length(R)){ for(k in 1:nrow(mpop[[e]])){ if(mpop[[e]][k,1]>0 ){ # only if there is a ID for an individual in this row, the individual exists and thus can collect resources if(CCF[i]<1){ mpop[[e]][k,18+i] <- mpop[[e]][k,18+i] * CCF[i] # for a single RType item count mpop[[e]][k,29+i] <- mpop[[e]][k,29+i] * CCF[i] # for a single RType Value RColTab[[e]][3,i] <- mpop[[e]][k,18+i] * CCF[i] + RColTab[[e]][3,i] # for control only } else{ # if no correction is needed then in RColTab[[e]][3,i] will be the same as in RColTab[[e]][2,i] RColTab[[e]][3,i] <- mpop[[e]][k,18+i] + RColTab[[e]][3,i] # for control only } }} # end of k individual row } # end of i Resource Type in R and end of loop #### pop$Sum #### # getting the overall ResourceValue intake per individual for(k in 1:nrow(mpop[[e]])){ if(mpop[[e]][k,1]>0){ # only if there is a ID for an individual in this row, the individual exist and thus can collect resources mpop[[e]][k,7]<-sum(mpop[[e]][k,30:40]) }} patchDensBefore[[e]][g]<-sum(mpop[[e]][,1] > 0) # checking patchsize before dispersal; for control only #### end of harvest before dispersal in generation [g] }# end loop of for(e in 1:length(patchID)){ # go through all patches in metapop("mpop") #### Dispersal #### for(e in 1:length(patchID)){ # go through all patches in metapop("mpop") for(k in 1:N){ # only the N residents are able to disperse. Individuals in rows > N are individuals which just immigrated from another patch if(mpop[[e]][k,1]>0 & mpop[[e]][k,5] == 1){ # only if there is a ID for an individual in this row, the individual exists and thus may disperse # and only if it is still alive dispTh<-runif(1) # random dispersal-threshold if(mpop[[e]][k,4] > dispTh){ # if dispersal tendency of indi is higher than dispTh then this indi will engage in dispersal mpop[[e]][k,76] <- "y" # specify that this indi dispersed #dispersal costs dispDeath<-runif(1) # random dispersal-death-threshold if(dC <= dispDeath){ # if dispersalCost (e.g.dC=0.1) are not higher than random threshold than indies will not die and disperse successfull dP <- sample(patchID,1) # random patch to disperse mpop[[dP]][N*e+k,] <- mpop[[e]][k,] # locate to new patch# N*e-N gives the first row in mpop[[dP]] which is reserved for immigrants of mpop[[e]] # + k indicates the row for individual[k] } # end of dispersal cost mpop[[e]][k,] <- 0 # delete from old patch; i.e. indi k is gone either to another patch or dead } # dispersal decision specifies if disoTh is higher than threshold mpop[[e]][k,76] <- "n" # specify that this indi did not dispersed }# end of dispersal loop for indi [k] }# end of dispersal loop for mpop[[e]] }# end loop of for(e in 1:length(patchID)){ # go through all patches in metapop("mpop") #### harvest after dispersal #### #### start loop for each Indi in patch mpop[[e]] AFTER dispersal #### for(e in 1:length(patchID)){ # go through all patches in metapop("mpop") patchDensAfter[[e]][g]<-sum(mpop[[e]][,1] > 0) # checking patchsize after dispersal # assign resident status for each individual for(k in 1:nrow(mpop[[e]])){ mpop[[e]][k,75] <- ifelse(mpop[[e]][k,74] > 0 & mpop[[e]][k,74]!= e, "n", 0) if(mpop[[e]][k,74] > 0 & mpop[[e]][k,74] == e){mpop[[e]][k,75] <- "y"} } for(k in 1:nrow(mpop[[e]])){ # start loop for each Indi in patch mpop[[e]] if(mpop[[e]][k,1]>0 & mpop[[e]][k,5]==1){ # only if there is a ID for an individual in this row, the individual exists and thus can collect resources ##### RCollect ##### RCollect <- 5*max(SLaD) * Env[[patchEnv[e]]] + 5*max(SLaD)*mpop[[e]][k,3] * Env[[patchEnv[e]]] #### calculating death aD#### Pevent<-rbinom(length(SLaD), size = 1, prob = PP*mpop[[e]][k,3]) # gives for each day of SL if a predator encounter happens; the higher E the higher the prob. becomes cSumP<-cumsum(Pevent) # gives cumulative predator encounters lethal<-(1/(1+(mpop[[e]][k,2]*PLF)*(cSumP-1))) * PL # gives the probabilities of lethality per encounter corrected for anti-predation learning lethality<-lethal*Pevent death<-runif(length(SLaD)) # random death-threshold dead<-lethality>death # specifies if lethality is higher than threshold if(is.finite(match(1,dead)/length(SLaD))){ mpop[[e]][k,5] <- 0 mpop[[e]][k,6] <- match(1,dead)/length(SLaD) mpop[[e]][k,78] <- match(1,dead)+length(SLbD) # number of days it lived during SLaD + number of days SLbD } else { mpop[[e]][k,6] <- 1 mpop[[e]][k,78] <- length(SLbD)+length(SLaD) } #### correcting resources aD #### # corrected for individuals E and OF per RType and by dP (time of death) for(i in 1:length(R)){ mpop[[e]][k,40+i] <- RCollect[i]*(1-mpop[[e]][k,3]*OF[i]) * mpop[[e]][k,6] # sum of collected items per RType corrected for individuals E and OF per RType and by dP (time of death) RColTab[[e]][1,i] <- RCollect[i]*(1-mpop[[e]][k,3]*OF[i]) * mpop[[e]][k,6] + RColTab[[e]][1,i] # for control only } ## correcting for efficiency depending on sum of collected items per RType, HT, and L[k] nR <- vector("double", length = length(R)) for(i in 1:length(R)){ # start loop correcting values and account for individual handling efficiency nR[i]<- mpop[[e]][k,40+i] # if(nR[i] > 0 && R[i] != 0){ # (R[i] != 0) without "R0" if(nR[i] > 0 && R[i] > 1){ # (R[i] > 1) for "R0" and "R1" no efficiency correction needed, since R1 has a HT=1 if(mpop[[e]][k,2] > 0){ # only if individuals have an L > 0 it will be calculated if they gain something from RTypes>R1 eff_rf <- 1-sqrt(RHT[i]-1)/(sqrt(1:nR[i] + mpop[[e]][k,7+i])*(mpop[[e]][k,2]*RLF)) # eff_rf = efficiency_resource_factor# sqrt(1:nR[i]+mpop[[e]][k,7+i]) # the "+mpop[[e]][k,7+i]" is the invdiviuals experience with this resource type from before dispersal eff_rf <- ifelse(eff_rf >= 0 , eff_rf, 0) mpop[[e]][k,51+i]<-sum(eff_rf) # gives number of R items exploited corrected for efficiency mpop[[e]][k,62+i]<-sum(eff_rf*Rval[i]) # gives sum of value for specific RType exploited corrected for efficiency RColTab[[e]][2,i] <- sum(eff_rf) + RColTab[[e]][2,i] # for control only } else { # if indi L==0, O will be assigned # This step should be unnecessary as there should be zeros assignd already mpop[[e]][k,51+i] <- 0 # } # end of correcting resource sum(resource values) for all > 0 count and not "R0" or "R1" } else { # start "correction" for "R1" = fill in column "pop$RCol1V" and "pop$RCol1Cor" if(R[i]==1){ mpop[[e]][k,51+i]<-nR[i] mpop[[e]][k,62+i]<-mpop[[e]][k,51+i]*Rval[i] RColTab[[e]][2,i] <- nR[i] + RColTab[[e]][2,i] # for control only }}# end of counting "R1" } } }# end of collecting R itmes for each individual #### Competition aD #### # in the environment of this population there should be around that many resources in total per RType Rmax <- 10 * max(SLaD) * N * Env[[patchEnv[e]]] RmaxRtype <-Rmax / CF ### Competition correction with the proxi number of items an indi exploited CCF <- ifelse(colSums(mpop[[e]][,52:62]) <= 0, 0, RmaxRtype/colSums(mpop[[e]][,52:62])) # for all RTypes CCF=CompetitionCorrectioFactor = vector with CCF for each RType RColTab[[e]][5,] <- CCF # for control only ### now for all RTypes for which CCF is < 1 (since when less itmes has been exploited than avaiable = CCF >1 and does not need to be corrected) for(i in 1:length(R)){ for(k in 1:nrow(mpop[[e]])){ if(mpop[[e]][k,1]>0){ # only if there is a ID for an individual in this row, the individual exist and thus can collect resources if(CCF[i]<1){ mpop[[e]][k,51+i] <- mpop[[e]][k,51+i] * CCF[i] # for a single RType item count mpop[[e]][k,62+i] <- mpop[[e]][k,62+i] * CCF[i] # for a single RType Value RColTab[[e]][3,i] <- mpop[[e]][k,18+i] * CCF[i] + RColTab[[e]][3,i] # for control only } else{ # if no correction is needed than in RColTab[[e]][3,i] will be the same as in RColTab[[e]][2,i] RColTab[[e]][3,i] <- mpop[[e]][k,18+i] + RColTab[[e]][3,i] # for control only } }} # end of k individual row } # end of i Resource Type in R and end of loop #### pop$Sum #### # getting the overall ResourceValue intake per individual for(k in 1:nrow(mpop[[e]])){ if(mpop[[e]][k,1]>0){ # only if there is a ID for an individual in this row, the individum exist and thus can collect resources mpop[[e]][k,7]<-mpop[[e]][k,7]+sum(mpop[[e]][k,63:73]) }} ## counting successful residents: could patch survive without emigration? suc_residents[[e]][g]<-sum(mpop[[e]][,7] > 0 & mpop[[e]][,74] == e) #### running GEEs #### if(g == 1 | g%%10==0){ # create Metapop of this generation Metapop1<-data.frame(mpop[[1]]) Metapop<-subset(Metapop1, Metapop1[,75]!="0" ) for(l in 2:length(patchID)){ Metapop1<-data.frame(mpop[[l]]) Metapop1<-subset(Metapop1, Metapop1[,75]!="0" ) Metapop<-rbind(Metapop1, Metapop) #assign("Metapop", Metapop, envir = .GlobalEnv) } # run GEE # order DF so id=orignPatch is in order (otherwise GEE will falsely count each new number in a row as a new patch) Metapop <- Metapop[order(Metapop$orignPatch),] m1 <- geeglm(Metapop$D ~ Metapop$E + Metapop$L, id = Metapop$orignPatch, data = Metapop, family = gaussian, corstr = "exchangeable") if(g == 1){ mOut[[1]] <- summary(m1)$coef mOut[[1]][,5] <- 1 } if(g%%10==0){ mOut[[(g/10)+1]] <- summary(m1)$coef mOut[[(g/10)+1]][5] <- g } } #### creating new generation #### newmpop <- mpop[[e]] # just to give the right names and structure for new dataframe newmpop[,] <- 0 # new generation will be empty before reproduction of parent generation if(sum(mpop[[e]][,7]) > 0){ # offspring will be generated only if at least one individual was able to collect resources # if no individual collected any resources no offspring is produced and the population gets extinct newmpop[1:N,] <- mpop[[e]][sample(nrow(mpop[[e]]), N , replace=T, prob = mpop[[e]][,7]*(1-mpop[[e]][,2]/a)),] # "reproduction depends on the sumS - L/a per individual newmpop[1:N,1] <- 1:N # give ID for N individuals. only N offspring will be allowed in total thus N = carrying capacity (if i want different K for patches, i need to specify K per patch and let sample k times) } # setting new data frame to start conditions newmpop[1:N,5] <- 1 # newmpop$dead <- 0 newmpop$sumS <- 0 # Number of Collected resources newmpop$RCol0 <- 0 newmpop$RCol1 <- 0 newmpop$RCol2 <- 0 newmpop$RCol3 <- 0 newmpop$RCol4 <- 0 newmpop$RCol5 <- 0 newmpop$RCol6 <- 0 newmpop$RCol7 <- 0 newmpop$RCol8 <- 0 newmpop$RCol9 <- 0 newmpop$RCol10 <- 0 # Number of Collected resources (corrected for efficiency and competition) newmpop$RCol0Cor = 0 newmpop$RCol1Cor = 0 newmpop$RCol2Cor = 0 newmpop$RCol3Cor = 0 newmpop$RCol4Cor = 0 newmpop$RCol5Cor = 0 newmpop$RCol6Cor = 0 newmpop$RCol7Cor = 0 newmpop$RCol8Cor = 0 newmpop$RCol9Cor = 0 newmpop$RCol10Cor = 0 # Number of Collected resources Values (corrected for efficiency and competition) newmpop$RCol0V = 0 newmpop$RCol1V = 0 newmpop$RCol2V = 0 newmpop$RCol3V = 0 newmpop$RCol4V = 0 newmpop$RCol5V = 0 newmpop$RCol6V = 0 newmpop$RCol7V = 0 newmpop$RCol8V = 0 newmpop$RCol9V = 0 newmpop$RCol10V = 0 # Resources Collected after Dispersal newmpop$RCol0aD <- 0 newmpop$RCol1aD <- 0 newmpop$RCol2aD <- 0 newmpop$RCol3aD <- 0 newmpop$RCol4aD <- 0 newmpop$RCol5aD <- 0 newmpop$RCol6aD <- 0 newmpop$RCol7aD <- 0 newmpop$RCol8aD <- 0 newmpop$RCol9aD <- 0 newmpop$RCol10aD <- 0 # Number of Collected resources (corrected for efficiency and competition) newmpop$RCol0CoraD = 0 newmpop$RCol1CoraD = 0 newmpop$RCol2CoraD = 0 newmpop$RCol3CoraD = 0 newmpop$RCol4CoraD = 0 newmpop$RCol5CoraD = 0 newmpop$RCol6CoraD = 0 newmpop$RCol7CoraD = 0 newmpop$RCol8CoraD = 0 newmpop$RCol9CoraD = 0 newmpop$RCol10CoraD = 0 # Number of Collected resources Values (corrected for efficiency and competition) newmpop$RCol0VaD = 0 newmpop$RCol1VaD = 0 newmpop$RCol2VaD = 0 newmpop$RCol3VaD = 0 newmpop$RCol4VaD = 0 newmpop$RCol5VaD = 0 newmpop$RCol6VaD = 0 newmpop$RCol7VaD = 0 newmpop$RCol8VaD = 0 newmpop$RCol9VaD = 0 newmpop$RCol10VaD = 0 newmpop$Disperser = 0 newmpop$Daysalive = 0 # in the begining it is yet not determined whether an indi will disperse or not # newmpop$orignPatch = e # will not be changed #### mutation #### # follows gauss distribution with mean = trait value of parent mutateL <- runif(nrow(mpop[[e]])) < qL ; #if(printIndi == 1){print(paste("mutateL",mutateL))} mutateE <- runif(nrow(mpop[[e]])) < qE ; #if(printIndi == 1){print(paste("mutateE",mutateE))} mutateD <- runif(nrow(mpop[[e]])) < qD ;# if(printIndi == 1){print(paste("mutatePicky", mutatePicky))} # mutation for trait L for(k in 1:N){ if(mutateL[[k]] == T & newmpop$ID[[k]]>0){ rm <-rnorm(1, mean = newmpop$L[[k]], sd = 0.1) # random mutation with mean of trait value of parent and sd of 0.1 if(rm <= 1 & rm >= 0){ # newmpop$L[[k]] <- rm } else if(rm > 1) { newmpop$L[[k]] <- 1 # max trait value } else{ newmpop$L[[k]] <- 0 # min trait value } } } # end of mutation L # mutation for trait E for(k in 1:N){ if(mutateE[[k]] == T & newmpop$ID[[k]]>0){ rm <-rnorm(1, mean = newmpop$E[[k]], sd = 0.1) # random mutation with mean of trait value of parent and sd of 0.1 if(rm <= 1 & rm >= 0){ # newmpop$E[[k]] <- rm } else if(rm > 1) { newmpop$E[[k]] <- 1 # max trait value } else{ newmpop$E[[k]] <- 0 # min trait value } } } # end of mutation E # mutation for trait Dispersal for(k in 1:N){ if(mutateD[[k]] == T & newmpop$ID[[k]]>0){ rm <-rnorm(1, mean = newmpop$D[[k]], sd = 0.1) # random mutation with mean of trait value of parent and sd of 0.1 if(rm <= 1 & rm >= 0){ # newmpop$D[[k]] <- rm } else if(rm > 1) { newmpop$D[[k]] <- 1 # max trait value } else{ newmpop$D[[k]] <- 0 # min trait value } } } # end of mutation D ### new pop ready to go ### if(g!=length(G)){ # overwrite old pop only until end of G so that in last generation the individuals result will not be overwritten mpop[[e]] <- newmpop } } # end of run through all pops # for(e in 1:length(patchID)){ #### random extinction patches #### if(Ext == "y"){ if(g < length(G)-1){ if(g%%EX==0){ ex <- sample(1:length(patchID), 1*EXn, replace = FALSE) # choose a patch randomly for(i in 1: length(ex)){ mpop[[ex[i]]][,] <- 0 # delete all individuals }}} } }# End Loop of all Generations #for(g in 1:length(G)) { #### assignments #### if(g==length(G) & e==length(patchID)){ # end of season assignments assign("mpop", mpop, envir = .GlobalEnv) assign("dfLE", dfLE, envir = .GlobalEnv) assign("meanL",meanL, envir = .GlobalEnv) assign("meanE",meanE, envir = .GlobalEnv) assign("meanD",meanD, envir = .GlobalEnv) assign("RColTab",RColTab, envir = .GlobalEnv) assign("patchDensBefore",patchDensBefore, envir = .GlobalEnv) assign("patchDensAfter",patchDensAfter, envir = .GlobalEnv) assign("suc_residents",suc_residents, envir = .GlobalEnv) assign("mOut",mOut, envir = .GlobalEnv) if(writeDF== 1){ saveRDS(dfLE, file = paste(path_toSaveto,"dfLE_", Cyc, sep = "")) # saves a list with a list for each pop saveRDS(mpop, file = paste(path_toSaveto,"mpop", Cyc, sep = "")) } #### meanPlot #### if(printPlot == 1){ meanPlot(N=N,G=G,patchEnv=patchEnv,saveplot=saveplot,path_to=path_to,SLbD=SLbD,SLaD=SLaD) } #### Tab_Summary #### Tab_Summary<-data.frame(matrix(0, nrow = 1, ncol = 11)) names(Tab_Summary)<-c("SL","N","G","Patch","meanL","meanE","meanD", "Estimate", "Std.err", "Wald", "Pr(>|W|)") Tab_Summary[,1]<-max(SLbD)*2 Tab_Summary[,2]<-N Tab_Summary[,3]<-max(G) Tab_Summary[,4]<-length(patchID) Tab_Summary[,5]<- mean(sapply(meanL, "[[",length(G))) Tab_Summary[,6]<- mean(sapply(meanE, "[[",length(G))) Tab_Summary[,7]<- mean(sapply(meanD, "[[",length(G))) Tab_Summary[,8:11]<- mOut[[length(mOut)]]["Metapop$L",1:4] assign("Tab_Summary", Tab_Summary, envir = .GlobalEnv) if(writeDF == 1){ write.csv(Tab_Summary, file = paste(path_toSaveto,"Tab_Summary", Cyc, sep = "")) } }# end of end-season assignments } # end of function