Exact inference for the risk ratio with an imperfect diagnostic test

Summary

The risk ratio quantifies the risk of disease in a study population relative to a reference population. Standard methods of estimation and testing assume a perfect diagnostic test having sensitivity and specificity of 100%. However, this assumption typically does not hold, and this may invalidate naive estimation and testing for the risk ratio. We propose procedures that control for sensitivity and specificity of the diagnostic test, given the risks are measured by proportions, as it is in cross-sectional studies or studies with fixed follow-up times. These procedures provide an exact unconditional test and confidence interval for the true risk ratio. The methods also cover the case when sensitivity and specificity differ in the two groups (differential misclassification). The resulting test and confidence interval may be useful in epidemiological studies as well as in clinical and vaccine trials. We illustrate the method with real-life examples which demonstrate that ignoring sensitivity and specificity of the diagnostic test may lead to considerable bias in the estimated risk ratio.

INTRODUCTION

The risk ratio or relative risk (RR) quantifies the risk in a study population (say, in those exposed to a risk factor) relative to a reference population. Naive inference assumes that the applied diagnostic procedure is perfect (its sensitivity and specificity is 100%), i.e. no misclassification occurs. Unfortunately, this assumption usually does not hold, and ignoring this may result in misleading conclusions. Our aim is to present an exact unconditional test and confidence interval for RR controlling for sensitivity and specificity of the diagnostic test.

We address the situation, in which

• • risks are quantified by the proportion of those having the condition (disease, recovery, etc.), which is typical in cross-sectional epidemiological studies but may also occur in other study designs;

• • two independent random samples are drawn from the two populations;

• • there is no classification error in the group definition (i.e. in the exposure);

• • sensitivity and specificity of the diagnostic procedure are known (rather than estimated in the same or in another study).

In this case RR is defined as the ratio of two proportions, RR = p ₂/p ₁, where p ₂ is the proportion of diseased in the study group and p ₁ is that in the reference group. Many authors use the term ‘prevalence ratio’ (PR) for this measure preserving the term ‘risk ratio’ for the incidence ratio [1–3], while others call it ‘prevalence risk ratio’ (PRR) [4, 5]. In the following text we will use the term risk ratio (RR). In cross-sectional studies and in therapeutic or vaccine trials with fixed-length follow-up this is the most natural measure to compare the groups. Guidelines for vaccine studies define vaccine efficacy as 1 − RR, where the reference group is placebo. Despite this, in many studies odds ratios (OR) are calculated just because logistic regression has become a standard analysis tool readily available in most statistical software systems, although using OR instead of RR is repeatedly criticized by statisticians [1, 6–8].

While the impact of misclassification on the results of statistical analyses has been studied since the 1950s in the biomedical as well as in the social sciences [9–16], no exact test or confidence interval for the true PR has been proposed. In the next section we describe the proposed procedures, then we present two applications, finally we summarize the properties of the method. R code for the procedures is available at www2.univet.hu/users/jreiczig/RR_SeSp.

METHODS

Let us denote for population i (i = 1, 2) the true prevalence by p _i , the sensitivity by Se _i , and the specificity by Sp _i . Then the probability of a positive diagnosis (also called observed or apparent prevalence) in the ith population is p _ia  = p _i ·Se _i  + (1 − p _i )·(1 − Sp _i ). This implies that taking independent samples of sizes n ₁ and n ₂ from the two populations, the number of test positives x ₁ and x ₂ follow the binomial distribution with parameters n ₁, p _1a and n ₂, p _2a . Thus, the relative frequencies x ₁/n ₁ and x ₂/n ₂ are estimates of p _1a and p _2a , therefore we will denote them by ${\hat p}_{1a}$ and ${\hat p}_{2a}$ . What we will make use of in the following text is that the parameters p _1a and p _2a are in a one-to-one relationship with the true prevalences p ₁ and p ₂, i.e. a hypothesis about p ₁ and p ₂ can be mapped onto a corresponding hypothesis about p _1a and p _2a , and tested using their estimates ${\hat p}_{1a}$ and ${\hat p}_{2a}$ . The general equation, which describes the relationship between the parameters p _i and p _ia , and allows for Se ₁≠Se ₂ and/or Sp ₁≠Sp ₂ is

(1) $$\eqalign{p_{2a} &= p_{1a}{\rm RR}\displaystyle{{(Se_2 + Sp_2 - 1)} \over {(Se_1 + Sp_1 - 1)}}\cr &\quad - {\rm RR} (1 - Sp_1)\displaystyle{{(Se_2 + Sp_2 - 1)} \over {(Se_1 + Sp_1 - 1)}} + (1 - Sp_2).}$$

If Se ₁ = Se ₂ and Sp ₁ = Sp ₂, which can often be assumed in real-life applications, the equation simplifies to

(2) $$p_{2a} = p_{1a}{\rm RR} + \left( {1 - Sp} \right)\left( {1 - {\rm RR}} \right),$$

where Sp denotes the common specificity; and for a perfect diagnostic test, i.e. for Se ₁ = Se ₂ = Sp ₁ = Sp ₂ = 1, it reduces to

(3) $$p_{2a} = p_{1a}{\rm RR}.$$

Solving the equations (1)–(3) for RR, the following expressions are obtained:

(4) $${\rm RR} = \,\displaystyle{{(\,p_{2a} + Sp_2 - 1)(Se_1 + Sp_1 - 1)} \over {(\,p_{1a} + Sp_1 - 1)(Se_2 + Sp_2 - 1)}},$$

(5) $${\rm RR} = \left( {p_{2a} + Sp - 1} \right)/\left( {p_{1a} + Sp - 1} \right),$$

(6) $${\rm RR} = p_{2a}/p_{1a}.$$

The point estimates for RR can be obtained by replacing p _1a and p _2a by ${\hat p}_{1a}$ and ${\hat p}_{2a}$ in equations (4)–(6). Note that the parameter space for the true prevalences (p ₁, p ₂) is the unit square, whereas that for (p _1a , p _2a ) is a rectangle within the unit square, namely [1 − Sp ₁, Se ₁] × [1 − Sp ₂, Se ₂]. Note also that the estimates ( ${\hat p}_{1a}$ , ${\hat p}_{2a}$ ) form a point in the sample space (which is actually a grid of points), rather than in the parameter space. In formula

$$({\hat p}_{1a},\,{\hat p}_{2a})\, \in \,\left\{ {0, 1/n_1, 2/n_1, \ldots, 1} \right\} \times \left\{ {0, 1/n_2, 2/n_2, \ldots, 1} \right\}.$$

Assume now that we want to test for H ₀: RR = RR₀, where RR = p ₂/p ₁ is the true risk ratio. This H ₀ is a composite hypothesis, corresponding to a line segment in the parameter space of the true prevalences p ₁ and p ₂, namely the set of points in the unit square satisfying the equation p ₂ = p ₁RR₀. If we map this onto the parameter space of the observed binomials p _1a and p _2a , it will form another line segment. Its position depends on RR₀ as well as on the sensitivities and specificities according to equations (1)–(3), but it is always located within the rectangle [1 − Sp ₁, Se ₁] × [1 − Sp ₂, Se ₂]. Figure 1 illustrates the position of this line segment for H ₀: RR = 2 depending on the sensitivities Se ₁, Se ₂, and specificities Sp ₁, Sp ₂.

Fig. 1. Lines corresponding to the hypothesis H ₀: RR = 2 in the two-dimensional space of the parameters p _1a and p _2a of the observed binomial variables. The position of the line depends on sensitivity (Se ₁ and Se ₂) and specificity (Sp ₁ and Sp ₂) of the test in the two populations.

Testing for H ₀ is equivalent to testing whether the parameters p _1a and p _2a of the observed independent binomial variables are located on the line segment corresponding to H ₀. As this is also a composite hypothesis, it can be tested applying the intersection-union principle [17], which means that a critical (or rejection) region for a composite H ₀ can be obtained by constructing an appropriate critical region for each element of H ₀, and taking the intersection of these regions. The steps of constructing this critical region follow the logic of Reiczigel et al. [18].

1. (1) For each simple hypothesis h ₀ ∈ H ₀, i.e. for each point (p _1a , p _2a ) on the line segment representing H ₀, we construct a critical region (rejection region) C _h0 in the sample space, consisting of those points, which have probability under h ₀ less than or equal to the probability of the observed point ( ${\hat p}_{1a}$ , ${\hat p}_{2a}$ ). In formula

   $$\eqalign{C_{h0} =\,& \{ \left( {i/n_1,j/n_2} \right)\!\!:i \in \left\{ {0, \ldots, n_1} \right\},j \in \cr & \left\{ {0, \ldots, n_2} \right\},P_{h0} \left( {i/n_1,j/n_2} \right)\, \les \,P_{h0} ({\hat p}_{1a},\,{\hat p}_{2a})\},}$$
   where P _h0() denotes the probability of a point (or a set of points) in the sample space, given h ₀ is true. C _h0 can be considered as the two-dimensional generalization of that proposed by Sterne [19] for a single binomial proportion. Let P _h0(C _h0) denote the probability of C _h0 under h ₀, and let M = max{P _h0(C _h0), h ₀ ∈ H ₀}.

2. (2) Next, for each h ₀ ∈ H ₀ we determine the subset S _h0 of the sample space consisting of the points with the smallest probability under h ₀ so that P _h0(S _h0)⩽M but adding any further point to S _h0 would result in P _h0(S _h0) > M. Let C _H0 denote the intersection of all these subsets, i.e. $C_{H0}=$ $\cap _{h_0 \in H_0}S_{h0}$ . It is easy to see that C _H0 contains ( ${\hat p}_{1a}$ , ${\hat p}_{2a}$ ) on its boundary.

3. (3) Finally, the P value is defined as the highest probability of C _H0 under H ₀ (i.e. for all simple hypotheses h ₀ in H ₀). In formula, the P value is equal to max{P _h0(C _H0), h ₀ ∈ H ₀}.

Figure 2 illustrates how the resulting critical region for H ₀: RR = 2 depends on sensitivities and specificities, given the two observed prevalences are ${\hat p}_{1a}$  = 0·575 (n ₁ = 40) and ${\hat p}_{2a}$  = 0·667 (n ₂ = 36).

Fig. 2. Critical regions in the sample space for H ₀: RR = 2 with observed proportions ${\hat p}_{1a}$  = 0·575 (n ₁ = 40) and ${\hat p}_{2a}$  = 0·667 (n ₂ = 36), depending on the sensitivities and specificities. Black dots form the critical region, the black square represents the observed data. The line shows the location of H ₀ in the parameter space.

Confidence intervals for the true RR can be constructed by inverting the above test. That is, lower and upper confidence limits to a given confidence level (1 − α) are defined as the smallest and largest true RR₀ not rejected by the test, i.e.

$$L = {\rm inf}\{ {{\rm RR}_0}\!\!:p_{{\rm RR}_0} \gt \alpha\} \,{\rm and}\,U = {\rm sup}\{ {{\rm RR}_0}\!\!:p_{{\rm RR}_0} \gt \alpha\}, $$

where p _RR0 denotes the P value from testing for H ₀: RR = RR₀. Computationally, L and U are determined by increasing RR₀ in small steps and performing the test. Step size may depend on the required precision. In our implementation of the algorithm the default is a multiplicative increment with step size 0·001, i.e. RR₀ is increased or decreased as RR_0,next = 1·001*RR₀ or 0·999*RR ₀. Figure 3 illustrates the procedure for observed proportions ${\hat p}_{1a}$  = 0·575 (n ₁  = 40), ${\hat p}_{2a}$  = 0·667 (n ₂  = 36), and Se ₁ = Se ₂ = 0·91, Sp ₁ = Sp ₂ = 0·8.

Fig. 3. Illustration of the confidence interval construction for observed proportions ${\hat p}_{1a}$  = 0·575 (n ₁ = 40), ${\hat p}_{2a}$  = 0·667 (n ₂ = 36), and Se ₁ = Se ₂ = 0·91, Sp ₁ = Sp ₂ = 0·8. The confidence limits represent the smallest and largest true risk ratio (RR) not rejected by the test.

One-tailed testing, i.e. H ₀: RR = RR₀ against H ₁: RR > RR₀ (or H ₁: RR < RR₀) is also possible, although there are different options to define this. Perhaps the simplest method is that one side of the line representing H ₀ (i.e. the intersection of the critical region with that half-plane) is removed from the critical region. One-sided confidence intervals (CI) can be derived by inverting this one-sided test.

APPLICATIONS

Example 1

Everhart et al. [20] studied the seroprevalence of Helicobacter pylori infection in adults in the United States. The analysis was carried out stratified by age and ethnic group. The infection status was determined by an IgG ELISA assay having 91% sensitivity and 96% specificity in all groups. For illustration we now compare the group of the youngest (20–29 years) and the oldest (⩾70 years), in which the observed seroprevalence was 16·7% and 56·9%, respectively. For these groups, the ratio of observed prevalences is 56·9/16·7 = 3·41 (95% CI 3·00–3·88), whereas the correction for sensitivity and specificity results in the true PR of 4·17 (95% CI 3·58–4·96).

Example 2

Suwancharoen et al. [21] conducted a serological survey of leptospirosis in livestock in Thailand using the microscopic agglutination test (MAT) to determine serostatus of the examined animals. Five animal species were included in the study: cattle, buffaloes, pigs, sheep and goats. Infection status of each species was measured by seroprevalence, and all other species were compared to cattle as the reference group by calculating the PRs (which are same as risk ratios). In this study the MAT test was assumed to be perfect, as this test is usually regarded to be the gold standard test. However, other studies found that the sensitivity of the MAT test is far below 100% [22–25]. Cumberland et al. [22] found the sensitivity to be 30% for first acute-phase specimens, 63% for second acute-phase specimens, and 76% for convalescent specimens. At the same time specificity was 99%, 98%, and 97%, respectively. Limmathurotsakul et al. [25] estimated the sensitivity of MAT by a Bayesian analysis and found it to be 49·8%. These findings indicate the need for an adjustment of the PR estimates.

If we consider the most optimistic scenario of these, i.e. Se = 76% and Sp = 97% [22], we obtain the adjusted PRs shown in Table 1. The difference between the unadjusted and adjusted PR is far from negligible, for example in case of sheep the adjusted ratio is less than half of the unadjusted one. In case of buffaloes, adjusted and unadjusted 95% CIs do not even overlap.

Table 1. Prevalence ratios reported in Suwancharoen et al. [ 21 ] and adjusted prevalence ratios assuming sensitivity = 0·76 and specificity = 0·97 (same in all groups), as reported by Cumberland et al. [ 22 ]

PR, Prevalence ratio; CI, confidence interval.

DISCUSSION

If sensitivity and/or specificity of the diagnostic test is <100%, the observed and true prevalence may differ from each other, influencing also estimation and testing of relative risk measured by the PR. We proposed an exact unconditional test and CIs for the true PR. The method can be applied even if the sensitivities and specificities differ in the two groups, for example if patients are diagnosed by different methods or some sort of differential misclassification occurs.

Taking sensitivity and specificity into account may either increase or decrease the P value compared to the one obtained without considering sensitivity and specificity. For instance, consider testing for H ₀: RR = 2 with observed prevalences ${\hat p}_{1a}$  = 0·48 (n ₁ = 50) and ${\hat p}_{2a}$  = 0·62 (n ₂ = 50). Assuming Se ₁ = Se ₂ = Sp ₁ = Sp ₂ = 1 results in P = 0·0251, for Se ₁ = Se ₂ = 0·8 and Sp ₁ = Sp ₂ = 1 the P value increases to P = 0·0474, whereas for Se ₁ = Se ₂ = 0·6 and Sp ₁ = Sp ₂ = 1 a smaller P value of 0·0155 is obtained.

It may occur that the observed data contradict the given sensitivity and specificity. Let us assume, for example, that a certain diagnostic test is known to have Se = 0·8 and Sp = 1 and using this test we observe 90 positives out of 100. This observation is very unlikely even if the true prevalence is 100%, since under these conditions the observed variable has a binomial distribution with n = 100, P = 0·8, and the probability that it is ⩾90 is as low as 0·0057. The same problem may arise if the number of positives is much less than the expected minimum assuming the given sensitivity and specificity. In such cases one should consider the possibility that sensitivity and/or specificity data are incorrect.

The proposed method can be further developed in several directions. First, one can take into account misclassification also in the group definition, i.e. in the exposure. Brenner et al. [26] investigated this for the incidence ratio in a cohort study. Going further, similar methods could be worked out for models with several predictors, of which the categorical ones may also be affected by misclassification. Some results are known on correcting the OR obtained from logistic regression [27, 28] but similar results still lack for the PR.

There has been a long debate whether the PR or the OR is more appropriate to quantify the risk in the study group relative to the control group in a study design in which both of them can be calculated [1, 7, 29, 30]. Savu et al. [8] stated that RR or PR is more intuitively interpretable than the OR. In spite of this, many studies report OR estimates, even if the design permits calculation of RR or PR. Petersen & Deddens [31] emphasized that in cross-sectional studies, in particular when the disease is not rare, it is preferable to use PR instead of OR. Anyway, it is worth noting that our method can quite easily be adapted to testing the OR or the risk difference (RD), as the hypotheses H ₀: OR = OR₀ and H ₀: RD = RD₀ also correspond to a subset of the parameter space [0,1] × [0,1].

Another direction of potential improvement of the proposed methods is to extend them to the case when sensitivity and specificity are not taken as known but estimated from other samples, which may increase the variance of the RR estimate. This is analogous to the problem solved by Lang & Reiczigel [32] for estimating prevalence.