Added Model M4 (DD & Demo Stoch).

Added piece "%>% sapply(min, 0.999)" in M3 and M4 to avoid survival parameter > 1.
Also added this in 'init_calib' function.

All models tested and work.

.Rbuildignore

@@ -6,3 +6,4 @@
 ^\./tests/test_build_Leslie\.R$
 ^junk\.R$
 ^run_shiny\.R$
+^draws_histog\.R$

NAMESPACE

@@ -3,10 +3,12 @@
 export(M1_noDD_noDemoStoch)
 export(M2_noDD_WithDemoStoch)
 export(M3_WithDD_noDemoStoch)
+export(M4_WithDD_WithDemoStoch)
 export(build_Leslie)
 export(calibrate_params)
 export(elicitation)
 export(get_metrics)
+export(infer_DD)
 export(init_calib)
 export(make_transparent)
 export(match_lam_delta)

R/calibrate_params.R

@@ -90,7 +90,7 @@ init_calib  <- function(s, f, lam0){
   nac = length(s)
 
   inits_vr <- c(s,f) + d
-  inits_vr <- c(tail(inits_vr, nac), head(inits_vr, nac))
+  inits_vr <- c(tail(inits_vr, nac), head(inits_vr, nac) %>% sapply(min, 0.999))
   inits <- inits_vr[inits_vr != 0]
   return(inits)
 } # End function

R/infer_DD.R

+#' Infer Densisty-dependence Parameters
+#'
+#' infer_DD calculates the missing value of a parameter (either K or rMAX, whichever is NULL) in a density-dependence relationship.
+#'
+#' @param lambda_MAX a strictly positive number. Maximum (theoretical) population growth rate.
+#' @param K a strictly positive number. Carrying capacity (= maximum size that the population can ever reach)
+#' @param theta a strictly positive number. Parameter defining the shape of the density-dependence relationship.
+#' The relationship is defined as : lambda <- 1 + rMAX*(1-(N/K)^theta)
+#' where rMAX = lambda_MAX - 1
+#' @param pop_size_current a strictly positive number. Current population size.
+#' @param pop_growth_current a strictly positive number. Current population growth rate.
+#'
+#' @return a list of length 2, containing K and lambda_MAX.
+#' @export
+#'
+#' @examples
+#'## Infer K
+#'infer_DD(lambda_MAX = 1.15, K = NULL, theta = 1, pop_size_current = 200, pop_growth_current = 1.08)
+#'## Infer lambda_MAX
+#'infer_DD(lambda_MAX = NULL, K = 2000, theta = 1, pop_size_current = 200, pop_growth_current = 1.08)
+#'
+#'
+infer_DD <- function(lambda_MAX = NULL, K = NULL, theta = 1, pop_size_current, pop_growth_current){
+
+  if(!is.null(lambda_MAX)){
+    # Infer K
+    rMAX = lambda_MAX - 1
+    r_a = pop_growth_current - 1
+    N_a = pop_size_current
+
+    K <- (N_a/((1-(r_a/rMAX))^(1/theta))) %>% round
+
+  }else{
+
+    if(!is.null(K)){
+      # Infer rMAX
+      N_a = pop_size_current
+      r_a = pop_growth_current - 1
+
+      rMAX <- r_a/((1-(N_a/K))^theta)
+      lambda_MAX = rMAX + 1
+
+    }else{
+
+      stop("Either lambda_MAX or K must be provided")
+
+    } # end if 2
+  } # end if 1
+
+  return(list(lambda_MAX = lambda_MAX, K = K))
+
+}  # END FUNCTION
+################################################################################
+
+

R/models.R

@@ -65,6 +65,7 @@ M1_noDD_noDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 #'
 M2_noDD_WithDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 
+  # Number of age classes
   nac = length(s)
 
   # Survivors (to "natural mortality" (s) and Wind Turbine Fatalities (1-h))
@@ -90,9 +91,9 @@ M2_noDD_WithDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 
 
 ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-##  Model 3: WITH Density-Dependence (noDD), No Demographic Stochasticity     ~~
+##  Model 3: WITH Density-Dependence (DD), No Demographic Stochasticity       ~~
 ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#' Demographic model without Density-Dependence and without Demographic Stochasticity
+#' Demographic model with Density-Dependence, but without Demographic Stochasticity
 #' In addition to natural survivals and fecundities, this demographic model includes
 #' a specific harvest / fatality parameter.
 #' NOTE : This is a POST-BREEDING demographic model.
@@ -114,7 +115,8 @@ M2_noDD_WithDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 #' f <- c(0, 0, 0.05, 0.55)
 #' N1 <- c(50, 60, 75, 100)
 #' h <- 0.05
-#' M1_noDD_noDemoStoch(N1, s, f, h)
+#' DD_params <- list(rMAX = 0.15, K = 1200, theta = 1)
+#' M3_WithDD_noDemoStoch(N1, s, f, h,DD_params = DD_params)
 #'
 M3_WithDD_noDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 
@@ -132,7 +134,7 @@ M3_WithDD_noDemoStoch <- function(N1, s, f, h, DD_params = NULL){
   d <- match_lam_delta(diff_rel_lam = diff_rel_lam, s=s, f=f)
   vr_Nt <- c(s,f) + d
 
-  s_Nt <- head(vr_Nt, length(s))
+  s_Nt <- head(vr_Nt, length(s)) %>% sapply(min, 0.999)
   f_Nt <- tail(vr_Nt, length(f))
 
   # Build the LESLIE matrix
@@ -145,3 +147,64 @@ M3_WithDD_noDemoStoch <- function(N1, s, f, h, DD_params = NULL){
 
 } # END FUNCTION
 ################################################################################
+
+
+
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+##  Model 4: With Density-Dependence (DD), WITH Demographic Stochasticity     ~~
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#' Demographic model with Density-Dependence and INCLUDING Demographic Stochasticity
+#' In addition to natural survivals and fecundities, this demographic model includes
+#' a specific harvest / fatality parameter.
+#' NOTE : This is a POST-BREEDING demographic model.
+#'
+#' @param N1 a vector of population sizes for each age class at time t1
+#' @param s a vector of survival probabilities for each age class
+#' @param f a vector of fecundity values for each age class
+#' @param h a number. The harvest or fatality rate
+#' @param DD_params density-dependence parameters. Not used in this model.
+#'
+#' @return a vector of population sizes for each age class at time t2
+#' @export
+#'
+#' @examples
+#' s <- c(0.5, 0.7, 0.8, 0.95)
+#' f <- c(0, 0, 0.05, 0.55)
+#' N1 <- c(50, 60, 75, 100)
+#' h <- 0.05
+#' DD_params <- list(rMAX = 0.15, K = 1200, theta = 1)
+#' M4_WithDD_WithDemoStoch(N1, s, f, h, DD_params = DD_params)
+#'
+M4_WithDD_WithDemoStoch <- function(N1, s, f, h, DD_params){
+
+  # Extract DD parameters from list
+  rMAX = DD_params$rMAX
+  K = DD_params$K
+  theta = DD_params$theta
+
+  # Apply density dependence effect
+  lam_Nt <- 1 + rMAX*(1-(sum(N1)/K)^theta)
+
+  # Calibrate vital rates to match lam_Nt
+  A <- build_Leslie(s = s, f = f)
+  diff_rel_lam <- (lam_Nt - lambda(A))/lambda(A)
+  d <- match_lam_delta(diff_rel_lam = diff_rel_lam, s=s, f=f)
+  vr_Nt <- c(s,f) + d
+
+  s_Nt <- head(vr_Nt, length(s)) %>% sapply(min, 0.999)
+  f_Nt <- tail(vr_Nt, length(f))
+
+  # Number of age classes
+  nac = length(s)
+
+  # Survivors (to "natural mortality" (s) and Wind Turbine Fatalities (1-h))
+  S <- rbinom(nac, N1, (1-h)*s_Nt)
+  N2 <- c(rep(0, nac-1), tail(S,1)) + c(0, head(S,-1))
+
+  # Births
+  N2[1] <- sum(rpois(nac, f_Nt*N2))
+
+  return(N2)
+
+} # END FUNCTION
+################################################################################

devtools_history.R

@@ -41,3 +41,8 @@ usethis::use_build_ignore("run_shiny.R")
 ## Put it on Git
 library(usethis)
 usethis::use_git()
+
+
+# Ignore file "run_shiny.R"
+usethis::use_build_ignore("draws_histog.R")
+

draws_histog.R

+
+##==============================================================================
+##                         Check lambda draws                                 ==
+##==============================================================================
+draws_histog <- function(draws, mu, se){
+
+  # Plot histogram
+  h <- hist(draws, breaks = length(draws)/10, border = 0)
+
+  # Theoretical Normal Curve
+  par(new=T)
+  curve(dnorm(x, mean=mu, sd=se), add=FALSE, lwd=3, col="darkblue",
+        xlim = c(min(draws), max(draws)), axes = FALSE, xlab = "", ylab = "")
+
+} # End function
+################################################################################
+
+
+# draws_histog(draws = run0$lambdas, mu = pop_growth_mean, se = pop_growth_se)

junk.R

@@ -7,6 +7,53 @@ s <- c(0.5, 0.7, 0.8, 0.95)
 f <- c(0, 0, 0.05, 0.55)
 N1 <- pop_vector(pop_size = 200, pop_size_type = "Npair", s, f)
 h <- 0.05
+DD_params <- list(rMAX = 0.15, K = 1200, theta = 1)
+
+
+
+# Infer K
+infer_DD(lambda_MAX = 1.15, K = NULL, theta = 1, pop_size_current = 200, pop_growth_current = 1.08)
+
+lambda_MAX = 1.15
+K = NULL
+theta = 1
+pop_size_current = 200
+pop_growth_current = 1.08
+
+
+# Infer lambda_MAX
+infer_DD(lambda_MAX = NULL, K = 2000, theta = 1, pop_size_current = 200, pop_growth_current = 1.08)
+
+lambda_MAX = NULL
+K = 429
+theta = 1
+pop_size_current = 200
+pop_growth_current = 1.08
+
+# Get an error
+infer_DD(lambda_MAX = NULL, K = NULL, theta = 1, pop_size_current = 200, pop_growth_current = 1.08)
+
+
+# Infer K
+rMAX = 0.15
+lam_a = 1.08
+theta = 1
+N_a = 200
+r_a = lam_a - 1 ; r_a
+
+K <- (N_a/((1-(r_a/rMAX))^(1/theta))) %>% round
+K
+
+# Infer rMAX
+K = 2000
+lam_a = 1.08
+theta = 1
+N_a = 500
+r_a = lam_a - 1 ; r_a
+
+rMAX <- r_a/((1-(N_a/K))^theta)
+rMAX
+
 
 
 
@@ -41,12 +88,74 @@ lambda(A_Nt)
 
 
 
+A <- build_Leslie(s = s, f = f)
+diff_rel_lam <- (1.5 - lambda(A))/lambda(A)
+diff_rel_lam
+d <- match_lam_delta(diff_rel_lam = diff_rel_lam, s=s, f=f)
+c(s,f) + d
+
+
+match_lam_delta(diff_rel_lam = 0.5, s=s, f=f)
+
+
 
 
 
 
 
 
+M <- fatalities
+N0 <- intial_pop_vector
+cv_env <- coeff_var_environ
+
+# Number of years
+nyr <- time_horzion
+
+# Number of age classes
+nac <- length(s)
+
+# Number of fatalities scenario
+nsc <- length(fatalities)
+
+# Initiate Pop Size (output) Array
+N <- array(NA, dim = c(nac, nyr, nsc), dimnames = list(paste0("age", 1:nac),
+                                                       paste0("year", 1:nyr),
+                                                       paste0("scenario", 1:nsc)
+))
+N[,1,] <-  N0
+
+## Loops over time (years)
+for(t in 2:nyr){
+
+  # Environmental Stochasticity
+  ss = 1 - rnorm(nac, mean = qlogis(1-s), sd = cv_env/(s)) %>% plogis        ## sample thru the mortality rate : 1 - s
+  ff = rnorm(nac, mean = log(f), sd = cv_env) %>% exp
+
+  # Fatalities : constant number (M) or constant rate (h)
+  if(fatal_constant == "M"){
+    h <- M/apply(N[,t-1,], 2, sum)
+  } else {
+    h <- M/apply(N[,1,], 2, sum)
+  }
+
+  # Sample a seed for RNG
+  seed <- runif(1, 0, 1e6)
+
+  ## Projection : apply the LESLIE matrix calculation forward
+  # Scenario 0
+  for(j in 1:nsc){
+    set.seed(seed)
+    N[,t,j] <- model_demo(N1 = N[,t-1,j], s = ss, f = ff, h = h[j], DD_params = DD_params)
+  } # j
+
+} # t
+
+N
+
+
+
+
+
 
 
 
@@ -74,6 +183,37 @@ coeff_var_environ = coeff_var_environ; fatal_constant = fatal_constant
 
 ##
 
+
+# Environmental Stochasticity
+ss = 1 - rnorm(nac, mean = qlogis(1-s), sd = cv_env/(s)) %>% plogis        ## sample thru the mortality rate : 1 - s
+ff = rnorm(nac, mean = log(f), sd = cv_env) %>% exp
+ss ; ff
+
+# Fatalities : constant number (M) or constant rate (h)
+if(fatal_constant == "M"){
+  h <- M/apply(N[,t-1,], 2, sum)
+} else {
+  h <- M/apply(N[,1,], 2, sum)
+}
+h
+
+# Sample a seed for RNG
+seed <- runif(1, 0, 1e6)
+
+## Projection : apply the LESLIE matrix calculation forward
+# Scenario 0
+#for(j in 1:nsc){
+j = 1
+  set.seed(seed)
+  N[,t,j] <- model_demo(N1 = N[,t-1,j], s = ss, f = ff, h = h[j], DD_params = DD_params)
+#} # j
+N
+
+model_demo(N1 = N[,t-1,j], s = ss, f = ff, h = h[j], DD_params = DD_params)
+
+
+t=2
+
 N1 = N[,t-1,j]
 s = ss
 f = ff

man/M3_WithDD_noDemoStoch.Rd

@@ -2,7 +2,7 @@
 % Please edit documentation in R/models.R
 \name{M3_WithDD_noDemoStoch}
 \alias{M3_WithDD_noDemoStoch}
-\title{Demographic model without Density-Dependence and without Demographic Stochasticity
+\title{Demographic model with Density-Dependence, but without Demographic Stochasticity
 In addition to natural survivals and fecundities, this demographic model includes
 a specific harvest / fatality parameter.
 NOTE : This is a POST-BREEDING demographic model.}
@@ -27,7 +27,7 @@ theta (shape of DD relationshp)}
 a vector of population sizes for each age class at time t2
 }
 \description{
-Demographic model without Density-Dependence and without Demographic Stochasticity
+Demographic model with Density-Dependence, but without Demographic Stochasticity
 In addition to natural survivals and fecundities, this demographic model includes
 a specific harvest / fatality parameter.
 NOTE : This is a POST-BREEDING demographic model.
@@ -37,6 +37,7 @@ s <- c(0.5, 0.7, 0.8, 0.95)
 f <- c(0, 0, 0.05, 0.55)
 N1 <- c(50, 60, 75, 100)
 h <- 0.05
-M1_noDD_noDemoStoch(N1, s, f, h)
+DD_params <- list(rMAX = 0.15, K = 1200, theta = 1)
+M3_WithDD_noDemoStoch(N1, s, f, h,DD_params = DD_params)
 
 }

man/M4_WithDD_WithDemoStoch.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/models.R
+\name{M4_WithDD_WithDemoStoch}
+\alias{M4_WithDD_WithDemoStoch}
+\title{Demographic model with Density-Dependence and INCLUDING Demographic Stochasticity
+In addition to natural survivals and fecundities, this demographic model includes
+a specific harvest / fatality parameter.
+NOTE : This is a POST-BREEDING demographic model.}
+\usage{
+M4_WithDD_WithDemoStoch(N1, s, f, h, DD_params)
+}
+\arguments{
+\item{N1}{a vector of population sizes for each age class at time t1}
+
+\item{s}{a vector of survival probabilities for each age class}
+
+\item{f}{a vector of fecundity values for each age class}
+
+\item{h}{a number. The harvest or fatality rate}
+
+\item{DD_params}{density-dependence parameters. Not used in this model.}
+}
+\value{
+a vector of population sizes for each age class at time t2
+}
+\description{
+Demographic model with Density-Dependence and INCLUDING Demographic Stochasticity
+In addition to natural survivals and fecundities, this demographic model includes
+a specific harvest / fatality parameter.
+NOTE : This is a POST-BREEDING demographic model.
+}
+\examples{
+s <- c(0.5, 0.7, 0.8, 0.95)
+f <- c(0, 0, 0.05, 0.55)
+N1 <- c(50, 60, 75, 100)
+h <- 0.05
+DD_params <- list(rMAX = 0.15, K = 1200, theta = 1)
+M4_WithDD_WithDemoStoch(N1, s, f, h, DD_params = DD_params)
+
+}

man/infer_DD.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/infer_DD.R
+\name{infer_DD}
+\alias{infer_DD}
+\title{Infer Densisty-dependence Parameters}
+\usage{
+infer_DD(
+  lambda_MAX = NULL,
+  K = NULL,
+  theta = 1,
+  pop_size_current,
+  pop_growth_current
+)
+}
+\arguments{
+\item{lambda_MAX}{a strictly positive number. Maximum (theoretical) population growth rate.}
+
+\item{K}{a strictly positive number. Carrying capacity (= maximum size that the population can ever reach)}
+
+\item{theta}{a strictly positive number. Parameter defining the shape of the density-dependence relationship.
+The relationship is defined as : lambda <- 1 + rMAX*(1-(N/K)^theta)
+where rMAX = lambda_MAX - 1}
+
+\item{pop_size_current}{a strictly positive number. Current population size.}
+
+\item{pop_growth_current}{a strictly positive number. Current population growth rate.}
+}
+\value{
+a list of length 2, containing K and lambda_MAX.
+}
+\description{
+infer_DD calculates the missing value of a parameter (either K or rMAX, whichever is NULL) in a density-dependence relationship.
+}
+\examples{
+## Infer K
+infer_DD(lambda_MAX = 1.15, K = NULL, theta = 1, pop_size_current = 200, pop_growth_current = 1.08)
+## Infer lambda_MAX
+infer_DD(lambda_MAX = NULL, K = 2000, theta = 1, pop_size_current = 200, pop_growth_current = 1.08)
+
+
+}

run_analysis.R

@@ -19,7 +19,7 @@ pop_growth_se = 0.05
 survivals <- c(0.5, 0.7, 0.8, 0.95)
 fecundities <- c(0, 0, 0.05, 0.55)
 
-model_demo = M3_WithDD_noDemoStoch #M2_noDD_WithDemoStoch #
+model_demo = M4_WithDD_WithDemoStoch #M3_WithDD_noDemoStoch # M2_noDD_WithDemoStoch #M1_noDD_noDemoStoch #
 time_horzion = 30
 coeff_var_environ = 0.10
 fatal_constant = "h"
@@ -51,42 +51,19 @@ run0 <- run_simul(nsim, cumuated_impacts,
                   model_demo, time_horzion, coeff_var_environ, fatal_constant)
 
 
-# saved time (ratio): 493/12
 # save(run0, file = "./data/run0.rda")
 names(run0)
-
 run0$time_run
+# saved time (ratio): 493/12
 
 N <- run0$N
 out <- get_metrics(N)
 out[time_horzion,"avg",]
 
-# Impact par parc
-# for(j in 2:length())
-#j=4
-#out[time_horzion, -2, j] - out[time_horzion, -2, j-1]
-
-
 plot_traj(N, xlab = "Annee", ylab = "Taille de population (totale)")
 
 plot_impact(N, xlab = "Annee", ylab = "Taille de population (totale)")
 
 
-##==============================================================================
-##                         Check lambda draws                                 ==
-##==============================================================================
-draws_histog <- function(draws, mu, se){
-
-  # Plot histogram
-  h <- hist(draws, breaks = length(draws)/10, border = 0)
-
-  # Theoretical Normal Curve
-  par(new=T)
-  curve(dnorm(x, mean=mu, sd=se), add=FALSE, lwd=3, col="darkblue",
-        xlim = c(min(draws), max(draws)), axes = FALSE, xlab = "", ylab = "")
-
-} # End function
-################################################################################
-
-
+source("draws_histog.R")
 draws_histog(draws = run0$lambdas, mu = pop_growth_mean, se = pop_growth_se)