Why is Stochastic Lambda > Deterministic?
Eric R. Scott
2022-04-21
1 TODO
Instead of deterministic IPM using different vital rates models, use the random effects models, but plug in
newdata = (log_size_prev = log_size, year = "newlevel")to get mean fitted values.Plot population trajectories for all models
Read Rees & Ellner 2009. Consider either fixed effect of year + kernel sampling OR random effect of year + parameter sampling. Maybe no added value for random effect of year + kernel sampling??
Rees, Mark, and Stephen P. Ellner. “Integral Projection Models for Populations in Temporally Varying Environments.” Ecological Monographs 79, no. 4 (2009): 575–94. https://doi.org/10.1890/08-1474.1.
Last compiled: 2022-05-03
2 Purpose
Figure out what the heck is going on.
3 Load Data
4 Pick two years, worst and best
Fit determinsitic IPMs to all years separately to identify a bad year and a good year. Then subset data to only contain those years
single_years <-
data_full %>%
filter(year>1999) %>%
group_by(year_fct) %>%
group_split()
#fit vital rates for FF
vit_list_yr_ff <-
single_years %>%
map(~dplyr::filter(.x, habitat == "1-ha")) %>%
map(~{
vit_list_det = c(
list(
vit_surv = surv_det(.x),
vit_size = size_det(.x),
vit_flwr = flwr_det(.x),
vit_size_sdlg = size_sdlg_det(.x),
vit_surv_sdlg = surv_sdlg_det(.x)
),
vit_other_ff)
})
# make IPMs for FF
ipms_yr_ff <-
map(vit_list_yr_ff, ~{
make_proto_ipm_det(.x, pop_vec = pop_vec_ff) %>%
make_ipm(iterations = 1000, #only needs 100 to converge
normalize_pop_size = FALSE, # to run as PVA
usr_funs = list(get_scat_params = get_scat_params))
})
# fit vital rates for CF
vit_list_yr_cf <-
single_years %>%
map(~dplyr::filter(.x, habitat == "CF")) %>%
map(~{
vit_list_det = c(
list(
vit_surv = surv_det(.x),
vit_size = size_det(.x),
vit_flwr = flwr_det(.x),
vit_size_sdlg = size_sdlg_det(.x),
vit_surv_sdlg = surv_sdlg_det(.x)
),
vit_other_cf)
})
#make IPMs for CF
ipms_yr_cf <-
map(vit_list_yr_cf, ~{
make_proto_ipm_det(.x, pop_vec = pop_vec_cf) %>%
make_ipm(iterations = 1000, #only needs 100 to converge
normalize_pop_size = FALSE, # to run as PVA
usr_funs = list(get_scat_params = get_scat_params))
})
#lambdas
lambdas_ff <- ipms_yr_ff %>% map_dbl(lambda)
lambdas_cf <- ipms_yr_cf %>% map_dbl(lambda)#worst year
(worst_year_ff <- c(2000:max(data_full$year))[which.min(lambdas_ff)]); min(lambdas_ff)## [1] 2003## [1] 0.9397404(worst_year_cf <- c(2000:max(data_full$year))[which.min(lambdas_cf)]); min(lambdas_cf)## [1] 2003## [1] 0.8810649#best year
(best_year_ff <- c(2000:max(data_full$year))[which.max(lambdas_ff)]); max(lambdas_ff)## [1] 2001## [1] 1.005335(best_year_cf <- c(2000:max(data_full$year))[which.max(lambdas_cf)]); max(lambdas_cf)## [1] 2002## [1] 1.027696plot(2000:max(data_full$year), lambdas_cf, type = "l", col = "red"); lines(x =2000:max(data_full$year), y = lambdas_ff, col = "blue", type = "l")
2003 was the worst year in both habitats. 2001 was a good year in both habitats (but not the best in CF)
#combine good and bad year
two_yrs <- data_full %>% filter(year %in% c(2003, 2001))
two_yrs_ff <- two_yrs %>% filter(habitat == "1-ha")
two_yrs_cf <- two_yrs %>% filter(habitat == "CF")5 Fit IPMs
5.1 Deterministic
Fit a deterministic IPM with just those two years of data
Fragments:
vit_list_det_ff = c(
list(
vit_surv = surv_det(two_yrs_ff),
vit_size = size_det(two_yrs_ff),
vit_flwr = flwr_det(two_yrs_ff),
vit_size_sdlg = size_sdlg_det(two_yrs_ff),
vit_surv_sdlg = surv_sdlg_det(two_yrs_ff)
),
vit_other_ff)
ipm_det_ff = make_proto_ipm_det(vit_list_det_ff, pop_vec = pop_vec_ff) %>%
make_ipm(iterations = 1000, #only needs 100 to converge
normalize_pop_size = FALSE, # to run as PVA
usr_funs = list(get_scat_params = get_scat_params))
# FF
lambda(ipm_det_ff)## lambda
## 0.9644405CF:
vit_list_det_cf = c(
list(
vit_surv = surv_det(two_yrs_cf),
vit_size = size_det(two_yrs_cf),
vit_flwr = flwr_det(two_yrs_cf),
vit_size_sdlg = size_sdlg_det(two_yrs_cf),
vit_surv_sdlg = surv_sdlg_det(two_yrs_cf)
),
vit_other_cf)
ipm_det_cf = make_proto_ipm_det(vit_list_det_cf, pop_vec = pop_vec_cf) %>%
make_ipm(iterations = 1000, #only needs 100 to converge
normalize_pop_size = FALSE, # to run as PVA
usr_funs = list(get_scat_params = get_scat_params))
lambda(ipm_det_cf)## lambda
## 0.95297915.2 Stochastic
The stochastic IPM requires re-writing the proto_ipm
because it’s now using different parameter sets for year. This now uses
vital rates models where there is a random effect of year fit such that
each year can have a different shape smoother (but not different
wiggliness).
## forests fragments
vit_list_stoch_ff = c(
list(
vit_surv = surv_raneff(two_yrs_ff),
vit_size = size_raneff(two_yrs_ff),
vit_flwr = flwr_raneff(two_yrs_ff),
vit_size_sdlg = size_sdlg_raneff(two_yrs_ff),
vit_surv_sdlg = surv_sdlg_raneff(two_yrs_ff)
),
vit_other_ff)
proto_ipm_stoch_ff <-
init_ipm(
sim_gen = "general",
di_dd = "di",
det_stoch = "stoch",
kern_param = "kern"
) %>%
# Define growth and survival kernel
define_kernel(
"P_yr", #_yr makes a different P kernel for each value of `yr` from par_set_indices
formula = s_z * G_z1z * d_log_size,
family = "CC",
# Probability of survival
s_z = predict(vit_surv,
newdata = tibble(log_size_prev = log_size_1, year_fct = yr),
type = 'response'),
# Growth
G_z1z = dt.scaled(
log_size_2,
df = size_nu,
mean = size_mu,
sd = size_sd
),
size_mu = predict(vit_size,
newdata = tibble(log_size_prev = log_size_1, year_fct = yr),
type = 'response'),
size_nu = get_scat_params(vit_size)["nu"],
size_sd = get_scat_params(vit_size)["sd"],
data_list = vit_list_stoch_ff,
states = list(c('log_size')),
uses_par_sets = TRUE,
par_set_indices = list(yr = c(2001, 2003)),
evict_cor = TRUE,
evict_fun = truncated_distributions("t.scaled", "G_z1z")
) %>%
# Define fecundity and recruitment kernel
define_kernel(
"go_sdlg_yr",
formula = p_f * fruits * seeds * g_e * d_log_size,
family = "CD",
# probability of flowering
p_f = predict(vit_flwr,
newdata = tibble(log_size_prev = log_size_1, year_fct = yr),
type = 'response'),
fruits = predict(vit_fruits,
newdata = tibble(log_size_prev = log_size_1),
type = 'response'),
seeds = exp(coef(vit_seeds)[1]), #mean seeds per fruit
g_e = vit_germ_est,
data_list = vit_list_stoch_ff,
states = list(c('log_size', 'sdlg')),
uses_par_sets = TRUE,
par_set_indices = list(yr = c(2001, 2003))
) %>%
define_kernel(
name = "stay_sdlg",
formula = 0,
family = "DD",
states = list(c("sdlg")),
evict_cor = FALSE
) %>%
# Enter into post-seedling
define_kernel(
name = "leave_sdlg_yr",
formula = s_sdlg * G_z2_sdlg,
family = "DC",
s_sdlg = predict(vit_surv_sdlg,
newdata = tibble(year_fct = yr),
type = 'response'), #intercept from intercept-only glm
G_z2_sdlg = dt.scaled(
x = log_size_2,
df = size_sdlg_nu,
mean = size_sdlg_mu,
sd = size_sdlg_sd
),
size_sdlg_mu = coef(vit_size_sdlg)[1], #intercept from intercept-only glm
size_sdlg_nu = get_scat_params(vit_size_sdlg)["nu"],
size_sdlg_sd = get_scat_params(vit_size_sdlg)["sd"],
data_list = vit_list_stoch_ff,
states = list(c("sdlg", "log_size")),
uses_par_sets = TRUE,
par_set_indices = list(yr = c(2001, 2003)),
evict_cor = TRUE,
evict_fun = truncated_distributions("t.scaled", "G_z2_sdlg")
) %>%
# define implementation with midpoint rule
define_impl(make_impl_args_list(
kernel_names = c("P_yr", "go_sdlg_yr", "stay_sdlg", "leave_sdlg_yr"),
int_rule = "midpoint",
state_start = c("log_size", "log_size", "sdlg", "sdlg"),
state_end = c("log_size", "sdlg", "sdlg", "log_size")
)) %>%
#define lower bound, upper bound, and number of meshpoints for log_size
define_domains(log_size = c(0, 8.018296, 100)) %>%
#arbitrary starting population state
define_pop_state(pop_vectors = pop_vec_ff)## forests fragments
vit_list_stoch_cf = c(
list(
vit_surv = surv_raneff(two_yrs_cf),
vit_size = size_raneff(two_yrs_cf),
vit_flwr = flwr_raneff(two_yrs_cf),
vit_size_sdlg = size_sdlg_raneff(two_yrs_cf),
vit_surv_sdlg = surv_sdlg_raneff(two_yrs_cf)
),
vit_other_cf)
proto_ipm_stoch_cf <-
init_ipm(
sim_gen = "general",
di_dd = "di",
det_stoch = "stoch",
kern_param = "kern"
) %>%
# Define growth and survival kernel
define_kernel(
"P_yr", #_yr makes a different P kernel for each value of `yr` from par_set_indices
formula = s_z * G_z1z * d_log_size,
family = "CC",
# Probability of survival
s_z = predict(vit_surv,
newdata = tibble(log_size_prev = log_size_1, year_fct = yr),
type = 'response'),
# Growth
G_z1z = dt.scaled(
log_size_2,
df = size_nu,
mean = size_mu,
sd = size_sd
),
size_mu = predict(vit_size,
newdata = tibble(log_size_prev = log_size_1, year_fct = yr),
type = 'response'),
size_nu = get_scat_params(vit_size)["nu"],
size_sd = get_scat_params(vit_size)["sd"],
data_list = vit_list_stoch_cf,
states = list(c('log_size')),
uses_par_sets = TRUE,
par_set_indices = list(yr = c(2001, 2003)),
evict_cor = TRUE,
evict_fun = truncated_distributions("t.scaled", "G_z1z")
) %>%
# Define fecundity and recruitment kernel
define_kernel(
"go_sdlg_yr",
formula = p_f * fruits * seeds * g_e * d_log_size,
family = "CD",
# probability of flowering
p_f = predict(vit_flwr,
newdata = tibble(log_size_prev = log_size_1, year_fct = yr),
type = 'response'),
fruits = predict(vit_fruits,
newdata = tibble(log_size_prev = log_size_1),
type = 'response'),
seeds = exp(coef(vit_seeds)[1]), #mean seeds per fruit
g_e = vit_germ_est,
data_list = vit_list_stoch_cf,
states = list(c('log_size', 'sdlg')),
uses_par_sets = TRUE,
par_set_indices = list(yr = c(2001, 2003))
) %>%
define_kernel(
name = "stay_sdlg",
formula = 0,
family = "DD",
states = list(c("sdlg")),
evict_cor = FALSE
) %>%
# Enter into post-seedling
define_kernel(
name = "leave_sdlg_yr",
formula = s_sdlg * G_z2_sdlg,
family = "DC",
s_sdlg = predict(vit_surv_sdlg,
newdata = tibble(year_fct = yr),
type = 'response'), #intercept from intercept-only glm
G_z2_sdlg = dt.scaled(
x = log_size_2,
df = size_sdlg_nu,
mean = size_sdlg_mu,
sd = size_sdlg_sd
),
size_sdlg_mu = coef(vit_size_sdlg)[1], #intercept from intercept-only glm
size_sdlg_nu = get_scat_params(vit_size_sdlg)["nu"],
size_sdlg_sd = get_scat_params(vit_size_sdlg)["sd"],
data_list = vit_list_stoch_cf,
states = list(c("sdlg", "log_size")),
uses_par_sets = TRUE,
par_set_indices = list(yr = c(2001, 2003)),
evict_cor = TRUE,
evict_fun = truncated_distributions("t.scaled", "G_z2_sdlg")
) %>%
# define implementation with midpoint rule
define_impl(make_impl_args_list(
kernel_names = c("P_yr", "go_sdlg_yr", "stay_sdlg", "leave_sdlg_yr"),
int_rule = "midpoint",
state_start = c("log_size", "log_size", "sdlg", "sdlg"),
state_end = c("log_size", "sdlg", "sdlg", "log_size")
)) %>%
#define lower bound, upper bound, and number of meshpoints for log_size
define_domains(log_size = c(0, 8.018296, 100)) %>%
#arbitrary starting population state
define_pop_state(pop_vectors = pop_vec_cf)Iterate IPM with 1000 iterations
year_seq = sample(c(2001, 2003), size = 1000, replace = TRUE)ipm_stoch_ff <-
proto_ipm_stoch_ff %>%
make_ipm(iterations = 1000,
normalize_pop_size = FALSE, # to run as PVA
kernel_seq = as.character(year_seq),
usr_funs = list(get_scat_params = get_scat_params))
ipm_stoch_cf <-
proto_ipm_stoch_cf %>%
make_ipm(iterations = 1000,
normalize_pop_size = FALSE, # to run as PVA
kernel_seq = as.character(year_seq),
usr_funs = list(get_scat_params = get_scat_params))6 Investigation
First, with just the worst and best years, is the stochastic lambda still higher than the deterministic?
all_lams <-
list(ff = ipm_stoch_ff, cf = ipm_stoch_cf) %>%
map_df(
~tibble(year = .x$env_seq,
lambda = lambda(.x, type_lambda = "all")[,1])%>%
#remove burnin
slice(101:n()),
.id = "habitat"
)tribble(~model, ~habitat, ~year, ~lambda,
"det", "ff", "both combined", lambda(ipm_det_ff),
"stoch", "ff", "both combined", lambda(ipm_stoch_ff, log = FALSE),
"det", "cf", "both combined", lambda(ipm_det_cf),
"stoch", "cf", "both combined", lambda(ipm_stoch_ff, log = FALSE),
"det", "ff", "2001", lambdas_ff[2],
"det", "ff", "2003", lambdas_ff[4],
"det", "cf", "2001", lambdas_cf[2],
"det", "cf", "2003", lambdas_cf[4]
) %>%
bind_rows(
all_lams %>%
group_by(habitat, year) %>%
summarize(lambda = mean(lambda)) %>%
mutate(year = paste(year, "mean"),
model = "stoch")
) %>% arrange(model, habitat)## `summarise()` has grouped output by 'habitat'. You can override using the
## `.groups` argument.## # A tibble: 12 × 4
## model habitat year lambda
## <chr> <chr> <chr> <dbl>
## 1 det cf both combined 0.953
## 2 det cf 2001 1.01
## 3 det cf 2003 0.881
## 4 det ff both combined 0.964
## 5 det ff 2001 1.01
## 6 det ff 2003 0.940
## 7 stoch cf both combined 0.965
## 8 stoch cf 2001 mean 0.973
## 9 stoch cf 2003 mean 0.952
## 10 stoch ff both combined 0.965
## 11 stoch ff 2001 mean 0.966
## 12 stoch ff 2003 mean 0.964Stochastic lambda is still slightly higher than deterministic
Within stochastic IPMs, the mean lambdas for individual transition years are consistent with good and bad years, but both are higher than deterministic lambda with both years combined!
The deterministic lambda for 2001 is a lot higher than the mean of the 2001 lambdas from the stochastic simulation and the deterministic lambda for 2003 is a lot lower than the mean of the 2003 lambdas from the stochastic simulation. This makes sense to me, as fitting year as a random effect rather than a fixed effect will shrink estimates toward a global mean.
How do they compare at the level of vital rates?
newdata <- expand.grid(
log_size_prev = seq(min(two_yrs$log_size_prev, na.rm = TRUE), max(two_yrs$log_size_prev, na.rm = TRUE), length.out = 50),
year_fct = as.character(c(2001, 2003))
)
safe_predict <- possibly(predict, NA)
#deterministic, just 2001
det_fitted_2001_ff <-
vit_list_yr_ff[[2]] %>%
map_df(~safe_predict(.x, newdata = newdata, type = "response")) %>%
bind_cols(newdata) %>%
select(-vit_fruits, -vit_seeds, -vit_germ_est) %>%
pivot_longer(starts_with("vit_"), names_to = "vital_rate", values_to = ".fitted", names_prefix = "vit_") %>%
mutate(year_fct = "det. 2001")
#deterministic, just 2003
det_fitted_2003_ff <-
vit_list_yr_ff[[4]] %>%
map_df(~safe_predict(.x, newdata = newdata, type = "response")) %>%
bind_cols(newdata) %>%
select(-vit_fruits, -vit_seeds, -vit_germ_est) %>%
pivot_longer(starts_with("vit_"), names_to = "vital_rate", values_to = ".fitted", names_prefix = "vit_") %>%
mutate(year_fct = "det. 2003")
#deterministic with both years combined
det_fitted_ff <-
vit_list_det_ff %>%
map_df(~safe_predict(.x, newdata = newdata, type = "response")) %>%
bind_cols(newdata) %>%
select(-vit_fruits, -vit_seeds, -vit_germ_est) %>%
pivot_longer(starts_with("vit_"), names_to = "vital_rate", values_to = ".fitted", names_prefix = "vit_") %>%
mutate(year_fct = "det. both years")
#stochastic
stoch_fitted_ff <-
vit_list_stoch_ff %>%
map_df(~safe_predict(.x, newdata= newdata, type = "response")) %>%
bind_cols(newdata) %>%
select(-vit_fruits, -vit_seeds, -vit_germ_est) %>%
pivot_longer(starts_with("vit_"), names_to = "vital_rate", values_to = ".fitted", names_prefix = "vit_") %>%
mutate(year_fct = paste("stoch.", year_fct))
bind_rows(stoch_fitted_ff, det_fitted_ff, det_fitted_2001_ff, det_fitted_2003_ff) %>%
ggplot(aes(x = log_size_prev, y = .fitted, color = year_fct, linetype = year_fct)) +
geom_line(alpha = 0.9) +
facet_wrap(~vital_rate, scales = "free") +
labs(title = "Forest Fragments")
newdata <- expand.grid(
log_size_prev = seq(min(two_yrs$log_size_prev, na.rm = TRUE), max(two_yrs$log_size_prev, na.rm = TRUE), length.out = 50),
year_fct = as.character(c(2001, 2003))
)
safe_predict <- possibly(predict, NA)
#deterministic, just 2001
det_fitted_2001_cf <-
vit_list_yr_cf[[2]] %>%
map_df(~safe_predict(.x, newdata = newdata, type = "response")) %>%
bind_cols(newdata) %>%
select(-vit_fruits, -vit_seeds, -vit_germ_est) %>%
pivot_longer(starts_with("vit_"), names_to = "vital_rate", values_to = ".fitted", names_prefix = "vit_") %>%
mutate(year_fct = "det. 2001")
#deterministic, just 2003
det_fitted_2003_cf <-
vit_list_yr_cf[[4]] %>%
map_df(~safe_predict(.x, newdata = newdata, type = "response")) %>%
bind_cols(newdata) %>%
select(-vit_fruits, -vit_seeds, -vit_germ_est) %>%
pivot_longer(starts_with("vit_"), names_to = "vital_rate", values_to = ".fitted", names_prefix = "vit_") %>%
mutate(year_fct = "det. 2003")
#deterministic with both years combined
det_fitted_cf <-
vit_list_det_cf %>%
map_df(~safe_predict(.x, newdata = newdata, type = "response")) %>%
bind_cols(newdata) %>%
select(-vit_fruits, -vit_seeds, -vit_germ_est) %>%
pivot_longer(starts_with("vit_"), names_to = "vital_rate", values_to = ".fitted", names_prefix = "vit_") %>%
mutate(year_fct = "det. both years")
#stochastic
stoch_fitted_cf <-
vit_list_stoch_cf %>%
map_df(~safe_predict(.x, newdata= newdata, type = "response")) %>%
bind_cols(newdata) %>%
select(-vit_fruits, -vit_seeds, -vit_germ_est) %>%
pivot_longer(starts_with("vit_"), names_to = "vital_rate", values_to = ".fitted", names_prefix = "vit_") %>%
mutate(year_fct = paste("stoch.", year_fct))
bind_rows(stoch_fitted_cf, det_fitted_cf, det_fitted_2001_cf, det_fitted_2003_cf) %>%
ggplot(aes(x = log_size_prev, y = .fitted, color = year_fct, linetype = year_fct)) +
geom_line(alpha = 0.9) +
facet_wrap(~vital_rate, scales = "free") +
labs(title = "Continuous Forest")
For the stochastic models, the two years are nearly indistinguishable for seedlings. I believe this is because of shrinkage—there is little information for either year so means are shrunk toward global average. You can see this more clearly for flowering probability where the order of lines at the larger end is 2001 deterministic (essentially a fixed effect), 2001 stochastic (random effect), both years combined, 2003 stochastic (random effect), and 2003 deterministic (fixed effect).