R/sim_surv.R

In ggcostoya/tpcurves2: T-P Curves

##* SIMULATE SURVIVAL

#' Simulate survival
#'
#' Simulates the survival of a population based on temperature performance
#'
#' @param TPD A population's thermal-performance data.
#'
#' @param Te The environmental temperature.
#'
#' @param Sp The relationship between survival and performance data.
#'
#' @param Degree The degree to fit the polynomial regressions relating Tm to P.
#'
#' @param Pmin The mimum performance to be considered.
#'
#' @return A new TP population dataset with only survivors from TPD adding to the ID, P and T columns the "Parent" column to identify parentage.
#'
#' @examples
#'
#' @export

sim_surv <- function(TPD, Te, Sp, Degree, Pmin){

  # add a paternity column if no parentage is provided.
  if(is.null(TPD$Parent)){TPD$Parent <- rep(NA, nrow(TPD))}

  # nest TPD (join Tm and P into a nested dataframe for each unique ind ID).
  TPD <- TPD %>% nest(TPD = c(Tm, P))

  # define the Survivors TPD, for now, empty rows.
  Survivors <- data.frame(ID = NA, Tm = NA, P = NA, Parent = NA)

  # start of the loop to determine surviving individual.
  for(i in 1:nrow(TPD)){

    # select an individual from the population TPD and unnest it's TPD.
    ind <- TPD[i,] %>% unnest(cols = c(TPD))

    # get the TPC from the individual's TPD using the provided Degree & Pmin arguments.
    TPCind <- get_tpc(Tm = ind$Tm, P = ind$P, Degree = Degree, Pmin = Pmin)

    # extract the individual's Pmax from the TPC.
    Pmax <- TPCind %>% filter(P == max(P, na.rm = T)) %>% select(P) %>% as.numeric()

    # extract individual's Pte by finding the Tm value in the TPC more closely matching Te.
    Pte <- TPCind %>% filter(Tm == Closest(Tm, Te)) %>% select(P) %>% as.numeric()

    # define the individual's relative performance to the maximum at Te. To keep Pr < 1 if Pte/Pmax > 1 then Pr = 1.
    Pr <- ifelse(Pte/Pmax < 1, Pte/Pmax, 1)

    # define the probability of survival using Pr and Sp. Sp !> 1.
    S <- Pr * Sp

    # determine if the individual survives.
    Survival <- rbinom(1, 1, S)

    # if the individual survives add its TPD to the Survivors population
    if(Survival == 1){Survivors <- rbind(Survivors, ind)} else {Survivors <- Survivors}

  }

  # tidy Survivors data and offer return.
  Survivors <- Survivors[-c(1),]
  return(Survivors)

}