Modelling the reduction of faecal egg count data

Description

Models the reduction in faecal egg counts data with a (un)paired (zero-inflated) Poisson-gamma model formulation using the Stan modelling language. It is computationally several-fold faster compare to conventional MCMC techniques. For the installation instruction of Stan, please read: Stan Installation.

Usage

1
2
3
4
5
fecr_stan(preFEC, postFEC,rawCounts = FALSE, preCF = 50,
  postCF = preCF, paired = TRUE, zeroInflation=TRUE, 
  muPrior, kappaPrior, deltaPrior, phiPrior, 
  nsamples = 12000, nburnin = 2000, thinning = 1, nchain = 1, 
  ncore = 1, adaptdelta = 0.9, verbose = FALSE)

Arguments

preFEC

vector of pre-treatment faecal egg counts

postFEC

vector of post-treatment faecal egg counts

rawCounts

logical. If true, preFEC and postFEC correspond to raw counts (as counted on the McMaster slide). Otherwise they correspond to calculated epgs (raw counts times correction factor). Defaults to FALSE.

preCF

correction factor(s) before treatment

postCF

correction factor(s) after treatment

paired

logical. If true, uses the model for the paired design. Otherwise uses the model for the unpaired design

zeroInflation

logical. If true, uses the model with zero-inflation. Otherwise uses the model without zero-inflation

muPrior

a list with hyper-prior information for the baseline mean parameter μ. The default prior is list(priorDist = "gamma",hyperpars=c(1,0.001)), i.e. a gamma distribution with shape 1 and rate 0.001, its 90% probability mass lies between 51 and 2996

kappaPrior

a list with hyper-prior information for the dispersion parameter κ. The default prior is list(priorDist = "gamma",hyperpars=c(1,0.7)), i.e. a gamma distribution with shape 1 and rate 0.7, its 90% probability mass lies between 0.1 and 4.3 with a median of 1

deltaPrior

a list with hyper-prior information for the reduction δ. The default prior is list(priorDist = "beta",hyperpars=c(1,1))

phiPrior

a list with hyper-prior information for the zero-inflation parameter φ,The default prior is list(priorDist = "beta",hyperpars=c(1,1))

nsamples

a positive integer specifying how many iterations for each chain (including burn-in samples)

nburnin

number of burn-in samples

thinning

thinning parameter, a positive integer specifying the period for saving samples

nchain

a positive integer specifying the number of chains

ncore

number of cores to use when executing the chains in parallel

adaptdelta

the target acceptance rate, a value between 0 and 1

verbose

logical. If true, prints progress and debugging information

Details

The first time each model with non-default priors is applied, it can take up to 20 seconds for stan to compile the model. Currently the function only support prior distributions with two parameters. For a complete list of supported priors and their parameterization, please consult the list of distributions in Stan.

Sometimes the function outputs informational message from Stan regarding the Metropolis proposal rejections, this is due to the sampler hit the boundary of the parameter space. For some variables, the boundary point is not supported in the distribution. This is not a concern if there are only a few such warnings.

Value

An object of S4 class stanfit representing the fitted results. For more information, please see the stanfit-class in rstan reference manual.

Prints out the posterior summary of fecr as the reduction, meanEPG.untreated as the mean faecal egg counts before treatment, and meanEPG.treated as the mean faecal egg counts after treatment. The posterior summary contains the mean, standard deviation (sd), 2.5%, 25%, 50%, 75% and 97.5% percentiles, the 95% highest posterior density interval (HPDLow95 and HPDHigh95) and the posterior mode. NOTE: we recommend to use the 95% HPD interval and the mode for further statistical analysis.

Author(s)

Craig Wang craig.wang@uzh.ch

See Also

simData2s for simulating faecal egg counts data with two samples

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
## Not run: 
## load the sample data as a vector
data(epgs)

## apply zero-infation model to the data vector
model<-fecr_stan(epgs[,1],epgs[,2],rawCounts=FALSE,preCF=10,paired=TRUE,zeroInflation=TRUE)
samples<-stan2mcmc(model$stan.samples)

## a demonstration
demo("fecm_stan", package = "eggCounts") 

## End(Not run)