R/sim.R

In ggcostoya/limon: Evolutionary Constraints among Thermal Performance Traits

###*  SIMULATION FUNCTIONS

## SIMULATE SURVIVAL

#' Simulate survival
#'
#' Simulates the survival of an individual (defined by a TPC) at a particular environmental temperature.
#'
#' @param tpc An individual's TPC
#' @param te The environmental T
#' @param spsp Survival probability scaling parameter. Equivalent to the probability of survival at Pmax. Corresponds to the "p" parameter of a binomial distribution of size 1.
#' @param pmin The percentage of Pmax corresponding to Pmin.
#'
#' @return 1/0 If the individual survived or not.
#'
#' @export

sim_surv <- function(tpc, te, spsp, pmin){

  pmax <- get_tpts(tpc, pmin) %>% filter(tpt == "pmax") %>% select(value) %>% as.numeric() # extract pmax

  p_t <- tpc %>% filter(t == Closest(t,te)) %>% select(p) %>% as.numeric() # get performance at environmental temperature

  p_r <- p_t/pmax # performance relative to pmax

  s <- p_r*spsp # actual probability of survival

  s <- ifelse(s > 1, 1, s) # keep s below 1

  s <- ifelse(is.na(s), 0, s) # if survival is NA set it to 0

  survival <- rbinom(1,1,s) # determine actual survival

  return(as.numeric(survival))

}

## SIMULATE REPRODUCTION

#' Simulate reproduction
#'
#' Simulates the reproduction of an individual (defined by a TPC) at a particular environmental temperature
#'
#' @param tpc The individuals' TPC.
#' @param te The environmental T
#' @param rpsp Reproduction probability scaling parameter. Equivalent to the probability of reproduction at Pmax. Corresponds to the "lambda" parameter of a poisson distribution.
#' @param pmin The percentage of Pmax coresponding to Pmin
#'
#' @return An integer indicating the number of offspring produced
#'
#' @export

sim_repr <- function(tpc, te, rpsp, pmin){

  pmax <- get_tpts(tpc, pmin) %>% filter(tpt == "pmax") %>% select(value) %>% as.numeric() # extract pmax from the tpc

  p_t <- tpc %>% filter(t == Closest(t,te)) %>% select(p) %>% as.numeric() # get the performance at environmental temperature

  p_r <- p_t/pmax # performance relative to pmax

  r <- p_r * rpsp # actual probability of reproduction

  r <- ifelse(is.na(r), 0, r) # if reproduction probability is NA set it to 0

  n_off <- rpois(1,r) # determine the number of offspring generated where r is the lambda of a Poisson distribution

  return(as.numeric(n_off))

}

### SIMULATE OFFSPRING GENERATION

#' Simulate offspring
#'
#' Simulates offspring produced by an individual defined by some thermal-performance traits (TPTs).
#'
#' @param tpts The individual's TPTs
#' @param n_off The number of offspring to be generated by the individual
#' @param mu The mutation probability (constant for all TPTs)
#' @param sdm The percentage of the original TPT value that will act as the standard deviation of the normal distribution of mean = 0 from which the amount of change occuring on a mutating trait will be sampled.
#' @param cmtx A genetic correlations matrix for the TPTs. This parameter is optional and CANNOT INCLUDE Pmin.
#' @param pmin The percentage of Pmax corresponding to Pmin.
#' @param samples The number of samples to be generated for each offspring.
#' @param error The error to be introduced into the offspring's TPD in units of standard deviation of a normal distribution with 0 as mean.
#'
#' @return A TPD for the offspring generated with id, t and p as columns.
#'
#' @export

sim_offspring <- function(tpts, n_off, mu, sdm, cmtx, pmin, samples, error){

  # preparing the holder data
  offspring_tpd <- tibble(id = character(), t = logical(), p = logical())

  # loop to repeat the process for all offspring
  for(i in 1:n_off){

    # take pmin from tpts
    tpts <- filter(tpts, tpt != "pmin")

    # determine in which traits does mutation happen
    where_mutation <- rbinom(nrow(tpts),1,mu)

    # determine the amount of change due mutation
    change_mutation <- map_dbl(as.numeric(tpts$value), function(x){rnorm(1,0,sd = abs(sdm*x))})*where_mutation

    # apply the change to the tpts

    # if no correlations matrix is provided
    if(missing(cmtx)){

      # determine offspring tpts
      offspring_tpts <- tpts %>% transmute(tpt = tpt, value = value + change_mutation)

      # if correlation matrix is provided
    }else{

      # changes due correlation between traits
      change_corr <- as.vector(change_mutation[1:6] %*% cmtx)

      # determine offspring tpts
      offspring_tpts <- tpts %>% transmute(tpt = tpt, value = value + change_corr)

    }

    # add pmin back to the offspring tpts
    add_pmin <- tibble(tpt = "pmin", value = pmin)
    offspring_tpts <- bind_rows(offspring_tpts, add_pmin)

    # generate tpd based on the tpts
    offspring <- gen_tpd(tpts = offspring_tpts,samples = samples, error = error)

    # add id to the new individual
    offspring<- bind_cols(tibble(id = rep(gen_id(), nrow(offspring))),offspring)

    # bind the offspring data
    offspring_tpd <- bind_rows(offspring_tpd, offspring)

  }

  return(offspring_tpd)

}


### * Simulate survival throughout a generation

#' Simulate survival throughout a generation
#'
#' Simulates an individual's survival for a period of time corresponding to a generation.
#'
#' @param tpc An individual's thermal-performance curve (TPC)
#' @param tseq A vector of temperature values which will also determine how many days conform a generation.
#' @param spsp Survival probability scaling parameter. Equivalent to the probability of survival at Pmax. Corresponds to the "p" parameter of a binomial distribution of size 1.
#' @param pmin The percentage of Pmax corresponding to Pmin.
#'
#' @return A vector of 0/1 for survival at each moment in time. If a 0 occurs, all subsequent values are 0.
#'
#' @examples
#'
#' @export

sim_surv_gen <- function(tpc, tseq, spsp, pmin){

  # holder vector of survival values
  survival <- rep(NA,length(tseq))

  # loop to determine survival
  for(i in 1:length(tseq)){survival[i] <- sim_surv(tpc, tseq[i], spsp = spsp, pmin = pmin)}

  # no revival
  for(i in 2:length(survival)){survival[i] <- ifelse(survival[i-1] == 0 & survival[i] == 1, 0, survival[i])}

  return(survival)

}

### * Simulate reproduction throughout a generation

#' Simulate reproduction throughout a generation.
#'
#' Simulates an individual's reproductive output throughout a generation.
#'
#' @param tpc An individual's thermal-performance curve (TPC)
#' @param tseq A vector of temperature values which will also determine how many days conform a generation.
#' @param surv A vector of survival values (0/1) in the format of the output of sim_surv_gen.
#' @param rpsp Reproduction probability scaling parameter. Equivalent to the probability of reproduction at Pmax. Corresponds to the "lambda" parameter of a poisson distribution.
#' @param pmin The percentage of Pmax corresponding to Pmin.
#'
#' @return A vector of integers determining offspring produced each day the individual survived.
#'
#' @examples
#'
#' @export

sim_repr_gen <- function(tpc, tseq, surv, rpsp, pmin){

  # subset only the number of days in which the individual survived
  days_alive <- surv[surv > 0]

  # subset temperature sequence for days in which individual was alive
  tseq <- tseq[1:length(days_alive)]

  # holder vector of reproduction values
  reproduction <- rep(NA, length(tseq))

  for(i in 1:length(tseq)){reproduction[i] <- sim_repr(tpc, tseq[i], rpsp = rpsp, pmin = pmin)}

  return(reproduction)

}

### * Simulate generation temperature sequence.

#' Simulate generation temperature sequence.
#'
#' Simulates a sequence of temperature values for the length of an entire generation
#'
#' @param mean_te Mean temperature
#' @param sd_te Standard deviation temperature
#' @param n_days Number of days conforming the generation.
#'
#' @return A sequence of temperature values
#'
#' @examples
#' tseq <- sim_gen_tseq(30,5,365)
#' plot(tseq, type = "l", xlab = "Day", ylab = "Temperature")
#'
#' @export

sim_gen_tseq <- function(mean_te, sd_te, n_days){

  tseq <- rnorm(n_days, mean = mean_te, sd = sd_te)

  return(tseq)
}

### * Simulate multiple generation temperature sequence

#' Simulate multiple temperature sequences
#'
#' Simulates a sequence of temperature values spaning over multiple generations
#'
#' @param mean_tes A vector of mean temperature values as long as the number of generations to generate.
#' @param sd_tes A vector of sd. values as long as the number of generations to simulate.
#' @param n_days The number of days per generation.
#'
#' @return A tibble of two columns, generation number and temperature "t"
#'
#' @examples
#' tseq <- sim_tseqs(mean_tes = c(30,31,32), sd_tes = c(0.5,1,1.5), n_days = 365)
#' plot(tseq$t, type = "l", xlab = "Day", ylab = "Temperature", lwd = 1.5)
#' abline(v = c(365,365*2), lty = 2, col = "lightgrey")
#'
#' @export

sim_tseqs <- function(mean_tes, sd_tes, n_days){

  # empty holder data frame
  tseqs <- tibble(gen = logical(), t = logical())

  # loop to generate temperature data
  for(i in 1:length(mean_tes)){

    # generate temperature sequence for the current generation & put it nicely in a tibble
    tseq_gen <- tibble(gen = rep(i,n_days), t = sim_gen_tseq(mean_te = mean_tes[i], sd_te = sd_tes[i], n_days = n_days))

    # attach it to existing dataset
    tseqs <- bind_rows(tseqs, tseq_gen)

  }

  return(tseqs)

}

### * Simulate a generation

#' Simulate a generation
#'
#' Simulates the survival and reproduction of a population over a generation.
#'
#' @param pop_tpd The thermal-performance dataset (TPD) of the original population.
#' @param tseq A sequence of temperature values generated by the sim_gen_tseq function
#' @param spsp Survival probability scaling parameter. Equivalent to the probability of survival at Pmax. Corresponds to the "p" parameter of a binomial distribution of size 1.
#' @param rpsp Reproduction probability scaling parameter. Equivalent to the probability of reproduction at Pmax. Corresponds to the "lambda" parameter of a poisson distribution.
#' @param sdm The percentage of the original TPT value that will act as the standard deviation of the normal distribution of mean = 0 from which the amount of change occuring on a mutating trait will be sampled.
#' @param cmtx A genetic correlations matrix for the TPTs. This parameter is optional and CANNOT INCLUDE Pmin.
#' @param pmin The percentage of Pmax corresponding to Pmin.
#' @param samples The number of samples to be generated for each offspring.
#' @param error The error to be introduced into the offspring's TPD in units of standard deviation of a normal distribution with 0 as mean.
#'
#' @return A TPD with the original and the new generation of individuals.
#'
#' @examples
#'
#' @export

sim_generation <- function(pop_tpd, tseq, spsp, rpsp, cmtx, mu, sdm, pmin, samples, error){

  # set the generation column if the starting popultion tpd does not have it
  if(!"gen" %in% colnames(pop_tpd)){pop_tpd <- as_tibble(cbind(gen = rep(0, nrow (pop_tpd)), pop_tpd))}

  # holder tibble for next generation
  tpdt1 <- tibble(gen = logical(), id = character(), t = logical(), p = logical())

  # select current generation
  tpdt0 <- pop_tpd %>% filter(gen == max(gen,na.rm = T)) %>% nest(tpd = c(t,p))

  # start loop to assign survival, reproduction and offspring for each individual
  for(i in 1:nrow(tpdt0)){

    # extract an individual data
    ind_tpd <- tpdt0 %>% slice(i) %>% select(tpd) %>% unnest(cols = c(tpd)) # tpd

    # determine the individual's tpc and tpts
    ind_tpc <- ind_tpd %>% fit_tpd() %>% gen_tpc()
    ind_tpts <- get_tpts(ind_tpc, pmin)

    # determine individual survival throughout the geneartion
    ind_surv <- sim_surv_gen(ind_tpc, tseq, spsp, pmin)

    # if the individual survived till the end of the generation include in next generation
    if(last(ind_surv) == 1){

      # select the individual's data for next generation and add 1 to the gen
      ind <- tpdt0 %>% slice(i) %>% mutate(gen = gen + 1) %>% unnest(cols = c(tpd))

      # bind the individual's data to the tpdt1 dataset
      tpdt1 <- suppressWarnings(bind_rows(tpdt1, ind))

    }

    # determine individual reproduction throughout the generation
    ind_repr <- sum(sim_repr_gen(ind_tpc, tseq, ind_surv, rpsp, pmin))

    # if individual reproduction is bigger than zero then generate offspring data
    if(ind_repr > 0){

      # if no correlations matrix is provided
      if(missing(cmtx)){

        # generate individual offspring
        ind_offspring <- sim_offspring(tpts = ind_tpts, n_off = ind_repr, mu = mu, sdm = sdm, pmin = pmin, samples = samples, error = error)

        # if the correlations matrix is provided
      }else{

        # generate individual offspring
        ind_offspring <- sim_offspring(tpts = ind_tpts, n_off = ind_repr, mu = mu, sdm = sdm, cmtx = cmtx, pmin = pmin, samples = samples, error = error)

      }

      # add generation data to the ind_offspring dataset
      ind_offspring <- as_tibble(cbind(gen = rep(mean(tpdt0$gen, na.rm = T) + 1, nrow(ind_offspring)), ind_offspring))

      # bind the individual's offspring tpd to th tpdt1 dataset
      tpdt1 <- suppressWarnings(bind_rows(tpdt1, ind_offspring))

    }

  }

  # bind tpdt1 to the simulation dataset (pop_tpd)
  pop_tpd <- bind_rows(pop_tpd, tpdt1)

  return(pop_tpd)

}

# Simulation

#' Simulation
#'
#' Simulates the evolution of a population based on thermal conditions
#'
#' @param pop_tpd The thermal-performance dataset (TPD) of the original population
#' @param tseqs A dataset of temperature values spanning multiple generations in the form of the output of the sim_tseqs function. It's length determines the number of generations.
#' @param spsp Survival probability scaling parameter. Equivalent to the probability of survival at Pmax. Corresponds to the "p" parameter of a binomial distribution of size 1.
#' @param rpsp Reproduction probability scaling parameter. Equivalent to the probability of reproduction at Pmax. Corresponds to the "lambda" parameter of a poisson distribution.
#' @param sdm The percentage of the original TPT value that will act as the standard deviation of the normal distribution of mean = 0 from which the amount of change occuring on a mutating trait will be sampled.
#' @param cmtx A genetic correlations matrix for the TPTs. This parameter is optional and CANNOT INCLUDE Pmin.
#' @param pmin The percentage of Pmax corresponding to Pmin.
#' @param samples The number of samples to be generated for each offspring.
#' @param error The error to be introduced into the offspring's TPD in units of standard deviation of a normal distribution with 0 as mean.
#'
#' @return A TPD with the original and all subsequent simulated generations of individuals
#'
#' @examples
#'
#' @export

simulation <- function(pop_tpd, tseqs, spsp, rpsp, cmtx, mu, sdm, pmin, samples, error){

  # determine number of generations from tseqs dataset
  n_gen <- max(tseqs$gen, na.rm = T)

  # loop to run the simulation each generation
  for(i in 1:n_gen){

    # select temperature data for that generation
    tseq_gen <- tseqs %>% filter(gen == i) %>% select(t) %>% as.data.frame() %>% .[,1]

    # run the simulation for that generation

    # if no correlations matrix is provided
    if(missing(cmtx)){

      # simulate generation
      new_pop_tpd <- sim_generation(pop_tpd = pop_tpd, tseq = tseq_gen, spsp = spsp, rpsp = rpsp, mu = mu, sdm = sdm, pmin = pmin, samples = samples, error = error)

      # if the correlations matrix is provided
    }else{

      # simulate generation
      new_pop_tpd <- sim_generation(pop_tpd = pop_tpd, tseq = tseq_gen, spsp = spsp, rpsp = rpsp, mu = mu, cmtx = ctmx, sdm = sdm, pmin = pmin, samples = samples, error = error)

    }

    # stop the loop if no rows are added (if none survived or reproduced)
    if(nrow(pop_tpd) == nrow(new_pop_tpd)){break}

    # reassign pop_tpd
    pop_tpd <- new_pop_tpd

  }

  return(pop_tpd)

}