Freecalc population sensitivity with imperfect sensitivity and specificity
Analyse the results of 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;
- Sample size tested;
- Number tested positive;
- 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 - cluster-level sensitivity) and Type II (1 - cluster-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; and
- The desired precision of results (number of digits to be displayed after the decimal point).
The results are displayed in terms of the null and alternate hypothesis:
- The probability of the null hypothesis is the probability of observing this many reactors or fewer, if the population was diseased at a level equal to or greater than the specified design prevalence. If this probability is small, we can conclude that t is very unlikely that the popultion is diseased. If the probability is large, then there is not enough evidence to conclude that the population is free from disease.
- The probability of the alternative hypothesis is the probability of observing this many reactors or more if the poplation was truly disease free. If this is small, then it is very unlikely that the population is free from disease. If it is large, then it is connsistent with there being no disease in the population.
- If both null and alternative probabilities are small, it suggests that the population is not free from disease, but the prevalence is less than the design prevalence specified.
- If both null and alternative probabilities are large, The sample size was too small to distinguish a population with the specified design prevalence from a disease-free population.
No results
No example available
No references available
######################################
# Program to do freecalc calculations
######################################
# 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:18])
# cat("
", a1)
N<- a1[1] # population size
# cat("
", N)
n<- a1[2] # sample size
n.pos<- a1[3] # observed number positive
se<- a1[4] # test sensitivity
sp<- a1[5] # test specificity
dis<- a1[6] # number/prevalence of diseased elemnets
prev.type<- ifelse(dis < 1, 1, 0) # 0 = number of diseased elements, 1 = prevalence
type1<- a1[7] # type I error
type2<- a1[8] # type II error
method<- a1[9] # calculation method 0 = hypergeo, 1 = binom approx, 2 = simple binom
pop.threshold<- a1[10] # Infinite population threshold
option<- 1 # calculation option - 0 = sample size, 1 = analyse data
digits<- a1[11] # number of digits in results
# email<- a0[22]
prec<- 10^-digits
prev<- ifelse(prev.type, dis, dis/N)
dis<- ifelse(prev.type, max(1, round(dis*N, 0)), dis)
method.lst<- c("Modified hypergeometric exact", "Modified binomial approximation", "Simple binomial (large population)")
heading<- c("FreeCalc sample size estimation", "FreeCalc data analysis")
# cat("
method:", method)
filename<- digest(Sys.time)
# data.file<- ifelse(test, paste(fpath, "ts_data_1.txt", sep = ""), a0[13]) #test name
tmp.path<- paste(fpath, "tmp/", sep = "")
tmp.file<- paste(fpath, "tmp/", filename, sep = "")
sinkfile<- paste(fpath, "tmp/", filename, ".txt", sep="") # fpath,
# cat("
", data.file)
inputs<- array(0, dim = c(11, 1))
rownames(inputs)<- c("Population size", "Sample size", "Number positive",
"Test sensitivity", "Test specificity", "Design prevalence",
"Diseased elements", "Analysis method", "Target Type I error",
"Target Type II error",
"Population threshold for infinite probability formula")
inputs[1:5,1]<- a1[1:5]
inputs[6,1]<- prev
inputs[7,1]<- dis
inputs[8,1]<- method.lst[method+1]
inputs[9:11,1]<- a1[c(7, 8, 10)]
# cat("
Test: 0")
#
# analyse data
if (method == 2 || N > pop.threshold) { # simple binomial
P1<- 1-sep.binom.imperfect(n, n.pos+1, se, sp, prev)
} else { # modified hypergeometric exact
P1<- 1-sep.freecalc(N,n,n.pos+1,se,sp,prev)
}
P2<- 1-sph.binom(n, n.pos, sp)
Px<- 1-sph.binom(n, n.pos+1, sp)
SpH<- ifelse(1-P2 < prec, paste("<", format(prec, nsmall = digits, scientific = F)), format(round(1-P2, digits), scientific = F))
SpHx<- ifelse(1-Px < prec, paste("<", format(prec, nsmall = digits, scientific = F)), format(round(1-Px, digits), scientific = F))
SeH<- ifelse(1-P1 < prec, paste("<", format(prec, nsmall = digits, scientific = F)), format(round(1-P1, digits), scientific = F))
P1.txt<- ifelse(P1 < prec, paste("<", format(prec, nsmall = digits, scientific = F)), format(round(P1, digits), scientific = F))
P2.txt<- ifelse(P2 < prec, paste("<", format(prec, nsmall = digits, scientific = F)), format(round(P2, digits), scientific = F))
conc<- c(paste("These results are not adequate to reject the null hypothesis and conclude that the population is free from disease. The sample size was too small to distinguish a population with prevalence of ", prev*100, "% from a disease-free population.", sep = ""),
paste("These results are adequate to reject the null hypothesis and conclude that the population is free from disease (at the expected minimum prevalence of ", prev*100, "%) at the ", SeH, " confidence level.", sep = ""),
paste("These results are not adequate to conclude that the population is free from disease (at the expected minimum prevalence of ", prev*100, "%). The confidence level is only ", SeH, ". We may conclude that the population is diseased at a confidence level of ", SpH, "", sep = ""),
paste("These results are adequate to reject the null hypothesis but not to accept the alternate hypothesis. We may conclude that the population is not free from disease, but that the prevalence is less than ", prev*100, "%", sep = ""))
result<- array("", dim = c(6, 1))
rownames(result)<- c("Null hypothesis:", "Alternative hypothesis:", "Conclusion:", "Method:", "Cluster-sensitivity (SeH)", "Cluster-specificity (SpH)")
result[1,1]<- paste("Probability of observing <= ", n.pos, " reactors in a sample of ",
n, " individuals from a population with a disease prevalence of ", round(prev*100, 1), "% = ", P1.txt, ".", sep = "")
result[2,1]<- paste("Probability of observing >= ", n.pos, " reactors in a sample of ",
n, " individuals from a disease free population = ", P2.txt, ".", sep = "")
result[3,1]<- ifelse(P1 > type1 && P2 > type2, conc[1], ifelse(P1 <= type1 && P2 <= type2, conc[4], ifelse(P1<= type1, conc[2], conc[3])))
result[4,1]<- ifelse(N > pop.threshold, method.lst[3], method.lst[method+1])
result[5,1]<- SeH
result[6,1]<- SpHx
seh.1<- 0
sph.1<- 0
# res1<- array(0, dim = c(n+1, 4))
# colnames(res1)<- c("Cut-point", "SeH", "SpH", "Success")
# res1[,1]<- 0:n
success<- 0
for (cp in 0:n) {
if (identical(method, 2) || N > pop.threshold) { # simple binomial
P1<- 1-sep.binom.imperfect(n, cp+1, se, sp, prev)
} else { # modified hypergeometric exact
P1<- 1-sep.freecalc(N,n,cp+1,se,sp,prev)
}
# P2<- SeHSimpleBin(n, cp, se, sp, 0)
Px<- 1-sph.binom(n, cp+1, sp)
# SpH<- ifelse(1-P2 < prec, paste("<", format(prec, nsmall = digits, scientific = F)), format(round(1-P2, digits), scientific = F))
seh.1<- 1 - P1 # ifelse(1-P1 < prec, paste("<", format(prec, nsmall = digits, scientific = F)), format(round(1-P1, digits), scientific = F))
sph.1<- 1 - Px # ifelse(1-Px < prec, paste("<", format(prec, nsmall = digits, scientific = F)), format(round(1-Px, digits), scientific = F))
if (Px < type2) {
if (P1 < type1) success<- 1
break
}
}
# tmp1<- analyse.data(N, n, n.pos, se, sp, dis, prev, prec, method, method.lst)
# result<- tmp1[[1]]
res1<- round(c(cp+1, seh.1, sph.1, success), digits)
res2<- array("", dim = c(3, 1))
res2[2:3]<- res1[2:3]
res2[1]<- paste(">=", res1[1], " positives required to consider herd positive.", sep = "")
colnames(res2)<- paste("Optimum cut-point, SeH and SpH for given population and sample sizes and target Type I and Type II errors.")
rownames(res2)<- c("Cut-point", "SeH", "SpH")
if (res1[4]) {
tmp2<- paste("Target herd sensitivity and target herd specificity achieved, assuming >=", res1[1], "reactor(s) required to consider herd positive." )
} else {
tmp2<- paste("Unable to achieve both target herd sensitivity and target herd specificity with given sample size. Maximum achieved herd sensitivity for herd specificity >=", 1 - type2, "is", res1[2], "assuming >=", res1[1], "reactor(s) required to consider herd positive." )
}
res2<- rbind(res2, "Conclusion" = tmp2)
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)
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"))
output<- paste(output, HTMLStream(res2, cellborder = 0, classfirstline = "mbg", classfirstcolumn = "mbg", classcellinside = "left_mar", cellalign = "left", align = "left"))
output<- paste(output, "Download excel file of results
\n")
output<- paste(output, "Detailed results
\n", sep = "")
# write output to file
file.name<- file(paste(tmp.file, "_result.html", sep = ""), open = "wt")
cat(output, file=file.name)
close(file.name) #)
# cat("
Binomial approximation incorrect")
cat(output)
