## ----eval = FALSE------------------------------------------------------------- # # my_ipm <- init_ipm(sim_gen = "simple", di_dd = "di", det_stoch = "det") # # ## ----eval = FALSE------------------------------------------------------------- # # my_ipm <- define_kernel( # proto_ipm = my_ipm, # name = "P", # # # The formula describes how the vital rates generate the sub-kernel. This is # # equivalent to equation 3 above. # # formula = s * g, # # # Next, we define the vital rate expressions for the P kernel (Eqs 4-6). # # Here, we create an expression for the inverse of the link function from # # our GLM for survival. In this case, it's the inverse logit (Eq 4). # # s = 1/(1 + exp(-(s_int + s_slope * dbh_1))), # # # Growth is defined by Eqs 5-6 above. There is no inverse link function required # # because we use a Gaussian GLM with an identity link function (i.e. no # # transformation). # # g = dnorm(dbh_2, g_mu, g_sd), # Eq 5 # g_mu = g_int + g_slope * dbh_1, # Eq 6 # # # The family describes the type of transition that kernel produces. See below. # # family = "CC", # # data_list = list( # s_int = 0.2, # coef(my_surv_mod)[1] # s_slope = 0.5, # coef(my_surv_mod)[2] # g_int = 0.1, # coef(my_grow_mod)[1] # g_slope = 1.033, # coef(my_grow_mod)[2] # g_sd = 2.2 # sd(resid(my_grow_mod)) # ), # # # states should be a list of the state variables that the kernel operates on # # states = list(c('dbh')), # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "g") # ) # ## ----eval = FALSE------------------------------------------------------------- # # # the %>% takes the result of the first operation and passes it as the first # # argument to the second function. # # my_ipm <- init_ipm(sim_gen = "simple", di_dd = "di", det_stoch = "det") %>% # define_kernel( # name = "P", # formula = s * g, # family = "CC", # s = 1/(1 + exp(-(s_int + s_slope * dbh_1))), # g = dnorm(dbh_2, g_mu, g_sd), # g_mu = g_int + g_slope * dbh_1, # # data_list = list( # s_int = -3, # s_slope = 0.5, # g_int = 0.1, # g_slope = 1.033, # g_sd = 2.2 # ), # states = list(c('dbh')), # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions(fun = "norm", # target = "g") # ) # ## ----eval = FALSE------------------------------------------------------------- # # my_ipm <- my_ipm %>% # define_kernel( # name = "F", # formula = r_r * r_s * r_d, # Eq 7 # family = "CC", # # # Eq 8. Again, we use the inverse logit transformation to compute pr(flowering) # r_r = 1/(1 + exp(-(r_r_int + r_r_slope * dbh_1))), # # # Eq 9. We exponentiate this because of the log link in our seed production model # r_s = exp(r_s_int + r_s_slope * dbh_1), # # # Eq 10. In this case, both the mean and standard deviation are constants # # r_d = dnorm(dbh_2, r_d_mu, r_d_sd), # data_list = list( # r_r_int = -5, # coef(my_flower_mod)[1] # r_r_slope = 0.1, # coef(my_flower_mod)[2] # r_s_int = -3, # coef(my_seed_mod)[1] # r_s_slope = 0.03, # coef(my_seed_mod)[2] # r_d_mu = 1.2, # mean(my_recr_data$size_2, na.rm = TRUE) # r_d_sd = 0.7 # sd(my_recr_data$size_2, na.rm = TRUE) # ), # states = list(c('dbh')), # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "r_d") # ) # ## ----eval = FALSE------------------------------------------------------------- # # my_ipm <- my_ipm %>% # # define_impl( # # # Create one list of lists with all the info # # list( # # # The components of the big list should be the sub-kernels. Each sub-kernel # # requires 3 pieces of information: the numerical integration rule, the # # state variable it modifies (state_start), and the state variable it # # generates (state_end). # # P = list(int_rule = "midpoint", # state_start = "dbh", # state_end = "dbh"), # F = list(int_rule = "midpoint", # state_start = "dbh", # state_end = "dbh") # ) # ) # ## ----eval = FALSE------------------------------------------------------------- # # # Alternative 1 - call make_impl_args_list() before beginning the IPM creation pipe # # impl_args <- make_impl_args_list( # kernel_names = c("P", "F"), # int_rule = rep("midpoint", 2), # state_start = rep("dbh", 2), # state_end = rep("dbh", 2) # ) # # my_ipm <- my_ipm %>% # define_impl(impl_args) # # my_ipm <- my_ipm %>% # # # Alternative 2, put the call to make_impl_args_list() inside of define_impl(). # # define_impl( # make_impl_args_list( # kernel_names = c("P", "F"), # int_rule = rep("midpoint", 2), # state_start = rep("dbh", 2), # state_end = rep("dbh", 2) # ) # ) # # ## ----eval = FALSE------------------------------------------------------------- # # my_ipm <- my_ipm %>% # define_domains( # dbh = c( # the name of the state variable. # 1, # the lower bound for the domain. L from Eq 1 # 30, # the upper bound for the domain. U from Eq 1 # 200 # the number of mesh points to use for integration # ) # ) # ## ----eval = FALSE------------------------------------------------------------- # # my_ipm <- my_ipm %>% # define_pop_state(n_dbh = rep(1/200, 200)) # ## ----eval = FALSE------------------------------------------------------------- # # my_ipm <- make_ipm(my_ipm, # iterations = 100) # # lambda_ipmr <- lambda(my_ipm) # repro_value <- left_ev(my_ipm) # stable_dist <- right_ev(my_ipm) # # # make_ipm_report works on either proto_ipms or made ipm objects. This # # requires the 'rmarkdown' package to run. # make_ipm_report(my_ipm$proto_ipm, title = "my_ipm_report", render_output = TRUE) # ## ----eval = FALSE------------------------------------------------------------- # # library(ipmr) # # data(iceplant_ex) # # # growth model # # grow_mod <- lm(log_size_next ~ log_size, data = iceplant_ex) # grow_sd <- sd(resid(grow_mod)) # # # # survival model # # surv_mod <- glm(survival ~ log_size, data = iceplant_ex, family = binomial()) # # # Pr(flowering) model # # repr_mod <- glm(repro ~ log_size, data = iceplant_ex, family = binomial()) # # # Number of flowers per plant model # # seed_mod <- glm(flower_n ~ log_size, data = iceplant_ex, family = poisson()) # # # New recruits have no size(t), but do have size(t + 1) # # recr_data <- subset(iceplant_ex, is.na(log_size)) # # recr_mu <- mean(recr_data$log_size_next) # recr_sd <- sd(recr_data$log_size_next) # # # This data set doesn't include information on germination and establishment. # # Thus, we'll compute the realized recruitment parameter as the number # # of observed recruits divided by the number of flowers produced in the prior # # year. # # recr_n <- length(recr_data$log_size_next) # # flow_n <- sum(iceplant_ex$flower_n, na.rm = TRUE) # # # Now, we bundle everything into a list as we did before. We can # # pass the model objects directly into the list, and do not need to extract # # coefficients. # # params <- list(recr_mu = recr_mu, # recr_sd = recr_sd, # grow_sd = grow_sd, # surv_mod = surv_mod, # grow_mod = grow_mod, # repr_mod = repr_mod, # seed_mod = seed_mod, # recr_n = recr_n, # flow_n = flow_n) # # # The lower and upper bounds for integration. Adding 20% on either end to minimize # # eviction. The minimum value is negative, so we multiply by 1.2, rather than 0.8. # # If it were positive, we would need to change that to 0.8. # # L <- min(c(iceplant_ex$log_size, iceplant_ex$log_size_next), na.rm = TRUE) * 1.2 # U <- max(c(iceplant_ex$log_size, iceplant_ex$log_size_next), na.rm = TRUE) * 1.2 # ## ----eval = FALSE------------------------------------------------------------- # # pred_ipm <- init_ipm(sim_gen = "simple", di_dd = "di", det_stoch = "det") %>% # define_kernel( # name = "P", # family = "CC", # formula = s * g, # # # Instead of the inverse logit transformation, we use predict() here. # # We have to be sure that the "newdata" argument of predict is correctly specified. # # This means matching the names used in the model itself (log_size) to the names # # we give the domains (sa_1). In this case, we use "sa" (short for surface area) as # # the new data to generate predictions for. # # s = predict(surv_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # g_mu = predict(grow_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # # # We specify the rest of the kernel the same way. # # g = dnorm(sa_2, g_mu, grow_sd), # states = list(c("sa")), # data_list = params, # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions(fun = "norm", # target = "g") # ) %>% # define_kernel( # name = "F", # family = "CC", # formula = r_p * r_s * r_d * r_r, # # # As above, we use predict(model_object). We make sure the names of the "newdata" # # match the names in the vital rate model formulas, and the values match the # # names of domains they use. # # r_p = predict(repr_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # r_s = predict(seed_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # r_r = recr_n / flow_n, # r_d = dnorm(sa_2, recr_mu, recr_sd), # states = list(c("sa")), # data_list = params, # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions(fun = "norm", # target = "r_d") # ) %>% # define_impl( # make_impl_args_list( # kernel_names = c("P", "F"), # int_rule = rep('midpoint', 2), # state_start = rep("sa", 2), # state_end = rep("sa", 2) # ) # ) %>% # define_domains( # sa = c(L, # U, # 100) # ) %>% # define_pop_state(n_sa = runif(100)) %>% # make_ipm(iterate = TRUE, # iterations = 200) # ## ----eval = FALSE------------------------------------------------------------- # # params <- list(recr_mu = recr_mu, # recr_sd = recr_sd, # grow_sd = grow_sd, # surv_mod = use_vr_model(surv_mod), # wrap the model # grow_mod = use_vr_model(grow_mod), # wrap the model # repr_mod = use_vr_model(repr_mod), # wrap the model # seed_mod = use_vr_model(seed_mod), # wrap the model # recr_n = recr_n, # flow_n = flow_n) # # # pred_ipm <- init_ipm(sim_gen = "simple", di_dd = "di", det_stoch = "det") %>% # define_kernel( # name = "P", # family = "CC", # formula = s * g, # s = predict(surv_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # g_mu = predict(grow_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # g = dnorm(sa_2, g_mu, grow_sd), # states = list(c("sa")), # data_list = params, # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "g") # ) %>% # define_kernel( # name = "F", # family = "CC", # formula = r_p * r_s * r_d * r_r, # r_p = predict(repr_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # r_s = predict(seed_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # r_r = recr_n / flow_n, # r_d = dnorm(sa_2, recr_mu, recr_sd), # states = list(c("sa")), # data_list = params, # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "r_d") # ) %>% # define_impl( # make_impl_args_list( # kernel_names = c("P", "F"), # int_rule = rep('midpoint', 2), # state_start = rep("sa", 2), # state_end = rep("sa", 2) # ) # ) %>% # define_domains( # sa = c(L, # U, # 100) # ) %>% # define_pop_state(n_sa = runif(100)) %>% # make_ipm(iterate = TRUE, # iterations = 200) # ## ----eval = FALSE------------------------------------------------------------- # # R_0 <- function(ipm) { # # Fm <- ipm$sub_kernels$F # Pm <- ipm$sub_kernels$P # # I <- diag(nrow(Pm)) # # N <- solve(I - Pm) # # R <- Fm %*% N # # out <- Re(eigen(R)$values[1]) # # return(out) # # } # # gen_T <- function(ipm) { # # R0 <- R_0(ipm) # lam <- lambda(ipm) %>% as.vector() # # out <- log(R0) / log(lam) # # return(out) # # } # # R_0(pred_ipm) # gen_T(pred_ipm) # ## ----eval = FALSE------------------------------------------------------------- # # library(ipmr) # # # Define some fixed parameters # # fixed_list <- list( # s_int = 1.03, # fixef(my_surv_mod)[1] - uses fixef because we now have a model with random effects # s_slope = 2.2, # fixef(my_surv_mod)[2] # g_int = 3.7, # fixef(my_grow_mod)[1] # g_slope = 0.92, # fixef(my_grow_mod)[2] # sd_g = 0.9, # sd(resid(my_grow_mod)) # r_r_int = 0.09, # coef(my_repro_mod)[1] - uses coef because there are no random effects in this model # r_r_slope = 0.05, # coef(my_repro_mod)[2] # r_s_int = 0.1, # fixef(my_flower_mod)[1] # r_s_slope = 0.005, # fixef(my_flower_mod)[2] # mu_rd = 9, # mean(my_recr_data$size_2, na.rm = TRUE) # sd_rd = 2 # sd(my_recr_data$size_2, na.rm = TRUE) # ) # ## ----eval = FALSE------------------------------------------------------------- # # # Now, simulate some random intercepts for growth (g_), survival (s_), # # and offspring production (r_s_). This part is for the purpose of the example. # # # First, we create vector of values that each random component can take. # # g_r_int <- rnorm(5, 0, 0.3) # unlist(ranef(my_grow_mod)) for an lme4 output # s_r_int <- rnorm(5, 0, 0.7) # unlist(ranef(my_surv_mod)) for an lme4 output # r_s_r_int <- rnorm(5, 0, 0.2) # unlist(ranef(my_flower_mod)) for an lme4 output # # # We'll call our hypothetical sites 1, 2, 3, 4, and 5. The "r" prefix is to # # remind us that these are random quantities. # # nms <- paste("r_", 1:5, sep = "") # # names(g_r_int) <- paste('g_', nms, sep = "") # names(s_r_int) <- paste('s_', nms, sep = "") # names(r_s_r_int) <- paste('r_s_', nms, sep = "") # # # Each set of parameters is converted to a named list. The names should match # # the variables referenced in each define_kernel() call. # # g_params <- as.list(g_r_int) # s_params <- as.list(s_r_int) # r_s_params <- as.list(r_s_r_int) # # # add them all together using c() # # all_params_list <- c(fixed_list, g_params, s_params, r_s_params) # ## ----eval = FALSE------------------------------------------------------------- # # my_ipm <- init_ipm(sim_gen = "simple", di_dd = "di", det_stoch = "det") %>% # define_kernel( # # # Our P kernels will vary from site to site, so we index it with "_site" # # name = 'P_site', # # # Similarly, our survival and growth functions will vary from site to site # # so these are also indexed # # formula = s_site * g_site, # family = "CC", # # # The linear predictor for the survival function can be split out # # into its own expression as well. This might help keep track of things. # # Survival is indexed by site as well. # # s_lin_site = s_int + s_r_site + s_slope * ht_1, # s_site = 1 / (1 + exp(-s_lin_site)), # # # Again, we modify the vital rate expression to include "_site". # # g_site = dnorm(ht_2, mean = mu_g_site, sd = sd_g), # mu_g_site = g_int + g_slope * ht_1 + g_r_site, # # data_list = all_params_list, # states = list(c('ht')), # # # Here, we tell ipmr that the model has some parameter sets, and # # provide a list describing the values the index can take. The values in # # par_set_indices are substituted for "site" everywhere in the model, except # # for the data list. This is why we had to make sure that the names there # # matched the levels we supply here. # # uses_par_sets = TRUE, # par_set_indices = list(site = 1:5), # # # We must also index the variables in the eviction function # # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "g_site") # # ) %>% # define_kernel( # # # The F kernel also varies from site to site # # name = "F_site", # formula = r_r * r_s_site * r_d, # family = "CC", # # # In this example, we didn't include a site level effect for probability # # of flowering, only seed production. Thus, this expression is NOT indexed. # # r_r_lin = r_r_int + r_r_slope * ht_1, # r_r = 1 / (1 + exp(- r_r_lin)), # # # We index the seed production expression with the site effect # # r_s_site = exp(r_s_int + r_s_r_site + r_s_slope * ht_1), # r_d = dnorm(ht_2, mean = mu_rd, sd = sd_rd), # data_list = all_params_list, # states = list(c('ht')), # # # As in the P kernel, we specify the values the index can have. # # uses_par_sets = TRUE, # par_set_indices = list(site = 1:5), # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "r_d") # ) %>% # define_impl( # make_impl_args_list( # # # The impl_args are also modified with the index # # kernel_names = c("P_site", "F_site"), # int_rule = rep("midpoint", 2), # state_start = rep("ht", 2), # state_end = rep("ht", 2) # ) # ) %>% # define_domains(ht = c(0.2, 40, 100)) %>% # # # We also append the suffix in define_pop_state(). THis will create a deterministic # # simulation for every "site" # # define_pop_state(n_ht_site = runif(100)) %>% # make_ipm(iterate = TRUE, # iterations = 100) # # lambda(my_ipm) # # ## ----eval = FALSE------------------------------------------------------------- # # library(ipmr) # # # Define some fixed parameters # # fixed_list <- list( # s_int = 1.03, # fixef(my_surv_mod)[1] - uses fixef because we now have a model with random effects # s_slope = 2.2, # fixef(my_surv_mod)[2] # g_int = 3.7, # fixef(my_grow_mod)[1] # g_slope = 0.92, # fixef(my_grow_mod)[2] # sd_g = 0.9, # sd(resid(my_grow_mod)) # r_r_int = 0.09, # coef(my_repro_mod)[1] - uses coef because there are no random effects in this model # r_r_slope = 0.05, # coef(my_repro_mod)[2] # r_s_int = 0.1, # fixef(my_flower_mod)[1] # r_s_slope = 0.005, # fixef(my_flower_mod)[2] # mu_fd = 9, # mean(my_recr_data$size_2, na.rm = TRUE) # sd_fd = 2 # sd(my_recr_data$size_2, na.rm = TRUE) # ) # ## ----eval = FALSE------------------------------------------------------------- # # Now, simulate some random intercepts for growth (g_), survival (s_), # # and offspring production (r_s_). This part is for the purpose of the example. # # # First, we create vector of values corresponding to # # g_r_int <- rnorm(5, 0, 0.3) # unlist(ranef(my_grow_mod)) for an lme4 output # s_r_int <- rnorm(5, 0, 0.7) # unlist(ranef(my_surv_mod)) for an lme4 output # r_s_r_int <- rnorm(5, 0, 0.2) # unlist(ranef(my_flower_mod)) for an lme4 output # # nms <- paste("r_", 1:5, sep = "") # # names(g_r_int) <- paste('g_', nms, sep = "") # names(s_r_int) <- paste('s_', nms, sep = "") # names(r_s_r_int) <- paste('r_s_', nms, sep = "") # # # Each set of parameters is converted to a named list. The names should match # # the variables referenced in each define_kernel() call. # # g_params <- as.list(g_r_int) # s_params <- as.list(s_r_int) # r_s_params <- as.list(r_s_r_int) # # # add them all together using c() # # all_params_list <- c(fixed_list, g_params, s_params, r_s_params) # ## ----eval = FALSE------------------------------------------------------------- # # inv_logit <- function(sv, int, slope) { # return( # 1/(1 + exp(-(int + slope * sv))) # ) # } # # # same as above, but handles the extra term from the random effect we simulated. # # inv_logit_r <- function(sv, int, slope, r_eff) { # return( # 1/(1 + exp(-(int + slope * sv + r_eff))) # ) # } # # pois_r <- function(sv, int, slope, r_eff) { # return( # exp( # int + slope * sv + r_eff # ) # ) # } # # my_funs <- list(inv_logit = inv_logit, # inv_logit_r = inv_logit_r, # pois_r = pois_r) ## ----eval = FALSE------------------------------------------------------------- # # my_ipm <- init_ipm(sim_gen = "simple", # di_dd = "di", # det_stoch = "stoch", # kern_param = "kern") %>% # define_kernel( # name = 'P_yr', # P becomes P_yr # formula = s_yr * g_yr, # g and s become g_yr and s_yr, respectively # family = "CC", # # # Note the usage of the inv_logit_r, which we defined in the block above. # # it is passed to make_ipm() # # s_yr = inv_logit_r(ht_1, s_int, s_slope, s_r_yr), # g_yr = dnorm(ht_2, mean = mu_g_yr, sd = sd_g), # mu_g_yr = g_int + g_slope * ht_1 + g_r_yr, # # # all_params_list contains the named parameters g_r_1, g_r_2, s_r_1, s_r_2, etc. # # This is the only level where the user is required to fully expand the name # # X par_set_indices combinations. # # data_list = all_params_list, # states = list(c('ht')), # uses_par_sets = TRUE, # par_set_indices = list(yr = 1:5), # evict_cor = TRUE, # # # reference to g_yr in evict_fun is also updated # # evict_fun = truncated_distributions("norm", "g_yr") # # ) %>% # define_kernel( # name = "F_yr", # Update the names as we did for the P kernel # formula = r_r * r_s_yr * r_d, # family = "CC", # r_r = inv_logit(ht_1, r_r_int, r_r_slope), # r_s_yr = pois_r(ht_1, r_s_int, r_s_slope, r_s_r_yr), # r_d = dnorm(ht_2, mean = mu_fd, sd = sd_fd), # data_list = all_params_list, # states = list(c('ht')), # uses_par_sets = TRUE, # par_set_indices = list(yr = 1:5), # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "r_d") # ) %>% # define_impl( # make_impl_args_list( # kernel_names = c("P_yr", "F_yr"), # int_rule = rep("midpoint", 2), # state_start = rep("ht", 2), # state_end = rep("ht", 2) # ) # ) %>% # define_domains(ht = c(0.2, 40, 100)) %>% # define_pop_state(n_ht = runif(100)) %>% # make_ipm(usr_funs = my_funs, # kernel_seq = sample(1:5, 100, replace = TRUE), # iterate = TRUE, # iterations = 100) # # # # The default for stochastic lambda is to compute mean(log(ipm$pop_state$lambda)). # # It also removes the first 10% of iterations to handle "burn in". The amount # # removed can be adjusted with the burn_in parameter. # # stoch_lambda <- lambda(my_ipm) # # # lambda(comp_method = 'pop_size', type = 'all') will compute the population # # growth rate for every time step as the sum(n_ht_t_1) / sum(n_ht_t). # # iteration_lambdas <- lambda(my_ipm, type_lambda = 'all') # # ## ----eval = FALSE------------------------------------------------------------- # # library(ipmr) # # # Define the fixed parameters in a list. The temperature and precipitation # # coefficients are defined as s_temp, s_precip, g_temp, and g_precip. # # constant_params <- list( # s_int = -10, # s_slope = 1.5, # s_precip = 0.00001, # s_temp = -0.003, # g_int = 0.2, # g_slope = 1.01, # g_sd = 1.2, # g_temp = -0.002, # g_precip = 0.004, # r_r_int = -3.2, # r_r_slope = 0.55, # r_s_int = -0.4, # r_s_slope = 0.5, # r_d_mu = 1.1, # r_d_sd = 0.1 # ) # # # Now, we create a set of environmental covariates. In this example, we use # # a normal distribution for temperature and a Gamma for precipitation. # # env_params <- list( # temp_mu = 8.9, # temp_sd = 1.2, # precip_shape = 1000, # precip_rate = 2 # ) # # # We define a wrapper function that samples from these distributions # # sample_env <- function(env_params) { # # # We generate one value for each covariate per iteration, and return it # # as a named list. We can reference the names in this list in vital rate # # expressions. # # temp_now <- rnorm(1, # env_params$temp_mu, # env_params$temp_sd) # # precip_now <- rgamma(1, # shape = env_params$precip_shape, # rate = env_params$precip_rate) # # out <- list(temp = temp_now, precip = precip_now) # # return(out) # # } # # # # Again, we can define our own functions and pass them into calls to make_ipm. This # # isn't strictly necessary, but can make the model code more readable/less error prone. # # inv_logit <- function(lin_term) { # 1/(1 + exp(-lin_term)) # } # # ## ----eval = FALSE------------------------------------------------------------- # # init_pop_vec <- runif(100) # # param_resamp_model <- init_ipm(sim_gen = "simple", # di_dd = "di", # det_stoch = "stoch", # kern_param = "param") %>% # define_kernel( # name = 'P', # formula = s * g, # family = 'CC', # # # Parameters created by define_env_state() can be referenced by name just like # # any other parameter in the model. # # g_mu = g_int + g_slope * surf_area_1 + g_temp * temp + g_precip * precip, # s_lin_p = s_int + s_slope * surf_area_1 + s_temp * temp + s_precip * precip, # s = inv_logit(s_lin_p), # g = dnorm(surf_area_2, g_mu, g_sd), # # # data_list = constant_params, # states = list(c('surf_area')), # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "g") # ) %>% # define_kernel( # name = 'F', # formula = r_r * r_s * r_d, # family = 'CC', # r_r_lin_p = r_r_int + r_r_slope * surf_area_1, # r_r = inv_logit(r_r_lin_p), # r_s = exp(r_s_int + r_s_slope * surf_area_1), # r_d = dnorm(surf_area_2, r_d_mu, r_d_sd), # data_list = constant_params, # states = list(c('surf_area')), # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "r_d") # ) %>% # define_impl( # make_impl_args_list( # kernel_names = c("P", "F"), # int_rule = rep('midpoint', 2), # state_start = rep('surf_area', 2), # state_end = rep('surf_area', 2) # ) # ) %>% # define_domains(surf_area = c(0, 10, 100)) # ## ----eval = FALSE------------------------------------------------------------- # # param_resamp_model <- param_resamp_model %>% # # define_env_state( # env_covs = sample_env(env_params), # data_list = list(env_params = env_params, # sample_env = sample_env) # ) %>% # define_pop_state( # pop_vectors = list( # n_surf_area = init_pop_vec # ) # ) %>% # make_ipm(usr_funs = list(inv_logit = inv_logit), # iterate = TRUE, # iterations = 10, # return_sub_kernels = TRUE) # # # By default, lambda computes stochastic lambda with stochastic models # # lambda(param_resamp_model) # # # We can get lambdas for each time step by adding type_lambda = "all" # # lambda(param_resamp_model, type_lambda = 'all') # # # If we want to see the actual draws that were used at each step of the # # model iteration, we can access these using the output's $env_seq slot. # # param_resamp_model$env_seq # ## ----eval = FALSE------------------------------------------------------------- # # param_resamp_model <- init_ipm(sim_gen = "simple", # di_dd = "di", # det_stoch = "stoch", # kern_param = "param") %>% # define_kernel( # name = 'P', # formula = s * g, # family = 'CC', # # # Parameters created by define_env_state() can be referenced by name just like # # any other parameter in the model. # # g_mu = g_int + g_slope * surf_area_1 + g_temp * temp + g_precip * precip, # s_lin_p = s_int + s_slope * surf_area_1 + s_temp * temp + s_precip * precip, # s = inv_logit(s_lin_p), # g = dnorm(surf_area_2, g_mu, g_sd), # # # data_list = constant_params, # states = list(c('surf_area')), # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "g") # ) %>% # define_kernel( # name = 'F', # formula = r_r * r_s * r_d, # family = 'CC', # r_r_lin_p = r_r_int + r_r_slope * surf_area_1, # r_r = inv_logit(r_r_lin_p), # r_s = exp(r_s_int + r_s_slope * surf_area_1), # r_d = dnorm(surf_area_2, r_d_mu, r_d_sd), # data_list = constant_params, # states = list(c('surf_area')), # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "r_d") # ) %>% # define_impl( # make_impl_args_list( # kernel_names = c("P", "F"), # int_rule = rep('midpoint', 2), # state_start = rep('surf_area', 2), # state_end = rep('surf_area', 2) # ) # ) %>% # define_domains(surf_area = c(0, 10, 100)) # # # ## ----eval = FALSE------------------------------------------------------------- # # env_states <- data.frame(precip = rgamma(10, shape = 1000, rate = 2), # temp = rnorm(10, mean = 8.9, sd = 1.2)) # ## ----eval = FALSE------------------------------------------------------------- # # sample_env <- function(env_states, iteration) { # # out <- as.list(env_states[iteration, ]) # names(out) <- names(env_states) # # return(out) # # } # ## ----eval = FALSE------------------------------------------------------------- # # param_resamp_model <- param_resamp_model %>% # define_env_state(env_params = sample_env(env_states, # iteration = t), # "t" indexes the current model iteration # data_list = list( # env_states = env_states, # sample_env = sample_env # )) %>% # define_pop_state( # pop_vectors = list( # n_surf_area = init_pop_vec # ) # ) %>% # make_ipm(usr_funs = list(inv_logit = inv_logit), # iterate = TRUE, # iterations = 10) ## ----eval = FALSE------------------------------------------------------------- # # # param_resamp_model <- init_ipm(sim_gen = "simple", # di_dd = "di", # det_stoch = "stoch", # kern_param = "param") %>% # define_kernel( # name = 'P', # formula = s * g, # family = 'CC', # # # Parameters created by define_env_state() can be referenced by name just like # # any other parameter in the model. # # g_mu = g_int + g_slope * surf_area_1 + g_temp * temp + g_precip * precip, # s_lin_p = s_int + s_slope * surf_area_1 + s_temp * temp + s_precip * precip, # s = inv_logit(s_lin_p), # g = dnorm(surf_area_2, g_mu, g_sd), # # # data_list = constant_params, # states = list(c('surf_area')), # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "g") # ) %>% # define_kernel( # name = 'F', # formula = r_r * r_s * r_d, # family = 'CC', # r_r_lin_p = r_r_int + r_r_slope * surf_area_1, # r_r = inv_logit(r_r_lin_p), # r_s = exp(r_s_int + r_s_slope * surf_area_1), # r_d = dnorm(surf_area_2, r_d_mu, r_d_sd), # data_list = constant_params, # states = list(c('surf_area')), # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions("norm", "r_d") # ) %>% # define_impl( # make_impl_args_list( # kernel_names = c("P", "F"), # int_rule = rep('midpoint', 2), # state_start = rep('surf_area', 2), # state_end = rep('surf_area', 2) # ) # ) %>% # define_domains(surf_area = c(0, 10, 100)) # # # ## ----eval = FALSE------------------------------------------------------------- # # env_sampler <- function(time, init_temp = 10, init_precip = 1500) { # # t_est <- init_temp + 0.2 * time + rnorm(1, mean = 0, sd = 0.2) # # p_est <- init_precip + -3.3 * time + rnorm(1, mean = 0, sd = 100) # # if(p_est <= 0) p_est <- 1 # # out <- list(temp = t_est, # precip = p_est) # # return(out) # # } # # ## ----eval = FALSE------------------------------------------------------------- # # param_resamp_ipm <- param_resamp_model %>% # define_env_state( # env_params = env_sampler(time = t, init_temp = 10, init_precip = 1500), # data_list = list() # ) %>% # define_pop_state( # n_surf_area = init_pop_vec # ) %>% # make_ipm( # iterations = 100, # usr_funs = list(env_sampler = env_sampler, # inv_logit = inv_logit) # ) # # lambda(param_resamp_ipm) # ## ----eval = FALSE------------------------------------------------------------- # # library(ipmr) # # data(iceplant_ex) # # grow_mod <- lm(log_size_next ~ log_size, data = iceplant_ex) # grow_sd <- sd(resid(grow_mod)) # # surv_mod <- glm(survival ~ log_size, data = iceplant_ex, family = binomial()) # repr_mod <- glm(repro ~ log_size, data = iceplant_ex, family = binomial()) # seed_mod <- glm(flower_n ~ log_size, data = iceplant_ex, family = poisson()) # recr_data <- subset(iceplant_ex, is.na(log_size)) # # recr_mu <- mean(recr_data$log_size_next) # recr_sd <- sd(recr_data$log_size_next) # recr_n <- length(recr_data$log_size_next) # flow_n <- sum(iceplant_ex$flower_n, na.rm = TRUE) # # params <- list(recr_mu = recr_mu, # recr_sd = recr_sd, # grow_sd = grow_sd, # surv_mod = surv_mod, # grow_mod = grow_mod, # repr_mod = repr_mod, # seed_mod = seed_mod, # recr_n = recr_n, # flow_n = flow_n) # # L <- min(c(iceplant_ex$log_size, iceplant_ex$log_size_next), na.rm = TRUE) * 1.2 # U <- max(c(iceplant_ex$log_size, iceplant_ex$log_size_next), na.rm = TRUE) * 1.2 # # obs_ipm <- init_ipm(sim_gen = "simple", # di_dd = "di", # det_stoch = "det") %>% # define_kernel( # name = "P", # family = "CC", # formula = s * g, # # # Instead of the inverse logit transformation, we use predict() here. # # We have to be sure that the "newdata" argument of predict is correctly specified. # # This means matching the names used in the model itself (log_size) to the names # # we give the domains. In this case, we use "sa" (short for surface area) as # # the new data to generate predictions for. # # s = predict(surv_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # g_mu = predict(grow_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # # # We specify the rest of the kernel the same way. # # g = dnorm(sa_2, g_mu, grow_sd), # states = list(c("sa")), # data_list = params, # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions(fun = "norm", # target = "g") # ) %>% # define_kernel( # name = "F", # family = "CC", # formula = r_p * r_s * r_d * r_r, # # # As above, we use predict(model_object). We make sure the names of the "newdata" # # match the names in the vital rate model formulas, and the values match the # # names of domains they use. # # r_p = predict(repr_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # r_s = predict(seed_mod, # newdata = data.frame(log_size = sa_1), # type = 'response'), # r_r = recr_n / flow_n, # r_d = dnorm(sa_2, recr_mu, recr_sd), # states = list(c("sa")), # data_list = params, # uses_par_sets = FALSE, # evict_cor = TRUE, # evict_fun = truncated_distributions(fun = "norm", # target = "r_d") # ) %>% # define_impl( # make_impl_args_list( # kernel_names = c("P", "F"), # int_rule = rep('midpoint', 2), # state_start = rep("sa", 2), # state_end = rep("sa", 2) # ) # ) %>% # define_domains( # sa = c(L, # U, # 100) # ) %>% # define_pop_state(n_sa = runif(100)) %>% # make_ipm(iterate = TRUE, # iterations = 100) # # lambda_obs <- lambda(obs_ipm) # ## ----eval = FALSE------------------------------------------------------------- # # all_lambdas <- numeric(50L) # # recr_data <- subset(iceplant_ex, is.na(log_size)) # # adults <- subset(iceplant_ex, !is.na(log_size)) # # use_proto <- obs_ipm$proto_ipm # ## ----eval = FALSE------------------------------------------------------------- # # # for(i in 1:50) { # # sample_ind <- seq(1, nrow(adults), by = 1) # # boot_ind <- sample(sample_ind, size = nrow(adults), replace = TRUE) # # boot_data <- rbind(adults[boot_ind, ], # recr_data) # # grow_mod <- lm(log_size_next ~ log_size, data = boot_data) # grow_sd <- sd(resid(grow_mod)) # # surv_mod <- glm(survival ~ log_size, data = boot_data, family = binomial()) # repr_mod <- glm(repro ~ log_size, data = boot_data, family = binomial()) # seed_mod <- glm(flower_n ~ log_size, data = boot_data, family = poisson()) # # recr_mu <- mean(recr_data$log_size_next) # recr_sd <- sd(recr_data$log_size_next) # recr_n <- length(recr_data$log_size_next) # flow_n <- sum(boot_data$flower_n, na.rm = TRUE) # # params <- list(recr_mu = recr_mu, # recr_sd = recr_sd, # grow_sd = grow_sd, # surv_mod = surv_mod, # grow_mod = grow_mod, # repr_mod = repr_mod, # seed_mod = seed_mod, # recr_n = recr_n, # flow_n = flow_n) # # # Insert the new vital rate models into the proto_ipm, then rebuild the IPM. # # parameters(use_proto) <- params # # boot_ipm <- use_proto %>% # make_ipm(iterate = TRUE, # iterations = 100) # # all_lambdas[i] <- lambda(boot_ipm) # # } # # # Plot the results # # hist(all_lambdas) # abline(v = lambda_obs, col = 'red', lwd = 2, lty = 2) # ## ----eval = FALSE------------------------------------------------------------- # # browseVignettes("ipmr") #