Nothing
R/fec.power.R
In bayescount: Power Calculations and Bayesian Analysis of Count Distributions and FECRT Data using MCMC
Defines functions
fec.power
fec.power.limits
Documented in
fec.power
fec.power.limits
fec.power.limits <- function(meanepg=200, g.faeces=3, sensitivity=1/25, replicates=1, animals=10, coeffvarrep=0.4, coeffvarind=0.3, coeffvargroup=0.7, true.sample=FALSE, lower.limit=NA, upper.limit=NA, iterations=100000, power = 0.95, confidence = 0.99, feedback=FALSE, forcesim=FALSE){
# Power is the required power, confidence is the confidence in this power when using simulations
if(power >= 1) stop("Required power must be < 1")
if(confidence >= 1) stop("Confidence must be < 1")
conf <- (1-confidence)/2
lci <- 0+conf
uci <- 1-conf
if(replicates < 1 | animals < 1) stop("Specified values for animals and replicates must be greater than 0")
if(!is.na(lower.limit) & !is.na(upper.limit)) stop("One or both of lower.limit or upper.limit must be non-fixed (NA)")
fix.lower <- !is.na(lower.limit)
fix.upper <- !is.na(upper.limit)
fixed.lower <- lower.limit
fixed.upper <- upper.limit
target <- power
if(feedback){
if(any(.Platform$GUI == c("AQUA", "Rgui"))){
#feedback <- FALSE
warning("Printing the progress of the function using the GUI versions of R may massively increase the time taken, I suggest setting feedback=FALSE or using a command line version of R instead")
}
##### CHECK ON WINDOWS / UNIX
}
# If we have 1 animal and true.sample=T then we can forget about cv.group and use nbpower
if(animals==1 & true.sample==TRUE) coeffvargroup <- 10^-10
if(coeffvargroup > (coeffvarrep+coeffvarind)/10) approximate <- FALSE else approximate <- TRUE
if(forcesim) approximate <- FALSE
if(animals==1 & true.sample==FALSE) warning("NOTE: The power calculated is for the population from which the animal is derived, to calculate the power for the true mean of this individual use true.sample=TRUE")
lowerl <- lower.limit
upperl <- upper.limit
start <- Sys.time()
if(!approximate){
if(true.sample){
out <- .C("poweranalysissamplefixed", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.integer(iterations), as.integer(feedback), numeric(iterations), PACKAGE="bayescount")
lo <- length(out)
meancounts <- out[[lo]]
}else{
out <- .C("poweranalysispopulationfixed", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.integer(iterations), as.integer(feedback), numeric(iterations), PACKAGE="bayescount")
lo <- length(out)
meancounts <- out[[lo]]
}
mcs <- meancounts
f <- function(limit){
prob <- limit
#limits <- HPDinterval(mcmc(mcs), prob=prob)[1,]
limits <- quantile(mcs, probs=c(0+((1-prob)/2), 1-((1-prob)/2)))
nin <- sum(mcs <= limits[2] & mcs >= limits[1])
nout <- length(mcs)-nin
med <- qbeta(0.5, nin+1, nout+1)
return(med)
}
if(fix.lower){
f <- function(limit){
#limits <- HPDinterval(mcmc(mcs), prob=prob)[1,]
nin <- sum(mcs <= limit & mcs >= fixed.lower)
nout <- length(mcs)-nin
med <- qbeta(0.5, nin+1, nout+1)
return(med)
}
}
if(fix.upper){
f <- function(limit){
#limits <- HPDinterval(mcmc(mcs), prob=prob)[1,]
nin <- sum(mcs <= fixed.upper & mcs >= -limit)
nout <- length(mcs)-nin
med <- qbeta(0.5, nin+1, nout+1)
return(med)
}
}
limits <- c(0,1)
if(fix.upper) limits <- c(-Inf,0)
if(fix.lower) limits <- c(0,Inf)
bsres <- binary.search(f, target, limits)
if(bsres$status!="OK"){
if(bsres$status!="Absolute value not possible but this is closest") stop(bsres$status)
}
if(fix.upper) limits <- c(-bsres$value, fixed.upper)
if(fix.lower) limits <- c(fixed.lower, bsres$value)
if(!fix.upper & !fix.lower) limits <- quantile(mcs, prob=c(0+((1-bsres$value)/2), 1-((1-bsres$value)/2)))
nin <- sum(mcs <= limits[2] & mcs >= limits[1])
nout <- length(mcs)-nin
power <- qbeta(c(lci,0.5,uci), nin+1, nout+1)
names(limits) <- c("lower.limit", "upper.limit")
names(power) <- c(paste("lower.", confidence*100, "%.estimate", sep=""), "median.estimate", paste("upper.", confidence*100, "%.estimate", sep=""))
return(list(limits=limits, power=power))
}else{
#Embed within wrapper function (that calculates cvw from rep, sample, g.faeces etc and works out if approximation to multiple animals can be used):
if(true.sample==FALSE | animals > 1){
# If we're using this with true.sample=FALSE then cvgroup must be so small compared to cvind and cvrep that we can just pretend that we're drawing all replicates from the same animal (with a slightly larger cv and animalmean=populationmean):
coeffvarind <- sqrt(coeffvarind^2 + coeffvargroup^2 + coeffvarind^2*coeffvargroup^2)
replicates <- replicates*animals
}
coeff.var <- sqrt(coeffvarind^2 + (coeffvarrep^2)/g.faeces + (coeffvarind^2 * (coeffvarrep^2)/g.faeces))
eggs.counted <- replicates*meanepg*sensitivity
shape <- replicates / coeff.var^2
f <- function(limit){
accuracy <- limit
upper.tol <- replicates*sensitivity*(meanepg+meanepg*accuracy)
lower.tol <- replicates*sensitivity*(meanepg-meanepg*accuracy)
return(pnbinom(upper.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) - pnbinom(lower.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) + suppressWarnings(dnbinom(lower.tol, size=shape, mu=eggs.counted)))
}
if(fix.lower){
f <- function(limit){
upper.tol <- replicates*sensitivity*(limit)
lower.tol <- replicates*sensitivity*(fixed.lower)
return(pnbinom(upper.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) - pnbinom(lower.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) + suppressWarnings(dnbinom(lower.tol, size=shape, mu=eggs.counted)))
}
}
if(fix.upper){
f <- function(limit){
upper.tol <- replicates*sensitivity*(fixed.upper)
lower.tol <- replicates*sensitivity*(-limit)
return(pnbinom(upper.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) - pnbinom(lower.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) + suppressWarnings(dnbinom(lower.tol, size=shape, mu=eggs.counted)))
}
}
limits <- c(0,Inf)
if(fix.upper) limits <- c(-Inf,0)
if(fix.lower) limits <- c(0,Inf)
bsres <- binary.search(f, target, limits)
if(bsres$status!="OK"){
if(bsres$status!="Absolute value not possible but this is closest") stop(bsres$status)
}
if(fix.upper) limits <- c(-bsres$value, fixed.upper)
if(fix.lower) limits <- c(fixed.lower, bsres$value)
if(!fix.upper & !fix.lower) limits <- c((meanepg-meanepg*bsres$value), (meanepg+meanepg*bsres$value))
power <- replicate(3, bsres$objective)
names(limits) <- c("lower.limit", "upper.limit")
names(power) <- c("absolute.value", "absolute.value", "absolute.value")
return(list(limits=limits, power=power))
}
}
fec.power <- function(meanepg=200, g.faeces=3, sensitivity=1/25, replicates=1, animals=10, coeffvarrep=0.4, coeffvarind=0.3, coeffvargroup=0.7, true.sample=FALSE, accuracy=0.1, lower.limit=meanepg*(1-accuracy), upper.limit=meanepg*(1+accuracy), maxiterations=1000000, precision=2, confidence = 0.99, feedback=FALSE, forcesim=FALSE){
if(confidence >= 1) stop("Confidence must be < 1")
conf <- (1-confidence)/2
lci <- 0+conf
uci <- 1-conf
if(feedback){
if(any(.Platform$GUI == c("AQUA", "Rgui"))){
#feedback <- FALSE
warning("Printing the progress of the function using the GUI versions of R may massively increase the time taken, I suggest setting feedback=FALSE or using a command line version of R instead")
}
##### CHECK ON WINDOWS / UNIX
}
if(replicates < 1 | animals < 1) stop("Specified values for animals and replicates must be greater than 0")
# If we have 1 animal and true.sample=T then we can forget about cv.group and use nbpower
if(animals==1 & true.sample==TRUE) coeffvargroup <- 10^-10
if(coeffvargroup > (coeffvarrep+coeffvarind)/10) approximate <- FALSE else approximate <- TRUE
if(forcesim) approximate <- FALSE
if(animals==1 & true.sample==FALSE) warning("NOTE: The power calculated is for the population from which the animal is derived, to calculate the power for the true mean of this individual use true.sample=TRUE")
lowerl <- lower.limit
upperl <- upper.limit
start <- Sys.time()
if(approximate){
if(true.sample==FALSE){
# If we're using this with true.sample=FALSE then cvgroup must be so small compared to cvind and cvrep that we can just pretend that we're drawing all replicates from the same animal (with a slightly larger cv and animalmean=populationmean):
coeffvarind <- sqrt(coeffvarind^2 + coeffvargroup^2 + coeffvarind^2*coeffvargroup^2)
replicates <- replicates*animals
}
# Produces the same results as poweranalysispopulation (only for animal=1 and/or cv.group<<cv.other)
coeff.var <- sqrt(coeffvarind^2 + (coeffvarrep^2)/g.faeces + (coeffvarind^2 * (coeffvarrep^2)/g.faeces))
eggs.counted <- replicates*meanepg*sensitivity
upper.tol <- replicates*sensitivity*(upperl)
lower.tol <- replicates*sensitivity*(lowerl)
shape <- replicates / coeff.var^2
nbpower <- pnbinom(upper.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) - pnbinom(lower.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) + suppressWarnings(dnbinom(lower.tol, size=shape, mu=eggs.counted))
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output = list(roundedci=replicate(3, round(nbpower, precision)), ci=replicate(3, nbpower))
}else{
if(true.sample){
if(is.na(precision)){
out <- .C("poweranalysissamplefixed", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.integer(maxiterations), as.integer(feedback), numeric(maxiterations), PACKAGE="bayescount")
lo <- length(out)
meancounts <- out[[lo]]
nin <- sum(meancounts >= lowerl & meancounts <= upperl)
nout <- length(meancounts)-nin
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output=list(roundedci=round(qbeta(c(lci, 0.5, uci), nin+1, nout+1), precision), ci=qbeta(c(lci, 0.5, uci), nin+1, nout+1), within=nin, without=nout, total=nin+nout)
}else{
out <- .C("poweranalysissample", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.numeric(lowerl), as.numeric(upperl), as.integer(maxiterations), as.integer(precision), as.numeric(lci), as.numeric(uci), as.integer(feedback), as.integer(0), as.integer(0), PACKAGE="bayescount")
lo <- length(out)
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output=list(roundedci = round(qbeta(c(lci, 0.5, uci), out[[lo-1]]+1, (out[[lo]]-out[[lo-1]])+1), precision), ci=qbeta(c(lci, 0.5, uci), out[[lo-1]]+1, (out[[lo]]-out[[lo-1]])+1), within=out[[lo-1]], without=out[[lo]]-out[[lo-1]], total=out[[lo]])
}
}else{
# coeffvarrep is the cv between subsamples taken from the same replicate and is effected by gfaeces (calculation done in C)
# coeffvarind is cv between replicates and is effected by replicates. This actually means it is EITHER the combined cv between samples and days within an animal OR just the cv between samples if taken from a single faecal pile (in which case the cv between days is actually part of the cv between animals and will consequently make things a bit wrong due to the true sample effect). Doesn't reflect subsamples (what I usually call replicates) as this is exactly the same as changing the sensitivity and g.faeces.
# coeffvargroup is cv between animals and is always calculated explicitly since there is no shortcut. If we have only 1 animal or coeffvargroup=0, then powernb works perfectly.
if(is.na(precision)){
out <- .C("poweranalysispopulationfixed", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.integer(maxiterations), as.integer(feedback), numeric(maxiterations), PACKAGE="bayescount")
lo <- length(out)
meancounts <- out[[lo]]
nin <- sum(meancounts >= lowerl & meancounts <= upperl)
nout <- length(meancounts)-nin
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output=list(roundedci=round(qbeta(c(lci, 0.5, uci), nin+1, nout+1), precision), ci=qbeta(c(lci, 0.5, uci), nin+1, nout+1), within=nin, without=nout, total=nin+nout)
}else{
out <- .C("poweranalysispopulation", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.numeric(lowerl), as.numeric(upperl), as.integer(maxiterations), as.integer(precision), as.numeric(lci), as.numeric(uci), as.integer(feedback), as.integer(0), as.integer(0), PACKAGE="bayescount")
lo <- length(out)
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output=list(roundedci = round(qbeta(c(lci, 0.5, uci), out[[lo-1]]+1, (out[[lo]]-out[[lo-1]])+1), precision), ci=qbeta(c(lci, 0.5, uci), out[[lo-1]]+1, (out[[lo]]-out[[lo-1]])+1), within=out[[lo-1]], without=out[[lo]]-out[[lo-1]], total=out[[lo]])
}
}
}
if(approximate){
names(output$roundedci) <- c("absolute.value", "absolute.value", "absolute.value")
names(output$ci) <- c("absolute.value", "absolute.value", "absolute.value")
}else{
names(output$roundedci) <- c(paste("lower.", confidence*100, "%", sep=""), "median", paste("upper.", confidence*100, "%", sep=""))
names(output$ci) <- c(paste("lower.", confidence*100, "%", sep=""), "median", paste("upper.", confidence*100, "%", sep=""))
}
return(output)
}
FEC.power.limits <- fec.power.limits
FEC.precision <- fec.power.limits
fec.precision <- fec.power.limits
count.precision <- fec.precision
FEC.power <- fec.power
count.power <- fec.power
Try the bayescount package in your browser
Any scripts or data that you put into this service are public.
bayescount documentation built on May 7, 2022, 5:05 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions fec.power fec.power.limits
Documented in fec.power fec.power.limits
fec.power.limits <- function(meanepg=200, g.faeces=3, sensitivity=1/25, replicates=1, animals=10, coeffvarrep=0.4, coeffvarind=0.3, coeffvargroup=0.7, true.sample=FALSE, lower.limit=NA, upper.limit=NA, iterations=100000, power = 0.95, confidence = 0.99, feedback=FALSE, forcesim=FALSE){
# Power is the required power, confidence is the confidence in this power when using simulations
if(power >= 1) stop("Required power must be < 1")
if(confidence >= 1) stop("Confidence must be < 1")
conf <- (1-confidence)/2
lci <- 0+conf
uci <- 1-conf
if(replicates < 1 | animals < 1) stop("Specified values for animals and replicates must be greater than 0")
if(!is.na(lower.limit) & !is.na(upper.limit)) stop("One or both of lower.limit or upper.limit must be non-fixed (NA)")
fix.lower <- !is.na(lower.limit)
fix.upper <- !is.na(upper.limit)
fixed.lower <- lower.limit
fixed.upper <- upper.limit
target <- power
if(feedback){
if(any(.Platform$GUI == c("AQUA", "Rgui"))){
#feedback <- FALSE
warning("Printing the progress of the function using the GUI versions of R may massively increase the time taken, I suggest setting feedback=FALSE or using a command line version of R instead")
}
##### CHECK ON WINDOWS / UNIX
}
# If we have 1 animal and true.sample=T then we can forget about cv.group and use nbpower
if(animals==1 & true.sample==TRUE) coeffvargroup <- 10^-10
if(coeffvargroup > (coeffvarrep+coeffvarind)/10) approximate <- FALSE else approximate <- TRUE
if(forcesim) approximate <- FALSE
if(animals==1 & true.sample==FALSE) warning("NOTE: The power calculated is for the population from which the animal is derived, to calculate the power for the true mean of this individual use true.sample=TRUE")
lowerl <- lower.limit
upperl <- upper.limit
start <- Sys.time()
if(!approximate){
if(true.sample){
out <- .C("poweranalysissamplefixed", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.integer(iterations), as.integer(feedback), numeric(iterations), PACKAGE="bayescount")
lo <- length(out)
meancounts <- out[[lo]]
}else{
out <- .C("poweranalysispopulationfixed", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.integer(iterations), as.integer(feedback), numeric(iterations), PACKAGE="bayescount")
lo <- length(out)
meancounts <- out[[lo]]
}
mcs <- meancounts
f <- function(limit){
prob <- limit
#limits <- HPDinterval(mcmc(mcs), prob=prob)[1,]
limits <- quantile(mcs, probs=c(0+((1-prob)/2), 1-((1-prob)/2)))
nin <- sum(mcs <= limits[2] & mcs >= limits[1])
nout <- length(mcs)-nin
med <- qbeta(0.5, nin+1, nout+1)
return(med)
}
if(fix.lower){
f <- function(limit){
#limits <- HPDinterval(mcmc(mcs), prob=prob)[1,]
nin <- sum(mcs <= limit & mcs >= fixed.lower)
nout <- length(mcs)-nin
med <- qbeta(0.5, nin+1, nout+1)
return(med)
}
}
if(fix.upper){
f <- function(limit){
#limits <- HPDinterval(mcmc(mcs), prob=prob)[1,]
nin <- sum(mcs <= fixed.upper & mcs >= -limit)
nout <- length(mcs)-nin
med <- qbeta(0.5, nin+1, nout+1)
return(med)
}
}
limits <- c(0,1)
if(fix.upper) limits <- c(-Inf,0)
if(fix.lower) limits <- c(0,Inf)
bsres <- binary.search(f, target, limits)
if(bsres$status!="OK"){
if(bsres$status!="Absolute value not possible but this is closest") stop(bsres$status)
}
if(fix.upper) limits <- c(-bsres$value, fixed.upper)
if(fix.lower) limits <- c(fixed.lower, bsres$value)
if(!fix.upper & !fix.lower) limits <- quantile(mcs, prob=c(0+((1-bsres$value)/2), 1-((1-bsres$value)/2)))
nin <- sum(mcs <= limits[2] & mcs >= limits[1])
nout <- length(mcs)-nin
power <- qbeta(c(lci,0.5,uci), nin+1, nout+1)
names(limits) <- c("lower.limit", "upper.limit")
names(power) <- c(paste("lower.", confidence*100, "%.estimate", sep=""), "median.estimate", paste("upper.", confidence*100, "%.estimate", sep=""))
return(list(limits=limits, power=power))
}else{
#Embed within wrapper function (that calculates cvw from rep, sample, g.faeces etc and works out if approximation to multiple animals can be used):
if(true.sample==FALSE | animals > 1){
# If we're using this with true.sample=FALSE then cvgroup must be so small compared to cvind and cvrep that we can just pretend that we're drawing all replicates from the same animal (with a slightly larger cv and animalmean=populationmean):
coeffvarind <- sqrt(coeffvarind^2 + coeffvargroup^2 + coeffvarind^2*coeffvargroup^2)
replicates <- replicates*animals
}
coeff.var <- sqrt(coeffvarind^2 + (coeffvarrep^2)/g.faeces + (coeffvarind^2 * (coeffvarrep^2)/g.faeces))
eggs.counted <- replicates*meanepg*sensitivity
shape <- replicates / coeff.var^2
f <- function(limit){
accuracy <- limit
upper.tol <- replicates*sensitivity*(meanepg+meanepg*accuracy)
lower.tol <- replicates*sensitivity*(meanepg-meanepg*accuracy)
return(pnbinom(upper.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) - pnbinom(lower.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) + suppressWarnings(dnbinom(lower.tol, size=shape, mu=eggs.counted)))
}
if(fix.lower){
f <- function(limit){
upper.tol <- replicates*sensitivity*(limit)
lower.tol <- replicates*sensitivity*(fixed.lower)
return(pnbinom(upper.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) - pnbinom(lower.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) + suppressWarnings(dnbinom(lower.tol, size=shape, mu=eggs.counted)))
}
}
if(fix.upper){
f <- function(limit){
upper.tol <- replicates*sensitivity*(fixed.upper)
lower.tol <- replicates*sensitivity*(-limit)
return(pnbinom(upper.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) - pnbinom(lower.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) + suppressWarnings(dnbinom(lower.tol, size=shape, mu=eggs.counted)))
}
}
limits <- c(0,Inf)
if(fix.upper) limits <- c(-Inf,0)
if(fix.lower) limits <- c(0,Inf)
bsres <- binary.search(f, target, limits)
if(bsres$status!="OK"){
if(bsres$status!="Absolute value not possible but this is closest") stop(bsres$status)
}
if(fix.upper) limits <- c(-bsres$value, fixed.upper)
if(fix.lower) limits <- c(fixed.lower, bsres$value)
if(!fix.upper & !fix.lower) limits <- c((meanepg-meanepg*bsres$value), (meanepg+meanepg*bsres$value))
power <- replicate(3, bsres$objective)
names(limits) <- c("lower.limit", "upper.limit")
names(power) <- c("absolute.value", "absolute.value", "absolute.value")
return(list(limits=limits, power=power))
}
}
fec.power <- function(meanepg=200, g.faeces=3, sensitivity=1/25, replicates=1, animals=10, coeffvarrep=0.4, coeffvarind=0.3, coeffvargroup=0.7, true.sample=FALSE, accuracy=0.1, lower.limit=meanepg*(1-accuracy), upper.limit=meanepg*(1+accuracy), maxiterations=1000000, precision=2, confidence = 0.99, feedback=FALSE, forcesim=FALSE){
if(confidence >= 1) stop("Confidence must be < 1")
conf <- (1-confidence)/2
lci <- 0+conf
uci <- 1-conf
if(feedback){
if(any(.Platform$GUI == c("AQUA", "Rgui"))){
#feedback <- FALSE
warning("Printing the progress of the function using the GUI versions of R may massively increase the time taken, I suggest setting feedback=FALSE or using a command line version of R instead")
}
##### CHECK ON WINDOWS / UNIX
}
if(replicates < 1 | animals < 1) stop("Specified values for animals and replicates must be greater than 0")
# If we have 1 animal and true.sample=T then we can forget about cv.group and use nbpower
if(animals==1 & true.sample==TRUE) coeffvargroup <- 10^-10
if(coeffvargroup > (coeffvarrep+coeffvarind)/10) approximate <- FALSE else approximate <- TRUE
if(forcesim) approximate <- FALSE
if(animals==1 & true.sample==FALSE) warning("NOTE: The power calculated is for the population from which the animal is derived, to calculate the power for the true mean of this individual use true.sample=TRUE")
lowerl <- lower.limit
upperl <- upper.limit
start <- Sys.time()
if(approximate){
if(true.sample==FALSE){
# If we're using this with true.sample=FALSE then cvgroup must be so small compared to cvind and cvrep that we can just pretend that we're drawing all replicates from the same animal (with a slightly larger cv and animalmean=populationmean):
coeffvarind <- sqrt(coeffvarind^2 + coeffvargroup^2 + coeffvarind^2*coeffvargroup^2)
replicates <- replicates*animals
}
# Produces the same results as poweranalysispopulation (only for animal=1 and/or cv.group<<cv.other)
coeff.var <- sqrt(coeffvarind^2 + (coeffvarrep^2)/g.faeces + (coeffvarind^2 * (coeffvarrep^2)/g.faeces))
eggs.counted <- replicates*meanepg*sensitivity
upper.tol <- replicates*sensitivity*(upperl)
lower.tol <- replicates*sensitivity*(lowerl)
shape <- replicates / coeff.var^2
nbpower <- pnbinom(upper.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) - pnbinom(lower.tol, size=shape, mu=eggs.counted, lower.tail=TRUE) + suppressWarnings(dnbinom(lower.tol, size=shape, mu=eggs.counted))
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output = list(roundedci=replicate(3, round(nbpower, precision)), ci=replicate(3, nbpower))
}else{
if(true.sample){
if(is.na(precision)){
out <- .C("poweranalysissamplefixed", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.integer(maxiterations), as.integer(feedback), numeric(maxiterations), PACKAGE="bayescount")
lo <- length(out)
meancounts <- out[[lo]]
nin <- sum(meancounts >= lowerl & meancounts <= upperl)
nout <- length(meancounts)-nin
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output=list(roundedci=round(qbeta(c(lci, 0.5, uci), nin+1, nout+1), precision), ci=qbeta(c(lci, 0.5, uci), nin+1, nout+1), within=nin, without=nout, total=nin+nout)
}else{
out <- .C("poweranalysissample", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.numeric(lowerl), as.numeric(upperl), as.integer(maxiterations), as.integer(precision), as.numeric(lci), as.numeric(uci), as.integer(feedback), as.integer(0), as.integer(0), PACKAGE="bayescount")
lo <- length(out)
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output=list(roundedci = round(qbeta(c(lci, 0.5, uci), out[[lo-1]]+1, (out[[lo]]-out[[lo-1]])+1), precision), ci=qbeta(c(lci, 0.5, uci), out[[lo-1]]+1, (out[[lo]]-out[[lo-1]])+1), within=out[[lo-1]], without=out[[lo]]-out[[lo-1]], total=out[[lo]])
}
}else{
# coeffvarrep is the cv between subsamples taken from the same replicate and is effected by gfaeces (calculation done in C)
# coeffvarind is cv between replicates and is effected by replicates. This actually means it is EITHER the combined cv between samples and days within an animal OR just the cv between samples if taken from a single faecal pile (in which case the cv between days is actually part of the cv between animals and will consequently make things a bit wrong due to the true sample effect). Doesn't reflect subsamples (what I usually call replicates) as this is exactly the same as changing the sensitivity and g.faeces.
# coeffvargroup is cv between animals and is always calculated explicitly since there is no shortcut. If we have only 1 animal or coeffvargroup=0, then powernb works perfectly.
if(is.na(precision)){
out <- .C("poweranalysispopulationfixed", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.integer(maxiterations), as.integer(feedback), numeric(maxiterations), PACKAGE="bayescount")
lo <- length(out)
meancounts <- out[[lo]]
nin <- sum(meancounts >= lowerl & meancounts <= upperl)
nout <- length(meancounts)-nin
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output=list(roundedci=round(qbeta(c(lci, 0.5, uci), nin+1, nout+1), precision), ci=qbeta(c(lci, 0.5, uci), nin+1, nout+1), within=nin, without=nout, total=nin+nout)
}else{
out <- .C("poweranalysispopulation", as.numeric(meanepg), as.numeric(g.faeces), as.numeric(sensitivity), as.integer(replicates), as.integer(animals), as.numeric(coeffvarrep), as.numeric(coeffvarind), as.numeric(coeffvargroup), as.numeric(lowerl), as.numeric(upperl), as.integer(maxiterations), as.integer(precision), as.numeric(lci), as.numeric(uci), as.integer(feedback), as.integer(0), as.integer(0), PACKAGE="bayescount")
lo <- length(out)
time <- timestring(start, Sys.time(), units="s", show.units=FALSE)
output=list(roundedci = round(qbeta(c(lci, 0.5, uci), out[[lo-1]]+1, (out[[lo]]-out[[lo-1]])+1), precision), ci=qbeta(c(lci, 0.5, uci), out[[lo-1]]+1, (out[[lo]]-out[[lo-1]])+1), within=out[[lo-1]], without=out[[lo]]-out[[lo-1]], total=out[[lo]])
}
}
}
if(approximate){
names(output$roundedci) <- c("absolute.value", "absolute.value", "absolute.value")
names(output$ci) <- c("absolute.value", "absolute.value", "absolute.value")
}else{
names(output$roundedci) <- c(paste("lower.", confidence*100, "%", sep=""), "median", paste("upper.", confidence*100, "%", sep=""))
names(output$ci) <- c(paste("lower.", confidence*100, "%", sep=""), "median", paste("upper.", confidence*100, "%", sep=""))
}
return(output)
}
FEC.power.limits <- fec.power.limits
FEC.precision <- fec.power.limits
fec.precision <- fec.power.limits
count.precision <- fec.precision
FEC.power <- fec.power
count.power <- fec.power
Try the bayescount package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.