FreeCalc: Calculate sample size for freedom testing with imperfect tests
Calculate the required sample size and cut-point for testing to demonstrate population freedom from disease using imperfect tests an allowing for small populations.
This utility uses the methods described by:
Cameron and Baldock (1998):
A new probability formula for surveys to substantiate freedom from disease.
Prev. Vet. Med. 34:1-17 and
Cameron (1999): Survey Toolbox for Livestock
Diseases - A practical manual and software package for active surveillance of livestock
diseases in developing countries. Australian Centre for International Agricultural
Research, Canberra, Australia.
These methods are also the same as those used in the
FreeCalc Program.
Inputs include:
- Size of the population sampled;
- Test sensitivity and specificity;
- Design prevalence (the hypothetical prevalence to be detected). Design prevalence can be specified as either a fixed number of elements from the population or a proportion of the population;
- Type I (1 - population-sensitivity) and Type II (1 - population-specificity) error values for determining whether to accept/reject the null or alternative hypothesis;
- Calculation method: hypergeometric (for small populations), or simple binomial (for large populations);
- The population size threshold, above which the simple binomial method is used regardless of which calculation method has been selected;
- The maximum upper limit for required sample size; and
- The desired precision of results (number of digits to be displayed after the decimal point).
The results are presented as:
- The minimum sample size and corresponding cut-point number of positives to achieve the specified type I and type II errors for the given population, design prevalence and test performance;
- achieved type I and Type II error levels and corresponding population-level sensitivities and specificities;
- A descriptive interpretation of the results; and
- an error message if the desired error levels cannot be achieved within the limits of populatuon and/or maximum sample size.
No results
No example available
No references available
######################################
# Program to do freecalc calculations for sample size
######################################
# uses RSurveillance package
x<- 1
rm(list = ls())
test<- ifelse(length(commandArgs()) == 2, TRUE, FALSE)
fpath<- ifelse(test, "webRootUrl", "rtoolsPath")
# cat("
Test:",length(commandArgs()))
# load header scripts
source(paste(fpath, "R/epi_head.R", sep = ""))
source(paste(fpath, "R/HTMLStream.R", sep = ""))
source(paste(fpath, "R/epitools_functions.r", sep = ""))
# extract command arguments 1 = se, 2 = sp, 3 = population size, 4 = number of diseased elements,
# 5 = sample size, 6 = number positive,
# 7 = method, 8 = type1 error, 9 = type2 error, 10 = max sample size, 11 = infinite population threshold
# method 1 = modified hypergeomertric, 2 = modified binomial, 3 = simple binomial
# cat("
", a0[20:22])
a1<- type.convert(a0[8:17])
# cat("
", a1)
N<- a1[1] # population size
se<- a1[2] # test sensitivity
sp<- a1[3] # test specificity
dis<- a1[4] # number/prevalence of diseased elemnets
prev.type<- ifelse(dis < 1, 1, 0) # 0 = number of diseased elements, 1 = prevalence
# n<- a1[5] # sample size
# n.pos<- a1[6] # observed number positive
type1<- a1[5] # type I error
type2<- a1[6] # type II error
method<- a1[7] # calculation method 0 = hypergeo, 1 = binom approx, 2 = simple binom
pop.threshold<- a1[8] # Infinite population threshold
max.ss<- a1[9] # maximum sample size
option<- 0 # calculation option - 0 = sample size, 1 = analyse data
digits<- a1[10] # number of digits in results
prec<- 10^-digits
prev<- ifelse(prev.type, dis, dis/N)
dis<- ifelse(prev.type, max(1, round(dis*N, 0)), dis)
method<- ifelse(N > pop.threshold, 2, method)
method.lst<- c("Modified hypergeometric exact", "Modified binomial approximation", "Simple binomial (large population)")
heading<- c("FreeCalc sample size estimation", "FreeCalc data analysis")
filename<- digest(Sys.time)
tmp.path<- paste(fpath, "tmp/", sep = "")
tmp.file<- paste(fpath, "tmp/", filename, sep = "")
sinkfile<- paste(fpath, "tmp/", filename, ".txt", sep="") # fpath,
# cat(46)
inputs<- array("", dim = c(11, 1))
inputs[1:3,1]<- a1[1:3]
inputs[4,1]<- prev
inputs[5,1]<- dis
inputs[6:7,1]<- a1[5:6]
inputs[8,1]<- method.lst[method+1]
inputs[9:11,1]<- a1[8:10]
rownames(inputs)<- c("Population size", "Test sensitivity", "Test specificity", "Design prevalence",
"Diseased units", "Target Type I error", "Target Type II error", "Analysis method",
"Population threshold for binomial default", "Maximum sample size", "Digits")
# cat("
method:", method)
# function to get cut-point for given sample size
get.cp<- function(n.cp, se, sp, type2) {
cp<- 0
SpH<- 0
while (SpH < 1-type2) {
cp<- cp+1 # probability of observed result from diseas-free popn
SpH<- sph.binom(n.cp, cp, sp)
}
return (list("cp" = cp, "SpH" = SpH))
} # end of get.cp
# function to do sample size calcuations
get.ss<- function(N, se, sp, prev, type1, type2, dis, method, method.lst, maxSS) {
N1<- min(N, maxSS)
brks<- c(50, 100, 1000, 5000, 10000, Inf)
steps<- c(5, 10, 50, 100, 200)
step<- steps[which(N1 < brks)[1]]
ss<- seq(0, N1, by = step)
ss[1]<- 1
if (length(ss) == 1) ss[2]<- N1
cp<- 0
SpH<- 0
SeH<- 0
P1<- 0
success<- F
for (s in 1:length(ss)) {
tmp<- get.cp(ss[s], se, sp, type2)
cp[s]<- tmp[[1]]
SpH[s]<- tmp[[2]]
if (identical(method, 2) || N > pop.threshold) { # simple binomial
P1[s]<- 1-sep.binom.imperfect(ss[s], cp[s], se, sp, prev)
} else { # modified hypergeometric exact
P1[s]<- 1-sep.freecalc(N,ss[s], cp[s],se,sp,prev)
}
SeH[s]<- 1-P1[s]
cp[s]<- cp[s]-1
if (P1[s] <= type1) {
success<- T
n1<- ss[s]
break
}
} # end of s loop
if (success) {
ss[(s+1):(s+step)]<- (ss[s-1]+1):(ss[s-1]+step)
for (x in 1:step) {
tmp<- get.cp(ss[s+x], se, sp, type2)
cp[s+x]<- tmp[[1]]
SpH[s+x]<- tmp[[2]]
if (identical(method, 2) || N > pop.threshold) { # simple binomial
P1[s+x]<- 1-sep.binom.imperfect(ss[s+x], cp[s+x], se, sp, prev)
} else { # modified hypergeometric exact
P1[s+x]<- 1-sep.freecalc(N,ss[s+x], cp[s+x],se,sp,prev)
}
SeH[s+x]<- 1-P1[s+x]
cp[s+x]<- cp[s+x]-1
if (P1[s+x] <= type1) {
success<- T
n1<- ss[s+x]
break
}
} # end of x loop
result<- array("", dim = c(8, 1))
rownames(result)<- c("Required sample size:", "Cut-point number of positives:",
"Type I error:", "Type II error:", "Population-level sensitivity:",
"Population-level specificity:", "Interpretation:", "Method:")
result[1,1]<- n1
result[2,1]<- cp[s+x]
result[3,1]<- round(1-SeH[s+x], digits)
result[4,1]<- round(1-SpH[s+x], digits)
result[5,1]<- round(SeH[s+x], digits)
result[6,1]<- round(SpH[s+x], digits)
result[7,1]<- paste("If a random sample of ", n1, " units is taken from a population of ",
N, " and ", cp[s+x], " or fewer reactors are found, the probability that the population is diseased at a prevalence of ",
prev, " is ", round(1-SeH[s+x], digits), ".", sep = "")
result[8,1]<- method.lst[method+1]
return(list(success, "result" = result, "n" = ss[1:(s+x)], "seh" = SeH[1:(s+x)], "sph" = SpH[1:(s+x)], "cp" = cp[1:(s+x)], P1[1:(s+x)], "sam.size" = (s+x)))
} else {
result<- array("", dim = c(1, 1))
result[1,1]<- ifelse(identical(N1, N), "Unable to achieve required accuracy by sampling every unit",
paste("Unable to achieve required accuracy by sampling every unit up to the maximum sample size of ", maxSS, ".", sep = ""))
return(list(success, "result" = result))
} # end of if/else
# return(list("seh"=SeH, "sph"=SpH, "cp"=cp, "sam.size"=n1, "success"=success))
} # end of get.ss function
# cat("
", N)
tmp<- get.ss(N, se, sp, prev, type1, type2, dis, method, method.lst, max.ss)
# cat("
", test)
success<- tmp[[1]]
result<- tmp[[2]]
if (success) {
ss<- tmp[[3]]
SeH<- round(tmp[[4]], digits)
SpH<- round(tmp[[5]], digits)
cp.list<- tmp[[6]]
P1<- round(tmp[[7]], digits)
n<- tmp[[8]]
cp<- cp.list[n]
result1<- cbind("n" = ss, SeH, SpH, "Cut-point" = cp.list, "Type 1 error" = P1, "Type 2 error" = 1-SpH)
}
tmp<- gsub("", "", result)
result.txt<- gsub("", "", tmp)
# cat("
Test:", method)
file.name<- file(paste(tmp.file, "_result.xls", sep = ""), open = "wt")
writeLines(c(heading[option+1], date(), ""), con = file.name)
writeLines("Inputs", con = file.name)
write.table(inputs, file = file.name, sep = "\t", append = T, col.names = F)
writeLines(c("", "Results"), con = file.name)
write.table(result.txt, file = file.name, sep = "\t", append = T, col.names = F, row.names = T)
if (success) {
writeLines(c("", "Detailed Results"), con = file.name)
write.table(result1, file = file.name, sep = "\t", append = T, col.names = F, row.names = T)
}
close(file.name)
# output results
d1<- paste(substr(date(), 1, 10), substr(date(), 20,24), " @", substr(date(), 11, 16))
# reformat headers
output<- paste("", heading[option+1], "
\n")
output<- paste(output, "", "Analysed: ", d1, "
\n", sep="")
output<- paste(output, "Inputs
\n")
output<- paste(output, HTMLStream(inputs, cellborder = 0, classfirstline = "mbg", classfirstcolumn = "mbg", classcellinside = "left_mar", cellalign = "center", align = "left"))
output<- paste(output, "Results
\n")
output<- paste(output, HTMLStream(result, cellborder = 0, classfirstline = "mbg", classfirstcolumn = "mbg", classcellinside = "left_mar", cellalign = "left", align = "left"))
if (success) {
output<- paste(output, "Download excel file of results
\n")
output<- paste(output, "Detailed results
\n", sep = "")
}
# email.results(email, output, fpath, filename, "freecalc2")
cat(output)
