fecr_stanExtra {eggCounts}R Documentation

Model the reduction of faecal egg count using custom models

Description

Models the reduction in faecal egg counts with custom model formulation using Stan modelling language (for advanced users).

Usage

fecr_stanExtra(preFEC, postFEC, rawCounts = FALSE, preCF = 50, postCF = preCF, 
  modelName = NULL, modelCode = NULL, modelFile = NULL, modelData = NULL,
  nsamples = 2000, nburnin = 1000, thinning = 1, nchain = 2, 
  ncore = 1, adaptDelta = 0.95, verbose = FALSE)

Arguments

preFEC

numeric vector. Pre-treatment faecal egg counts. Not required if modelData is supplied.

postFEC

numeric vector. Post-treatment faecal egg counts. Not required if modelData is supplied.

rawCounts

logical. If TRUE, preFEC and postFEC correspond to raw counts (as counted on equipment). Otherwise they correspond to calculated epgs (raw counts times correction factor). Defaults to FALSE. Not required if modelCode or modelFile is supplied.

preCF

a positive integer or a vector of positive integers. Pre-treatment correction factor(s). Not required if modelCode or modelFile is supplied.

postCF

a positive integer or a vector of positive integers. Post-treatment correction factor(s). Not required if modelCode or modelFile is supplied.

modelName

string. One of four availale models ("Po", "UPo", "ZIPo", "ZIUPo") from eggCountsExtra package, which corresponds to outlier-adjusted version of paired, unpaired, paired with zero inflation and unpaired with zero inflation models. Not required if modelCode or modelFile is supplied.

modelCode

stan model code. Not required when modelName or modelFile is supplied.

modelFile

stan model file with file extension ‘*.stan’. Not required when modelName or modelCode is supplied.

modelData

stan data list. A named list or environment providing the data for the model, or a character vector for all the names of objects in the current enviroment used as data. Not required when modelName is supplied.

nsamples

a positive integer. Number of samples for each chain (including burn-in samples).

nburnin

a positive integer. Number of burn-in samples.

thinning

a positive integer. Thinning parameter, i.e. the period for saving samples.

nchain

a positive integer. Number of chains.

ncore

a positive integer. Number of cores to use when executing the chains in parallel.

adaptDelta

numeric. The target acceptance rate, a numeric value between 0 and 1.

verbose

logical. If TRUE, prints progress and debugging information.

Details

If modelName is one of c("Po", "UPo", "ZIPo", "ZIUPo"), then outlier-adjusted models are used.

The first time each model is applied, it can take up to 20 seconds for Stan to compile the model.

The default number of samples per chain is 2000, with 1000 burn-in samples. Normally this is sufficient in Stan. If the chains do not converge, one should tune the MCMC parameters until convergence is reached to ensure reliable results.

Value

Prints out the posterior summary of FECR as the reduction, meanEPG.untreated as the mean pre-treatment epg, and meanEPG.treated as the mean after-treatment epg. The posterior summary contains the mean, standard deviation (sd), 2.5%, 50% and 97.5% percentiles, the 95% highest posterior density interval (HPDLow95 and HPDHigh95) and the posterior mode.

The returned value is a list that consists of:

stan.model

an object of class stanmodel-class that was used

stan.samples

an object of S4 class stanfit representing the fitted results

posterior.summary

a data.frame that is the same as the printed posterior summary. Not available for custom models.

Author(s)

Craig Wang

Examples

## Not run: 
library(eggCountsExtra)
data(epgs) ## load sample data

## apply paired model with outliers 
model1 <- fecr_stanExtra(epgs$before, epgs$after, rawCounts=FALSE, 
         preCF=10, modelName = "Po")
samples <- stan2mcmc(model1$stan.samples)
fecr_probs(model1$stan.samples, threshold = 0.99)

## apply a simple custom model
code <- "data{
  int J; // number of animals
  int y_before[J]; // after treatment McMaster count
  int y_after[J]; // before treatment McMaster count
}
parameters{
  real<lower=0> mu;
  real<lower=0,upper=1> delta;
}
model{
  mu ~ gamma(1,0.001);
  delta ~ beta(1,1);
  y_before ~ poisson(mu);
  y_after ~ poisson(mu*delta);
}"

dat <- list(J = nrow(epgs), y_before = epgs$before,
            y_after = epgs$after)
model2 <- fecr_stanExtra(modelCode = code, modelData = dat)

## End(Not run)

[Package eggCounts version 2.3-2 Index]