Recap Hui-Walter

• Model assumptions
  • The population is divided into two or more populations with different prevalences in which two or more tests are evaluated
  • Se and Sp are the same in all populations
  • The tests are conditionally independent given the disease status.

Criticism of the Hui-Walter paradigm assumption (1)

• The assumption (1) of distinct prevalences is necessary for the Hui-Walter model because otherwise, the data can be collapsed into a single 2x2 table with only three degrees of freedom for estimation.

• The smaller the difference between disease prevalences, the larger are the posterior credible intervals, indicating a loss in precision.

• The smallest difference in prevalence assessed by simulation was 10%. In case of rare diseases, it might be difficult to find populations with prevalences higher than 10%.

Criticism of the Hui-Walter paradigm assumption (2)

• In reality, the assumption of constant Se and Sp is questionable as already illustrated in the figure. Here the same diagnostic tests have been applied in the same laboratory to experimentally and naturally infected cats.

Criticism of the Hui-Walter paradigm assumption (2)

• If assumption (2) is not satisfied, the accuracies would differ between two populations, this would add four additional parameters to be estimated, while there are only three additional degrees of freedom.

• Se and Sp are assumed to vary with external factors.

• Se, for example, especially when detecting an infectious agent, may depend on the prevalence and the stage of disease.

• The occurrence of a so-called spectrum bias contradicts this assumption.

Criticism of the Hui-Walter paradigm assumption (3)

• Assumption (3) demanding conditional independence was the first to be questioned by Vacek (1985).

• Not accounting for potential conditional dependence may lead to misleading, biased estimates with a positive correlation leading to an over-estimation and a negative correlation to an under-estimation of the test Se and Sp.

Condtional independencies

• Test are considered conditionally independent if the probability of getting a given test result on one test does not depend on the result from the other test, given the disease status of the individual.

• Conditional independence implies that given that an animal is diseased (or not) the probability \(P\) of positive (or negative) outcomes for T, the test results of the first test, is the same - regardless of the known outcome for the second test, T.

Conditional independence

\[P(T\textsubscript{1}\textsuperscript{+},T\textsubscript{2}\textsuperscript{+} \vert D\textsuperscript{+}) = P(T\textsubscript{1}\textsuperscript{+} \vert D\textsuperscript{+}) \times P(T\textsubscript{2}\textsuperscript{+} \vert D\textsuperscript{+})]\]

Conditional dependence

\[P(T\textsubscript{1}\textsuperscript{+},T\textsubscript{2}\textsuperscript{+} \vert D\textsuperscript{+}) \neq P(T\textsubscript{1}\textsuperscript{+} \vert D\textsuperscript{+}) \times P(T\textsubscript{2}\textsuperscript{+} \vert D\textsuperscript{+})\]

Conditional dependencies

• Conditional dependence, in contrast, implies that \[P(T\textsubscript{1}\textsuperscript{+} \vert T\textsubscript{2}\textsuperscript{+}) \neq P(T\textsubscript{1}\textsuperscript{+} \vert T\textsubscript{2}\textsuperscript{--})\] and / or \[P(T\textsubscript{1}\textsuperscript{--} \vert T\textsubscript{2}\textsuperscript{--}) \neq P(T\textsubscript{1}\textsuperscript{--} \vert T\textsubscript{2}\textsuperscript{+})\]

Conditional (in)dependencies Interpretation

• Seen from a biological perspective, conditional dependency between two diagnostic tests could occur if both tests are based on the same biological principle.

• For example, the Sps of two ELISAs might be conditionally dependent because they are both affected by the same cross-reacting agent. Another example would be two PCRs utilising fecal material which might contain substances potentially inhibiting the PCR reaction.

Conditional dependencies

• Obviously, conditional dependencies or covariances are additional parameters to be estimated, which in the frequentist situation (without any constraints put on the parameters) would lead to a non-identifiable problem (over-parameterisation).
• Whereas under the assumption of conditional independence at least three tests per sample allowing to estimate seven parameters are needed, under the assumption of conditional dependence 15 parameters need to estimated thus leading to non-identifiability.
• It is of course vital, that the parameters of a latent class model are identifiable to obtain meaningful estimates.

Conditional dependencies

Conditional dependencies

• For example for two diagnostic tests named T and T the probabilities of the four different options of binary test results (+ +, + –, – +, – –) including also two conditional dependencies
  • \(covs12\), the covariance between the sensitivities of test 1 and 2
  • \(covc12\), the covariance between specificities of test 1 and 2) could be modelled as follows:

Example of a COVID-19 data set

Dealing with correlation

It helps to consider the data simulation as a (simplified) biological process (where my parameters are not representative of real life!).

2 pops - 3 tests

# The probability of infection with COVID in two populations:
prevalence <- c(0.01,0.05)
# The probability of shedding COVID in the nose conditional on infection:
nose_shedding <- 0.8
# The probability of shedding COVID in the throat conditional on infection:
throat_shedding <- 0.8
# The probability of detecting virus with the antigen test:
antigen_detection <- 0.75
# The probability of detecting virus with the PCR test:
pcr_detection <- 0.999
# The probability of random cross-reaction with the antigen test:
antigen_crossreact <- 0.05
# The probability of random cross-reaction with the PCR test:
pcr_crossreact <- 0.01

Simulating latent states:

N <- 20000
Populations <- length(prevalence)

covid_data <- tibble(Population = sample(seq_len(Populations), N, replace=TRUE)) %>%
  ## True infection status:
  mutate(Status = rbinom(N, 1, prevalence[Population])) %>%
  ## Nose shedding status:
  mutate(Nose = Status * rbinom(N, 1, nose_shedding)) %>%
  ## Throat shedding status:
  mutate(Throat = Status * rbinom(N, 1, throat_shedding))

Simulating test results:

covid_data <- covid_data %>%
  ## The nose swab antigen test may be false or true positive:
  mutate(NoseAG = case_when(
    Nose == 1 ~ rbinom(N, 1, antigen_detection),
    Nose == 0 ~ rbinom(N, 1, antigen_crossreact)
  )) %>%
  ## The throat swab antigen test may be false or true positive:
  mutate(ThroatAG = case_when(
    Throat == 1 ~ rbinom(N, 1, antigen_detection),
    Throat == 0 ~ rbinom(N, 1, antigen_crossreact)
  )) %>%
  ## The PCR test may be false or true positive:
  mutate(ThroatPCR = case_when(
    Throat == 1 ~ rbinom(N, 1, pcr_detection),
    Throat == 0 ~ rbinom(N, 1, pcr_crossreact)
  ))

The overall sensitivity of the tests can be calculated as follows:

covid_sensitivity <- c(
  # Nose antigen:
  nose_shedding*antigen_detection + (1-nose_shedding)*antigen_crossreact,
  # Throat antigen:
  throat_shedding*antigen_detection + (1-throat_shedding)*antigen_crossreact,
  # Throat PCR:
  throat_shedding*pcr_detection + (1-throat_shedding)*pcr_crossreact
)
covid_sensitivity
## [1] 0.6100 0.6100 0.8012

The overall specificity of the tests is more straightforward:

covid_specificity <- c(
  # Nose antigen:
  1 - antigen_crossreact,
  # Throat antigen:
  1 - antigen_crossreact,
  # Throat PCR:
  1 - pcr_crossreact
)
covid_specificity
## [1] 0.95 0.95 0.99

Model specification

prob[1,p] <-  prev[p] * ((1-se[1])*(1-se[2])*(1-se[3]) 
                         +covse12 +covse13 +covse23) +
              (1-prev[p]) * (sp[1]*sp[2]*sp[3] 
                             +covsp12 +covsp13 +covsp23)

prob[2,p] <- prev[p] * (se[1]*(1-se[2])*(1-se[3]) 
                           -covse12 -covse13 +covse23) +
               (1-prev[p]) * ((1-sp[1])*sp[2]*sp[3] 
                              -covsp12 -covsp13 +covsp23)

## snip ##
        
# Covariance in sensitivity between tests 1 and 2:
covse12 ~ dunif( (se[1]-1)*(1-se[2]) , 
                     min(se[1],se[2]) - se[1]*se[2] )
# Covariance in specificity between tests 1 and 2:
covsp12 ~ dunif( (sp[1]-1)*(1-sp[2]) , 
                     min(sp[1],sp[2]) - sp[1]*sp[2] )

Template Hui-Walter

The model code and data format for an arbitrary number of populations (and tests) can be determined automatically using the template_huiwalter function from the runjas package:

template_huiwalter(
  covid_data %>% select(Population, NoseAG, ThroatAG, ThroatPCR), 
  outfile = 'covidmodel.txt', covariance=TRUE)

This generates self-contained model/data/initial values etc

model{

    ## Observation layer:

    # Complete observations (N=20000):
    for(p in 1:Populations){
        Tally_RRR[1:8,p] ~ dmulti(prob_RRR[1:8,p], N_RRR[p])

        prob_RRR[1:8,p] <- se_prob[1:8,p] + sp_prob[1:8,p]
    }


    ## Observation probabilities:

    for(p in 1:Populations){

        # Probability of observing NoseAG- ThroatAG- ThroatPCR- from a true positive::
        se_prob[1,p] <- prev[p] * ((1-se[1])*(1-se[2])*(1-se[3]) +covse12 +covse13 +covse23)
        # Probability of observing NoseAG- ThroatAG- ThroatPCR- from a true negative::
        sp_prob[1,p] <- (1-prev[p]) * (sp[1]*sp[2]*sp[3] +covsp12 +covsp13 +covsp23)

        # Probability of observing NoseAG+ ThroatAG- ThroatPCR- from a true positive::
        se_prob[2,p] <- prev[p] * (se[1]*(1-se[2])*(1-se[3]) -covse12 -covse13 +covse23)
        # Probability of observing NoseAG+ ThroatAG- ThroatPCR- from a true negative::
        sp_prob[2,p] <- (1-prev[p]) * ((1-sp[1])*sp[2]*sp[3] -covsp12 -covsp13 +covsp23)

        # Probability of observing NoseAG- ThroatAG+ ThroatPCR- from a true positive::
        se_prob[3,p] <- prev[p] * ((1-se[1])*se[2]*(1-se[3]) -covse12 +covse13 -covse23)
        # Probability of observing NoseAG- ThroatAG+ ThroatPCR- from a true negative::
        sp_prob[3,p] <- (1-prev[p]) * (sp[1]*(1-sp[2])*sp[3] -covsp12 +covsp13 -covsp23)

        # Probability of observing NoseAG+ ThroatAG+ ThroatPCR- from a true positive::
        se_prob[4,p] <- prev[p] * (se[1]*se[2]*(1-se[3]) +covse12 -covse13 -covse23)
        # Probability of observing NoseAG+ ThroatAG+ ThroatPCR- from a true negative::
        sp_prob[4,p] <- (1-prev[p]) * ((1-sp[1])*(1-sp[2])*sp[3] +covsp12 -covsp13 -covsp23)

        # Probability of observing NoseAG- ThroatAG- ThroatPCR+ from a true positive::
        se_prob[5,p] <- prev[p] * ((1-se[1])*(1-se[2])*se[3] +covse12 -covse13 -covse23)
        # Probability of observing NoseAG- ThroatAG- ThroatPCR+ from a true negative::
        sp_prob[5,p] <- (1-prev[p]) * (sp[1]*sp[2]*(1-sp[3]) +covsp12 -covsp13 -covsp23)

        # Probability of observing NoseAG+ ThroatAG- ThroatPCR+ from a true positive::
        se_prob[6,p] <- prev[p] * (se[1]*(1-se[2])*se[3] -covse12 +covse13 -covse23)
        # Probability of observing NoseAG+ ThroatAG- ThroatPCR+ from a true negative::
        sp_prob[6,p] <- (1-prev[p]) * ((1-sp[1])*sp[2]*(1-sp[3]) -covsp12 +covsp13 -covsp23)

        # Probability of observing NoseAG- ThroatAG+ ThroatPCR+ from a true positive::
        se_prob[7,p] <- prev[p] * ((1-se[1])*se[2]*se[3] -covse12 -covse13 +covse23)
        # Probability of observing NoseAG- ThroatAG+ ThroatPCR+ from a true negative::
        sp_prob[7,p] <- (1-prev[p]) * (sp[1]*(1-sp[2])*(1-sp[3]) -covsp12 -covsp13 +covsp23)

        # Probability of observing NoseAG+ ThroatAG+ ThroatPCR+ from a true positive::
        se_prob[8,p] <- prev[p] * (se[1]*se[2]*se[3] +covse12 +covse13 +covse23)
        # Probability of observing NoseAG+ ThroatAG+ ThroatPCR+ from a true negative::
        sp_prob[8,p] <- (1-prev[p]) * ((1-sp[1])*(1-sp[2])*(1-sp[3]) +covsp12 +covsp13 +covsp23)

    }


    ## Priors:

    # Prevalence in population 1:
    prev[1] ~ dbeta(1,1)

    # Prevalence in population 2:
    prev[2] ~ dbeta(1,1)


    # Sensitivity of NoseAG test:
    se[1] ~ dbeta(1,1)T(1-sp[1], )
    # Specificity of NoseAG test:
    sp[1] ~ dbeta(1,1)

    # Sensitivity of ThroatAG test:
    se[2] ~ dbeta(1,1)T(1-sp[2], )
    # Specificity of ThroatAG test:
    sp[2] ~ dbeta(1,1)

    # Sensitivity of ThroatPCR test:
    se[3] ~ dbeta(1,1)T(1-sp[3], )
    # Specificity of ThroatPCR test:
    sp[3] ~ dbeta(1,1)


    # Covariance in sensitivity between NoseAG and ThroatAG tests:
    covse12 ~ dunif( (se[1]-1)*(1-se[2]) , min(se[1],se[2]) - se[1]*se[2] )  ## if the sensitivity of these tests may be correlated
    # covse12 <- 0  ## if the sensitivity of these tests can be assumed to be independent
    # Calculated relative to the min/max for ease of interpretation:
    corse12 <- ifelse(covse12 < 0, -covse12 / ((se[1]-1)*(1-se[2])), covse12 / (min(se[1],se[2]) - se[1]*se[2]))

    # Covariance in specificity between NoseAG and ThroatAG tests:
    covsp12 ~ dunif( (sp[1]-1)*(1-sp[2]) , min(sp[1],sp[2]) - sp[1]*sp[2] )  ## if the specificity of these tests may be correlated
    # covsp12 <- 0  ## if the specificity of these tests can be assumed to be independent
    # Calculated relative to the min/max for ease of interpretation:
    corsp12 <- ifelse(covsp12 < 0, -covsp12 / ((sp[1]-1)*(1-sp[2])), covsp12 / (min(sp[1],sp[2]) - sp[1]*sp[2]))

    # Covariance in sensitivity between NoseAG and ThroatPCR tests:
    covse13 ~ dunif( (se[1]-1)*(1-se[3]) , min(se[1],se[3]) - se[1]*se[3] )  ## if the sensitivity of these tests may be correlated
    # covse13 <- 0  ## if the sensitivity of these tests can be assumed to be independent
    # Calculated relative to the min/max for ease of interpretation:
    corse13 <- ifelse(covse13 < 0, -covse13 / ((se[1]-1)*(1-se[3])), covse13 / (min(se[1],se[3]) - se[1]*se[3]))

    # Covariance in specificity between NoseAG and ThroatPCR tests:
    covsp13 ~ dunif( (sp[1]-1)*(1-sp[3]) , min(sp[1],sp[3]) - sp[1]*sp[3] )  ## if the specificity of these tests may be correlated
    # covsp13 <- 0  ## if the specificity of these tests can be assumed to be independent
    # Calculated relative to the min/max for ease of interpretation:
    corsp13 <- ifelse(covsp13 < 0, -covsp13 / ((sp[1]-1)*(1-sp[3])), covsp13 / (min(sp[1],sp[3]) - sp[1]*sp[3]))

    # Covariance in sensitivity between ThroatAG and ThroatPCR tests:

    covse23 ~ dunif( (se[2]-1)*(1-se[3]) , min(se[2],se[3]) - se[2]*se[3] )  ## if the sensitivity of these tests may be correlated
    # covse23 <- 0  ## if the sensitivity of these tests can be assumed to be independent
    # Calculated relative to the min/max for ease of interpretation:
    corse23 <- ifelse(covse23 < 0, -covse23 / ((se[2]-1)*(1-se[3])), covse23 / (min(se[2],se[3]) - se[2]*se[3]))

    # Covariance in specificity between ThroatAG and ThroatPCR tests:
    covsp23 ~ dunif( (sp[2]-1)*(1-sp[3]) , min(sp[2],sp[3]) - sp[2]*sp[3] )  ## if the specificity of these tests may be correlated
    # covsp23 <- 0  ## if the specificity of these tests can be assumed to be independent
    # Calculated relative to the min/max for ease of interpretation:
    corsp23 <- ifelse(covsp23 < 0, -covsp23 / ((sp[2]-1)*(1-sp[3])), covsp23 / (min(sp[2],sp[3]) - sp[2]*sp[3]))

}

#monitor# se, sp, prev, covse12, corse12, covsp12, corsp12, covse13, corse13, covsp13, corsp13, covse23, corse23, covsp23, corsp23

## Inits:
inits{
"se" <- c(0.5, 0.99, 0.5)
"sp" <- c(0.99, 0.75, 0.99)
"prev" <- c(0.05, 0.95)
"covse12" <- 0
"covse13" <- 0
"covse23" <- 0
"covsp12" <- 0
"covsp13" <- 0
"covsp23" <- 0
}
inits{
"se" <- c(0.99, 0.5, 0.99)
"sp" <- c(0.75, 0.99, 0.75)
"prev" <- c(0.95, 0.05)
"covse12" <- 0
"covse13" <- 0
"covse23" <- 0
"covsp12" <- 0
"covsp13" <- 0
"covsp23" <- 0
}

## Data:
data{
"Populations" <- 2
"N_RRR" <- c(9823, 10177)
"Tally_RRR" <- structure(c(8693, 491, 451, 25, 91, 21, 20, 31, 8604, 539, 470, 27, 135, 80, 134, 188), .Dim = c(8, 2))
}

And can be run directly from R:

results <- run.jags('covidmodel.txt')
results

Template Hui-Walter

• Modifying priors must still be done directly in the model file
  • Same for adding .RNG.seed and the deviance monitor
• The model needs to be re-generated if the data changes
  • But remember that your modified priors will be reset
• There must be a single column for the population (as a factor), and all of the other columns (either factor, logical or numeric) are interpreted as being test results

Activating covariance terms

Find the lines for the covariances that we want to activate (i.e. the two Throat tests):

You will also need to uncomment out the relevant initial values for BOTH chains:

Session Summary

• Correlation terms add complexity to the model in terms of:
  • Opportunity to make a coding mistake
  • Reduced identifiability