simulating a fake dataset

simulate_one_cardoso.R

simulate_one_cardoso <- function(.cardoso_island = cardoso_island){
  
  spp_names <- unique(.cardoso_island$morphospecies)
  nspp <- length(spp_names)
  bromeliad_names <- unique(.cardoso_island$Bromeliad)
  nbrom <- length(bromeliad_names)
  
  # simulate from the parameters:
  b_intercept <- rnorm(1, 1, 0.2)
  b_max_water_c <- rnorm(1, 0.8, 0.2)
  b_BS_mean_scale <- rnorm(1, -0.4, 0.2)
  b_max_water_c_by_BS_mean_scale <- rnorm(1, 0, 0.2)
  # b_mean_lg_max <- rnorm(1, 0, 2)
  # b_mean_lg_max_by_max_water_c <- rnorm(1, 0, 2)
  b_zi_intercept <- rnorm(1, -1, 0.2)
  b_zi_max_water_c <- rnorm(1, 0.3, 0.3)
  b_zi_BS_mean_scale <- rnorm(1, -0.2, 0.2)
  b_zi_max_water_c_by_BS_mean_scale <- rnorm(1, 0, 0.2)
  # b_zi_mean_lg_max <- rnorm(1, 0, 2)
  # b_zi_mean_lg_max_by_max_water_c <- rnorm(1, 0, 2)
  
  # simulate from hyperpriors:
  
  # simulate the correlations
  species_corrs <- rethinking::rlkjcorr(1, 8)
  bromeliad_corrs <- rethinking::rlkjcorr(1, 2)
  
  # simulate the standard deviations
  
  species_sd <- rexp(8, 5)
  bromeliad_sd <- rexp(2, 4)
  
  # now matrix multiply to get covariance matrix
  species_Sigma <- diag(species_sd) %*% species_corrs %*% diag(species_sd)
  
  bromeliad_Sigma <- diag(bromeliad_sd) %*% bromeliad_corrs %*% diag(bromeliad_sd)
  
  # use MASS_mvnorm to simulate group effects. these are centered on 0 because
  # they are going to modify the fixed effects in the model for each group
  # (right??)
  
  species_values <- MASS::mvrnorm(nspp, mu = rep(0, 8), Sigma = species_Sigma)
  
  bromeliad_values <- MASS::mvrnorm(nbrom, mu = rep(0, 2), Sigma = bromeliad_Sigma)
  
  # add names to the random effects
  ranef_names <- c(
    "spp_intercept",
    "spp_BS_mean_scale", 
    "spp_max_water_c",
    "spp_max_water_c_by_BS_mean_scale",
    "spp_zi_intercept",
    "spp_zi_BS_mean_scale", 
    "spp_zi_max_water_c",
    "spp_zi_max_water_c_by_BS_mean_scale"
  )
  
  species_effects <- species_values %>% 
    as.data.frame() %>% 
    set_names(ranef_names) %>% 
    mutate(morphospecies = spp_names)
  
  bromeliad_effects <- bromeliad_values %>% 
    as.data.frame %>% 
    set_names(c("brom_intercept", "brom_zi_intercept")) %>% 
    mutate(Bromeliad = bromeliad_names)
  
  # combine with bromeliad and body size data
  
  cardoso_simulated <- .cardoso_island %>% 
    select(morphospecies, Bromeliad, approx_biomass_mg_c, max_vol_log_c) %>% 
    # add in every constant
    mutate(
      b_intercept = b_intercept,
      b_max_water_c = b_max_water_c,
      b_BS_mean_scale = b_BS_mean_scale,
      b_max_water_c_by_BS_mean_scale = b_max_water_c_by_BS_mean_scale,
      b_zi_intercept = b_zi_intercept,
      b_zi_max_water_c = b_zi_max_water_c,
      b_zi_BS_mean_scale = b_zi_BS_mean_scale,
      b_zi_max_water_c_by_BS_mean_scale = b_zi_max_water_c_by_BS_mean_scale
    ) %>% 
    # add in the groups
    left_join(species_effects, by = "morphospecies") %>% 
    left_join(bromeliad_effects, by = "Bromeliad") %>% 
    # calculate the average effect of each part:
    mutate(
      # zero-inflated part:
      zi = 
        b_zi_intercept + 
        b_zi_max_water_c * max_vol_log_c + 
        b_zi_BS_mean_scale * approx_biomass_mg_c + 
        b_zi_max_water_c_by_BS_mean_scale * max_vol_log_c * approx_biomass_mg_c + 
        # groups
        spp_zi_intercept + 
        spp_zi_BS_mean_scale * approx_biomass_mg_c + 
        spp_zi_max_water_c * max_vol_log_c + 
        spp_zi_max_water_c_by_BS_mean_scale * max_vol_log_c * approx_biomass_mg_c,
      count = 
        b_intercept + 
        b_max_water_c * max_vol_log_c + 
        b_BS_mean_scale * approx_biomass_mg_c + 
        b_max_water_c_by_BS_mean_scale * max_vol_log_c * approx_biomass_mg_c + 
        # groups
        spp_intercept + 
        spp_BS_mean_scale * approx_biomass_mg_c + 
        spp_max_water_c * max_vol_log_c + 
        spp_max_water_c_by_BS_mean_scale * max_vol_log_c * approx_biomass_mg_c,
    ) %>%
    # zi is on the logit scale, and count is on the log scale. I need to simulate the liklihood.. how do i do that?
    mutate(
      # does it even colonize??
      p_nonzero = rethinking::inv_logit(zi),
      present = rbinom(2450, 1, p_nonzero),
      # similarly, we transform the counts onto the observation scale then simulate them
      possible_abd = rpois(2450, lambda = exp(count)),
      abundance = present * possible_abd
    )
}