Tutorial: statistical methods for the meta-analysis of diagnostic test accuracy studies

Motivating example

Worldwide, sepsis and its sequelae still remain a frequent cause of acute illness and death in patients with community and nosocomial acquired infections [14]. Sepsis may be seen as systemic inflammatory response due to infection. However, a gold standard for the proof of infection is missing. Depending on prior antibiotic therapy, bacteremia is found only in approximately 30% of patients with sepsis. Furthermore, early clinical signs of sepsis, like fever, tachycardia, and leucocytosis, are unspecific and overlap with signs also seen in a multitude of systemic inflammatory response syndromes (SIRS) in the absence of infection, especially in surgical patients. Other signs, such as arterial hypotension, thrombocytopenia, or elevated lactate levels indicate, too late, the progression to organ dysfunction. Thus, delay in diagnosis and treatment of sepsis causes increased mortality.

In sepsis numerous humoral and cellular systems are activated, followed by a release of a multitude of mediators and other molecules that mediate the host response to infection. Several potential diagnostic indicators measured in the bloodstream have been evaluated for their clinical ability to assess the diagnosis and severity of sepsis. One of these, the 116 amino acid polypeptide procalcitonin (PCT) is frequently used when it comes to identify bacterial infections.

In this tutorial we will use Procalcitonin as an example for a meta-analysis of DTA studies using data from [15]. This is a meta-analysis of Procalcitonin (PCT) for diagnosis of sepsis in critically ill patients. Data sources were Medline, Embase, ISI Web of Knowledge, the Cochrane Library, Scopus, BioMed Central, and Science Direct, from inception to Feb 21, 2012, and reference lists of identified primary studies. Articles written in English, German, or French that investigated Procalcitonin for differentiation of septic patients – those with sepsis, severe sepsis, or septic shock – from those with a systemic inflammatory response syndrome of non-infectious origin were included. Excluded studies were studies of healthy people, patients without probable infection, and children younger than 28 days. Two independent investigators extracted patient and study characteristics, discrepancies were resolved by consensus. The search returned 3,487 reports, of which 31 fulfilled the inclusion criteria, accounting for 3,244 patients. Table 1 shows PCT data for diagnosis of sepsis which were extracted from the 31 studies.

Forest plots for sensitivity and specificity

One way to perform a diagnostic meta-analysis is to analyze sensitivity and specificity separately as those are key parameters when evaluating the performance of a binary diagnostic test [16]. This requires knowledge of a reference or gold standard which denotes the disease status D. The potential outcomes of a 2 × 2 table showing the disease status D in the columns and test results T in the rows are shown in Table 2. For a detailed description and examples see e.g. Schlattmann [17].

In a diagnostic meta-analysis we have for each individual study (i=1,…,k) study specific sensitivities (true positive rate, TPR) with Se i ˆ = TP i TP i + FN i . Here TP _i are the true positives and FN _i are the false negatives according to a gold standard. Looking at the study from Ahmadinejad (2009) we find a sensitivity equal to Se ˆ = 63 63 + 8 = 0.887 . Likewise, for each study specificity (true negative rate, TNR) is given by Sp ˆ i = TN i TN i + FP i . Here TN _i are the true negatives and FP _i are the false positives. Again, for the study by Ahmadinejad specificity is given by spec = 38 38 + 11 = 0.776 .

For a graphical presentation for each study sensitivity and specificity are calculated together with a 95% confidence interval and displayed in a so called forest plot. There are several ways to construct a confidence interval for a binomial proportion with different statistical properties [18, 19]. Figure 1 shows a forest plot of sensitivity of PCT on the left hand side and of specificity on the right hand side. In this plot set we see considerable heterogeneity for sensitivity ranging from 0.415 to 0.969. Likewise for specificity we find heterogenous results which range from 0.526 to 1.000.

Figure 1:

Univariate forest plots for sensitivity and specificity in Procalcitonin studies [15].

Fixed and random effects models

Statistically speaking sensitivity and specificity are proportions and can be treated as such in a meta-analysis [20].

Standard fixed effects models for meta-analyses can be applied [21], [22], [23]. One approach is using log transformed odds of sensitivity and specificity (logit transform). Odds are defined as odds = p 1 − p where p is a probability. If we have a fair coin the odds for head are 0.5 1 − 0.5 = 1 meaning that head and tail are equally likely. A logit transform of sensitivity is given as logit ( Se ) = log ( Se 1 − Se ) , where log denotes the natural logarithm.

Summary estimates of logit-transformed sensitivity and logit-transformed specificity, respectively are obtained as a weighted average of the respective logit-transformed proportions of the individual studies. Weights are given by the inverse of the respective study specific variances. This has the disadvantage that in the case of zero entries undefined log odds occur. Thus, in the past years there has been a lively discussion how to avoid undefined log odds [24], [25], [26] by adding e.g. 0.5 to each cell of the study specific 2 × 2 table in case of zero cells.

To avoid this, we apply logistic regression models potentially with random effects aka generalized linear mixed models. That is, we assume that sensitivity and specificity respectively follow a binomial distribution. Thus, for each study:

A common effect logistic regression model for sensitivity has the form

This is a generalized linear model with binomial errors, linear predictor β ₀ and logistic link function. The left hand side shows the natural logarithm of the odds of sensitivity. The unknown parameter β ₀ can be estimated using maximum likelihood using numerous statistical software packages such as R [13]. Also, this is a so called common effect model, since it assumes the overall sensitivity in each study is identical and given by

For the PCT data an application of two univariate common effects models for sensitivity and specificity yields the results presented in Table 3. Overall, assuming a common effect we find a sensitivity equal to 0.735 with 95% CI (0.715, 0.755) and a specificity equal to 0.747 with 95% CI (0.723, 0.769).

A common effect model assumes that the underlying true sensitivity is the same in each study. The overall variation and, therefore, the confidence intervals will reflect only random variation within each study but not any potential heterogeneity between the studies. Of course, the same applies for specificity.

Whether pooling of the data in this way is appropriate should be decided after investigating the heterogeneity of the study results. If the results vary substantially, no fixed effects pooled estimator should be presented [27]. As a result only estimators e.g. for selected subgroups should be calculated. The previous remark notwithstanding, a fixed effects meta-analysis is always valuable, since it tests the null-hypothesis that diagnostic accuracy was identical in all trials [28]. If the null-hypothesis is rejected then the alternative may be asserted that at least one study differs.

One way to address heterogeneity is the calculation of Cochran’s Q-statistic and the I ² measure. This describes the percentage of the variability in effect estimates that is due to heterogeneity rather than sampling error (chance). For sensitivity we find I ˆ 2 =69.8% and for specificity I ˆ 2 =67.4%. Likewise, the test for heterogeneity turns out to be statistically significant (sensitivity Q=99.4, df=30, p<0.001, specificity Q=92.06, df=30, p<0.001).

Thus, the investigation of heterogeneity between studies is a main task in each meta-analysis [29]. Here a common effect model is not appropriate. Alternatively, a random effects model which incorporates variation between studies should be considered.

A random effects logistic regression model has then the form

Again, this is a generalized linear model with binomial errors, linear predictor β ₀ and logistic link function. Additionally, we assume variability between studies given by the study specific departure b _i from the overall intercept β ₀. For the b _i a normal distribution with expectation zero and heterogeneity variance τ ² is assumed. The latter indicates variability between studies, i.e. heterogeneity. Both unknown parameters again can be estimated using maximum likelihood. Table 4 shows the result for the PCT data.

Overall, assuming a random effects model we find in Table 4 a sensitivity equal to 0.768 with 95% CI (0.720, 0.810) with heterogeneity variance τ ˆ 2 =0.36. For specificity, we find a value equal to 0.793 with 95% CI (0.743, 0.836) and heterogeneity variance τ ˆ 2 =0.379. Thus, we find substantial heterogeneity between studies.

Overall this approach seems to provide useful results in terms of sensitivity and specificity as e.g. investigated by Simel and Bossuyt [30]. However, we do not have any information on the correlation between sensitivity and specificity and the magnitude of the overall diagnostic performance.

Diagnostic odds ratio (DOR)

So far, we have considered sensitivity and specificity as a pair for each study. There have been many attempts to merge the results of a diagnostic study into one single measure. One proposal is the diagnostic odd ratio (DOR) [11, 31]

This is the ratio of the odds of a positive test result for a person with the disease divided by the odds for a positive test result for a healthy person. The value of a DOR ranges from 0 to infinity, where higher values indicate better discriminatory test performance. The synthesis of diagnostic odds ratios is straightforward and follows standard meta-analysis methods. Summary estimates of diagnostic odds ratios are obtained as a weighted average of the respective log transformed DORs of the individual studies. The weights are given by the inverse of the respective study specific variances.

First, investigating heterogeneity between studies we find substantial heterogeneity (Q=89.00, df=30, p<0.001, I ˆ 2 =66.3%. Thus, we apply a random effects mode with an overall DOR=11.698, 95% CI (8.301, 16.486). Thus, we see discriminatory potential of Procalcitonin for the diagnosis of sepsis.

Apart from challenges in interpreting diagnostic odds ratios, a disadvantage is that it is impossible to weight the true positive and false positive rates separately. Likewise, it is impossible to distinguish between tests with high sensitivity and low specificity and tests with low sensitivity and high specificity. Furthermore no direct investigation of the correlation between sensitivity and specificity is possible. Thus, bivariate models are preferable and introduced in Section 3.2.

Publication bias

Publication bias is a major form of bias in any meta-analysis. That is, if the studies that are included in a review have results that systematically differ from relevant studies that are missed, then the findings will be compromised by publication bias. Thus, researchers are advised to perform a thorough literature search and to investigate publication bias. Following Deeks et al. [32] we present the effective sample size funnel plot together with the associated regression test of asymmetry. The effective sample size plot (Figure 2) takes the DOR on the x-axis of the plot and 1 / ESS on the y-axis. ESS stands for effective sample size and is proportional to 1 n 1 + 1 n 2 . This test is based on the regression of log(DOR) against 1 / ESS , weighting by ESS.

Figure 2:

Diagnostic odds ratio against 1 / ESS , where ESS stands for effective sample size.

Unfortunately, for our example there is publication bias present (Test result: t=4.11, df=29, p-value=0.0003). More details are shown in Section 4.3. As result in a first step a repeated literature search would take place.

Plots of sensitivity and specificity in the summary receiver operator curve (sROC) space

Procalcitonin is a continuous diagnostic marker. Until now we have assumed that we are dealing with a binary diagnostic test. A frequent cut-off value equals 0.5 g/L. Values larger or equal than 0.5 g/L indicate a positive test and smaller values indicate a negative test result and thus we have transformed the continuous marker Procalcitonin into a binary test. Obviously, other cut-off values could be used. For example we could apply a cut off value ≥2.0 g/L. As a result increasing the cut-off value from 0.5 g/L to 2.0 g/L will lead to a decreased sensitivity and an increased specificity. This idea is depicted in Figure 3.

Figure 3:

Illustration of the cut-off value problem for Procalcitonin. Variation of the cut-off value c leads to an increased specificity and decreased sensitivity if c is moved to the right, and vice versa if c is moved to the left.

Descriptive statistics of the Procalcitonin data applied in Schlattmann [17] find a median PCT value equal to 0.2 g/L with a minimum equal to 0.01 g/L and a maximum of 200 g/L. Obviously, we could use any value between minimum and maximum as a cut off value and calculate the corresponding sensitivity and specificity.

This is done when we create a receiver operator curve (ROC) [33] which is obtained by calculating the sensitivity and specificity of every observed data value and plotting sensitivity against 1-specificity. A test that perfectly discriminates between the two groups would yield a “curve” that coincided with the left and top sides of the plot since we would not have any false negative (FN) or false positive (FP) values. A useless test would give a straight line from the bottom left corner to the top right. This implies that a true positive and a false positive test result are equally likely.

The performance of the test can be assessed by using the area under the receiver operating characteristic curve (AUC). This area may be interpreted as the probability that a random person with the disease has a higher value of the measurement than a random person without the disease. A perfect test would have an AUC=1 and a useless test has an AUC=0.5. This is shown in Figure 4.

Figure 4:

Procalcitonin as biomarker for sepsis: Receiver operator curve and area under the curve. The point shows sensitivity and specificity for a cut-off value of 0.5 g/L.

In diagnostic meta-analyses often only a single cut-off value for a specific study is provided. Hence not the study specific ROC curve is available but only the corresponding TP; FN, FP and TN as shown in Table 1, where e.g. Ahmadinejad applies a cut-off value of 0.5 g/L.

In order to display variation between studies due to different cut-off values plots in ROC space may be constructed. Here, a simple scatterplot of sensitivity vs. 1-specificity of each study is useful. Additional information showing also the variability within a study is shown in a cross-hair plot [34] which shows 1-specificity (false positive rate) vs. sensitivity together with the respective study specific 95% confidence intervals.

In Figure 5 the scatterplot shows variation in cut-off points as well in accuracy. Looking at the cross-hair plot on the right side we see also high variability of sensitivities and false positive rates indicating considerable heterogeneity.

Figure 5:

Scatterplot and cross-hair plot in ROC-space for the Procalcitonin data.

Univariate meta-analyses provide single estimates of sensitivity and specificity. Here, we might be interested in a joint pair together with a confidence region. Also we saw that heterogeneity is common in DTA studies. One reason is variation in cut-off points used in the individual studies. Another reason might be due to differences in the respective patient populations. Thus, we might be interested in a prediction region which shows where future studies might fall. Finally, the construction of a summary ROC curve across studies (sROC) might be of interest. These aims can be reached using an appropriate model, that is a bivariate statistical model.

Bivariate generalized linear mixed modes

The logistic models used so far have the disadvantage to ignore the bivariate structure of the data. Thus, frequently a bivariate linear random effects model is used for a DTA meta-analysis which was introduced by Reitsma et al. [35]. This model uses logit-transformed sensitivity and logit-transformed specificity simultaneously. Here it is assumed that the true logit-transformed sensitivities of the individual studies follow a normal distribution with a common mean value and between-study variability as in the univariate random effects model. Variation between studies can be attributed to unobserved heterogeneity due to e.g. heterogeneous study populations. Likewise, for the true logit-transformed specificities a normal distribution with a common mean value and between-study variability is assumed.

Now, this model introduces potential correlation between the true logit-transformed sensitivity and specificity within studies by assuming a bivariate normal distribution for the random effects. Besides variability between studies in the true underlying sensitivities and specificities, there is also variation due to sampling. Studies differ in size and thus in variation. Thus, on the second level of the model study specific variances of logit-transformed sensitivity and specificity are incorporated in order to take sampling variability into account.

As a result of this bivariate model approach summary estimates for sensitivity and specificity are obtained. In addition, based on the model’s assumption of bivariate normality an sROC curve can then be constructed from the parameter estimates of the model. Performing a bivariate linear random effects model for meta-analysis of diagnostic accuracy can be done using the ‘reitsma’ function implemented in the freely available R-package ‘mada’ [36].

However, to synthesize data, an exact binomial rendition [37] of the linear bivariate mixed-effects regression model developed by van Houwelingen et al. [38] for meta-analysis of treatment trials, modified for synthesis of diagnostic test data builds an alternative. As in the linear mixed effects model the correlation between sensitivity and specificity is taken care of. Furthermore, in contrast to a logit transformation no ad hoc continuity correction to avoid zero cells in the 2 × 2 table is required. Thus, this model is preferable as shown in simulation studies [39] and empirical comparisons [40]. Hence, in the following we concentrate on this bivariate logistic regression model with random effects (bivariate GLMM).

Since we present our results in ROC space we make a slight shift of presentation. We now model the false positive rate, i.e. 1-specificity. As in the case of univariate models we assume a binomial distribution for sensitivity and 1-specificity respectively. Hence the binomial distribution depicts within study variability of the i=1,…,k studies:

A bivariate random effects logistic regression model has then the form

Between study variability is addressed using a bivariate normal distribution with

Here Σ denotes the covariance matrix of the bivariate random effects distribution, where σ _μ ² denotes the between study variability of sensitivity on the logit scale. Likewise σ _ν ² denotes the between study variability of 1-specificity on the logit scale, whereas ρ denotes the correlation between sensitivity and 1-specificity. Estimation can again be done using maximum likelihood with general statistical software such as R as shown in Section 4.4.

For our example we obtain the following results shown in Table 5.

Based on the bivariate mixed effects logistic regression model we obtain an overall sensitivity equal to 0.767 and an overall specificity equal to 0.792. In terms of heterogeneity we find a variance between studies for sensitivity on the logit scale σ _μ ² equal to 0.357 and likewise for 1-specificity a variance σ _ν ² equal to 0.384. Importantly, we find a positive correlation, which implies a negative correlation=−0.23 between sensitivity and specificity. Only in this case the construction of a sROC curve is recommended [41].

Summary receiver operator curve (sROC curve)

According to item 21 of the PRISMA statement for DTA meta-analyses, [42], test accuracy, including variability should be reported. This includes summary results as well as confidence and prediction intervals respectively.

One way to address diagnostic test accuracy is to estimate the receiver operator curve based on the available data from the different studies. There are several methods available for sROC curve construction [43]. Here we apply the regression line of logit transformed sensitivity η based on logit transformed 1-specifcicity ξ. That is

(6) η = β 0 + ρ σ μ σ ν σ v 2 ( ξ − β 1 )

When transformed to the ROC space we obtain the sROC curve indicating the median sensitivity for a specific false positive rate. Figure 6 shows the sROC curve, the joint estimate of sensitivity and 1-specificity together with a 95% confidence and prediction region. This prediction region indicates the extent of statistical heterogeneity by depicting a region within which, assuming the model is correct, we have 95% confidence that the true sensitivity and specificity of a future study will take place. Obviously, for Procalcitonin we find substantial heterogeneity.

Figure 6:

Procalcitonin as biomarker for sepsis: sROC curve (solid line). The point shows the joint estimate of sensitivity and 1-specificity together with a 95% confidence (dashed line) region and 95% prediction region (dotted line).

When evaluating the diagnostic performance of a biomarker the area under the curve is of interest. To restrict the computation of the AUC to the observed false positive rates leads to the partial area under the curve (pAUC). This summary index is considered to be more practically relevant than the area under the entire ROC curve (AUC) because it avoids extrapolation. For the data at hand we obtain a pAUC equal to 0.629 and for completeness an AUC=0.799 indicating helpful diagnostic performance.

Importing and manipulating data

Make sure you are working in the right directory. Please give the path to your directory where you save the file containing the data. For example:

The data from our example are read from an Excel .csv file and stored under the name ’PCT’. The command ‘read.csv2’ reads Excel files in .csv format. First comes the name of the file. Then ‘header=T’ implies that the first line contains the variable names.

The object ‘PCT’ contains the data and can be modified. For example the data column TP’ contains the true positives as explained in Table 2.

The command ’attach’ provides access to the individual elements of the data object ‘PCT’.

In a first step we create a new variable called ‘n1’ which is a new column in our data set. To do this the syntax ‘PCT$n1’ is applied. Important, by using ‘PCT$n1’ a new column ‘n1’ is added to the dataframe ‘PCT’. This variable contains the total number of diseased persons per study and is given as the sum of true positives TP and false negatives FN. In a similar way we create the variable ‘n2’, i.e. the total number of healthy individuals. The command ‘head’ shows the first six lines of the dataframe ‘PCT’.

The symbol ‘#’ indicates a comment which will not be executed by the program.

Two univariate meta-analyses

Diagnostic odds ratio (DOR) and publication bias

For the calculation of the DOR and the assessment of publication bias we use again the library ‘meta’ with the function ‘metabin’. Necessary arguments are the true positives TP, the number of diseased n ₁, the false positives FP and the number of healthy subjects n ₂.

Based on the stored results in ‘m.dor’ we show the results with ‘summary(m.dor)’ and construct the funnel plot shown in Figure 2.

Next we perform the regression test for publication bias:

Clearly, we find publication bias. In real life this needs to be investigated further, e.g. by a repeated literature search.

Bivariate meta-analysis

Bivariate logistic regression model with random effects

In order to use the bivariate logistic regression model with random effects we first need to transpose the data from ‘wide’ to ‘long’ format. Furthermore, we need new variables indicating disease status called ‘healthy’ and ‘diseased’. Also we need new outcome variables “positive” for positive test results and “negative” vice versa. We can use the function ‘reshape’ where we create a new dataframe under the name ‘long’. Next, the new variable ‘healthy’ as ‘1-diseased’ is created and the data are sorted by study ‘id’.

Looking at the first six lines of the dataframe ‘long’ with the command ‘head(long)’ gives:

The data in long form contains the necessary information for the calculation of the bivariate random effects logistic model. Next, we use the function ‘glmer’ from the library ‘lme4’. Now, our dependent variable is the combination of positive test results and negative test results in a matrix of the respective columns of the dataframe ‘long’ using.

Now, apply the covariate “healthy” for healthy subjects coded ‘1’ if true and ‘0’ otherwise. This covariate quantifies the mean false positive rate on the logit scale. The covariate ‘diseased’ is coded ‘1’ for diseased subjects and ‘0’ otherwise and quantifies the mean true positive rate on the logit scale (i.e. sensitivity).

We do not want an intercept, thus our formula is ‘∼0+healthy + diseased’ for the fixed effects. For bivariate random effects we use ‘+(0+healthy + diseased| Study)’. Finally, we assume a binomial distribution thus family=binomial). Hence the following code is applied and the result is stored as ‘pct.glmm2’.

Using the command ‘summary’ with ‘pct.glmm2’ as an argument shows the results. Let’s start with the covariance matrix of the random effects:

We are interested in variability between studies, thus ‘Groups name’ refers to ‘Study’. The variance equal to 0.3565 associated to ‘diseased’ is the variability between studies for the true positive rates σ _μ ² in Eq (5). Likewise the variance equal to 0.3838 refers to the variability between studies σ _ν ² for false positives rates. Finally, ‘rho’ denotes the correlation ρ in Eq (5).

Next we look at the fixed effects:

The estimate of the covariate ‘diseased’ is an estimate of β ₀ equal to 1.1892. Likewise the estimate −1.3395 for the covariate healthy is an estimate of β ₁ which denotes the false positive rate on the logit scale.

In order obtain the estimates on the original scale we use again the command ‘lsmeans’

The first line of the output with the covariates ‘diseased’ equal to 1 and healthy equal to ‘0’ refers to sensitivity which is equal to 0.767. Likewise, the second line shows the false positive rate equal to 0.208. This completes the results shown in Table 5.