R/plotCobweb.R

In tessington/quantecol: Introduction to Quantitative Ecology

plotCobweb <- function(modeloutput, modelname, modelpars, viewcode = FALSE) {
  #' Cobweb plots for logistic, Ricker, Gompertz, or theta-logistic models.
  #'
  #' @param modeloutput an array of population sizes at time "t" generated by a discrete time model.  Equivalent to "output" used in Section 14.1.2
  #' @param modelname character string indicating model used: "logistic", "ricker", "gompertz", "theta-logistic"
  #' @param modelpars a list containing model parameters.  For logistic, this will be list(r = r, K = K).  For ricker, use list(alphal = alphal, K = K).  For gompertz, use list(b = b, K = K).  For theta-logistic, use list(r = r, K = K, theta = theta)
  #' @param viewcode TRUE or FALSE (default) indicating whether to print the function code
  #'
  #'
  #' @return a cobweb plot in the figure window
  #' @export
  #'
  #' @examples
  #'
  #' #Logistic Equation:
  #' r <- 1.5
  #' K <- 100
  #' N0 <- 10
  #' tmax <- 50
  #' years <- 0:tmax
  #' output <- rep(x = NA, times = length(years))
  #' output[1] <- N0
  #' n.loop <- length(years)
  #' for (i in 2:n.loop) {
  #'   n.t <- output[i-1]
  #'   n.t.plus.1 <- n.t + r * n.t * (1 - n.t / K)
  #'   output[i] <- n.t.plus.1
  #' }
  #' plotCobweb(modeloutput = output, modelpars = list(r = r, K = K), modelname = "logistic")
  #'
  #'
  #'
  # Error checking
  if(!modelname  %in% c("logistic", "ricker", "gompertz", "theta-logistic")) {
    stop("Incorrect modelname specification.  Must be either logistic, ricker, gompertz or theta-logistic")
  }

  if(is.null(modelpars$K)) stop("parameter K not provided in modelpars list")

  # more error checking for parameters
  if (modelname == "logistic" | modelname == "theta-logistic") {
    if(is.null(modelpars$r)) stop("parameter r not provied in modelpars list")
  }
  if (modelname == "gompertz") {
    if(is.null(modelpars$b)) stop("parameter b not provided in modelpars list")
    if(modelpars$b>0) stop("parameter b must be negative")
  }
  if (modelname =="ricker") {
    if(is.null(modelpars$alphal)) stop("parameter alphal not provided in modelpars list")
  }
  if (modelname == "theta-logistic") {
    if(is.null(modelpars$theta)) stop("parameter theta not provided in modelpars list")
  }
  if (length(modeloutput)<2) stop("modeloutput must be an array with at least two elements")

  # End of error checking

  # extract carrying capacity from modelpars list
  K <- modelpars$K
  # make a list of N at time t
  nt<- seq(from = 0, to = 2.0 * K, length.out = 100)
  # do this if model is "logistic"
  if(modelname == "logistic") {

    # extract r from modelpars list
    r <- modelpars$r
    # calculate Nt+1 from Nt
    ntplus1 <- nt + r * nt* (1 - nt / K)
  }
  if(modelname == "ricker") {
    # extract alpha from modelpars list
    alphal <- modelpars$alphal
    # calculate Nt+1 from Nt
    ntplus1<- nt* exp(alphal * (1 - nt / K))
  }
  if(modelname == "gompertz") {
    # extract b from modelpars list
    b = modelpars$b
    nt[1] <- 0.001 # not defined for nt = 0
    # calculate Nt+1 from Nt
    ntplus1 <-nt* (nt/ K)^b
  }
  if(modelname == "theta-logistic") {
    # extract r from modelpars list
  r <- modelpars$r
  theta <- modelpars$theta
  # calculate Nt+1 from Nt
  ntplus1<- nt + r * nt * (1 - (nt/K)^theta)
  }

  # plot recursive function
  plot(x = nt, y = ntplus1,
       type = "l",
       lwd = 3,
       xlab = expression(paste("N"["t"], sep = "")),
       ylab = expression(paste("N"["t"+1], sep = "")),
       ylim = c(0, 2.0 * K),
       xlim <- c(0, 2.0 * K),
       las =1,
       cex.lab = 1.25,
       cex.axis = 1.25
  )
  # Add Nt = Nt+1 line (Replacement line)
  abline (a = 0,
          b = 1,
          col = "gray",
          lwd = 3)

  # plot vertical lines
  n.loops <- length(modeloutput) -1
  for (i in 1:n.loops){
    xs <- rep(x = modeloutput[i], times = 2)
    ys <- c(modeloutput[i], modeloutput[i+1])
    lines(xs, ys,
          lwd = 1
    )
  }
  # plot horizontal lines
  for (i in 1:n.loops) {
    xs <- c(modeloutput[i], modeloutput[i+1])
    ys <- rep (x = modeloutput[i+1], times = 2)
    lines(xs, ys,
          lwd = 1
    )
  }
  if(viewcode) cat(' # extract carrying capacity from modelpars list
  K <- modelpars$K
  # make a list of N at time t
  nt<- seq(from = 0, to = 1.5 * K, length.out = 100)
  # do this if model is "logistic"
  if(modelname == "logistic") {
    # extract r from modelpars list
    r <- modelpars$r
    # calculate Nt+1 from Nt
    ntplus1 <- nt + r * nt* (1 - nt / K)
  }
  # do this if model is "ricker"
  if(modelname == "ricker") {
    # extract alpha from modelpars list
    alphal <- modelpars$alphal
    # calculate Nt+1 from Nt
    ntplus1<- nt* exp(alphal * (1 - nt / K))
  }
  # do this if model is "gompertz"
  if(modelname == "gompertz") {
    # extract b from modelpars list
    b = modelpars$b
    nt[1] <- 0.001 # not defined for nt = 0
    # calculate Nt+1 from Nt
    ntplus1 <-nt* (nt/ K)^b
  }
  # do this if model is "theta-logistic
  if(modelname == "theta-logistic") {
    # extract r from modelpars list
  r <- modelpars$r
  theta <- modelpars$theta
  # calculate Nt+1 from Nt
  ntplus1<- nt + r * nt * (1 - (nt/K)^theta)
  }

  # plot recursive function
  plot(x = nt, y = ntplus1,
       type = "l",
       lwd = 3,
       xlab = expression(paste("N"["t"], sep = "")),
       ylab = expression(paste("N"["t"+1], sep = "")),
       ylim = c(0, 1.5 * K),
       xlim <- c(0, 1.5 * K),
       las =1,
       cex.lab = 1.25,
       cex.axis = 1.25
  )
  # Add Nt = Nt+1 line (Replacement line)
  abline (a = 0,
          b = 1,
          col = "gray",
          lwd = 3)

  # plot vertical lines
  n.loops <- length(modeloutput) -1
  for (i in 1:n.loops){
    xs <- rep(x = modeloutput[i], times = 2)
    ys <- c(modeloutput[i], modeloutput[i+1])
    lines(xs, ys,
          lwd = 1
    )
  }
  # plot horizontal lines
  for (i in 1:n.loops) {
    xs <- c(modeloutput[i], modeloutput[i+1])
    ys <- rep (x = modeloutput[i+1], times = 2)
    lines(xs, ys,
          lwd = 1
    )
  }
', sep = "\n")
  }