Skip to content

Instantly share code, notes, and snippets.

@aammd
Created May 27, 2018 13:23
Show Gist options
  • Save aammd/21057acfbfd2bc4c31d04051cf931acc to your computer and use it in GitHub Desktop.
Save aammd/21057acfbfd2bc4c31d04051cf931acc to your computer and use it in GitHub Desktop.
simulating a fake dataset
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
)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment