library(lifelihood)
#> Loading required package: tidyverse
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ dplyr 1.2.0 ✔ readr 2.2.0
#> ✔ forcats 1.0.1 ✔ stringr 1.6.0
#> ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
#> ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
#> ✔ purrr 1.2.1
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag() masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidyverse)
df <- datapierrick |>
mutate(
par = as.factor(par),
spore = as.factor(spore),
block = rep(1:2, each = nrow(datapierrick) / 2)
)
generate_clutch_vector <- function(N) {
return(paste(
"pon",
rep(c("start", "end", "size"), N),
rep(1:N, each = 3),
sep = "_"
))
}
clutchs <- generate_clutch_vector(28)
lifelihoodData <- as_lifelihoodData(
df = df,
sex = "sex",
sex_start = "sex_start",
sex_end = "sex_end",
maturity_start = "mat_start",
maturity_end = "mat_end",
clutchs = clutchs,
block = "block",
death_start = "death_start",
death_end = "death_end",
covariates = c("par", "spore"),
model_specs = c("wei", "lgn", "wei")
)
results <- lifelihood(
lifelihoodData,
path_config = use_test_config("config_pierrick")
)
#> [1] "/Users/runner/work/_temp/Library/lifelihood/bin/lifelihood-macos /Users/runner/work/Lifelihood/Lifelihood/lifelihood_7289_8302_5270_8561/temp_file_data_lifelihood.txt /Users/runner/work/Lifelihood/Lifelihood/lifelihood_7289_8302_5270_8561/temp_param_range_path.txt 0 25 FALSE 0 FALSE 0 7289 8302 5270 8561 10 20 1000 0.3 NULL 2 2 50 1 1 0.001"
#> Warning in check_estimation(results): Estimation of 'expt_death' is close to
#> the maximum bound: expt_death~=323.958894644731 (bound=322.38). Consider
#> increasing maximum bound.
summary(results)
#>
#> === LIFELIHOOD RESULTS ===
#>
#> Sample size: 550
#>
#> --- Model Fit ---
#> Log-likelihood: -32234.677
#> AIC: 64517.4
#> BIC: 64620.8
#>
#> --- Key Parameters ---
#>
#> Mortality:
#> expt_death (Intercept) -0.909 (0.000)
#> expt_death eff_expt_death_par_1 -0.084 (0.000)
#> expt_death eff_expt_death_par_2 -0.010 (0.000)
#> expt_death eff_expt_death_spore_1 -0.753 (0.000)
#> expt_death eff_expt_death_spore_2 -0.036 (0.000)
#> expt_death eff_expt_death_spore_3 -0.418 (0.000)
#> expt_death eff_expt_death_par_1:spore_1 0.008 (0.000)
#> expt_death eff_expt_death_par_2:spore_1 3.575 (0.000)
#> expt_death eff_expt_death_par_1:spore_2 -0.634 (0.000)
#> expt_death eff_expt_death_par_2:spore_2 -1.121 (0.000)
#> expt_death eff_expt_death_par_1:spore_3 9.736 (0.000)
#> expt_death eff_expt_death_par_2:spore_3 -0.382 (0.000)
#> survival_param2 (Intercept) -4.877 (0.000)
#>
#> Maturity:
#> expt_maturity (Intercept) -1.518 (0.000)
#> expt_maturity eff_expt_maturity_par_1 0.148 (0.000)
#> expt_maturity eff_expt_maturity_par_2 0.081 (0.000)
#> maturity_param2 (Intercept) -7.405 (0.000)
#>
#> Reproduction:
#> expt_reproduction (Intercept) -1.789 (0.000)
#> expt_reproduction eff_expt_reproduction_par_1 -1.696 (0.000)
#> expt_reproduction eff_expt_reproduction_par_2 -1.103 (0.000)
#> reproduction_param2 (Intercept) -1.122 (0.000)
#> reproduction_param2 eff_reproduction_param2_par_1 -0.439 (0.000)
#> reproduction_param2 eff_reproduction_param2_par_2 -0.800 (0.000)
#> n_offspring (Intercept) -2.543 (0.000)
#>
#> --- Convergence ---
#> All parameters within bounds
#>
#> ======================lifelihood offers a lifelihood::simulate_life_history() function that lets you simulate new observations based on estimates you made.
Fitting
First, we need to fit the model with lifelihood::lifelihood():
Default simulations
By default, lifelihood will simulate all life history events (maturity, reproduction, and death):
simulate_life_history(results) |> head()
#> # A tibble: 6 × 149
#> mortality maturity clutch_1 n_offspring_clutch_1 clutch_2 n_offspring_clutch_2
#> <dbl> <dbl> <dbl> <int> <dbl> <int>
#> 1 149. 12.8 15.2 3 20.6 3
#> 2 98.7 12.0 14.7 8 22.7 11
#> 3 84.2 12.6 22.3 3 24.0 2
#> 4 135. 12.7 18.1 5 24.1 3
#> 5 129. 12.3 15.3 3 20.0 5
#> 6 67.8 13.4 20.7 4 22.3 2
#> # ℹ 143 more variables: clutch_3 <dbl>, n_offspring_clutch_3 <int>,
#> # clutch_4 <dbl>, n_offspring_clutch_4 <int>, clutch_5 <dbl>,
#> # n_offspring_clutch_5 <int>, clutch_6 <dbl>, n_offspring_clutch_6 <int>,
#> # clutch_7 <dbl>, n_offspring_clutch_7 <int>, clutch_8 <dbl>,
#> # n_offspring_clutch_8 <int>, clutch_9 <dbl>, n_offspring_clutch_9 <int>,
#> # clutch_10 <dbl>, n_offspring_clutch_10 <int>, clutch_11 <dbl>,
#> # n_offspring_clutch_11 <int>, clutch_12 <dbl>, …But you can specify which event you want:
simulate_life_history(results, event = "maturity") |> head()
#> # A tibble: 6 × 1
#> maturity
#> <dbl>
#> 1 14.4
#> 2 12.4
#> 3 12.3
#> 4 14.5
#> 5 12.0
#> 6 13.0Simulations with visit masks
lifelihood lets you specify visit masks that are used to simulate data that more closely reflects how the original data was measured by adding constraints to the interval dates.
If you scroll to the top, you’ll see that we passed a column name from our dataframe to the block argument. This column represents the block to which each individual belongs. Since we provided this value, we just need to add use_censoring = TRUE to include censoring time intervals in the simulation:
results |>
simulate_life_history(event = "maturity", use_censoring = TRUE) |>
head()
#> # A tibble: 6 × 3
#> maturity maturity_start maturity_end
#> <dbl> <int> <int>
#> 1 12.5 12 13
#> 2 12.3 12 13
#> 3 13.5 13 14
#> 4 11.9 11 12
#> 5 13.1 13 14
#> 6 13.3 13 14Right now, visit masks are “deduced” from the original dataset, but you can provide your own visit masks with the visits argument.
It must be a dataframe with 2 columns: block (the same name as passed to lifelihood::as_lifelihoodData() in the block argument) and exactly visit. For each block, visit corresponds to the ages at which the events of individuals were recorded.
Let’s define this visits dataframe:
Now we can pass this to the simulate_life_history() function:
results |>
simulate_life_history(
event = "maturity",
use_censoring = TRUE,
visits = visits
) |>
head()
#> # A tibble: 6 × 3
#> maturity maturity_start maturity_end
#> <dbl> <int> <int>
#> 1 12.9 11 13
#> 2 12.6 11 13
#> 3 13.4 13 15
#> 4 12.3 11 13
#> 5 13.0 11 13
#> 6 12.8 11 13