fecr_stan.Rd
Models the reduction in faecal egg counts with Bayesian hierarchical models. See Details for a list of model choices.
fecr_stan(preFEC, postFEC, rawCounts = FALSE, preCF = 50, postCF = preCF, paired = TRUE, indEfficacy = TRUE, zeroInflation = FALSE, muPrior, kappaPrior, deltaPrior, phiPrior, deltakappaPrior, nsamples = 2000, nburnin = 1000, thinning = 1, nchain = 2, ncore = 1, adaptDelta = 0.95, saveAll = FALSE, verbose = FALSE)
preFEC | numeric vector. Pre-treatment faecal egg counts. |
---|---|
postFEC | numeric vector. Post-treatment faecal egg counts. |
rawCounts | logical. If TRUE, |
preCF | a positive integer or a vector of positive integers. Pre-treatment correction factor(s). |
postCF | a positive integer or a vector of positive integers. Post-treatment correction factor(s). |
paired | logical. If TRUE, uses the model for the paired design. Otherwise uses the model for the unpaired design |
indEfficacy | logical. If TRUE, uses the paired model allowing for individual efficacy. Only use in combination with |
zeroInflation | logical. If TRUE, uses the model with zero-inflation. Otherwise uses the model without zero-inflation. |
muPrior | a named list. Prior for the group mean epg parameter \(\mu\). The default prior is |
kappaPrior | a named list. Prior for the group dispersion parameter \(\kappa\). The default prior is |
deltaPrior | a named list. Prior for the reduction parameter \(\delta\). The default prior is |
phiPrior | a named list. Prior for the zero-inflation parameter \(\phi\). The default prior is |
deltakappaPrior | a named list. Prior information for the shape parameter of reduction \(\delta_\kappa\). The default prior is |
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. |
saveAll | logical. If TRUE, posterior samples for all parameters are saved in the |
verbose | logical. If TRUE, prints progress and debugging information. |
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.
NOTE: Based on our simulation studies, we recommend to use (2.5%, 97.5%) as the 95% credible interval and the median as summary statistics of reduction for the individual efficacy model. For all other models, we recommend to use the 95% HPD interval and the mode.
The returned value is a list that consists of:
an object of S4 class stanfit
representing the fitted results
a data.frame that is the same as the printed posterior summary
unpaired without zero-inflation: set paired = FALSE
, indEfficacy = FALSE
, zeroInflation = FALSE
unpaired with zero-inflation: set paired = FALSE
, indEfficacy = FALSE
, zeroInflation = TRUE
paired without zero-inflation: set paired = TRUE
, indEfficacy = FALSE
, zeroInflation = FALSE
paired with zero-inflation: set paired = TRUE
, indEfficacy = FALSE
, zeroInflation = TRUE
paired with individual efficacy: set paired = TRUE
, indEfficacy = TRUE
, zeroInflation = FALSE
Consider using non-default prior for \(\delta\) when,
there is on average an increase in egg counts after treatment
there are divergent-sample warnings
there are non-convergence warnings
Two examples of useful non-default priors include:
list(priorDist = "normal", hyperpars = c(1, 5))
for stablizing the reduction parameter without being informative.
list(priorDist = "beta", hyperpars = c(0, 5))
for allowing up to 4-fold increase of egg count after treatment.
The first time each model with non-default priors is applied, it can take up to 20 seconds 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.
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.
Individual efficacy models: Craig Wang, Paul R. Torgerson, Ray M. Kaplan, Melissa M. George, Reinhard Furrer. (2018) Modelling anthelmintic resistance by extending eggCounts package to allow individual efficacy, International Journal for Parasitology: Drugs and Drug Resistance, Volume 8, Pages 386-393. https://doi.org/10.1016/j.ijpddr.2018.07.003
Zero-inflation models: Craig Wang, Paul R. Torgerson, Johan Hoglund, Reinhard Furrer. (2017) Zero-inflated hierarchical models for faecal egg counts to assess anthelmintic efficacy, Veterinary Parasitology, Volume 235, Pages 20-28. http://dx.doi.org/10.1016/j.vetpar.2016.12.007
Other models: Paul R. Torgerson, Michaela Paul, Reinhard Furrer. (2014) Evaluating faecal egg count reduction using a specifically designed package 'eggCounts' in R and a user friendly web interface, International Journal for Parasitology, Volume 44, Pages 299-303. http://dx.doi.org/10.1016/j.ijpara.2014.01.005
simData2s
for simulating faecal egg counts data with two samples
## load sample data data(epgs) ## apply individual efficacy model to the data vectors model <- fecr_stan(epgs$before, epgs$after, rawCounts = FALSE, preCF = 50, paired = TRUE, indEfficacy = TRUE)#> #> SAMPLING FOR MODEL 'indefficacy' NOW (CHAIN 1). #> Chain 1: Gradient evaluation took 3e-05 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds. #> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup) #> Chain 1: Iteration: 500 / 2000 [ 25%] (Warmup) #> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup) #> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling) #> Chain 1: Iteration: 1500 / 2000 [ 75%] (Sampling) #> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling) #> Chain 1: Elapsed Time: 1.04448 seconds (Warm-up) #> Chain 1: 0.972609 seconds (Sampling) #> Chain 1: 2.01709 seconds (Total) #> #> SAMPLING FOR MODEL 'indefficacy' NOW (CHAIN 2). #> Chain 2: Gradient evaluation took 5.8e-05 seconds #> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds. #> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup) #> Chain 2: Iteration: 500 / 2000 [ 25%] (Warmup) #> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup) #> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling) #> Chain 2: Iteration: 1500 / 2000 [ 75%] (Sampling) #> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling) #> Chain 2: Elapsed Time: 0.926934 seconds (Warm-up) #> Chain 2: 0.70542 seconds (Sampling) #> Chain 2: 1.63235 seconds (Total) #> Model: Bayesian model without zero-inflation for paired design allowing individual efficacy #> Number of Samples: 2000 #> Warm-up samples: 1000 #> Thinning: 1 #> Number of Chains 2 #> mean sd 2.5% 50% 97.5% HPDLow95 #> FECR 0.8802 0.0893 0.6609 0.8960 0.9988 0.7063 #> meanEPG.untreated 1415.8815 701.7183 511.4374 1271.0191 3161.4554 391.9358 #> meanEPG.treated 172.3279 169.2501 1.5226 126.2742 658.0284 0.0000 #> mode HPDHigh95 #> FECR 0.9076 1.0000 #> meanEPG.untreated 986.5780 2843.4982 #> meanEPG.treated 88.2204 511.9502 #> #> NOTE: there is no evidence of non-convergence since all parameters have potential scale reduction factors (Brooks and Gelman, 1998) less than 1.1.## convert to MCMC object and inspect the summary samples <- stan2mcmc(model$stan.samples) summary(samples)#> #> Iterations = 1:2000 #> Thinning interval = 1 #> Number of chains = 1 #> Sample size per chain = 2000 #> #> 1. Empirical mean and standard deviation for each variable, #> plus standard error of the mean: #> #> Mean SD Naive SE Time-series SE #> FECR 0.8802 0.08929 0.001996 0.001996 #> meanEPG.untreated 1415.8815 701.71834 15.690899 15.690899 #> meanEPG.treated 172.3279 169.25012 3.784548 3.950002 #> kappa 0.2374 0.08820 0.001972 0.001887 #> delta_mu 0.4208 0.25356 0.005670 0.005992 #> delta_shape 0.5455 0.49689 0.011111 0.011111 #> #> 2. Quantiles for each variable: #> #> 2.5% 25% 50% 75% 97.5% #> FECR 0.66087 0.8322 0.8960 0.9474 0.9988 #> meanEPG.untreated 511.43741 922.0963 1271.0191 1732.2854 3161.4554 #> meanEPG.treated 1.52260 59.0009 126.2742 231.7798 658.0284 #> kappa 0.10280 0.1725 0.2254 0.2880 0.4361 #> delta_mu 0.08409 0.2021 0.3633 0.6072 0.9478 #> delta_shape 0.08564 0.2306 0.3841 0.6712 2.0404 #>