Pathfinder to Rstan

Interfaces
5

Thnaks @andrjohns , I was of the impression that if you delete the draw files, then the fit$draws() command fails to produce them.

If this is not the case, it would simplify things a lot, for the use of cmdstanR in R packages.

Interestingly on some models of mine, vb get’s closer to HMC than Pathfinder. Below I am plotting the point estimates. (tagging @avehtari in case you find this of interest)

image
image1356×1004 44.8 KB

I am using sccomp to call all three fitting methods. Below the code

devtools::install_github(stemangiola/sccomp@cmdstanr)
library(sccomp)
path  =
  seurat_obj |>
  sccomp_estimate(
  formula_composition = ~ continuous_covariate * type ,
  formula_variability = ~ 1,
  sample, cell_group,
  cores = 1,
  mcmc_seed = 42,
  max_sampling_iterations = 1000, inference_method = "pathfinder"
  )
var  =
  seurat_obj |>
  sccomp_estimate(
  formula_composition = ~ continuous_covariate * type ,
  formula_variability = ~ 1,
  sample, cell_group,
  cores = 1,
  mcmc_seed = 42,
  max_sampling_iterations = 1000, inference_method = "variational"
  )
hmc  =
  seurat_obj |>
  sccomp_estimate(
  formula_composition = ~ continuous_covariate * type ,
  formula_variability = ~ 1,
  sample, cell_group,
  cores = 1,
  mcmc_seed = 42,
  max_sampling_iterations = 1000, inference_method = "hmc"
  )


tibble(
hmc = hmc |> pull(c_effect),
vb = var |> pull(c_effect),
path = path |> pull(c_effect)
) |> GGally::ggpairs()

Here is how I am calling those fitting methods

mod$sample(
        data = data_for_model ,
        chains = chains,
        parallel_chains = chains,
        threads_per_chain = 1,
        iter_warmup = warmup_samples,
        iter_sampling = as.integer(output_samples /chains),
        save_warmup = FALSE,
        init = init,
        output_dir = output_directory
      )

model$variational(
          data = data,
          output_samples = output_samples,
          iter = iter,
          tol_rel_obj = 0.01,
          output_dir = output_dir,
          init = init,
          ...
        )

 model$pathfinder(
          data = data,
          tol_rel_obj = tol_rel_obj,
          output_dir = output_dir,
          init = init,
          num_paths=10,
          ...
        )