Required packages
library(epiR)
library(readxl)
library(dplyr)
library(ggplot2)
library(tidyr)
library(ggbreak)
library(patchwork)
# Analytical sensitivity
analytical_val <- as.table(matrix(c(33,0,1,4), nrow = 2, byrow = TRUE))
rval <- epi.tests(analytical_val, conf.level = 0.95)
print(rval)
## Outcome + Outcome - Total
## Test + 33 0 33
## Test - 1 4 5
## Total 34 4 38
##
## Point estimates and 95% CIs:
## --------------------------------------------------------------
## Apparent prevalence * 0.87 (0.72, 0.96)
## True prevalence * 0.89 (0.75, 0.97)
## Sensitivity * 0.97 (0.85, 1.00)
## Specificity * 1.00 (0.40, 1.00)
## Positive predictive value * 1.00 (0.89, 1.00)
## Negative predictive value * 0.80 (0.28, 0.99)
## Positive likelihood ratio Inf (NaN, Inf)
## Negative likelihood ratio 0.03 (0.00, 0.20)
## False T+ proportion for true D- * 0.00 (0.00, 0.60)
## False T- proportion for true D+ * 0.03 (0.00, 0.15)
## False T+ proportion for T+ * 0.00 (0.00, 0.11)
## False T- proportion for T- * 0.20 (0.01, 0.72)
## Correctly classified proportion * 0.97 (0.86, 1.00)
## --------------------------------------------------------------
## * Exact CIs
Loading dataset from clinical validation
data <- read_excel("SupFile10. qPCR and LAMP TcHSP70 results for clinical samples copy.xlsx")
## New names:
## • `` -> `...13`
## • `` -> `...14`
Number of total tests
## [1] 109
Number of different patients
data %>%
distinct(Sample_ID) %>%
nrow()
## [1] 100
Print samples with more than one sample tested
dup <- data[duplicated(data$Sample_ID), ]
table(dup$Sample_ID)
##
## 10032 114 18 8664 9297 9340 9362 9399
## 1 1 2 1 1 1 1 1
Summary by test
table(pivot_data$Method, pivot_data$Test_result)
##
## NEGATIVE POSITIVE
## LAMP 92 17
## qPCR 87 22
Note: true positive=17, true negative 87, false positive= 0, false
negative= 5
Statistics by test
by_test <- as.table(matrix(c(17,0,5,87), nrow = 2, byrow = TRUE))
rval <- epi.tests(by_test, conf.level = 0.95)
print(rval)
## Outcome + Outcome - Total
## Test + 17 0 17
## Test - 5 87 92
## Total 22 87 109
##
## Point estimates and 95% CIs:
## --------------------------------------------------------------
## Apparent prevalence * 0.16 (0.09, 0.24)
## True prevalence * 0.20 (0.13, 0.29)
## Sensitivity * 0.77 (0.55, 0.92)
## Specificity * 1.00 (0.96, 1.00)
## Positive predictive value * 1.00 (0.80, 1.00)
## Negative predictive value * 0.95 (0.88, 0.98)
## Positive likelihood ratio Inf (NaN, Inf)
## Negative likelihood ratio 0.23 (0.11, 0.49)
## False T+ proportion for true D- * 0.00 (0.00, 0.04)
## False T- proportion for true D+ * 0.23 (0.08, 0.45)
## False T+ proportion for T+ * 0.00 (0.00, 0.20)
## False T- proportion for T- * 0.05 (0.02, 0.12)
## Correctly classified proportion * 0.95 (0.90, 0.98)
## --------------------------------------------------------------
## * Exact CIs
kappa_by_test<-epi.kappa(by_test, method = "cohen", alternative = c("greater"), conf.level = 0.95)
print(kappa_by_test)
## $prop.agree
## obs exp
## 1 0.9541284 0.7051595
##
## $pindex
## est se lower upper
## 1 -0.6422018 0.0518233 -0.7437736 -0.54063
##
## $bindex
## est se lower upper
## 1 -0.04587156 0.0518233 -0.1474434 0.05570024
##
## $pabak
## est lower upper
## 1 0.9082569 0.7923862 0.9698802
##
## $kappa
## est se lower upper
## 1 0.8444191 0.06796335 0.7112133 0.9776248
##
## $z
## test.statistic p.value
## 1 12.42462 9.606401e-36
##
## $mcnemar
## test.statistic df p.value
## 1 5 1 0.02534732
Summary per patient
Positive patients by LAMP in at least one sample
LAMP <- pivot_data %>%
filter(Method == "LAMP", Test_result == "POSITIVE") %>%
select(Sample_ID) %>%
distinct() %>% # Use distinct instead of unique for better readability
pull(Sample_ID)
LAMP
## [1] "9516" "9519" "9435" "9557" "366" "9778" "10032" "9364" "10084"
## [10] "114" "9399" "10049" "10070" "18"
Positive patients by qPCR in at least one sample
qPCR <- pivot_data %>%
filter(Method == "qPCR", Test_result == "POSITIVE") %>%
select(Sample_ID) %>%
distinct() %>%
pull(Sample_ID)
qPCR
## [1] "9297" "9516" "9519" "9435" "9557" "366" "9778" "10032" "9364"
## [10] "10084" "114" "9399" "10049" "10070" "18" "439"
Patients that were pisitive by qPCR but not by LAMP
difference<-setdiff(qPCR, LAMP)
difference
## [1] "9297" "439"
Statistics by patient
by_patient <- as.table(matrix(c(14,0,2, 84), nrow = 2, byrow = TRUE))
rval <- epi.tests(by_patient, conf.level = 0.95)
print(rval)
## Outcome + Outcome - Total
## Test + 14 0 14
## Test - 2 84 86
## Total 16 84 100
##
## Point estimates and 95% CIs:
## --------------------------------------------------------------
## Apparent prevalence * 0.14 (0.08, 0.22)
## True prevalence * 0.16 (0.09, 0.25)
## Sensitivity * 0.88 (0.62, 0.98)
## Specificity * 1.00 (0.96, 1.00)
## Positive predictive value * 1.00 (0.77, 1.00)
## Negative predictive value * 0.98 (0.92, 1.00)
## Positive likelihood ratio Inf (NaN, Inf)
## Negative likelihood ratio 0.12 (0.03, 0.46)
## False T+ proportion for true D- * 0.00 (0.00, 0.04)
## False T- proportion for true D+ * 0.12 (0.02, 0.38)
## False T+ proportion for T+ * 0.00 (0.00, 0.23)
## False T- proportion for T- * 0.02 (0.00, 0.08)
## Correctly classified proportion * 0.98 (0.93, 1.00)
## --------------------------------------------------------------
## * Exact CIs
kappa_by_test<-epi.kappa(by_patient, method = "cohen", alternative = c("greater"), conf.level = 0.95)
print(kappa_by_test)
## $prop.agree
## obs exp
## 1 0.98 0.7448
##
## $pindex
## est se lower upper
## 1 -0.7 0.05047772 -0.7989345 -0.6010655
##
## $bindex
## est se lower upper
## 1 -0.02 0.05047772 -0.1189345 0.07893451
##
## $pabak
## est lower upper
## 1 0.96 0.8592321 0.9951373
##
## $kappa
## est se lower upper
## 1 0.9216301 0.05485893 0.8141086 1.029152
##
## $z
## test.statistic p.value
## 1 16.8 1.22022e-63
##
## $mcnemar
## test.statistic df p.value
## 1 2 1 0.1572992
Correlation between LAMP and parasitic load estimated based on a
qPCR callibration curve
data_pos<-data %>% filter(qPCR=="POSITIVE")
p<-ggplot(data_pos, aes(x= qPCR, y=Parasitic_load, color=LAMP)) +
geom_jitter(width = 0.3, height = 1)+
scale_color_manual(values = c("#c87137ff", "#53c675ff")) +
scale_y_break(c(500000, 950000), scales = 0.5)+
scale_y_break(c(65000, 300000), scales = 0.5)
p
ggsave(p, file = "PCR_LAMP.svg", width = 3, height = 2.5, units = "in")