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)Customize parameter boundaries and estimation warning
Source:vignettes/customize-parameter-boundaries-and-estimation-warning.qmd
- What is the required data format to work with lifelihood?
- Setting up the configuration file
- How to use the lifelihood package
Load libraries
Data preparation
Load the dataset from .csv file:
df |> head()
#> type geno sex_start sex_end sex mat_start mat_end clutch_start1 clutch_end1
#> 1 1 0 0 1000 0 0 1000 NA NA
#> 2 2 0 0 1000 0 2 3 4 3
#> 3 0 1 0 1000 0 3 4 2 4
#> 4 0 1 0 1000 0 3 4 2 4
#> 5 0 1 0 1000 0 3 4 2 4
#> 6 1 0 0 1000 0 0 1000 NA NA
#> clutch_size1 clutch_start2 clutch_end2 clutch_size2 death_start death_end
#> 1 NA NA NA NA 0.1 1
#> 2 4 2 4 5 9.0 10
#> 3 5 NA NA NA 5.0 6
#> 4 5 NA NA NA 5.0 6
#> 5 5 NA NA NA 5.0 6
#> 6 NA NA NA NA 0.1 1Prepare input parameters for the lifelihood() function:
# name of the columns of the clutchs into a single vector
clutchs <- c(
"clutch_start1",
"clutch_end1",
"clutch_size1",
"clutch_start2",
"clutch_end2",
"clutch_size2"
)Note: If you have a large number of clutches, it is easier to generate this vector programmatically, particularly if your dataset contains a large number of clutches. See the Generate clutch names vignette.
Create the lifelihoodData object
dataLFH <- as_lifelihoodData(
df = df,
sex = "sex",
sex_start = "sex_start",
sex_end = "sex_end",
maturity_start = "mat_start",
maturity_end = "mat_end",
clutchs = clutchs,
death_start = "death_start",
death_end = "death_end",
covariates = c("geno", "type"),
model_specs = c("gam", "lgn", "wei")
)Get the results
Let’s run the analysis with default parameters.
results <- lifelihood(
lifelihoodData = dataLFH,
path_config = use_test_config("config")
)
#> [1] "/Users/runner/work/_temp/Library/lifelihood/bin/lifelihood-macos /Users/runner/work/Lifelihood/Lifelihood/lifelihood_3450_719_2309_3733/temp_file_data_lifelihood.txt /Users/runner/work/Lifelihood/Lifelihood/lifelihood_3450_719_2309_3733/temp_param_range_path.txt 0 25 FALSE 0 FALSE 0 3450 719 2309 3733 10 20 1000 0.3 NULL 2 2 50 1 1 0.001"Warning
What is it?
When runnning lifelihood() function, you might encounter the following warning:
## Warning in check_estimation(results_lifelihood = results): Estimation of
## 'increase_death_hazard' is close to the maximum bound:
## increase_death_hazard≃9.9864992332278. Consider increasing maximum bound.This warning indicates that the estimation of the increase_death_hazard parameter is close to the maximum bound.
This is not an error, but a signal that the estimation of the increase_death_hazard parameter is approaching its upper limit. During the optimization process, the algorithm may have been constrained by the current parameter boundaries, potentially affecting the accuracy of the estimation.
To address this warning, you can consider increasing the maximum bound of the increase_death_hazard parameter. This will give the optimization algorithm more flexibility to find the best estimate.
Customize parameter boundaries
You can get the parameter boundaries with the default_bounds_df() function and by passing the lifelihoodData object:
bounds_df <- default_bounds_df(dataLFH)
bounds_df
#> param min max
#> 1 expt_death 0.001 40
#> 2 survival_param2 0.05 500
#> 3 ratio_expt_death 0.01 100
#> 4 prob_death 1e-05 0.99999
#> 5 sex_ratio 1e-05 0.99999
#> 6 expt_maturity 0.001 8
#> 7 maturity_param2 0.005 600
#> 8 ratio_expt_maturity 0.01 100
#> 9 expt_reproduction 0.001 10
#> 10 reproduction_param2 0.0025 10
#> 11 n_offspring 1 50
#> 12 increase_death_hazard 1e-05 10
#> 13 tof_reduction_rate 1e-07 10
#> 14 increase_tof_n_offspring 1e-07 10
#> 15 lin_decrease_hazard -20 20
#> 16 quad_senescence -20 20
#> 17 quad_decrease_hazard -10 10
#> 18 quad_change_n_offspring -10 10
#> 19 tof_n_offspring -10 10
#> 20 fitness 0.001 1000Since 10 seems to not be high enough, let’s try with 80:
bounds_df[bounds_df$name == "increase_death_hazard", "max"] <- 80Once it’s changed, you just have to call lifelihood() again with the param_bounds_df argument:
results <- lifelihood(
lifelihoodData = dataLFH,
path_config = use_test_config("config"),
param_bounds_df = bounds_df
)
#> [1] "/Users/runner/work/_temp/Library/lifelihood/bin/lifelihood-macos /Users/runner/work/Lifelihood/Lifelihood/lifelihood_6413_92_8212_8973/temp_file_data_lifelihood.txt /Users/runner/work/Lifelihood/Lifelihood/lifelihood_6413_92_8212_8973/temp_param_range_path.txt 0 25 FALSE 0 FALSE 0 6413 92 8212 8973 10 20 1000 0.3 NULL 2 2 50 1 1 0.001"Now we don’t get any warning!