Difference between revisions of "Population growth"

library(OpasnetUtils)

objects.latest(Op_enXXXX) # growth2

growth2(100, 30, 0.1, 0.05, 10)

# This is code Op_enXXXX/initiate on  page [[Population growth]]

library(OpasnetUtils)

growth2 <-function(
  N,    # size of population in the beginning
  tmax, # end of follow-up time
  r,    # growth rate (fraction per time unit)
  sdr,  # standard deviation of r
  it    # number of iterations
){
  popsize <- data.frame(time=c(0:tmax)) # the data.frame for the final data with tmax rows
  for(i in 1:it){   # doing all this it times
    pop.size <- c(N, rep(NA, tmax))
    for(t in 1:tmax){                   # this is the same bit as before
      nt <- ifelse(t==1, N, nt1)
      nt1 <- nt*exp(rnorm(1, mean=r, sd=sdr))
      pop.size[t+1] <- nt1
    }
    popsize <- cbind(popsize, pop.size) # adding the data of each iteration to the same data frame
  }
  popsize$mean.size <- rowMeans(popsize[,c(2:it)]) # mean sizes for times from iterations
  quantiles <- data.frame()
  # with a loop calculating quantiles for each row, aka. each time.
  for(p in 1:tmax+1){
    quantiles <- rbind(quantiles, quantile(popsize[p,c(2:it+1)], probs=c(0.025,0.975)))
  }
  quantiles <- stack(quantiles) #flipping the quantile
  quantiles[c((tmax+1):(2*tmax)),"values"] <- rev(quantiles[c((tmax+1):(2*tmax)),"values"])
  quantiles$time <- c(1:tmax, tmax:1) # adding the time as x-axis points
  
  plot(x=popsize$time, y=popsize$mean.size, type="n",
       xlab = "Time", ylab="Population size", main="Population growth over time",
       ylim = c(0, max(quantiles$value))) #y lim so that the upper quantile is also visible
  polygon(quantiles$time, quantiles$values, col="palegreen") # quantile lines
  lines(popsize$time, popsize$mean.size) # pop size line on top
}

objects.store(growth2)
cat("Function growth2 stored.\n")