Vignette for simulation analysis
08 March, 2024
Setup
library(lubridate)
library(terra)
library(raster)
library(tidyverse)
library(tidyterra)
library(conflicted)
library(ggplot2)
library(glmmTMB)
library(amt) # dev version needed: remotes::install_github("jmsigner/amt")
library(cowplot)
library(MASS)
library(mixedSSA)
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")
set.seed(11102023)
n <- 30 ### Number of individuals to simulate Introduction
We simulate trajectories following Integrated step-selection
analysis (ISSA) (Equation 1) for a population with
individual-specific, known parameter values to demonstrate the
application of mixed ISSAs that allow for individual variability in
movement parameters (Figure 1). Individual tracks were simulated using
the amt package. We also compare the estimates from the
mixed ISSA to those obtained by fitting a separate model to true
individual-specific parameters from each individual. This comparison
highlights the benefits of using random effects in ISSA; namely; the
ability to borrow information across like individuals, which leads to
improved estimates of individual-specific parameters.
Raster Generation
We generated a categorical raster with two classes (Habitat A, Habitat B) where the individual selects for the habitat B and moves differently in the two habitats (Section 3.1).
tc2 <- raster(nrows=1000, ncols=1000, xmn=-500, xmx=500, ymn=-500,ymx=500)
tc2[] <- runif(1000000, 0,50) ## assigning simulated values to the raster
# producing a new raster with moving window of size 11x13 and calculating mean values for the neighborhood
tc2 <- raster::focal(tc2, w = matrix(1, 11, 13), mean)
tc2 <- raster::focal(tc2, w = matrix(1, 11, 13), mean)
# we did the focal operation twice to produce a smoothed out edges of the interspersed habitats
reclass_df <-matrix(c(12,24.5, 1, 24.5,38, 0), ncol=3, byrow=T)
tc_reclass <- raster::reclassify(tc2, reclass_df)
r <- rast(tc_reclass)
names(r) <- "x1"
r <- as.factor(r)
cls <- data.frame(id=c(0,1), x1=c("Habitat A", "Habitat B"))
levels(r) <- cls
plot(r)
#### Initializing
sig <-matrix(c(0.25,0.15,0.15, 0.15, 0.2, 0.15, 0.15,0.15, 0.25), ncol =3, byrow = T)
param.value <-data.frame(id = NA, shape = NA, scale = NA, kappa = NA)
ssf_sim <-list()Simulating tracks using the amt functions
The simulation of movement tracks in amt requires three
steps: 1. Specifying parameters governing the movements for the
individual using the make_issf_model function
2. Specifying a redistribution kernel for generating movements based on
an initial location and bearing
3. Generating paths from the simulated redistribution kernel specified
for a user-provided number of steps
Warning: The simulation takes a while to run (a few hours to days) and requires amt version 0.2.2. Please check the version and compatibility before the simulation run.
for(i in 1:n){
start <- make_start(c(runif(1,-20,20), runif(1,-20,20)), ta_ = 0)
param.value[i,1] <- i
param.value[i,2:4] <- mvrnorm(n=1, mu=c(3,1,4), Sigma = sig)
## simulating initial parameters using multivariate normal distribution
m <- make_issf_model(coefs = c(x1_end = 0.05,
"sl_" = 0, "log(sl_)" = 0, "cos(ta_)" = 0,
"sl_:x1_start" = 0.05, "log(sl_):x1_start" = 0.05,"cos(ta_):x1_start" = 0.05),
sl = make_gamma_distr(shape = param.value[i,2],
scale = param.value[i,3]
),
ta = make_vonmises_distr(kappa = param.value[i,4]
))
rdk <- redistribution_kernel(m, start = start, map = r, as.rast= FALSE,
landscape= "continuous")
ssf_sim[[i]] <- simulate_path(rdk, n.steps = floor(rnorm(1,2000,250))) ### simulate individual trajectories
ssf_sim[[i]]$id <-i
}head(param.value)## id shape scale kappa
## 1 1 3.755938 1.147982 4.667033
## 2 2 2.771625 1.019327 4.009808
## 3 3 3.403455 1.361682 4.521502
## 4 4 3.136668 1.301838 3.739852
## 5 5 3.087484 1.214830 3.736779
## 6 6 1.968182 0.482264 3.317193ssf_simdat <-do.call(rbind, ssf_sim) ### binding individual tracks into one datasetMaking tracks from the simulated trajectories:
ssf_simtr <- ssf_simdat |> make_track(x_, y_, t_ , id=id)Visualization of the tracks:
ssf_simtr_sample <- ssf_simtr |> filter(id %in% c(1:4))
ggplot() + geom_spatraster(data= r, show.legend = FALSE) +
geom_path(data= ssf_simtr_sample, aes(x1_, y1_), col= "black")+
facet_wrap(~ id, nrow =2)+ scale_fill_terrain_d() + theme_bw()
Individual Model fitting
Generating random steps from the tracks
ssfdat_sim_track_ind <- ssf_simtr %>% nest(data = -id) %>%
mutate( data = map(data, ~ .x %>%
steps() %>%
random_steps() %>%
extract_covariates(r, where = "both"))) %>%
unnest(cols = data)glimpse(ssfdat_sim_track_ind)## Rows: 691,933
## Columns: 14
## $ id <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 1…
## $ x1_ <dbl> -3.0809728, -3.0809728, -3.0809728, -3.0809728, -3.0809728, -…
## $ x2_ <dbl> -2.8117514, -1.9647963, -2.8893406, -1.8360494, -2.8686341, -…
## $ y1_ <dbl> -11.82080, -11.82080, -11.82080, -11.82080, -11.82080, -11.82…
## $ y2_ <dbl> -11.77324, -12.83395, -11.46173, -12.02318, -12.05308, -11.65…
## $ sl_ <dbl> 0.2733908, 1.5074243, 0.4070025, 1.2612660, 0.3147076, 1.8526…
## $ ta_ <dbl> 0.0883989278, -0.8235224352, 0.9940861008, -0.2476242936, -0.…
## $ t1_ <dttm> 2023-10-11 16:59:22, 2023-10-11 16:59:22, 2023-10-11 16:59:2…
## $ t2_ <dttm> 2023-10-11 17:59:22, 2023-10-11 17:59:22, 2023-10-11 17:59:2…
## $ dt_ <drtn> 1 hours, 1 hours, 1 hours, 1 hours, 1 hours, 1 hours, 1 hour…
## $ case_ <lgl> TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
## $ step_id_ <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4…
## $ x1_start <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ x1_end <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…### Fitting Individual model
ssffits <- ssfdat_sim_track_ind %>% nest(data=-id) %>%
mutate(mod = map(data, function(x) (try(fit_issf(case_ ~ sl_ + log(sl_)+ cos(ta_)+ x1_end +
sl_:x1_start + log(sl_):x1_start+ cos(ta_):x1_start+
strata(step_id_), data = x)))))
### Updating individual modelWe estimate the tentative individual movement parameters for each individual and then update the movement parameters using the update_gamma and update_vonmises functions of the amt package. In this vignette, we will be focusing on the movement parameters when the individual is in the reference habitat class (Habitat A).
for(i in 1:n){
indv_dat <- subset(ssfdat_sim_track_ind, ssfdat_sim_track_ind$id == i)
sl_fit <- amt::fit_distr(indv_dat$sl_[indv_dat$case_=="TRUE"], "gamma", na.rm = T)
ta_fit <- amt::fit_distr(indv_dat$ta_[indv_dat$case_=="TRUE"], "vonmises", na.rm = T)
tab <- broom::tidy(ssffits$mod[[i]]$model)
param.value[i, 5:6] <- as.data.frame(update_gamma(dist= sl_fit,
beta_sl = tab$estimate[tab$term == "sl_"] ,
beta_log_sl = tab$estimate[tab$term == "log(sl_)"])$params)[1,]
param.value[i, 7:8] <- as.data.frame(update_vonmises(dist= ta_fit,
beta_cos_ta = tab$estimate[tab$term == "cos(ta_)"])$params)[1,]
}
names(param.value)[5:8] <-c("shape.indv", "scale.indv", "kappa.indv", "mu.indv")
head(param.value)## id shape scale kappa shape.indv scale.indv kappa.indv mu.indv
## 1 1 3.755938 1.147982 4.667033 3.627582 1.2611523 4.744660 0
## 2 2 2.771625 1.019327 4.009808 2.764352 1.1300422 3.883596 0
## 3 3 3.403455 1.361682 4.521502 3.819950 1.1367329 4.808811 0
## 4 4 3.136668 1.301838 3.739852 3.243939 1.2286589 4.286453 0
## 5 5 3.087484 1.214830 3.736779 3.375905 1.0370127 3.746843 0
## 6 6 1.968182 0.482264 3.317193 1.854183 0.5113274 3.313954 0Mixed-effect Model fitting
Generating random steps from the tracks with one tentative distribution
The step to generate random steps is similar to the individual model but we are using one tentative distribution to generate the random steps for all individuals.
ssfdat_sim_track <- ssf_simtr %>% nest(data = -id) %>%
mutate( data = map(data, ~ .x %>%
steps())) %>%
unnest(cols = data) %>%
random_steps() %>%
extract_covariates(r, where = "both")
ssfdat_sim_track$step_id <- with(ssfdat_sim_track, paste(id, step_id_, sep = "_"))glimpse(ssfdat_sim_track)## Rows: 692,252
## Columns: 15
## $ id <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 1…
## $ x1_ <dbl> -3.0809728, -3.0809728, -3.0809728, -3.0809728, -3.0809728, -…
## $ x2_ <dbl> -2.8117514, -1.9392788, 0.8229708, -2.1675012, -2.6620170, 2.…
## $ y1_ <dbl> -11.82080, -11.82080, -11.82080, -11.82080, -11.82080, -11.82…
## $ y2_ <dbl> -11.773236, -11.105575, -8.735172, -11.683578, -11.812017, -1…
## $ sl_ <dbl> 0.2733908, 1.3472238, 4.9761305, 0.9237207, 0.4190479, 6.1793…
## $ ta_ <dbl> 0.08839893, 0.47317822, 0.58238216, 0.06263479, -0.06551059, …
## $ t1_ <dttm> 2023-10-11 16:59:22, 2023-10-11 16:59:22, 2023-10-11 16:59:2…
## $ t2_ <dttm> 2023-10-11 17:59:22, 2023-10-11 17:59:22, 2023-10-11 17:59:2…
## $ dt_ <drtn> 1 hours, 1 hours, 1 hours, 1 hours, 1 hours, 1 hours, 1 hour…
## $ case_ <lgl> TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
## $ step_id_ <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4…
## $ x1_start <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ x1_end <dbl> 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ step_id <chr> "11_3", "11_3", "11_3", "11_3", "11_3", "11_3", "11_3", "11_3…#saveRDS(ssfdat_sim_track, "Mixed_simdata_final.rds")Fitting the mixed-effect model using glmmTMB package
following Muff et al. (2020):
The code repository can be found at https://conservancy.umn.edu/handle/11299/204737
ssf_sim_model <- glmmTMB(case_ ~ sl_ + log(sl_)+ cos(ta_)+ x1_end +
sl_:x1_start + log(sl_):x1_start + cos(ta_):x1_start +
(1|step_id)+ (0+ sl_ +(log(sl_)) + cos(ta_)|id) +(0+ x1_end|id),
REML = TRUE, family=poisson(), data =ssfdat_sim_track, doFit=FALSE)
## we have to manually fix the variance of the intercept first and then fix the
## standard deviation of the first random term, which is the (1|step_id) component
## in the above model equation that is fixed as 1e3 here.
ssf_sim_model$parameters$theta[1] <- log(1e3)
## We need to tell glmmTMB not to change the variance by setting it to NA:
nvarparm7 <- length(ssf_sim_model$parameters$theta)
ssf_sim_model$mapArg <- list(theta=factor(c(NA, 1:(nvarparm7-1))))
## Then fit the model and look at the results:
snapper.ssf_sim<- glmmTMB:::fitTMB(ssf_sim_model)
#saveRDS(snapper.ssf_sim, "population_model_final.rds")summary(snapper.ssf_sim)## Family: poisson ( log )
## Formula:
## case_ ~ sl_ + log(sl_) + cos(ta_) + x1_end + sl_:x1_start + log(sl_):x1_start +
## cos(ta_):x1_start + (1 | step_id) + (0 + sl_ + (log(sl_)) +
## cos(ta_) | id) + (0 + x1_end | id)
## Data: ssfdat_sim_track
##
## AIC BIC logLik deviance df.resid
## 1000141.0 1000309.3 -500055.5 1000111.0 549431
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## step_id (Intercept) 1.000e+06 1.000e+03
## id sl_ 1.180e+00 1.086e+00
## log(sl_) 2.313e-01 4.810e-01 0.41
## cos(ta_) 2.216e-01 4.708e-01 0.58 0.46
## id.1 x1_end 1.206e-07 3.473e-04
## Number of obs: 549438, groups: step_id, 49956; id, 30
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.76059 4.47429 -0.393 0.69396
## sl_ -0.81092 0.19890 -4.077 4.56e-05 ***
## log(sl_) 1.36604 0.09160 14.913 < 2e-16 ***
## cos(ta_) 0.21756 0.09230 2.357 0.01842 *
## x1_end 0.01806 0.02331 0.775 0.43834
## sl_:x1_start 0.03094 0.01172 2.641 0.00827 **
## log(sl_):x1_start 0.02124 0.02623 0.810 0.41803
## cos(ta_):x1_start -0.08419 0.04691 -1.795 0.07271 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1After fitting the mixed-effect model using the glmmTMB package, we get the model output into the mixedSSA package for updating the movement parameters. First, we update the parameters in the step-length distribution and assume they follow a gamma distribution.
updated_sl_parameters <- mixedSSA::update_dist(model = snapper.ssf_sim,
dist_name = "gamma", # name of the step-length distribution assumed
beta_sl = "sl_", # the name of the step lengths coefficient in our model
beta_log_sl = "log(sl_)", # the name of the log(sl) coefficient in our model
interaction_var_name = "x1_start", # name of the habitat variable layer
random_effects_var_name = "id", # for individual update
quantiles = FALSE # for continuous interaction variables
)
# Extract the updated movement parameters for the reference habitat class and
# not include the first row of tentative movement parameters
param.value[,9:10] <- updated_sl_parameters@updated_parameters[-1,5:6] Then, we update the parameters in the turn-angle distribution assuming a von Mises distribution.
updated_ta_parameters <- mixedSSA::update_dist(model = snapper.ssf_sim,
dist_name = "vonmises", # name of the turn-angle distribution assumed
beta_cos_ta = "cos(ta_)", # the name of the cos(ta) coefficient in our model
interaction_var_name = "x1_start", # name of the habitat variable layer
random_effects_var_name = "id", # for individual update
quantiles = FALSE # for continuous interaction variables
)
# Extract the updated movement parameters for the reference habitat class and
# not include the first row of tentative movement parameters
param.value[,11:12] <- updated_ta_parameters@updated_parameters[-1,4:5]
names(param.value)[9:12] <-c("shape.all", "scale.all", "kappa.all", "mu.all") head(param.value)## id shape scale kappa shape.indv scale.indv kappa.indv mu.indv
## 1 1 3.755938 1.147982 4.667033 3.549135 1.2386194 4.880021 0
## 2 2 2.771625 1.019327 4.009808 2.803302 1.0278879 3.903209 0
## 3 3 3.403455 1.361682 4.521502 3.883147 1.1499781 4.656608 0
## 4 4 3.136668 1.301838 3.739852 3.182499 1.2884090 4.060282 0
## 5 5 3.087484 1.214830 3.736779 3.283903 1.1088030 3.620669 0
## 6 6 1.968182 0.482264 3.317193 1.854183 0.5113274 3.313954 0
## shape.all scale.all kappa.all mu.all
## 1 3.710071 1.1337195 4.694428 0
## 2 2.783951 1.0142899 3.974995 0
## 3 3.568963 1.2646067 4.555025 0
## 4 3.221643 1.2140585 4.026107 0
## 5 3.242791 1.1269102 3.799401 0
## 6 1.924404 0.4980171 3.403450 0Plotting for model comparison
We calculated the deviations between each estimated parameter and the parameter value used to simulate the data. We compared the deviation of individual and mixed-effect model to evaluate their relative performance.
param.value1 <- param.value |> mutate(indv_shape_bias = shape - shape.indv,
mixed_shape_bias = shape - shape.all,
indv_scale_bias = scale - scale.indv,
mixed_scale_bias = scale - scale.all,
indv_kappa_bias = kappa - kappa.indv,
mixed_kappa_bias = kappa - kappa.all)
shp_bias <- param.value1 |>
pivot_longer(cols = c(indv_shape_bias, mixed_shape_bias), names_to = "shapediff", values_to = "shape_value")|>
ggplot()+ geom_boxplot(aes(y= shape_value, x = shapediff, fill = shapediff), outlier.shape = NA)+
geom_jitter(aes(y= shape_value, x = shapediff), size = 0.6)+
scale_x_discrete(labels=c('Individual\nModel', 'Mixed-effect\nModel'))+
scale_fill_manual(breaks = c("indv_shape_bias", "mixed_shape_bias"),
values = c("indv_shape_bias" = "#ab1489",
"mixed_shape_bias"= "#69b3a2"))+
labs(y= "Shape Parameter deviation", x= "Model type")+ theme_bw()+
theme(legend.position = "none",
axis.title = element_text(size=16))
scl_bias <- param.value1 |>
pivot_longer(cols = c(indv_scale_bias, mixed_scale_bias), names_to = "scalediff", values_to = "scale_value")|>
ggplot()+ geom_boxplot(aes(y= scale_value, x = scalediff, fill = scalediff), outlier.shape = NA)+
geom_jitter(aes(y= scale_value, x = scalediff), size = 0.6)+
scale_x_discrete(labels=c('Individual\nModel', 'Mixed-effect\nModel'))+
scale_fill_manual(breaks = c("indv_scale_bias", "mixed_scale_bias"),
values = c("indv_scale_bias" = "#ab1489",
"mixed_scale_bias"= "#69b3a2"))+
labs(y= "Scale Parameter deviation", x= "Model type")+ theme_bw()+
theme(legend.position = "none",
axis.title = element_text(size=16))
kap_bias <- param.value1 |>
pivot_longer(cols = c(indv_kappa_bias, mixed_kappa_bias), names_to = "kappadiff", values_to = "kappa_value") |>
ggplot()+ geom_boxplot(aes(y= kappa_value, x = kappadiff, fill = kappadiff), outlier.shape = NA)+
geom_jitter(aes(y= kappa_value, x = kappadiff), size = 0.6)+
scale_x_discrete(labels=c('Individual\nModel', 'Mixed-effect\nModel'))+
scale_fill_manual(breaks = c("indv_kappa_bias", "mixed_kappa_bias"),
values = c("indv_kappa_bias" = "#ab1489",
"mixed_kappa_bias"= "#69b3a2"))+
labs(y= "Kappa Parameter deviation", x= "Model type")+ theme_bw()+
theme(legend.position = "none",
axis.title = element_text(size=16))
plot_grid(shp_bias, scl_bias, kap_bias, ncol=3, labels = c("A", "B", "C"))
sessionInfo()## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_India.utf8 LC_CTYPE=English_India.utf8
## [3] LC_MONETARY=English_India.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_India.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] mixedSSA_1.0.0 MASS_7.3-57 cowplot_1.1.1 amt_0.2.1.0
## [5] glmmTMB_1.1.7 conflicted_1.1.0 tidyterra_0.4.0 forcats_0.5.1
## [9] stringr_1.5.0 dplyr_1.0.9 purrr_0.3.4 readr_2.1.2
## [13] tidyr_1.2.0 tibble_3.2.1 ggplot2_3.4.4 tidyverse_1.3.2
## [17] raster_3.5-21 sp_1.5-0 terra_1.7-39 lubridate_1.8.0
##
## loaded via a namespace (and not attached):
## [1] assertive.base_0.0-9 TH.data_1.0-10
## [3] googledrive_2.0.0 minqa_1.2.4
## [5] colorspace_2.1-0 ellipsis_0.3.2
## [7] class_7.3-20 estimability_1.4.1
## [9] fs_1.5.2 rstudioapi_0.13
## [11] proxy_0.4-27 circular_0.4-95
## [13] farver_2.1.1 assertive.sets_0.0-3
## [15] hash_2.2.6.3 fansi_1.0.4
## [17] mvtnorm_1.1-3 xml2_1.3.3
## [19] assertive.data.uk_0.0-2 codetools_0.2-18
## [21] splines_4.2.1 cachem_1.0.6
## [23] knitr_1.42 jsonlite_1.8.8
## [25] assertive_0.3-6 nloptr_1.2.2.2
## [27] assertive.data.us_0.0-2 broom_1.0.0
## [29] dbplyr_2.2.1 compiler_4.2.1
## [31] httr_1.4.3 emmeans_1.8.0
## [33] backports_1.4.1 assertthat_0.2.1
## [35] Matrix_1.5-4.1 fastmap_1.1.0
## [37] gargle_1.2.0 cli_3.6.0
## [39] htmltools_0.5.5 tools_4.2.1
## [41] coda_0.19-4 gtable_0.3.3
## [43] glue_1.6.2 Rcpp_1.0.11
## [45] cellranger_1.1.0 jquerylib_0.1.4
## [47] vctrs_0.6.5 nlme_3.1-159
## [49] assertive.models_0.0-2 assertive.files_0.0-2
## [51] assertive.datetimes_0.0-3 xfun_0.38
## [53] rbibutils_2.2.16 lme4_1.1-30
## [55] rvest_1.0.2 lifecycle_1.0.4
## [57] googlesheets4_1.0.0 zoo_1.8-10
## [59] scales_1.2.1 hms_1.1.1
## [61] sandwich_3.0-0 TMB_1.9.4
## [63] assertive.matrices_0.0-2 assertive.strings_0.0-3
## [65] yaml_2.3.5 memoise_2.0.1
## [67] sass_0.4.5 stringi_1.7.8
## [69] highr_0.10 checkmate_2.1.0
## [71] e1071_1.7-11 boot_1.3-28
## [73] Rdpack_2.6 rlang_1.1.2
## [75] pkgconfig_2.0.3 assertive.data_0.0-3
## [77] evaluate_0.23 lattice_0.20-45
## [79] sf_1.0-13 labeling_0.4.2
## [81] assertive.properties_0.0-5 tidyselect_1.2.0
## [83] assertive.code_0.0-4 magrittr_2.0.3
## [85] R6_2.5.1 generics_0.1.3
## [87] multcomp_1.4-15 DBI_1.1.3
## [89] pillar_1.9.0 haven_2.5.0
## [91] withr_2.5.0 berryFunctions_1.22.0
## [93] fitdistrplus_1.1-8 assertive.numbers_0.0-2
## [95] units_0.8-0 abind_1.4-5
## [97] survival_3.3-1 modelr_0.1.8
## [99] crayon_1.5.1 assertive.types_0.0-3
## [101] KernSmooth_2.23-20 utf8_1.2.3
## [103] tzdb_0.3.0 rmarkdown_2.21
## [105] grid_4.2.1 readxl_1.4.0
## [107] data.table_1.14.2 reprex_2.0.1
## [109] digest_0.6.33 classInt_0.4-7
## [111] xtable_1.8-4 numDeriv_2016.8-1.1
## [113] munsell_0.5.0 assertive.reflection_0.0-5
## [115] bslib_0.4.2