inst/doc/simulations.R

In alien: Estimate Invasive and Alien Species (IAS) Introduction Rates

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(alien)

## -----------------------------------------------------------------------------
list_tsl <- list(
  b0  = c(-1,      # mean introduction rate of 0.3 sp/y
          0,       # mean introduction rate of  1  sp/y
          1),      # mean introduction rate of 2.7 sp/y
  b1  = c(-0.02,   # annual decrease of  5%
          -0.01,   # annual decrease of 10%
          -0.005,  # annual decrease of 20%
          0,       # constant introduction rate
          0.005,   # annual increase of  5% 
          0.01,    # annual increase of 10% 
          0.02),   # annual increase of 20% 
  pi0 = c(-4,      # annual observation probability of 0.017,   ~55 years delay
          -2,      # annual observation probability of 0.119,   ~8 years delay
          0,       # annual observation probability of   0.5,   ~2 years delay
          2,       # annual observation probability of   0.88,  ~0.13 years delay
          4),      # annual observation probability of   0.98,  ~0.02 year delay
  pi1 = c( 0,      # annual increase of     0% in probability of detection
           0.5,   # annual increase of  0.02% in probability of detection
           1,  # annual increase of   0.2% in probability of detection
           2), # effect of population size on detection probability
  pi2 = c( 0,      # no effect of population size on detection probability
           1e-7,   # normal effect of population size on detection probability
           0.001), # effect of population size on detection probability
  length = seq(30,200,15)
)

df_tsl <- expand.grid(list_tsl)

est_names <- c("b0_est", "b1_est", "pi0_est", "pi1_est", "pi2_est")
ses_names <- c("b0_se", "b1_se", "pi0_se", "pi1_se", "pi2_se")

## ----subsetting data----------------------------------------------------------
set.seed(3333) # Set a seed for replication
df_tsl <- df_tsl |> dplyr::sample_n(1)

## ----theme setup, echo = FALSE------------------------------------------------
ggplot2::theme_set(
  ggplot2::theme_bw()+
    ggplot2::theme(axis.title = ggplot2::element_text(size = 20),
                   axis.text = ggplot2::element_text(size = 18),
                   panel.grid = ggplot2::element_blank())
)

## ----show lambda, fig.width = 8, fig.height= 4.5, fig.align='center'----------
params <- as.numeric(df_tsl[1:5]) # parameters
names(params) <- c("b0","b1","pi0","pi1", "pi2") # set names
time_span <- as.numeric(df_tsl[6]) # time series length.

sampling_lambda <- alien:::calculate_lambda(mu = ~ time,
                                            pi = ~ time,
                                            data = data.frame(time = seq_len(time_span)),
                                            beta = params[c("b0","b1")],
                                            gamma = params[c("pi0","pi1")],
                                            growth_param = params["pi2"],
                                            type = "exponential")
ggplot2::ggplot()+
  ggplot2::aes(x = seq_len(time_span), y = sampling_lambda)+
  ggplot2::geom_line(linewidth = 1.1) + 
  ggplot2::labs(x = "Time", y = expression("\U03BB"[t]))


## ----simulation---------------------------------------------------------------
list <- apply(df_tsl,
                MARGIN = 1, FUN = function(row) {


                  # Defining variables and parameters:
                  params <- as.numeric(row[1:5])
                  names(params) <- c("b0","b1","pi0","pi1", "pi2")
                  time_span <- row[6] # time series length.

                  # Create a faux sampling vector based on the time span to be used for our sampling_proxy
                  sampling_vector <- scale(1:time_span) |> as.numeric()

                  # Create lambda vector
                  sampling_lambda <- alien:::calculate_lambda(mu = ~ time,
                                            pi = ~ time,
                                            data = data.frame(time = seq_len(time_span)),
                                            beta = params[c("b0","b1")],
                                            gamma = params[c("pi0","pi1")],
                                            growth_param = params["pi2"],
                                            type = "exponential")
                  
                  number_of_repetitions <- 1
                  
                  result <- lapply(seq_len(number_of_repetitions), function(iteration) {

                    # Creating an introduction
                    sampling_data <- rpois(time_span, sampling_lambda)

                    # model fitting start
                    naive_model <- glm(formula = sampling_data ~ seq_len(time_span), family = "poisson")

                    guess_all <- c(b0 = 0, b1 = 0, pi0 = 0, pi1 = 0, pi2 = 0)

                    snc_model <- alien::snc(y = sampling_lambda, type = "exponential",
                                            mu = ~ time,
                                            pi = ~ time, growth = T,
                                            init = guess_all,
                                            data = data.frame(time = seq_len(time_span)),
                                            control = list(maxit = 10000))


                    constant_detection_model <- alien::snc(y = sampling_lambda, type = "exponential",
                                                           mu = ~ time,
                                                           pi = ~ 1, growth = F,
                                                           init = guess_all[1:3],
                                                           data = data.frame(time = seq_len(time_span)),
                                                           control = list(maxit = 10000))


                    sampling_proxy_model <- alien::snc(y = sampling_lambda, type = "exponential",
                                                       mu = ~ time,
                                                       pi = ~ sampling_vector, growth = F,
                                                       init = guess_all[1:4],
                                                       data = data.frame(time = seq_len(time_span), sampling_vector = sampling_vector),
                                                       control = list(maxit = 10000))

                    # model fitting end
                    
                    # Summarize results start
                    naive_par <-  coefficients(summary(naive_model))[,1] |>  dplyr::bind_rows() |> `colnames<-`(est_names[1:2])
                    naive_se <-  coefficients(summary(naive_model))[,2] |> dplyr::bind_rows() |> `colnames<-`(ses_names[1:2])
                    naive_convergence <- tibble::tibble(converged = naive_model$converged)

                    naive_row <- dplyr::bind_cols(naive_par, naive_se, naive_convergence) |>
                      tibble::add_column(method = "naive")

                    list_of_rows <- lapply(list(snc_model, constant_detection_model, sampling_proxy_model), function(est) {

                      method <- dplyr::case_when(
                        identical(est, snc_model) ~ "snc",
                        identical(est, constant_detection_model) ~ "snc_3p",
                        identical(est, sampling_proxy_model) ~ "snc_vec"
                      )

                      if (any(class(est) == "try-error")){
                        est_row <- rep(NA,11)
                        names(est_row) <- c("b0_est", "b1_est", "pi0_est",
                                            "pi1_est", "pi2_est", "b0_se", "b1_se",
                                            "pi0_se", "pi1_se", "pi2_se", "converged")
                        est_row <- dplyr::bind_rows(est_row) |> tibble::add_column(method = method)
                      } else {

                        par <- est$coefficients$Estimate |> `names<-`(est_names[1:nrow(est$coefficients)]) |> dplyr::bind_rows()
                        ses <- est$coefficients$`Std.Err` |> `names<-`(ses_names[1:nrow(est$coefficients)])  |> dplyr::bind_rows()
                        convergence <- tibble::tibble(converged = est$convergence == 0)
                        est_row <- dplyr::bind_cols(par, ses, convergence) |> tibble::add_column(method = method)
                      }
                    }
                    )

                    params <- dplyr::bind_rows(naive_row,
                                        dplyr::bind_rows(list_of_rows)) |>
                      tibble::add_column(iteration = iteration)

                    # summarize results end
                    return(list(params = params,
                                full = list(simulation = sampling_data,
                                            estimates = snc_model,
                                            estimates_3p = constant_detection_model,
                                            estimates_vec = sampling_proxy_model,
                                            perfect_detection = coefficients(summary(naive_model))
                                )
                    )
                    )
                  })

                  simulation_estimates <- dplyr::bind_rows(purrr::map(result, 1))
                  full <- purrr::map(result, 2)

                  return(list(sim_params = row,
                              estimates = simulation_estimates,
                              full = full))
                })


## ----view results-------------------------------------------------------------

purrr::map(list, 2) |> dplyr::bind_rows() # for the summary of each iteration


## ----fig.width = 7------------------------------------------------------------
t <- seq(0, 4*pi, length.out = 100)
a <- 3
b <- 0.6
set.seed(100) # This is done to maintain reproducability within time series lengths
c.norm <- rnorm(100)
amp <- 2

sampling_vector <- 3 * a*sin(b*t)+c.norm*amp + 0.3 * 1:100

sampling_vector <- scale(sampling_vector) |> as.numeric()


ggplot2::ggplot()+
  ggplot2::aes(x = 1:100, y = sampling_vector)+
  ggplot2::geom_line(linewidth = 1.1) + 
  ggplot2::labs(x = "Time", y = "Sampling proxy value")

## ----lambda_with_proxy, fig.width = 8, fig.height= 4.5, fig.align='center'----
params["pi1"] <- 0.2
time_span <- 100
lambda_with_proxy <- alien:::calculate_lambda(mu = ~ time,
                                            pi = ~ sampling_vector,
                                            data = data.frame(time = seq_len(time_span),
                                                              sampling_vector = sampling_vector),
                                            beta = params[c("b0","b1")],
                                            gamma = params[c("pi0","pi1")],
                                            growth_param = params["pi2"],
                                            type = "exponential")

ggplot2::ggplot()+
  ggplot2::aes(x = 1:100, y = lambda_with_proxy)+
  ggplot2::geom_line(linewidth = 1.1) + 
  ggplot2::labs(x = "Time", y = expression("\U03BB"[t]))