Acknowledgements

Author list for this study:

Anastasia N. Tellis1, Sam M. Rowe2, Ronald Coilparampil3, Cheryl Jenkins3, Andrew Dart2, Ruth N. Zadoks2 and Katrina L. Bosward2, 4

Author affiliations: 1The University of Sydney, School of Life and Environmental Sciences, Faculty of Science, Camperdown, NSW, 2006 Australia 2The University of Sydney, Sydney School of Veterinary Science, Faculty of Science, Camden, NSW, 2570 Australia 3NSW Department of Primary Industries, Elizabeth Macarthur Agricultural Institute, Menangle, NSW, 2568 Australia 4Corresponding author

pre code {
 white-space: 
}
library(knitr)
library(readr)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.0.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
data <- read_csv("AlpacaDxData.csv")
## Parsed with column specification:
## cols(
##   id = col_double(),
##   state = col_character(),
##   CFT = col_double(),
##   IFA = col_double(),
##   ELISA = col_double(),
##   SP = col_double()
## )

Import data

Show first 10 rows of dataset

head(data, n=10)



Inspect Data

#summarytools::dfSummary(Baseline, style='grid')
print(summarytools::dfSummary(data,valid.col=FALSE, graph.magnif=0.8, style="grid"), method = "render")

Data Frame Summary

data

Dimensions: 131 x 6
Duplicates: 0
No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 id [numeric] Mean (sd) : 66 (38) min < med < max: 1 < 66 < 131 IQR (CV) : 65 (0.6) 131 distinct values 0 (0%)
2 state [character] 1. NSW 2. SA 3. VIC
41(31.3%)
18(13.7%)
72(55.0%)
 0 (0%)
3 CFT [numeric] Min : 0 Mean : 0.1 Max : 1
0:84(89.4%)
1:10(10.6%)
 37 (28.24%)
4 IFA [numeric] Min : 0 Mean : 0.1 Max : 1
0:123(93.9%)
1:8(6.1%)
 0 (0%)
5 ELISA [numeric] 1 distinct value
0:131(100.0%)
 0 (0%)
6 SP [numeric] Mean (sd) : 2.3 (3.4) min < med < max: -3.8 < 1.4 < 16.7 IQR (CV) : 2.5 (1.5) 64 distinct values 0 (0%)

Generated by summarytools 0.9.6 (R version 4.0.0)
2020-09-10


## Create data.subset This dataset removes the inconclusive results from the CFT (n=37) and will be used to make comparisons between the 3 tests at once.

data.subset <- data %>% filter(!is.na(CFT))

Create 2 x 2 tables

CFT vs IFA

library(janitor)
data.subset %>% tabyl(CFT,IFA)

CFT vs ELISA

data.subset %>% tabyl(CFT,ELISA)

IFA vs ELISA

data.subset %>% tabyl(IFA,ELISA)

Pairwise agreement between tests (Kappa) using data.subset

CFT vs IFA

library(psych)
library(broom)
CFT.IFA <- cohen.kappa(cbind(data.subset$CFT,data.subset$IFA)) %>% tidy %>% select(estimate,conf.low,conf.high) %>% slice(n()) 
cohen.kappa(cbind(data.subset$CFT,data.subset$IFA))
## Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)
## 
## Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
##                   lower estimate upper
## unweighted kappa -0.017     0.29   0.6
## weighted kappa   -0.017     0.29   0.6
## 
##  Number of subjects = 94

CFT vs ELISA

CFT.ELISA <- cohen.kappa(cbind(data.subset$CFT,data.subset$ELISA)) %>% tidy %>% select(estimate,conf.low,conf.high) %>% slice(n()) 
cohen.kappa(cbind(data.subset$CFT,data.subset$ELISA))
## Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)
## 
## Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
##                  lower estimate upper
## unweighted kappa     0        0     0
## weighted kappa       0        0     0
## 
##  Number of subjects = 94

CFT vs ELISA

IFA.ELISA <- cohen.kappa(cbind(data.subset$IFA,data.subset$ELISA)) %>% tidy %>% select(estimate,conf.low,conf.high) %>% slice(n()) 
cohen.kappa(cbind(data.subset$IFA,data.subset$ELISA))
## Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)
## 
## Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
##                     lower estimate   upper
## unweighted kappa -5.9e-08        0 5.9e-08
## weighted kappa   -5.9e-08        0 5.9e-08
## 
##  Number of subjects = 94

Table of agreement

Table1 <- rbind(CFT.IFA,CFT.ELISA,IFA.ELISA)
Table1$Comparison <- c("CFT vs IFA","CFT vs ELISA","IFA vs ELISA")
Table1$Estimate <- paste0(round(Table1$estimate,2)," (",round(Table1$conf.low,2),", ",round(Table1$conf.high,2),")")
Table1 <- Table1 %>% select(Comparison,Estimate)
Table1

Using Bayesian LCM for estimating test characteristics for tests

library(runjags)
library(rjags)
library(lattice)
library(mcmcplots)
library(epiR)

Prior distributions

Prevalence of infection

Very little is known about the seroprevalence of Q fever in camelids. A survey of 167 camels in Northern Iran by Prioz et al (2015) found that 29% tested positive with a commercial ELISA test kit.

Prevalence: Uniform ranging from 0 to 30%

curve(dunif(x,0.0001,0.3),col = "blue", xlab = "Prevalence")

Complement fixation test (CFT)

Muleme et al (2016) found SE (0.3) and SP (0.97) in goats in Victoria. Horigan et al (2011) found SE (0.36) and SP (0.98) in sheep, goats and cattle. Kittelberger et al (2009) found a SE of 0.68 when using known seropositive ruminants.

Sensitivity: Mode = 0.45, and a 97.5th percentile of 0.8 Specificity: Mode = 0.97, and a 2.5th percentile of 0.95

library(ggplot2)
#Se
epi.betabuster(mode=0.45,conf=0.95,greaterthan=F,x=0.8)[c(1,2)]
## $shape1
## [1] 2.521
## 
## $shape2
## [1] 2.859
curve(dbeta(x,2.52,2.86),col = "blue", xlab = "Sensitivity")

#Sp
epi.betabuster(mode=0.97,conf=0.95,greaterthan=T,x=0.95)[c(1,2)]
## $shape1
## [1] 100
## 
## $shape2
## [1] 4.061856
curve(dbeta(x,100,4.06),col = "blue", xlab = "Specificity")

N = 100000

CFTprior <- data.frame(y = c(rbeta(N,2.52,2.86),rbeta(N,100,4.06)),
           Distribution = rep(c("Sensitivity","Specificity"),each=N))

CFT <- ggplot(data = CFTprior, aes(x = y, group = Distribution)) + geom_density(col = NA, alpha = 0.5,aes(fill = Distribution, color = Distribution)) + theme_classic()+ theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank(),axis.text.x=element_blank(),axis.title.x=element_blank(),legend.position = "none",text=element_text(family="Times")) + ggtitle("Complement fixation test") + xlim(0,1)

CFT

IFA

Muleme et al (2016) found SE (0.95) and SP (0.925) in goats in Victoria. Wood et al (2019) found SE (0.74) and SP (0.98) in cattle.

Sensitivity: Mode = 0.87, and a 2.5th percentile of 0.50. Specificity: Mode = 0.95, and a 2.5th percentile of 0.90. I chose this because IFA are generally highly specific.

#Se
epi.betabuster(mode=0.87,conf=0.95,greaterthan=T,x=0.5)[c(1,2)]
## $shape1
## [1] 5.866
## 
## $shape2
## [1] 1.727103
curve(dbeta(x,5.87,1.73),col = "blue", xlab = "Sensitivity")

#Sp
epi.betabuster(mode=0.95,conf=0.95,greaterthan=T,x=0.9)[c(1,2)]
## $shape1
## [1] 99.698
## 
## $shape2
## [1] 6.194632
curve(dbeta(x,99.7,6.2),col = "blue", xlab = "Specificity")

IFAprior <- data.frame(y = c(rbeta(N,5.87,1.73),rbeta(N,99.7,6.2)),
           Distribution = rep(c("Sensitivity","Specificity"),each=N))

IFA <- ggplot(data = IFAprior, aes(x = y, group = Distribution)) + geom_density(col = NA, alpha = 0.5,aes(fill = Distribution, color = Distribution)) + theme_classic()+ theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank(),axis.text.x=element_blank(),axis.title.x=element_blank(),legend.position = "none",text=element_text(family="Times")) + ggtitle("Immunofluorescence assay")+ xlim(0,1)
IFA

ELISA

Studies have evaluate various ELISAs for Q Fever with SE estimates ranging from 0.55 to 0.88. Specificity has ranged from 0.96 to 0.98 (Horrigan et al, 2011, Muleme et al, 2016 and Wood et al, 2019). Furthermore, our experience using this assay indicates that SE may be very low, as very few cases test positive

Sensitivitiy: Mode = 0.7, 2.5th percentile = 0.2 Specificity: Mode = 0.97, 2.5th percentile = 0.90.

#Se
epi.betabuster(mode=0.7,conf=0.95,greaterthan=T,x=0.2)[c(1,2)]
## $shape1
## [1] 2.256
## 
## $shape2
## [1] 1.538286
curve(dbeta(x,2.26,1.54),col = "blue", xlab = "Sensitivity")

#curve(dunif(x,0.0001,0.88),col = "blue", xlab = "Sensitivity")

#Sp
epi.betabuster(mode=0.97,conf=0.95,greaterthan=T,x=0.9)[c(1,2)]
## $shape1
## [1] 53.581
## 
## $shape2
## [1] 2.626216
curve(dbeta(x,53.6,2.6),col = "blue", xlab = "Specificity")

ELISAprior <- data.frame(y = c(rbeta(N,2.26,1.54),rbeta(N,53.6,2.6)),
           Distribution = rep(c("Sensitivity","Specificity"),each=N))

ELISA <- ggplot(data = ELISAprior, aes(x = y, group = Distribution)) + geom_density(col = NA, alpha = 0.5,aes(fill = Distribution, color = Distribution)) + theme_classic()+ theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank(),axis.title.x=element_blank(),legend.position = "bottom",text=element_text(family="Times")) + ggtitle("ELISA") + xlim(0,1)

ELISA

Summary of priors

png("plot.png", width=500*10, height=720*10, res=100*15)
cowplot::plot_grid(CFT,IFA,ELISA, ncol=1, axis="lr", align="v",rel_heights = c(1,1,1.6))
dev.off()
## quartz_off_screen 
##                 2
cowplot::plot_grid(CFT,IFA,ELISA, ncol=1, axis="lr", align="v",rel_heights = c(1,1,1.6))

References:

Horigan, Mark W., et al. “Q fever diagnosis in domestic ruminants: comparison between complement fixation and commercial enzyme-linked immunosorbent assays.” Journal of veterinary diagnostic investigation 23.5 (2011): 924-931.

Kittelberger, R., et al. “Comparison of the Q-fever complement fixation test and two commercial enzyme-linked immunosorbent assays for the detection of serum antibodies against Coxiella burnetii (Q-fever) in ruminants: Recommendations for use of serological tests on imported animals in New Zealand.” New Zealand veterinary journal 57.5 (2009): 262-268.

Muleme, Michael, et al. “Bayesian validation of the indirect immunofluorescence assay and its superiority to the enzyme-linked immunosorbent assay and the complement fixation test for detecting antibodies against Coxiella burnetii in goat serum.” Clin. Vaccine Immunol. 23.6 (2016): 507-514.

Pirouz, Hossein Janati, et al. “Seroepidemiology of Q fever in one-humped camel population in northeast Iran.” Tropical animal health and production 47.7 (2015): 1293-1298.

Wood, Caitlin, et al. “Validation of an indirect immunofluorescence assay (IFA) for the detection of IgG antibodies against Coxiella burnetii in bovine serum.” Preventive veterinary medicine 169 (2019): 104698.

Run BLC model



Model 1: One population, 2 tests

library(runjags)
library(rjags)
library(lattice)
library(mcmcplots)
library(epiR)
library(dplyr)

#Model 2: One population, 3 tests -----

#Test 1 (T1) = Immunofluorescence Assay
#Test 2 (T2) = Commercial ELISA test

#Priors
Prior <- "
#Prevalence of IMI
pi ~ dunif(0,.3)

#Reference test
SeT1 ~ dbeta(5.87,1.73)
SpT1 ~ dbeta(99.7,6.2)
SeT2 ~ dbeta(2.26,1.54)
SpT2 ~ dbeta(53.6,2.62)"

#Define data function 
# ++ / +- / -+ / --
table(data$IFA,data$ELISA)
##    
##       0
##   0 123
##   1   8
pop <- matrix(c(123,8,0,0),1,4)
ns <- sum(pop)
dd<-list(n = ns, pop = pop)
modelData <- dump.format(dd)


#Define model function
mmString <- function(Priors){
  paste0(
    "model{
pop[1,1:4] ~ dmulti(cell[1:4],n)

cell[1] <- pi*((1-SeT1)*(1-SeT2)) + (1-pi)*(SpT1*SpT2) # p[0,0] n = 123
cell[2] <- pi*(SeT1*(1-SeT2)) + (1-pi)*((1-SpT1)*SpT2) # p[1,0] n = 8
cell[3] <- pi*((1-SeT1)*SeT2) + (1-pi)*(SpT1*(1-SpT2)) # p[0,1] n = 0
cell[4] <- pi*(SeT1*SeT2) + (1-pi)*((1-SpT1)*(1-SpT2)) # p[1,1] n = 0

",Priors,"

}")
}

#Create final model string
modelString <- mmString(Prior)


#Run models ----

#Run
print(system.time(
  x <- run.jags( model = modelString,
                 monitor = c("SeT1","SeT2","SpT1","SpT2","pi"),
                 data = modelData,
                 burnin = 5000,
                 sample = 15000,
                 thin = 5, 
                 n.chains=3, 
                 inits = NA
  )
))
## Compiling rjags model...
## Calling the simulation using the rjags method...
## Adapting the model for 1000 iterations...
## Burning in the model for 5000 iterations...
## Running the model for 75000 iterations...
## Simulation complete
## Calculating summary statistics...
## Calculating the Gelman-Rubin statistic for 5 variables....
## Finished running the simulation
##    user  system elapsed 
##   2.123   0.222   2.394
#Reports
out.coda <- list(   result = as.mcmc.list(x[[1]]),
                  model = modelString,
                  data = modelData)

out.coda$result %>% summary
## 
## Iterations = 6001:80996
## Thinning interval = 5 
## Number of chains = 3 
## Sample size per chain = 15000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean       SD  Naive SE Time-series SE
## SeT1 0.76901 0.144340 6.804e-04      6.951e-04
## SeT2 0.48037 0.239138 1.127e-03      1.201e-03
## SpT1 0.94400 0.016009 7.547e-05      7.731e-05
## SpT2 0.98589 0.008663 4.084e-05      4.229e-05
## pi   0.01802 0.019092 9.000e-05      1.174e-04
## 
## 2. Quantiles for each variable:
## 
##           2.5%      25%     50%     75%   97.5%
## SeT1 0.4340009 0.680003 0.79224 0.88281 0.97337
## SeT2 0.0802081 0.285611 0.46894 0.66994 0.92480
## SpT1 0.9094048 0.933901 0.94511 0.95522 0.97234
## SpT2 0.9645278 0.981354 0.98757 0.99230 0.99755
## pi   0.0004115 0.004679 0.01185 0.02477 0.06994
xyplot(out.coda$result)

densityplot(out.coda$result)

acfplot(out.coda$result)

mcmcplot(out.coda$result,style="plain")
## 
                                                                                
Preparing plots for SeT1.  20% complete.
## 
                                                                                
Preparing plots for SeT2.  40% complete.
## 
                                                                                
Preparing plots for SpT1.  60% complete.
## 
                                                                                
Preparing plots for SpT2.  80% complete.
## 
                                                                                
Preparing plots for pi.  100% complete.

Model 2: One population, 3 tests

#Model 1: One population, 3 tests -----

#Test 1 (T1) = Complement fixation test
#Test 2 (T2) = Immunofluorescence Assay
#Test 3 (T3) = Commercial ELISA test

Prior <- "
#Prevalence of IMI
pi ~ dunif(0,.3)

#Reference test
SeT1 ~ dbeta(2.52,2.86)
SpT1 ~ dbeta(100,4.06)
SeT2 ~ dbeta(5.87,1.73)
SpT2 ~ dbeta(99.7,6.2)
SeT3 ~ dbeta(2.26,1.54)
SpT3 ~ dbeta(53.6,2.62)"


#Define model function
mmString <- function(Priors){
  paste0(
    "model{
pop[1, 1:8] ~ dmulti(p[1:8], n)
p[1] <- pi * (1-SeT1) * (1-SeT2) * (1-SeT3) + (1-pi) * (SpT1) * (SpT2) * (SpT3)    #p[0,0,0] n = 80
p[2] <- pi * (SeT1) * (1-SeT2) * (1-SeT3) + (1-pi) * (1-SpT1) * (SpT2) * (SpT3)    #p[1,0,0] n =  7
p[3] <- pi * (1-SeT1) * (SeT2) * (1-SeT3) + (1-pi) * (SpT1) * (1-SpT2) * (SpT3)    #p[0,1,0] n =  4
p[4] <- pi * (SeT1) * (SeT2) * (1-SeT3) + (1-pi) * (1-SpT1) * (1-SpT2) * (SpT3)    #p[1,1,0] n =  3
p[5] <- pi * (1-SeT1) * (1-SeT2) * (SeT3) + (1-pi) * (SpT1) * (SpT2) * (1-SpT3)    #p[0,0,1] n =  0
p[6] <- pi * (SeT1) * (1-SeT2) * (SeT3) + (1-pi) * (1-SpT1) * (SpT2) * (1-SpT3)    #p[1,0,1] n =  0
p[7] <- pi * (1-SeT1) * (SeT2) * (SeT3) + (1-pi) * (SpT1) * (1-SpT2) * (1-SpT3)    #p[0,1,1] n =  0
p[8] <- pi * (SeT1) * (SeT2) * (SeT3) + (1-pi) * (1-SpT1) * (1-SpT2) * (1-SpT3)    #p[1,1,1] n =  0

",Priors,"

}")
}


#Data
#matrix(as.numeric(table(data.subset$CFT,data.subset$IFA,data.subset$ELISA)),1,8)
pop <- matrix(c(80,7,4,3,0,0,0,0),1,8)
ns <- sum(pop)
dd<-list(n = ns, pop = pop)
modelData <- dump.format(dd)


#Create final model string
modelString <- mmString(Prior)


#Run models ----

#Algorithm outcome ----

#Run
print(system.time(
  x <- run.jags( model = modelString,
                 monitor = c("SeT1","SeT2","SeT3","SpT1","SpT2","SpT3","pi"),
                 data = modelData,
                 burnin = 5000,
                 sample = 15000,
                 thin = 5, 
                 n.chains=3, 
                 inits = NA
  )
))
## Compiling rjags model...
## Calling the simulation using the rjags method...
## Adapting the model for 1000 iterations...
## Burning in the model for 5000 iterations...
## Running the model for 75000 iterations...
## Simulation complete
## Calculating summary statistics...
## Calculating the Gelman-Rubin statistic for 7 variables....
## Finished running the simulation
##    user  system elapsed 
##   4.730   0.449   5.340
#Reports
out.coda <- list(   result = as.mcmc.list(x[[1]]),
                  model = modelString,
                  data = modelData)

out.coda$result %>% summary
## 
## Iterations = 6001:80996
## Thinning interval = 5 
## Number of chains = 3 
## Sample size per chain = 15000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean      SD  Naive SE Time-series SE
## SeT1 0.57467 0.17322 8.165e-04      8.209e-04
## SeT2 0.75748 0.14750 6.953e-04      7.330e-04
## SeT3 0.29482 0.18400 8.674e-04      9.191e-04
## SpT1 0.94453 0.01847 8.708e-05      8.969e-05
## SpT2 0.94950 0.01718 8.100e-05      8.265e-05
## SpT3 0.98196 0.01092 5.147e-05      5.311e-05
## pi   0.05988 0.03860 1.820e-04      2.027e-04
## 
## 2. Quantiles for each variable:
## 
##          2.5%    25%     50%     75%  97.5%
## SeT1 0.222930 0.4543 0.58221 0.70276 0.8816
## SeT2 0.425262 0.6636 0.78015 0.87394 0.9724
## SeT3 0.043138 0.1547 0.25840 0.39965 0.7487
## SpT1 0.904155 0.9330 0.94595 0.95765 0.9767
## SpT2 0.910962 0.9391 0.95127 0.96179 0.9781
## SpT3 0.955212 0.9761 0.98406 0.99008 0.9968
## pi   0.005549 0.0317 0.05312 0.08053 0.1534
xyplot(out.coda$result)

densityplot(out.coda$result)

acfplot(out.coda$result)

mcmcplot(out.coda$result,style="plain")
## 
                                                                                
Preparing plots for SeT1.  14% complete.
## 
                                                                                
Preparing plots for SeT2.  29% complete.
## 
                                                                                
Preparing plots for SeT3.  43% complete.
## 
                                                                                
Preparing plots for SpT1.  57% complete.
## 
                                                                                
Preparing plots for SpT2.  71% complete.
## 
                                                                                
Preparing plots for SpT3.  86% complete.
## 
                                                                                
Preparing plots for pi.  100% complete.

NPV and PPV

#Create NPV plot -----
library(ggplot2)
valplotfunc <- function(SE,SP,Testname){
d <- matrix(c(seq(0,1,0.005))) %>% as.data.frame
d$P <- d$V1
d$V1 <- NULL
d$Test <- Testname
d$SE <- SE
d$SP <- SP

d$NPV <- (d$SP*(1-d$P)) / ((1-d$SE)*d$P + d$SP*(1-d$P))
d$PPV <- d$SE*d$P / (d$SE*d$P + (1 - d$SP) * (1-d$P))
d$FP <- d$P*(1-d$SE)
d}

#All IMI
a <- valplotfunc(0.58,0.95,"Complement fixation test")
b <- valplotfunc(0.78,0.95,"Immunofluorescence Assay")
c <- valplotfunc(0.26,0.98,"ELISA")

e <- rbind(a,b,c)

#windowsFonts(Times=windowsFont("Times New Roman"))

#NPV without legend
NPV <- ggplot(e) + coord_cartesian(ylim = (c(0,1))) + aes(x=P, y=NPV, group=Test, colour=Test) +
  labs(y = "Negative predictive value", x = "Seroprevalence") + geom_line(aes(colour=Test,linetype=Test),size=1.1) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1),axis.text=element_text(size=12,family="Times"),axis.title=element_text(size=14,face="bold",family="Times"),panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "grey"),legend.position="None",legend.text = element_text(colour="black", size=12,face="bold",family="Times")) + scale_color_brewer(palette="Set1")
NPV

PPV <- ggplot(e) + coord_cartesian(ylim = (c(0,1))) + aes(x=P, y=PPV, group=Test, colour=Test) +
  labs(y = "Positive predictive value", x = "Seroprevalence") + geom_line(aes(colour=Test,linetype=Test),size=1.1) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1),axis.text=element_text(size=12,family="Times"),axis.title=element_text(size=14,face="bold",family="Times"),panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "grey"),legend.position="noneee",legend.text = element_text(colour="black", size=12,face="bold",family="Times")) + scale_color_brewer(palette="Set1") 
PPV

NPVleg <- ggplot(e) + coord_cartesian(ylim = (c(0,1))) + aes(x=P, y=NPV, group=Test, colour=Test) +
  labs(y = "Negative predictive value", x = "Seroprevalence") + geom_line(aes(colour=Test,linetype=Test),size=1.1) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1),axis.text=element_text(size=12,family="Times"),axis.title=element_text(size=14,face="bold",family="Times"),panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "grey"),legend.position="bottom",legend.text = element_text(colour="black", size=12,face="bold",family="Times")) + scale_color_brewer(palette="Set1")
NPVleg

#Print without legend vertical
tiff("Vertical.tiff", units="in", width=3.5, height=7, res=300)
cowplot::plot_grid(PPV,NPV,nrow=2)
dev.off()
## quartz_off_screen 
##                 2
#Print legend
tiff("legend.tiff", units="in", width=9, height=9, res=300)
cowplot::plot_grid(NPVleg,nrow=2)
dev.off()
## quartz_off_screen 
##                 2