CNVRG (version 0.1)

varHMC: @encoding UTF-8 Perform Hamiltonian Monte Carlo sampling

Description

This function uses a compiled Dirichlet-multinomial model and performs Hamiltonian Monte Carlo sampling of posteriors using 'Stan'. After sampling it is important to check convergence. Use the summary function and shinystan to do this. If you use this function then credit 'Stan' and 'RStan' along with this package. Warning: data must be input in the correct organized format or this function will not provide accurate results. See vignette if you are unsure how to organize data. Warning: depending upon size of data to be analyzed this function can take a very long time to run.

Usage

varHMC(
  countData,
  starts,
  ends,
  algorithm = "NUTS",
  chains = 2,
  burn = 500,
  samples = 1000,
  thinning_rate = 2,
  cores = 1,
  params_to_save = c("pi", "p")
)

Arguments

countData

A matrix or data frame of counts, must only include numeric data. Data should be arranged so that the first n rows correspond to one treatment group and the next n rows correspond with the next treatment group, and so on. The row indices for the first and last sample in these groups are fed into this function via 'starts' and 'ends'.

starts

A vector defining the indices that correspond to the first sample in each treatment group.

ends

A vector defining the indices that correspond to the last sample in each treatment group.

algorithm

The algorithm to use when sampling. Either 'NUTS' or 'HMC' or 'Fixed_param'. If unsure, then pick NUTS. This is "No U-turn sampling". The abbreviation is from 'Stan'.

chains

The number of chains to run.

burn

The warm-up or 'burn-in' time.

samples

How many samples from the posterior to save.

thinning_rate

Thinning rate to use during sampling.

cores

The number of cores to use.

params_to_save

The parameters from which to save samples. Can be 'p', 'pi', 'theta'.

Value

A fitted 'Stan' object that includes the samples from the parameters designated.

Examples

# NOT RUN {
com_demo <-matrix(0, nrow = 10, ncol = 10)
com_demo[1:5,] <- c(rep(3,5), rep(7,5)) #Alternates 3 and 7
com_demo[6:10,] <- c(rep(7,5), rep(3,5)) #Reverses alternation
names(com_demo) <- rep("name", 10)
#These are toy data, many more samples, multiple chains, and a longer burn
#are likely advisable for real data.
fitstan_HMC <- varHMC(com_demo,starts = c(1,6),
ends=c(5,10),
chains = 1,
burn = 100,
samples = 150,
thinning_rate = 2)
# }