R/sim_repr.R

In ggcostoya/tpcurves2: T-P Curves

##* SIMULATE REPRODUCTION

#' Simulate reproduction
#'
#' Simulates the reproduction of a population
#'
#' @param TPD A population's thermal performance data.
#'
#' @param Te The environmental temperature
#'
#' @param Rp The relationship between performance and reproduction.
#'
#' @param Mu Mutation rate
#'
#' @param Mum Mean percentage of change on a TPT due to mutation
#'
#' @param Musd SD of the percentage of change on a TPT due to mutation
#'
#' @param Degree The degree to fit the polynomial regressions relating Tm to P.
#'
#' @param Pmin The mimum performance to be considered.
#'
#' @param Samples Samples to be taken from each new individual.
#'
#' @param Error Amount of error to be introduced in offspring samples.
#'
#' @return A new TP population dataset adding to the ID, P and T columns the "Parent" column to identify parents.
#'
#' @examples
#'
#' @export

sim_repr <- function(TPD, Te, Rp, Mu, Mum, Musd, Degree, Pmin, Samples, Error){

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

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

  # start of the loop to generate new individuals.
  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 individual's TPTs
    ind_tpts <- get_tpts(Tm = ind$Tm, P = ind$P, Degree = Degree, Pmin = Pmin)

    # get the individual's id
    ind_id <- as.character(ind$ID[1])

    # get the individual's TPC
    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)) %>% select(P) %>% as.numeric()

    # extract the 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)

    # determine the probability of survival using Pr and Rp. Rp can > 1
    R <- Pr * Rp

    # round down R and use binomial to sample rest to determine number of offspring
    N_off <- floor(R) + rbinom(1, 1, R - floor(R))

    # if N_off is >= 1 then generate the offspring
    if(N_off >=1){

      # loop for each of the individual's offspring
      for(k in 1:N_off){

        # get the base offspring's TPTs, first as a copy of the parents'
        off_tpts <- ind_tpts

        # determine TPTs in which mutation occurs
        mut_where <- rbinom(nrow(off_tpts), 1, Mu)

        # loop to implement mutations on the traits
        for(t in 1:nrow(off_tpts)){

          # set new trait value after if mutation occurs
          off_tpts$value[t] <- off_tpts$value[t] + mut_where[t]*rnorm(1, mean = off_tpts$value[t]*Mum, sd = off_tpts$value[t]*Musd)

        }

        # generate offspring's TPD from TPTs
        off_tpd <- gen_tpd(Topt = off_tpts$value[1], CTmax = off_tpts$value[2], CTmin = off_tpts$value[3],
                           Pmax = off_tpts$value[5], Pmin = Pmin, Error = Error, Samples = Samples, Degree = Degree )

        # generate offspring id, add parent identity and TPD
        off_tpd <- cbind(ID = rep(gen_id(), nrow(off_tpd)), off_tpd, Parent = rep(ind_id, nrow(off_tpd)))

        # add new offspring's TPD to the Offsrping population
        Offspring <- rbind(Offspring, off_tpd)

      }

    }

  }

  # tidy OFffspring data and offer return
  Offspring <- Offspring[-c(1),]
  return(Offspring)

}