Title: | Hierarchical Modelling of Faecal Egg Counts |
---|---|
Description: | An implementation of Bayesian hierarchical models for faecal egg count data to assess anthelmintic efficacy. Bayesian inference is done via MCMC sampling using 'Stan' <https://mc-stan.org/>. |
Authors: | Craig Wang [aut, cre] , Michaela Paul [aut], Tea Isler [ctb], Reinhard Furrer [ctb] , Trustees of Columbia University [cph] (src/init.cpp, tools/make_cc.R, R/stanmodels.R) |
Maintainer: | Craig Wang <[email protected]> |
License: | GPL (>= 3) |
Version: | 2.4 |
Built: | 2024-10-28 07:12:10 UTC |
Source: | CRAN |
This package implements Bayesian hierarchical models for the analysis of faecal egg count data. Bayesian inference is done via efficient MCMC sampling using Stan. Additional (experimental) models are available externally for handling FECs with potential outliers or bi-modality. The models are in eggCountsExtra package hosted on Github.
Package: | eggCounts |
Type: | Package |
Version: | 2.4 |
Date: | 2023-10-14 |
License: | GPL (>= 3) |
LazyLoad: | yes |
Stan is a probabilistic programming language for specifying Bayesian hierarchical models. It is computationally faster compared to conventional MCMC techniques. For the installation instruction and other information about Stan, please read here.
Craig Wang [email protected]
Michaela Paul
## Not run: ## Citations citation('eggCounts') ## History of changes file.show(system.file("NEWS", package = "eggCounts")) ## Demonstration demo("fecm_stan", package = "eggCounts") ## Install eggCountsExtra devtools::install_github("CraigWangUZH/eggCountsExtra") ## End(Not run)
## Not run: ## Citations citation('eggCounts') ## History of changes file.show(system.file("NEWS", package = "eggCounts")) ## Demonstration demo("fecm_stan", package = "eggCounts") ## Install eggCountsExtra devtools::install_github("CraigWangUZH/eggCountsExtra") ## End(Not run)
This is an example dataset containing 14 eggs per gram (epg) values in sheep before and after anthelmintic treatment of benzimidazole. The correction factor of the diagnostic technique was 50.
data(epgs)
data(epgs)
A data.frame containing 14 observations.
Craig Wang, Paul R. Torgerson, Johan Hoglund, Reinhard Furrer, Zero-inflated hierarchical models for faecal egg counts to assess anthelmintic efficacy, Veterinary Parasitology, Volume 235, 2017, Pages 20-28.
Models the mean of faecal egg counts with Bayesian hierarchical models. See Details for a list of model choices.
fec_stan(fec, rawCounts = FALSE, CF = 50, zeroInflation = TRUE, muPrior, kappaPrior, phiPrior, nsamples = 2000, nburnin = 1000, thinning = 1, nchain = 2, ncore = 1, adaptDelta = 0.95, saveAll = FALSE, verbose = FALSE)
fec_stan(fec, rawCounts = FALSE, CF = 50, zeroInflation = TRUE, muPrior, kappaPrior, phiPrior, nsamples = 2000, nburnin = 1000, thinning = 1, nchain = 2, ncore = 1, adaptDelta = 0.95, saveAll = FALSE, verbose = FALSE)
fec |
numeric vector. Faecal egg counts. |
rawCounts |
logical. If TRUE, |
CF |
a positive integer or a vector of positive integers. Correction factor(s). |
zeroInflation |
logical. If true, uses the model with zero-inflation. Otherwise uses the model without zero-inflation |
muPrior |
named list. Prior for the group mean epg parameter |
kappaPrior |
named list. Prior for the group dispersion parameter |
phiPrior |
named list. Prior for the zero-inflation parameter |
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. |
without zero-inflation: set zeroInflation = FALSE
with zero-inflation: set zeroInflation = TRUE
Note that this function only models the mean of egg counts, see fecr_stan()
for modelling the reduction.
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 User Guide.
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.
Prints out summary of meanEPG
as the posterior mean 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: we recommend to use the 95% HPD interval and the mode for further statistical analysis.
The returned value is a list that consists of:
stan.samples |
an object of S4 class |
posterior.summary |
a data.frame that is the same as the printed posterior summary |
Craig Wang
simData1s
for simulating faecal egg count data with one sample
## load the sample data data(epgs) ## apply zero-infation model model <- fec_stan(epgs$before, rawCounts = FALSE, CF = 50)
## load the sample data data(epgs) ## apply zero-infation model model <- fec_stan(epgs$before, rawCounts = FALSE, CF = 50)
Computes probability of the reduction parameter's marginal posterior density relative to a threshold.
fecr_probs(stanFit, threshold = 0.95, lessthan = TRUE, plot = TRUE, xlab, ylab, main, verbose = TRUE, ...)
fecr_probs(stanFit, threshold = 0.95, lessthan = TRUE, plot = TRUE, xlab, ylab, main, verbose = TRUE, ...)
stanFit |
|
threshold |
numeric. The default threshold is 0.95 (95%). |
lessthan |
logical. If TRUE, the probability less than the threshold is computed. Otherwise greater or equal to the threshold is computed. Default is TRUE. |
plot |
logical. If TRUE, the posterior density of the reduction is plotted with region less than the threshold shaded. |
xlab , ylab , main
|
strings. Arguments for plotting. Only used if |
verbose |
logical. If TRUE, a statement with computed probability is printed. |
... |
additional plotting arguments |
Returns a numeric value indicating the probability in percentage.
Craig Wang
## load sample data data(epgs) ## apply individual efficacy model to the data vectors model <- fecr_stan(epgs$before, epgs$after, rawCounts = FALSE, preCF = 10, paired = TRUE, indEfficacy = TRUE) fecr_probs(model$stan.samples)
## load sample data data(epgs) ## apply individual efficacy model to the data vectors model <- fecr_stan(epgs$before, epgs$after, rawCounts = FALSE, preCF = 10, paired = TRUE, indEfficacy = TRUE) fecr_probs(model$stan.samples)
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)
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 |
kappaPrior |
a named list. Prior for the group dispersion parameter |
deltaPrior |
a named list. Prior for the reduction parameter |
phiPrior |
a named list. Prior for the zero-inflation parameter |
deltakappaPrior |
a named list. Prior information for the shape parameter of reduction |
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. |
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 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 User Guide.
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.
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:
stan.samples |
an object of S4 class |
posterior.summary |
a data.frame that is the same as the printed posterior summary |
Craig Wang
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. <doi: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. <doi: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. <doi: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) ## convert to MCMC object and inspect the summary samples <- stan2mcmc(model$stan.samples) summary(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) ## convert to MCMC object and inspect the summary samples <- stan2mcmc(model$stan.samples) summary(samples)
Models the reduction in faecal egg counts with custom model formulation using Stan modelling language (for advanced users).
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)
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)
preFEC |
numeric vector. Pre-treatment faecal egg counts. Not required if |
postFEC |
numeric vector. Post-treatment faecal egg counts. Not required if |
rawCounts |
logical. If TRUE, |
preCF |
a positive integer or a vector of positive integers. Pre-treatment correction factor(s). Not required if |
postCF |
a positive integer or a vector of positive integers. Post-treatment correction factor(s). Not required if |
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 |
stan model code. Not required when |
modelFile |
stan model file with file extension ‘*.stan’. Not required when |
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 |
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. |
If modelName
is one of c("Po", "UPo", "ZIPo", "ZIUPo"), then outlier-adjusted models are used.
In paired models, outliers are those counts with postFEC > preFEC. Outlier weights are assigned as the inverse of postFEC/preFEC,
In unpaired models, outliers are those counts with postFEC greater than the 95th percentile of a Poisson distribution, where the Poisson mean is computed based on the mean of postFEC excluding postFEC > Q3 + 1.5*IQR. Q3 is the 75th percentile and IQR is the interquartile range. The lowest outlier weight is assigned as 0.01, and other outliers assigned proportionally.
In both cases, non-outliers are assigned with outlier weight = 1.
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.
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 |
stan.samples |
an object of S4 class |
posterior.summary |
a data.frame that is the same as the printed posterior summary. Not available for custom models. |
Craig Wang
## 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)
## 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)
Models the reduction in faecal egg counts with a simple Bayesian model formulation. The model is for paired design only, and it assumes Poisson distribution for the observed egg counts.
fecr_stanSimple(preFEC, postFEC, rawCounts = FALSE, preCF = 50, postCF = preCF, muPrior, deltaPrior, nsamples = 2000, nburnin = 1000, thinning = 1, nchain = 2, ncore = 1, adaptDelta = 0.95, saveAll = FALSE, verbose = FALSE)
fecr_stanSimple(preFEC, postFEC, rawCounts = FALSE, preCF = 50, postCF = preCF, muPrior, deltaPrior, 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 |
positive integer or vector of positive integers. Pre-treatment correction factor(s). |
postCF |
positive integer or vector of positive integers. Post-treatment correction factor(s). |
muPrior |
named list. Prior for the group mean epg parameter |
deltaPrior |
named list. Prior for the reduction parameter |
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. |
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 User Guide.
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.
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: we recommend to use the 95% HPD interval and the mode for further statistical analysis.
The returned value is a list that consists of:
stan.samples |
an object of S4 class |
posterior.summary |
A data.frame that is the same as the printed posterior summary |
Tea Isler
Craig Wang
simData2s
for simulating faecal egg counts data with two samples
## load sample data data(epgs) ## apply paired model with individual efficacy model <- fecr_stanSimple(epgs$before, epgs$after, rawCounts = FALSE, preCF = 10) samples <- stan2mcmc(model$stan.samples)
## load sample data data(epgs) ## apply paired model with individual efficacy model <- fecr_stanSimple(epgs$before, epgs$after, rawCounts = FALSE, preCF = 10) samples <- stan2mcmc(model$stan.samples)
Computes the standard Faecal Egg Count Reduction Test together with approximate confidence interval according to the WAAVP guidelines (Coles et al., 1992, 2006). The function also returns bootstrap confidence intervals.
fecrtCI(epg1, epg2, paired = FALSE, alpha = 0.05, R = 1999)
fecrtCI(epg1, epg2, paired = FALSE, alpha = 0.05, R = 1999)
epg1 |
numeric vector. Faecal egg counts in untreated animals. |
epg2 |
numeric vector. Faecal egg counts in treated animals. |
paired |
logical. If TRUE, indicates samples are paired. Otherwise samples are unpaired. |
alpha |
numeric. Confidence level of the interval. |
R |
numeric. Number of bootstrap replicates. |
A list with
estimate |
estimated percentage reduction in mean epg |
bootCI |
bootstrap confidence interval |
approxCI |
approximate confidence interval |
Michaela Paul
Coles GC, Bauer C, Borgsteede FHM, Geerts S, Klei TR, Taylor MA and Waller, PJ (1992). World Association for the Advancement of Veterinary Parasitology (WAAVP) methods for the detection of anthelmintic resistance in nematodes of veterinary importance, Veterinary Parasitology, 44:35-44.
Coles GC, Jackson F, Pomroy WE, Prichard RK, von Samson-Himmelstjerna G, Silvestre A, Taylor MA and Vercruysse J (2006). The detection of anthelmintic resistance in nematodes of veterinary importance, Veterinary Parasitology, 136:167-185.
data(epgs) fecrtCI(epgs$before, epgs$after, paired = TRUE)
data(epgs) fecrtCI(epgs$before, epgs$after, paired = TRUE)
Compute the shape parameters from a Beta distribution for based on some prior belief.
getPrior_delta(lower, upper, p = 0.7, mode, conc, plot = TRUE)
getPrior_delta(lower, upper, p = 0.7, mode, conc, plot = TRUE)
lower , upper , p
|
numeric. Prior belief about the reduction. There is |
mode , conc
|
numeric. Prior belief about the reduction. The mode and concentration parameters of a beta distribution. Higher concentration indicates smaller variance. Not used if |
plot |
logical. If TRUE, the prior distribution is plotted after parameters are found. |
The multiroot
function from rootSolve package is used to compute the parameters.
Returns Beta prior parameters for and the printed argument to use in a
fecr_stan()
or a fec_stan()
function call.
Tea Isler
Craig Wang
# there is 80% probability that the reduction is between 60% and 90% getPrior_delta(lower = 0.6, upper = 0.9, p = 0.8)
# there is 80% probability that the reduction is between 60% and 90% getPrior_delta(lower = 0.6, upper = 0.9, p = 0.8)
Compute the shape and rate parameters from a Gamma distribution for based on some prior belief about its cumulative distribution function.
getPrior_mu(x, px, y, py, s1 = 1, s2 = 0.001, plot = TRUE)
getPrior_mu(x, px, y, py, s1 = 1, s2 = 0.001, plot = TRUE)
x , px , y , py
|
numeric. Threshold of some prior belief about true epg. There is |
s1 , s2
|
numeric. Starting values. |
plot |
logical. If TRUE, the prior distribution is plotted after parameters are found. |
multiroot
function from rootSolve package is used to compute the parameters.
Returns Gamma prior parameters for and the printed argument to use in a
fecr_stan()
or a fec_stan()
function call.
Tea Isler
Craig Wang
# there is 30% probability that the mean epg is less than 200 # and 80% probability that the mean epg is less than 500 getPrior_mu(x = 200, px = 0.3, y = 500, py = 0.8)
# there is 30% probability that the mean epg is less than 200 # and 80% probability that the mean epg is less than 500 getPrior_mu(x = 200, px = 0.3, y = 500, py = 0.8)
Plot egg count data to reflect changes between before and after treatment.
plotCounts(data, paired = TRUE, points = TRUE, points.method = "jitter", xlabel = "", ylabel = "Faecal egg counts [epg]", ...)
plotCounts(data, paired = TRUE, points = TRUE, points.method = "jitter", xlabel = "", ylabel = "Faecal egg counts [epg]", ...)
data |
a data.frame with two columns, the first column is before treatment counts, the second column is after treatment counts. |
paired |
logical. If TRUE, uses the plot for the paired design. Otherwise uses the plot for the unpaired design. |
points |
logical. If TRUE, add individual points for unpaired plot. Not used if |
points.method |
string. It is used to separate coincident points if |
xlabel |
string. Label of x-axis. |
ylabel |
String. Label of y-axis. |
... |
Additional arguments for function |
For paired data, a xyplot is used. For unpiared data, a grouped boxplot is used.
A plot is drawn.
Craig Wang
data(epgs) plotCounts(epgs[,c("before","after")], paired = TRUE)
data(epgs) plotCounts(epgs[,c("before","after")], paired = TRUE)
Simulates (zero-inflated) egg count data
simData1s(n = 10, mean = 500, kappa = 0.5, phi = 1, f = 50, rounding = TRUE, seed = NULL)
simData1s(n = 10, mean = 500, kappa = 0.5, phi = 1, f = 50, rounding = TRUE, seed = NULL)
n |
positive integer. Sample size. |
mean |
numeric. True number of eggs per gram (epg). |
kappa |
numeric. Overdispersion parameter, |
phi |
numeric. Prevalence, i.e. proportion of infected animals, between 0 and 1. |
f |
positive integer. Correction factor of the egg counting technique, either an integer or a vector of integers with length |
rounding |
logical. If TRUE, the Poisson mean for the raw counts is rounded. The rounding applies since the mean epg is frequently reported as an integer value. For more information, see Details. |
seed |
integer. Random seed. |
In the simulation of raw (master
) counts, it follows a Poisson distribution with some mean. The mean is frequently rounded down if it has a very low value and rounding = TRUE
, hence there expects to be a bias overall when < 150. Set
rounding = FALSE
for not to have any bias in the simulated counts.
A data.frame with three columns, namely the observed epg (obs
), actual number of eggs counted (master
) and
true epg in the sample (true
).
Craig Wang
Michaela Paul
fec_stan
for analyzing faecal egg count data with one sample
fec <- simData1s(n = 10, mean = 500, kappa = 0.5, phi = 0.7)
fec <- simData1s(n = 10, mean = 500, kappa = 0.5, phi = 0.7)
Generates two samples of (zero-inflated) egg count data
simData2s(n = 10, preMean = 500, delta = 0.1, kappa = 0.5, deltaShape = NULL, phiPre = 1, phiPost = phiPre, f = 50, paired = TRUE, rounding = TRUE, seed = NULL)
simData2s(n = 10, preMean = 500, delta = 0.1, kappa = 0.5, deltaShape = NULL, phiPre = 1, phiPost = phiPre, f = 50, paired = TRUE, rounding = TRUE, seed = NULL)
n |
positive integer. Sample size. |
preMean |
numeric. True pre-treatment epg. |
delta |
numeric. Proportion of epg left after treatment, between 0 and 1. 1 - |
kappa |
numeric. Overdispersion parameter, |
deltaShape |
numeric. Shape parameter for the distribution of reductions. If NULL, the same reduction is applied to the latent true epg of each animal. |
phiPre |
numeric. Pre-treatment prevalence (i.e. proportion of infected animals), between 0 and 1. |
phiPost |
numeric. Post-treatment prevalence, between 0 and 1. |
f |
integer or vector of integers. Correction factor of the egg counting technique |
paired |
logical. If TRUE, paired samples are simulated. Otherwise unpaired samples are simulated. |
rounding |
logical. If TRUE, the Poisson mean for the raw counts is rounded. The rounding applies since the mean epg is frequently reported as an integer value. For more information, see Details. |
seed |
an integer that will be used in a call to set.seed before simulation. If NULL, a random seed is allocated. |
In the simulation of raw (master
) counts, it follows a Poisson distribution with some mean. The mean is frequently rounded down if it has a very low value and rounding = TRUE
, there expects to be some bias in the mean reduction when < 150 and
< 0.1. Set
rounding = FALSE
for not to have any bias.
A data.frame with six columns, namely the observed epg (obs
),
actual number of eggs counted (master
) and true epg in the sample (true
) for both pre- and post- treatment.
Craig Wang
Michaela Paul
fecr_stan
for analyzing faecal egg count data with two samples
fec <- simData2s(n = 10, preMean = 500, delta = 0.1, kappa = 0.5) ## show the bias when the true reduction should be 95% fec <- simData2s(n = 1e5, preMean = 150, delta = 0.05, kappa = 0.5, seed = 1) 1 - mean(fec$masterPost)/mean(fec$masterPre) ## without bias fec <- simData2s(n = 1e5, preMean = 150, delta = 0.05, kappa = 0.5, seed = 1, rounding = FALSE) 1 - mean(fec$masterPost)/mean(fec$masterPre)
fec <- simData2s(n = 10, preMean = 500, delta = 0.1, kappa = 0.5) ## show the bias when the true reduction should be 95% fec <- simData2s(n = 1e5, preMean = 150, delta = 0.05, kappa = 0.5, seed = 1) 1 - mean(fec$masterPost)/mean(fec$masterPre) ## without bias fec <- simData2s(n = 1e5, preMean = 150, delta = 0.05, kappa = 0.5, seed = 1, rounding = FALSE) 1 - mean(fec$masterPost)/mean(fec$masterPre)
Converts a stanfit
object into a mcmc
object for easier analysis.
stan2mcmc(stanFit)
stan2mcmc(stanFit)
stanFit |
a |
The output can be analyzed as a mcmc
object with the functions from the coda
package. NOTE: The resulting MCMC object does not contain warm-up samples and is already thinned.
A MCMC object with a list of relevant parameters.
Craig Wang
data(epgs) ## apply zero-infation model for the paired design model <- fecr_stan(epgs$before, epgs$after, rawCounts = FALSE, indEfficacy = FALSE, preCF = 10, paired = TRUE, zeroInflation = TRUE) samples <- stan2mcmc(model$stan.samples) summary(samples) plot(samples)
data(epgs) ## apply zero-infation model for the paired design model <- fecr_stan(epgs$before, epgs$after, rawCounts = FALSE, indEfficacy = FALSE, preCF = 10, paired = TRUE, zeroInflation = TRUE) samples <- stan2mcmc(model$stan.samples) summary(samples) plot(samples)