Revision Update

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: AutoPop
 Title: Autopolyploid population simulations
-Version: 1.0.0.0000
+Version: 1.0.0.1
 Authors@R: 
     c(person(given = "Michelle", family ="Gaynor", middle = "L", email = "[email protected]", role = c("aut", "cre"),
         comment = c(ORCID = "0000-0002-3912-6079")), 
@@ -18,7 +18,7 @@ Description: R-based joint-dynamic population simulation for
 License: GPL (>= 3)
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 Suggests: 
     testthat (>= 3.0.0)
 Imports:

NAMESPACE

@@ -1,8 +1,10 @@
 # Generated by roxygen2: do not edit by hand
 
 importFrom(MASS,ginv)
+importFrom(MASS,mvrnorm)
 importFrom(stats,dhyper)
 importFrom(stats,na.omit)
+importFrom(stats,pnorm)
 importFrom(stats,rbeta)
 importFrom(stats,rbinom)
 importFrom(stats,rmultinom)

R/alphabeta.calc.R

@@ -1,8 +1,9 @@
 #' Mature survival - Calculates alpha and beta from mean and variance.
 #'
-#' @description Calculates the alpha and beta parameters based on mu, the mean probability of mature survival, and the calculated variance.
+#' @description Calculates the alpha and beta parameters based on mu, here the mean probability of mature survival,
+#' and s calculated variance.
 #'
-#' @param mu Mean probability of mature survival. Must be a single integer between 0 and 1.
+#' @param mu Sample mean. Must be a single integer between 0 and 1.
 #' @param var Derived variance based on var.option. Must be a single integer between 0 and 1.
 #'
 #' @returns List with the parameters alpha and beta.

R/env.copula.R

@@ -0,0 +1,32 @@
+#' Environmental Copula -  Calculates correlation probability.
+#'
+#' @description Calculates the probability coefficient which allow cytotype and stage correlations
+#'
+#'
+#' @param rho Correlation coefficient. Must be a single integer between 0 and 1.
+#' @inheritParams gen.iter.f.choosy
+#'
+#' @returns A list of two matrix: X and U. Here, X is the matrix which was simulated from
+#' a multivariate normal distribution. The U matrix is simply the X matrix transformed into a uniform distribution.
+#' Each matrix has 3 columns, one for each cytotype. The length of the matrix matches the number of generations provided.
+#'
+#' @importFrom MASS mvrnorm
+#' @importFrom stats pnorm
+#'
+
+
+env.copula <- function(rho, generations){
+  # Create covariance matrix
+  n <- 3 # Three cytotypes
+  mean <- rep(0,n)
+  # covariance matrix for the cytotypes
+  SIGMA <- ((1-rho)*diag(n) + rho*matrix(1, nrow = n, ncol = n))
+  # Simulate from a multivariate normal distribution
+  X <- MASS::mvrnorm(n = generations, mu = mean, Sigma = SIGMA)
+  # Transform to uniform using pnorm, which is the cdf for a normal distribution
+  U <- pnorm(X)
+  out <- list(X, U)
+  return(out)
+}
+
+

R/form.autopop.R

@@ -26,10 +26,10 @@
 #' * C2: relative abundance of all diploids (ie. sum2x/sum).
 #' * C3: relative abundance of all triploids (ie. sum3x/sum).
 #' * C4: relative abundance of all tetraploids (ie. sum4x/sum)
-#'
+
 
 form.autopop <- function(popvect, generations){
-    plot.pop <- as.data.frame(do.call(rbind, popvect))
+    plot.pop <- as.data.frame(do.call(rbind, popvect))[,1:12]
     plot.pop$gen <- 1:(generations + 1)
     plot.pop[is.na(plot.pop)] = 0
 

R/gen.iter.f.choosy.R

@@ -19,6 +19,7 @@
 #' @inheritParams one.iter.f.choosy
 #' @inheritParams cytotype_repro_mate
 #' @inheritParams form.autopop
+#' @inheritParams env.copula
 #'
 #' @returns A single data frame as defined by `form.autopop()`. Each row is a generation. The columns are as follows,
 #' * V1: number of immature diploids.
@@ -47,42 +48,69 @@
 #' * C2: relative abundance of all diploids (sum2x/sum).
 #' * C3: relative abundance of all triploids (sum3x/sum).
 #' * C4: relative abundance of all tetraploids (sum4x/sum)
-#'
-#'
-
+#' * i1e: diploid immature survival probability during t - 1.
+#' * i2e: triploid immature survival probability during t - 1.
+#' * i3e: tetraploid immature survival probability during t - 1.
+#' * m1e: diploid mature survival probability during t - 1.
+#' * m2e: triploid mature survival probability during t - 1.
+#' * m3e: tetraploid mature survival probability during t - 1.
 
 
 #source("R/one-iter-f-choosy.R")
 #source("R/format-iter-choosy.R")
 #source("R/mature.surv.calc.R")
 
-gen.iter.f.choosy <- function(generations, init.pop, env.ci, aii.vec,
+gen.iter.f.choosy <- function(generations, init.pop, env.ci, env.sigma, aii.vec,
                               as.matur, as.msurv, d, gnum.base,
-                              b, cc, s, mc, density.type = "all",
-                              mate.lazy = FALSE){
+                              b, cc, s, mc,
+                              gam.density.type = "all", is.density.type = "all",  rho, mate.lazy = FALSE,
+                              rj = c(1,1),  yj = c(1,1),  env.immature.survival = FALSE){
 
   ## Set up first generation
   popvect <- list()
-  first.vec <- c(0, 0, 0, init.pop, 0, 0, 0, 0, 0, 0, 0, 0)
+  first.vec <- c(0, 0, 0, init.pop, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
   popvect[[1]] <- first.vec
-  as.msurv.e <- list()
-  as.msurv.e[[1]] <- as.msurv
-
 
+  as.msurv.e <- list()
+  as.msurv.e[[1]] <- c(0,0,0)
+  Xi.Ui.vec <- env.copula(rho = rho, generations = generations)
+  Xi.vec <- Xi.Ui.vec[[1]]
+  Ui.vec <- Xi.Ui.vec[[2]]
+  if(sum(as.msurv) != 0){
+  msurv.alpha.beta <- surv.shape(env.ci = env.ci, raw.means = as.msurv)
+  }
   ## Loop for set number of generations #########################################
   for(i in 1:generations){
-    # Environmental. Stochasticity
-    m1e <- mature.surv.calc(env.ci = env.ci, as.msurv = as.msurv.e[[1]][1])
-    m2e <- mature.surv.calc(env.ci = env.ci, as.msurv = as.msurv.e[[1]][2])
-    m3e <- mature.surv.calc(env.ci = env.ci, as.msurv = as.msurv.e[[1]][3])
+    if(sum(as.msurv) != 0){
+    # Mature Survival - Environmental. Stochasticity
+    m1e <- qbeta(Ui.vec[i,1], msurv.alpha.beta[1,1], msurv.alpha.beta[2,1])
+    m2e <- qbeta(Ui.vec[i,2], msurv.alpha.beta[1,2], msurv.alpha.beta[2,2])
+    m3e <- qbeta(Ui.vec[i,3], msurv.alpha.beta[1,3], msurv.alpha.beta[2,3])
     as.msurv.e[[i + 1]] <- c(m1e, m2e, m3e)
+    } else{
+      as.msurv.e[[i + 1]] <- c(0, 0, 0)
+    }
     ct.vec <- popvect[[i]]
-    popvect[[i + 1]] <- one.iter.f.choosy(ct.vec, aii.vec = aii.vec , as.matur = as.matur,
+    popvect[[i + 1]] <- one.iter.f.choosy(ct.vec, aii.vec = aii.vec, as.matur = as.matur,
                                           as.msurv.e.set = as.msurv.e[[i + 1]], d = d, gnum.base = gnum.base,
-                                          b = b, cc = cc, s = s, mc = mc, density.type, mate.lazy = mate.lazy)
+                                          b = b, cc = cc, s = s, mc = mc, mate.lazy = mate.lazy,
+                                          env.sigma = env.sigma, Xi.is = Xi.vec[i,],  gam.density.type = gam.density.type,
+                                          is.density.type = is.density.type,
+                                          rj = rj,
+                                          yj = yj, env.immature.survival = env.immature.survival)
   }
 
 
   onedf <- form.autopop(popvect, generations)
+  survival_rates <- as.data.frame(do.call(rbind, popvect))[1:(generations+1),13:15]
+  for(i in 1:generations){
+    onedf$i1e[i] <- survival_rates$V13[i]
+    onedf$i2e[i] <- survival_rates$V14[i]
+    onedf$i3e[i] <- survival_rates$V15[i]
+    onedf$m1e[i] <- as.msurv.e[[i]][1]
+    onedf$m2e[i] <- as.msurv.e[[i]][2]
+    onedf$m3e[i] <- as.msurv.e[[i]][3]
+  }
+
   return(onedf)
 }

R/mu.options.R

@@ -0,0 +1,30 @@
+#' Max mean allowed
+#'
+#' Given an env.ci value, we identified allowed means between 0 and 1.
+#'
+#' @param env.ci  Proportion of environmental variance used to define mature survival rate per generation
+#'   with a beta distribution. This number must be in between 0 and 1, but cannot be equal to 0 or 1.
+#'
+#' @return All values allowed and not allowed.
+#'
+
+mu.options <- function(env.ci){
+  p.mu <- seq(0.001,1, by = 0.0001)
+  notallowed <- c()
+  allowed <- c()
+  for(i in 1:length(p.mu)){
+
+    mu <- p.mu[i]
+    var <- (env.ci*(mu*(1-mu)))
+    less <- (mu*(mu-1))
+
+    if((var > less) == FALSE){
+      notallowed <- c(notallowed, mu)
+    }else{
+      allowed <- c(allowed, mu)
+    }
+  }
+ out <-  list(not_allowed = notallowed, allowed = allowed)
+ return(out)
+}
+

R/one.iter.f.choosy.R

@@ -9,50 +9,65 @@
 #'
 #' @param ct.vec Population composition at time t.
 #' Indicates the sum of each type at time t, `ct.vec = c(2x_immature, 3x_immature, 4x_immature, 2x_mature, 3x_mature, 4x_mature)`
-#' @param aii.vec The survival probability of an immature individual for each cytotype. Must be a list of three integers between 0 and 1. For example, `aii.vec = c(0.5, 0.3, 0.5)`.
-#' @param as.matur The probability of maturation from an immature stage to the mature stage for each cytotype. Must be a list of three integers between 0 and 1. For example, `as.matur = c(0.5, 0.3, 0.5)`.
-#' @param as.msurv.e.set The survival probability of a mature individual for each cytotype. Must be a list of three integers between 0 and 1. For example, `as.msurv.e.set = c(0.5, 0.3, 0.5)`.
+#' @param aii.vec The survival probability of an immature individual for each cytotype. Must be a list of three integers between 0 and 1. For example, `aii.vec = c(0.0005, 0.005, 0.0005)`.
+#' @param as.matur The probability of maturation from an immature stage to the mature stage for each cytotype. Must be a list of three integers between 0 and 1. For example, `as.matur = c(0.60, 0.06, 0.60)`.
+#' @param as.msurv.e.set The survival probability of a mature individual for each cytotype. Must be a list of three integers between 0 and 1. These are defined based on user supplied
+#'  mean probability of mature survival (`as.msurv`), proportion of environmental variance (`env.ci`),
 #'  This list is defined within the `gen.iter.f.choosy()` function.
 #' @param d Strength of density dependency on gamete production for each cytotype. Must be a list of three integers between 0 and 1. For example, `d = c(0.001, 0.009, 0.001)`.
 #' @param gnum.base  Mean number of gametes per individual per cytotype. Must be a list of three numeric values. For example, `gnum.base = c(100, 100, 100)`.
 #' @param mate.lazy Default = FALSE, this prevents selfing from occurring during outcrossing. However, this increases the computational time by 31x!
-#' @inheritParams cytotype_repro_mate
-#' @param density.type Default = "all", this sets the density at time t as all individuals at time t.
+#' @param gam.density.type Default = "all", this sets the density dependence type for number of gametes produced at time t as all individuals at time t.
+#' Alternatively, "like-cytotype" sets the density at time t for each cytotype based on only the total immature and mature individuals
+#' of that cytotype at time t.
+#' @param rj List of two indicates the stage dependent density dependent impact of immature and mature individuals on  number of gametes produced. Default is c(1,1).
+#' @param is.density.type Default = "all", this sets the density dependence type for immature surivival at time t as all individuals at time t.
 #' Alternatively, "like-cytotype" sets the density at time t for each cytotype based on only the total immature and mature individuals
 #' of that cytotype at time t.
-#' @returns List of 9, with 1:6 representing the number of individuals of each cytotype
-#'   and at both stages at time t + 1. Items 7:9 are the number of gametes sampled for
-#'   2x, 3x, and 4x individuals at time t.
-#
+#' @param yj List of two indicates the stage dependent density dependent impact of immature and mature individuals on the probability of immature survival. Default is c(1,1).
+#' @param env.sigma Sigma value to define environmental variance corresponding to immature survival rate per generation.
+#' @param Xi.is Sample from a multivariate normal distribution. Allows immature and mature survival probabilities be correlated,
+#' as well as allowing for a correlation among cytotypes based on a supplied `rho` value.
+#' @param env.immature.survival Default = FALSE. When FALSE, the mean immature survival rate is equal to the exponential of the inverse
+#' immature survival rate (`aii`) times the sum of immature individuals and mature individuals of the cytotype indicated by `is.density.type`,
+#' which may be modified by `yj`. When this variable equals TRUE, the previous mean, or the immature determinate survival rate, is transformed
+#' by log(mu/(1-mu)) and the Johnson SB distribution is sampled given the `env.sigma` and `Xi.is` (defined in `gen.iter.f.choosy()`).
+#' @inheritParams cytotype_repro_mate
+#' @returns List of 15, with 1:6 representing the number of individuals of each cytotype
+#'   and at both stages at time t + 1. Items 7:9 are the number of offspring generated for
+#'   2x, 3x, and 4x individuals at time t. Items 10:12 are the number of gametes generated from
+#'   2x, 3x, and 4x individuals at time t. Items 12:15 are the immature survival probabilities
+#'   at time t.
+#'
 #'
 #' @importFrom stats na.omit rbinom
 #'
 #'
 #'
 
 # One iteration function
-one.iter.f.choosy <- function(ct.vec, aii.vec,
-                              as.matur, as.msurv.e.set,
-                              d, gnum.base,
-                              s, b,
-                              cc, mc,
-                              density.type = "all", mate.lazy = FALSE){
+one.iter.f.choosy <- function(ct.vec, aii.vec, as.matur, as.msurv.e.set,
+                              d, gnum.base, s, b,
+                              cc, mc, gam.density.type, is.density.type,
+                              mate.lazy = FALSE, rj = c(1,1), yj = c(1,1), env.sigma,
+                              Xi.is, env.immature.survival){
 
     # Empty vector for next generation
     ctp1.vec <- rep(0, 6)
 
     # Calculate current generations gnum.vec
     ## gnum.vec = gnum.base times the exponential of the product of strength of density dependency
     ## times the sum of all individuals at time t
-    if(density.type == "all"){
-      gnum.vec <- c((gnum.base[1]*exp(-d[1]*(sum(ct.vec[1:6])))),
-                    (gnum.base[2]*exp(-d[2]*(sum(ct.vec[1:6])))),
-                    (gnum.base[3]*exp(-d[3]*(sum(ct.vec[1:6])))))
+    if(gam.density.type == "all"){
+      gam.density <- ((rj[1]*sum(ct.vec[1:3]))+(rj[2]*sum(ct.vec[4:6])))
+      gnum.vec <- c((gnum.base[1]*exp(-d[1]*gam.density)),
+                    (gnum.base[2]*exp(-d[2]*gam.density)),
+                    (gnum.base[3]*exp(-d[3]*gam.density)))
       gnum.vec[is.na(gnum.vec)] <- 0
-    } else if(density.type == "like-cytotype"){
-      gnum.vec <- c((gnum.base[1]*exp(-d[1]*(sum(c(ct.vec[1], ct.vec[4]))))),
-                    (gnum.base[2]*exp(-d[2]*(sum(c(ct.vec[2], ct.vec[5] ))))),
-                    (gnum.base[3]*exp(-d[3]*(sum(c(ct.vec[3], ct.vec[6]))))))
+    } else if(gam.density.type == "like-cytotype"){
+      gnum.vec <- c((gnum.base[1]*exp(-d[1]*((rj[1]*ct.vec[1]) + (rj[2]*ct.vec[4])))),
+                    (gnum.base[2]*exp(-d[2]*((rj[1]*ct.vec[2]) + (rj[2]*ct.vec[5])))),
+                    (gnum.base[3]*exp(-d[3]*((rj[1]*ct.vec[3]) + (rj[2]*ct.vec[6])))))
       gnum.vec[is.na(gnum.vec)] <- 0
     } else print("Density type unavailable. This option only allows density.type to be set as all or like-cytotype.")
 
@@ -68,10 +83,37 @@ one.iter.f.choosy <- function(ct.vec, aii.vec,
     trip.gamps <- off[2]; if(is.na(trip.gamps)){trip.gamps <- 0}
     tetr.gamps <- off[3]; if(is.na(tetr.gamps)){tetr.gamps <- 0}
 
+    # Setup Immature Survival Probability
+    if(is.density.type == "all"){
+      is.density <- ((yj[1]*sum(ct.vec[1:3]))+(yj[2]*sum(ct.vec[4:6])))
+      immature.det.surv <- c((exp(-aii.vec[1]*is.density)),
+                                (exp(-aii.vec[2]*is.density)),
+                                (exp(-aii.vec[3]*is.density)))
+    } else if(is.density.type == "like-cytotype"){
+      immature.det.surv <- c((exp(-aii.vec[1]*((yj[1]*ct.vec[1]) + (yj[2]*ct.vec[4])))),
+                                (exp(-aii.vec[2]*((yj[1]*ct.vec[2]) + (yj[2]*ct.vec[5])))),
+                                (exp(-aii.vec[3]*((yj[1]*ct.vec[3]) + (yj[2]*ct.vec[6])))))
+    } else print("Density type unavailable. This option only allows density.type to be set as all or like-cytotype.")
+
+    immature.det.surv[is.na(immature.det.surv)] <- 0
+    if(env.immature.survival == TRUE){
+    immature.mu<- c()
+    as.isurv <- c()
+    for(z in 1:3){
+      immature.mu[z] <- log(immature.det.surv[z]/(1-immature.det.surv[z]))
+      as.isurv[z] <- (exp(immature.mu[z] + (env.sigma*Xi.is[z]))/(1 + exp(immature.mu[z] + (env.sigma*Xi.is[z]))))
+    }
+    as.isurv[is.na(as.isurv)] <- 0
+    } else{
+      as.isurv <- immature.det.surv
+    }
+
+
+
     # Diploids - two steps: 1. Maturation+Survival and 2. Gamete production (defined above)
     ## Step 1:
     ### Apply survival and maturation
-    Ndip.imm <-  rbinom(n = 1, size = ct.vec[1], prob = exp(-aii.vec[1]*sum(ct.vec[1], ct.vec[4]))) #density
+    Ndip.imm <-  rbinom(n = 1, size = ct.vec[1], prob = as.isurv[1]) #density
     Ndip.mat <-  rbinom(n = 1, size = ct.vec[1], prob = as.matur[1])
     Ndip.surv <- rbinom(n = 1, size = ct.vec[4], prob = as.msurv.e.set[1])
 
@@ -84,7 +126,7 @@ one.iter.f.choosy <- function(ct.vec, aii.vec,
     # Triploids
     ## Step 1:
     ### Apply survival and maturation
-    Ntrip.imm <-  rbinom(n = 1, size = ct.vec[2], prob = exp(-aii.vec[2]*sum(ct.vec[2], ct.vec[5])))
+    Ntrip.imm <-  rbinom(n = 1, size = ct.vec[2], prob = as.isurv[2])
     Ntrip.mat <-  rbinom(n = 1, size = ct.vec[2], prob = as.matur[2])
     Ntrip.surv <- rbinom(n = 1, size = ct.vec[5], prob = as.msurv.e.set[2])
 
@@ -97,7 +139,7 @@ one.iter.f.choosy <- function(ct.vec, aii.vec,
     # Tetraploids
     ## Step 1:
     ### Apply survival and maturation
-    Ntetr.imm <-  rbinom(n = 1, size = ct.vec[3], prob = exp(-aii.vec[3]*sum(ct.vec[3], ct.vec[6])))
+    Ntetr.imm <-  rbinom(n = 1, size = ct.vec[3], prob = as.isurv[3])
     Ntetr.mat <-  rbinom(n = 1, size = ct.vec[3], prob = as.matur[3])
     Ntetr.surv <- rbinom(n = 1, size = ct.vec[6], prob = as.msurv.e.set[3])
 
@@ -106,8 +148,9 @@ one.iter.f.choosy <- function(ct.vec, aii.vec,
     ctp1.vec[6] <-  Ntetr.mat + Ntetr.surv
     ctp1.vec[3] <-  tetr.gamps + Ntetr.imm
 
+    ctp1.vec[is.na(ctp1.vec)] <- 0
     # Combined and Return
-    ctp1.vec <- c(ctp1.vec, off[1:6])
+    ctp1.vec <- c(ctp1.vec, off[1:6], as.isurv)
     return(ctp1.vec)
 }
 

R/surv.shape.R

@@ -0,0 +1,24 @@
+#' Survival shape -  Calculates the alpha and beta value for survival beta distributions
+#'
+#' @description Calculates the shape parameters needed to sample a beta distribution
+#' for both immature survival and mature survival
+#'
+#'
+#' @param raw.means Mean probability of mature survival or immature survival. Must be a list of three integer between 0 and 1.
+#' @param  env.ci Proportion of environmental variance used to define mature survival rate per generation.
+#'   Must be a single integer greater than or equal to 0 and less than 1.
+#'
+#' @returns Matrix of alpha and beta parameters where matrix[1,i] is the alpha value for the ith mean, and matrix[2,i]
+#' is the beta value for the ith mean.
+#'
+#' @importFrom stats na.omit rbeta
+#'
+
+
+surv.shape <- function(env.ci, raw.means){
+  var <- var.option(env.ci = env.ci, mu = raw.means)
+  alpha.beta <- sapply(1:length(raw.means), function(item) alphabeta.calc(raw.means[item], var[item]))
+  return(alpha.beta)
+}
+
+

R/var.option.R

@@ -14,5 +14,6 @@
 
 var.option <- function(env.ci, mu){
   var <- (env.ci*mu*(1-mu))
+  #var <- max(var, .Machine$double.xmin)
   return(var)
 }