library(MortalMove)
library(cmdstanr)
#> This is cmdstanr version 0.9.0.9000
#> - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
#> - Use set_cmdstan_path() to set the path to CmdStan
#> - Use install_cmdstan() to install CmdStan
set_cmdstan_path("YOUR_CMDSTAN_PATH") # Set your cmdstan path here
#> Warning: Path not set. Can't find directory: YOUR_CMDSTAN_PATH
Simulating data
sim <- simulate_data(n_animals = 100, n_fixes = 200, n_dead = 20, n_knots = 25)
sim[["raw_data"]] <- NULL
sim[["ind_cell_effect"]] <- 0
Fit stan model using cmdstan
model_file <- system.file("stan/mortality_model.stan", package = "MortalMove")
mod <- cmdstan_model(model_file, exe_file = system.file("stan/mortality_model.exe", package = "MortalMove"))
fit <- mod$sample(
data = sim,
chains = 2,
parallel_chains = 2,
iter_warmup = 500,
iter_sampling = 1500,
seed = 123
)
Check diagnostics
diagnostic_fit(fit, plot = TRUE, observed_y = rowSums(sim$time_step, na.rm = TRUE), delta = sim$delta)
#> This is posterior version 1.6.1
#>
#> Attaching package: 'posterior'
#> The following object is masked from 'package:FNN':
#>
#> entropy
#> The following objects are masked from 'package:stats':
#>
#> mad, sd, var
#> The following objects are masked from 'package:base':
#>
#> %in%, match
#> This is bayesplot version 1.13.0
#> - Online documentation and vignettes at mc-stan.org/bayesplot
#> - bayesplot theme set to bayesplot::theme_default()
#> * Does _not_ affect other ggplot2 plots
#> * See ?bayesplot_theme_set for details on theme setting
#>
#> Attaching package: 'bayesplot'
#> The following object is masked from 'package:posterior':
#>
#> rhat
#> This is loo version 2.8.0
#> - Online documentation and vignettes at mc-stan.org/loo
#> - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
#> # A tibble: 5 × 5
#> variable mean sd q5 q95
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 beta[1] 0.143 0.0508 0.0624 0.230
#> 2 beta[2] -0.0228 0.421 -0.734 0.651
#> 3 alpha[1] -0.0146 0.00362 -0.0205 -0.00865
#> 4 alpha[2] 0.176 0.0828 0.0424 0.317
#> 5 llambda -2.27 1.61 -4.92 0.329
#>
#> --- Rhat and ESS Summary ---
#> # A tibble: 5 × 4
#> variable rhat ess_bulk ess_tail
#> <chr> <dbl> <dbl> <dbl>
#> 1 beta[1] 1.00 2065. 1799.
#> 2 beta[2] 1.00 1725. 1577.
#> 3 alpha[1] 1.000 2202. 1880.
#> 4 alpha[2] 1.00 2153. 1706.
#> 5 llambda 1.00 2160. 1913.
#>
#> Computed from 3000 by 100 log-likelihood matrix.
#>
#> Estimate SE
#> elpd_loo -139.7 25.6
#> p_loo 2.8 0.6
#> looic 279.4 51.3
#> ------
#> MCSE of elpd_loo is 0.0.
#> MCSE and ESS estimates assume independent draws (r_eff=1).
#>
#> All Pareto k estimates are good (k < 0.7).
#> See help('pareto-k-diagnostic') for details.
#>
#> Posterior predictive p-value summary:
#> variable p_value
#> 1 Failed 0.3000
#> 2 Censored 0.7625
bayesplot::mcmc_trace(fit$draws(c("llambda","alpha","beta")))
Fit with spatial autocorrelation
# Fit with spatial autocorrelation
sim[["ind_cell_effect"]] <- 1
model_file <- system.file("stan/mortality_model.stan", package = "MortalMove")
mod <- cmdstan_model(model_file, exe_file = system.file("stan/mortality_model.exe", package = "MortalMove"))
# run fewer iterations due to time constraint
fit_spat <- mod$sample(
data = sim,
chains = 2,
parallel_chains = 2,
iter_warmup = 50,
iter_sampling = 150,
seed = 123
)
diagnostic_fit(fit_spat, plot = TRUE, observed_y = rowSums(sim$time_step, na.rm = TRUE), delta = sim$delta)
#> # A tibble: 5 × 5
#> variable mean sd q5 q95
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 beta[1] 0.139 0.0564 0.0536 0.243
#> 2 beta[2] 0.0363 0.424 -0.608 0.728
#> 3 alpha[1] -0.0185 0.00333 -0.0233 -0.0126
#> 4 alpha[2] 0.214 0.118 0.0100 0.403
#> 5 llambda -0.881 1.22 -3.18 0.828
#>
#> --- Rhat and ESS Summary ---
#> # A tibble: 5 × 4
#> variable rhat ess_bulk ess_tail
#> <chr> <dbl> <dbl> <dbl>
#> 1 beta[1] 1.02 389. 138.
#> 2 beta[2] 1.12 14.5 29.0
#> 3 alpha[1] 1.25 6.49 19.5
#> 4 alpha[2] 1.04 30.7 80.4
#> 5 llambda 1.75 3.33 15.0
#>
#> Computed from 300 by 100 log-likelihood matrix.
#>
#> Estimate SE
#> elpd_loo -137.5 25.4
#> p_loo 8.5 1.8
#> looic 274.9 50.9
#> ------
#> MCSE of elpd_loo is 0.2.
#> MCSE and ESS estimates assume independent draws (r_eff=1).
#>
#> All Pareto k estimates are good (k < 0.6).
#> See help('pareto-k-diagnostic') for details.
#>
#> Posterior predictive p-value summary:
#> variable p_value
#> 1 Failed 0.3000
#> 2 Censored 0.7625
bayesplot::mcmc_trace(fit_spat$draws(c("llambda","alpha","beta")))
