Compare two tests

This utility can be used to compare the performance of two tests evaluated in the same population. Kappa and proportional agreement for positive and negative results are calculated, as well as McNemar's Chi-squared value. Corresponding p-values are also calculated for both kappa and Chi-squared values. For infected populations, sensitivity values are calculated for both tests, as well as sensitivity covariances. Similarly, if the population is specified as uninfected, specificity values and specificity covariances are calculated.

See also:

Fleiss JL, 1981. Statistical methods for rates and proportions, 2nd Ed. John Wiley and Sons, New York. Chapter 13.
Cicchetti DV, Feinstein AR, 1990. High agreement but low kappa: II. Resolving the paradoxes. J. Clin. Epidemiol. 43: 551-558, for details of kappa(pos) and kappa(neg); and
Gardner IA, Stryhn H, Lind P, Collins MT, 2000. Conditional dependence between tests affects the diagnosis and surveillance of animal diseases. Prev. Vet. Med. 45: 107-122 for calculation and application of test covariances.

No results

No example available

No references available

				#############################################################
# Program to compare 2 tests
#############################################################


# check version and load header script
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<- a0[8:15]

# cat(a0)
# cat(a1)
digits<- 4

# initialise variables
test1<- a1[1]
test2<- a1[2]
test.data<- type.convert(a1[3:6])
conf<- type.convert(a1[7])
z.conf<- qnorm(1-(1-conf)/2)
names(a1)<- c("Test 1", "Test 2", "Cell a (++)", "Cell b (-+)", "Cell c (+-)", "Cell d (--)", "Confidence level")
n<- sum(test.data)
t1.pos<- test.data[1] + test.data[3]
t1.neg<- test.data[2] + test.data[4]
t2.pos<- test.data[1] + test.data[2]
t2.neg<- test.data[4] + test.data[3]
Pop.status<- type.convert(a1[8])          # status of source population 1 = unknown/mixed, 2 = all infected, 3 = all uninfected
p<- matrix(test.data, nrow = 2, byrow = T)/n

# table of inputs
inputs<- matrix(type.convert(a1[3:6]), nrow = 2, ncol = 2, byrow = T)
inputs<- addmargins(inputs)
colnames(inputs)<- paste(a1[1], c("+ve", "-ve", "Total"))
rownames(inputs)<- paste(a1[2], c("+ve", "-ve", "Total"))

# Calculate kappa values
obsAgree<- (test.data[1] + test.data[4])/n
chanceAgree.pos<- (t2.pos/n)*(t1.pos/n)
chanceAgree.neg<- (t2.neg/n)*(t1.neg/n)
chanceAgree.total<-  chanceAgree.pos + chanceAgree.neg
obsAgree.excess<- obsAgree - chanceAgree.total
maxAgree.excess<- 1 - chanceAgree.total
Kappa<- obsAgree.excess/maxAgree.excess
k.pos<- 2*test.data[1]/(2*test.data[1] + test.data[2] + test.data[3])
k.neg<- 2*test.data[4]/(2*test.data[4] + test.data[2] + test.data[3])
McNemar<- ((abs(test.data[2] - test.data[3]) - 1)^2)/(test.data[2] + test.data[3])
p.val.McNemar<- 1 - pchisq(McNemar, 1)
abs.diff.propn<- (test.data[2] - test.data[3])/n
se.abs.diff<- (sqrt((test.data[1] + test.data[4])*(test.data[2] + test.data[3]) + 4*test.data[2] * test.data[3]))/(n*sqrt(n))
rel.diff.propn<- (test.data[2] - test.data[3])/(test.data[2] + test.data[4])
se.rel.diff<- (1/(t1.neg^2))*(sqrt((n - test.data[1])*(test.data[2]*test.data[3] + t1.neg + t2.neg) - test.data[2]*test.data[3]*test.data[4]))
p1<- t2.pos/n
p.1<- t1.pos/n
p2<- t2.neg/n
p.2<- t1.neg/n
temp1<- 1/((1 - chanceAgree.total)*sqrt(n))
temp2<- sqrt(chanceAgree.total + chanceAgree.total^2 - (p1*p.1*(p1 + p.1) + p2*p.2*(p2 + p.2)))
se.kappa<- temp1*temp2
z.kappa<- abs(Kappa/se.kappa)
p.kappa<- 1 - pnorm(z.kappa)
A<- 0
for (i in 1:2) A<- A + p[i,i]*(1 - (sum(p[i,]) + sum(p[,i])) * (1 - Kappa))^2 
B<- (1 - Kappa)^2 * (p[1, 2]*(sum(p[,1]) + sum(p[2,]))^2 + p[2, 1]*(sum(p[,2]) + sum(p[1,]))^2) 
C<- (Kappa - chanceAgree.total*(1 - Kappa))^2
se.kappa.obs<- sqrt(A+B-C)/((1 - chanceAgree.total)*sqrt(n))


lc.kappa<- Kappa - z.conf * se.kappa.obs
uc.kappa<- Kappa + z.conf * se.kappa.obs
lc.abs.diff<- abs.diff.propn - z.conf*se.abs.diff
uc.abs.diff<- abs.diff.propn + z.conf*se.abs.diff
lc.rel.diff<- rel.diff.propn - z.conf*se.rel.diff
uc.rel.diff<- rel.diff.propn + z.conf*se.rel.diff

results<- array("", dim = c(14, 1))
results[1]<- round(Kappa, digits)
results[2]<- round(se.kappa, digits)
results[3]<- round(z.kappa, 2)
results[4]<- format(round(p.kappa, digits), scientific = F)
results[5]<- round(k.pos, digits)
results[6]<- round(k.neg, digits)
results[7]<- round(obsAgree, digits)
results[8]<- round(McNemar, digits)
results[9]<- round(p.val.McNemar, 2)
results[10]<- round(abs.diff.propn, digits)
results[11]<- round(rel.diff.propn, digits)
results[12]<- round(se.kappa.obs, digits)
results[13]<- round(lc.kappa, digits)
results[14]<- round(uc.kappa, digits)
# results[14]<- lc.abs.diff
# results[15]<- uc.abs.diff
# results[16]<- lc.rel.diff
# results[17]<- uc.rel.diff

dim(results)<- c(length(results), 1)
colnames(results)<- "Value"
rownames(results)<- c("Kappa", "SE for kappa = 0", "Z(kappa)", "p(kappa) - one-tailed", "Proportion positive agreement", 
                  "Proportion negative agreement",
                  "Overall proportion agreement", "McNemar's Chi sq", "p(Chi sq)", "Absolute diff. in proportions",
                  "Relative diff. in proportions", "SE for non-zero kappa", 
                  paste("Kappa lower ", conf*100, "% limit", sep = ""), 
                  paste("Kappa upper ", conf*100, "% limit", sep = ""))
#                 "Absolute diff. lower 95% limit", "Absolute diff. upper 95% limit",
#                 "Relative diff. lower 95% limit", "Relative diff. upper 95% limit")

#results<- round(results, digits)


######################################
# For infected populations only
######################################

# if the saqmple population are all infected!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Se.1<- t1.pos/n                           # sensitivity of test 1
Se.2<- t2.pos/n                           # sensitivity of test 1
p.2pos<- test.data[1]/n                   # proportion positive to both tests
Se.cov<- p.2pos - Se.1*Se.2               # covariance of sensitivities
Se.cov.min<- max(-1*(1 - Se.1)*(1 - Se.2), -1*Se.1*Se.2)
Se.cov.max<- min(Se.1*(1 - Se.2), Se.2*(1 - Se.1))
Se.cov.propn<- ifelse(Se.cov.max == 0, NA, Se.cov/Se.cov.max)

Se.results<- array(0, dim = c(7, 1))
Se.results[1]<- Se.1
Se.results[2]<- Se.2
Se.results[3]<- p.2pos
Se.results[4]<- Se.cov
Se.results[5]<- Se.cov.min
Se.results[6]<- Se.cov.max
Se.results[7]<- Se.cov.propn

rownames(Se.results)<- c(paste(test1, "sensitivity"), paste(test2, "sensitivity"), 
                       "Proportion ++", "Sensitivity covariance", "Min. sensitivity covariance",
                        "Max. sensitivity covariance", "Se covariance proportion")
Se.results<- round(Se.results, digits)


######################################
# For uninfected populations only
######################################

# if the sample population are all uninfected!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Sp.1<- t1.neg/n                           # specificity of test 1
Sp.2<- t2.neg/n                           # specificity of test 1
p.2neg<- test.data[4]/n                   # proportion negative to both tests
Sp.cov<- p.2neg - Sp.1*Sp.2               # covariance of sensitivities
Sp.cov.min<- max(-1*(1 - Sp.1)*(1 - Sp.2), -1*Sp.1*Sp.2)
Sp.cov.max<- min(Sp.1*(1 - Sp.2), Sp.2*(1 - Sp.1))
Sp.cov.propn<- ifelse(Sp.cov.max == 0, NA, Sp.cov/Sp.cov.max)

Sp.results<- array(0, dim = c(7, 1))
Sp.results[1]<- Sp.1
Sp.results[2]<- Sp.2
Sp.results[3]<- p.2neg
Sp.results[4]<- Sp.cov
Sp.results[5]<- Sp.cov.min
Sp.results[6]<- Sp.cov.max
Sp.results[7]<- Sp.cov.propn

rownames(Sp.results)<- c(paste(test1, "specificity"), paste(test2, "specificity"),
                       "Proportion --", "Specificity covariance", "Min. specificity covariance",
                        "Max. specificity covariance", "Sp covariance proportion")
Sp.results<- round(Sp.results, digits)


# write to html and file
# cat(test)
heading<- "Compare two tests"
subheadings<- "Test comparison statistics"
filename<- digest(Sys.time)
tmp.file<- paste(fpath, "tmp/", filename, sep = "")
result<- list(results)
if (Pop.status == 2) {
    result[[2]]<- Se.results
    subheadings[2]<- "Test sensitivities and covariance assuming an infected population"
} else if (Pop.status == 3) {
    result[[2]]<- Sp.results
    subheadings[2]<- "Test specificities and covariance assuming an uninfected population"
}
result.txt<- paste("The table below summarises relevant parameters for the comparison of two tests.
", 
      "Population status = ", ifelse(Pop.status == 1, "Unknown/Mixed", ifelse(Pop.status == 2, "All Infected", "All Uninfected")), "", sep = "")
output<- html.output(heading, subheadings, inputs, results = result, graphs = "", graph.headings = "", show.inputs = T, show.graphs = F, tmp.file, result.txt = result.txt) 
write.html(output, tmp.file)
cat(output)