Population sensitivity with imperfect test sensitivity and specificity

This utility estimates the probability of detecting disease herd or cluster-sensitivity in a large (infinite) population, if it is present at the specified design prevalence, assuming a test of known sensitivity and specificity and for a variable cut-point number fo positives to determine the test result. These analyses use the method from Martin et el. (1992) (Prev Vet Med, 14:33-43) and the binomial distribution function, assuming known test sensitivity and test specificity and a variable cut-point number of positives to declare a population infected (i.e. a variable (non-zero) number of positive positives can be allowed and still be recognised as free). The population is classified as infected if the number of positive results is equal to or greater than the cut-point.

Inputs are the sample size tested, test sensitivity, test specificity, the design (target) prevalence and the cut-point number of positives.

Outputs are:

the cluster-level sensitivity and specificity for the given sample size, test sensitivity, test specificity, design prevalence and cut-point number of positives; and
tables and graphs of cluster-level sensitivity and specificity values for a range of cut-point, sample size and design prevalence values.

No results

No example available

No references available

				######################################
# Program to calculate herd/flock level sensitivity and specificity
######################################
# uses RSurveillance

rm(list = ls())
test<- ifelse(length(commandArgs()) < 3, TRUE, FALSE)
fpath<- ifelse(test, "webRootUrl", "rtoolsPath")

# 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
    a1<- type.convert(a0[8:12])
# cat(a1)
digits<- 4
Prevalence<- a1[1]
Sens<- a1[2]
Spec<- a1[3]
# names(a1)<- c("Prevalence", "Sensitivity", "Specificity", "Sample Size")
s.size<- a1[4]               # number of sample sizes excluding user input
CP<- a1[5]
filename<- digest(Sys.time)
graphfile<- paste(fpath, "tmp/", filename, "a.png", sep="")
graph2<- paste(fpath, "tmp/", filename, "b.png", sep="")
graph3<- paste(fpath, "tmp/", filename, "c.png", sep="")
sinkfile<- paste(fpath, "tmp/", filename, ".txt", sep="")

# check for valid cut-point
if (CP > s.size) {
    cat("Cut-point number of reactors exceeds sample size. Please try again.")
    quit()
}


sink(sinkfile)

# table of inputs
inputs<- array("", dim = c(length(a1), 1))
inputs[1:length(a1), 1]<- a1[1:length(a1)]
rownames(inputs)<- c("Design prevalence", "Test sensitivity", "Test specificity", "Sample size", "Cut-point number of positives")

sample.size<- c(5, 10, 15, 20, 30, 50, 75, 100, 200, 300, 500)
Prev.list<- c(0.01, 0.02, 0.05, 0.1, 0.2, 0.3)
CP.list<- c(1, 2, 3, 4, 5, 10, 15)
s.sizes<- length(sample.size)
n.Prev<- length(Prev.list)
n.CP<- length(CP.list)

# P.Pos<- Prev.list * Sens + (1 - Prev.list)*(1 - Spec)

# SeH/SpH by cut-point
results<- array(0, dim = c(n.Prev + 1, n.CP))
rownames(results)<- c(paste("SeH - Prevalence = ", Prev.list, sep = ""), "SpH")
colnames(results)<- paste("CP = ", CP.list, sep = "")
for (c in 1:n.CP) {
  results[-nrow(results), c]<- sep.binom.imperfect(s.size, c=CP.list[c]-1, se=Sens, sp=Spec, pstar=Prev.list)  
  results[n.Prev+1, c]<- sph.binom(s.size, c=CP.list[c]-1, Spec) 
}
results<- round(results, digits)

# SeH/SpH by sample size
results1<- array(0, dim = c(n.Prev + 1, s.sizes))
rownames(results1)<- rownames(results)
colnames(results1)<- paste("Sample size = ", sample.size, sep = "")
for (c in 1:s.sizes) {
  results1[-nrow(results), c]<- sep.binom.imperfect(sample.size[c], c=CP, se=Sens, sp=Spec, pstar=Prev.list)  
  results1[n.Prev+1, c]<- sph.binom(sample.size[c], c=CP, Spec) 
}
results1<- round(results1, digits)

# actual SeH/SpH
SeH<- round(sep.binom.imperfect(s.size, c=CP, se=Sens, sp=Spec, pstar=Prevalence), digits)  
SpH<- round(sph.binom(s.size, c=CP, Spec), digits)


# graph results   by cut-point
    Title<-c("Cluster-level sensitivity and specificity")
    line.colours<- c("darkblue", "red", "darkgreen", "magenta", "brown", "purple", "orange")
    OpenGraphOutput(graphfile, pointsize = 12, ht = 6, wd = 8)
    plot(x = CP.list, y = results[1,], type="l", xlab = "Cut-point reactors", ylim = c(0, 1),
                xlim = c(0, max(CP.list)), ylab="SeH/SpH", main=Title, col = line.colours[1])
    legend.txt<- rownames(results)
    for (i in 2:(n.Prev + 1)) {
        lines(CP.list, results[i,], type="l",  col = line.colours[i])
#        legend.txt[i]<- rownames(results)[i]
    }
    legend(max(CP.list), 0.5, legend.txt, xjust = 1, yjust = 0.5,
        col = line.colours, lty = c(1, 1), plot = TRUE, cex = 0.7)
    CloseGraphOutput("B")

# graph SeH by prevalence
    Title<-c("Cluster-level sensitivity")
    OpenGraphOutput(graph2, pointsize = 12, ht = 6, wd = 8)
    plot(x = Prev.list, y = results[1:n.Prev,1], type="l", xlab = "Prevalence", ylim = c(0, 1),
                xlim = c(0, max(Prev.list)), ylab="SeH", main=Title, col = line.colours[1])
    legend.txt<- colnames(results)
    for (i in 2:n.CP) {
        lines(Prev.list, results[1:n.Prev,i], type="l",  col = line.colours[i])
#        legend.txt[i]<- rownames(results)[i]
    }
    legend(max(Prev.list), 0.5, legend.txt, xjust = 1, yjust = 0.5,
        col = line.colours, lty = c(1, 1), plot = TRUE, cex = 0.7)
    CloseGraphOutput("B")

# graph results by sample size
    Title<-c("Cluster-level sensitivity and specificity")
    line.colours<- c("darkblue", "red", "darkgreen", "magenta", "brown", "purple", "orange")
    OpenGraphOutput(graph3, pointsize = 12, ht = 6, wd = 8)
    plot(x = sample.size, y = results1[1,], type="l", xlab = "Sample size", ylim = c(0, 1),
                xlim = c(0, max(sample.size)), ylab="SeH/SpH", main=Title, col = line.colours[1])
    legend.txt<- rownames(results1)
    for (i in 2:(n.Prev + 1)) {
        lines(sample.size, results1[i,], type="l",  col = line.colours[i])
#        legend.txt[i]<- rownames(results)[i]
    }
    legend(max(sample.size), 0.5, legend.txt, xjust = 1, yjust = 0.5,
        col = line.colours, lty = c(1, 1), plot = TRUE, cex = 0.7)
    CloseGraphOutput("B")

sink()

# write to html and file
heading<- "Cluster-level sensitivity for testing with imperfect specificity and variable cut-points"
subheadings<- c("SeH and SpH for varying prevalence and cut-point", "SeH and SpH for varying prevalence and sample size at specified cut-point")
graph.txt<- c("Plots of SeH and SpH", "", "")
tmp.file<- paste(fpath, "tmp/", filename, sep = "")
result.txt<- c(paste("The table and graphs below summarise the cluster-sensitivity and cluster-specificity for various cut-points and prevalence levels, for specified sample size = ", s.size, ", test sensitivity = ", Sens, " and specificity = ", Spec, ".", 
      "For an assumed prevalence in the source population of ", Prevalence, ", and a sample size of ", s.size, ", the cluster-level sensitivity is ", SeH, ".
 For an uninfected population, a test specificity of ", Spec, " and a sample size of ", s.size, ", the cluster-level specificity is ", SpH, ".", sep = ""),
      paste("The table and graphs below summarise the cluster-sensitivity and cluster-specificity for various sample sizes and prevalence levels, for cut-point number of positives = ", CP, ", test sensitivity = ", Sens, " and specificity = ", Spec, ".", sep = ""))
output<- html.output(heading, subheadings, inputs, results = list(results, results1), graphs = c(graphfile, graph2, graph3), graph.headings = graph.txt, show.inputs = T, show.graphs = T, tmp.file, result.txt = result.txt) 
write.html(output, tmp.file)
cat(output)