# Illicit item detection in X-ray images for security applications

1st Georgios Batsis  Department of Informatics and Telematics
Harokopio University of Athens
Athens, Greece
<EMAIL_ADDRESS>2nd Ioannis Mademlis  Department of Informatics and Telematics
Harokopio University of Athens
Athens, Greece
<EMAIL_ADDRESS>3rd Georgios Th. Papadopoulos  Department of Informatics and
Telematics
Harokopio University of Athens
Athens, Greece
<EMAIL_ADDRESS>
###### Abstract

Automated detection of contraband items in X-ray images can significantly
increase public safety, by enhancing the productivity and alleviating the
mental load of security officers in airports, subways, customs/post offices,
etc. The large volume and high throughput of passengers, mailed parcels, etc.,
during rush hours make it a Big Data analysis task. Modern computer vision
algorithms relying on Deep Neural Networks (DNNs) have proven capable of
undertaking this task even under resource-constrained and embedded execution
scenarios, e.g., as is the case with fast, single-stage, anchor-based object
detectors. This paper proposes a two-fold improvement of such algorithms for
the X-ray analysis domain, introducing two complementary novelties. Firstly,
more efficient anchors are obtained by hierarchical clustering the sizes of
the ground-truth training set bounding boxes; thus, the resulting anchors
follow a natural hierarchy aligned with the semantic structure of the data.
Secondly, the default Non-Maximum Suppression (NMS) algorithm at the end of
the object detection pipeline is modified to better handle occluded object
detection and to reduce the number of false predictions, by inserting the
Efficient Intersection over Union (E-IoU) metric into the Weighted Cluster NMS
method. E-IoU provides more discriminative geometrical correlations between
the candidate bounding boxes/Regions-of-Interest (RoIs). The proposed method
is implemented on a common single-stage object detector (YOLOv5) and its
experimental evaluation on a relevant public dataset indicates significant
accuracy gains over both the baseline and competing approaches. This
highlights the potential of Big Data analysis in enhancing public safety.

###### Index Terms:

Illicit item detection, X-ray image analysis, Deep Neural Networks, Object
detection, Non-Maximum Suppression

##  I Introduction

Detecting contraband items using X-ray scanning of luggage, parcels, etc. is a
crucial requirement for ensuring public security (e.g. preventing terrorist
attacks, fighting smuggling of illegal goods, etc.). X-rays are
electromagnetic waves with wavelengths shorter than thee visible light, able
to penetrate most materials; X-ray scanners exploit this fundamental property
to screen items, such as luggage or packages (e.g., in airports, post/customs
offices, etc.). Human operators are able to detect a wide range of potential
threats, such as explosives, weapons, or sharp objects, using high-resolution
images generated by scanning machines [1]. However, fully manual screening has
important shortcomings: the quality of the scan image can be influenced by
several factors, such as occluded objects, cluttered environment or certain
material properties of the scanned items [2], while heavy traffic during rush
hours may mentally overload human security officers. Thus, illicit items may
be missed, due to the need for “the line to keep moving” or because of
perceptual limitations. The high volume and high throughput of X-ray scans in
such scenarios render manual screening ineffective and demand automated Big
Data analysis solutions.

Efficient automated X-ray image analysis/screening for automatic illicit item
detection is nowadays possible thanks to the advances of computer vision and
machine learning. Deep Neural Networks (DNNs) have proven to be remarkably
capable in supporting human operators for similar tasks, thus greatly
increasing their productivity and reducing the possibility of mistakes. Both
whole-image recognition and object detection methods have been proposed for
illicit/contraband item detection in X-ray images, based on DNNs. While the
former ones simply classify an entire image and assign it an overall class
label, algorithms of the latter type identify Regions-of-Interest (RoIs),
i.e., bounding boxes that localize (in 2D pixel coordinates) specific objects
visible in an input image. While there have been significant advancements in
object detection algorithms over the last few decades, achieving sufficient
performance in real-world scenarios continues to be a challenge [3]. The
majority of the proposed methods incorporates mechanisms designed to handle
domain-specific aspects (e.g., high occlusions, very cluttered backgrounds,
large class imbalance, etc.). Additionally, due to the typically cluttered
background of X-ray scan images of luggage, mailed parcels, etc., Non-Maximum
Suppression (NMS) is also particularly important for security applications.
NMS is the final refinement step incorporated to almost every visual object
detection framework, assigned the duty of merging/filtering any spatially
overlapping detected RoIs which correspond to a single visible object [4] [5].

Regarding image recognition, the method of [6] addressed the issue of limited
training data by employing a pretrained CNN and fine-tuning it in the X-ray
domain, while the method of [7] introduced a module named Class-balanced
Hierarchical Refinement (CHR) to enhance the prediction capacity of the CNN
under extreme class imbalance. This is an important issue in automated X-ray
screening, since negative images (where no illicit item is present) are
typically significantly more than the positive ones, with this fact reflected
in the relevant available datasets. CHR is separately evaluated on top of
three different CNNs: Res-Net101 [8], Inception-v3 [9] and Dense-Net [10].

Common DNNs for object detection have also been evaluated with regard to their
discrimination capacity and transferability between different X-ray scanners
[11]; examples include Faster R-CNN [12], Mask R-CNN [13] and RetinaNet [14].
However, modifying fast, anchor-based, single-stage object detectors such as
Single Shot MultiBox Detector (SSD) [15] or You Only Look Once (YOLO) [16] is
the most common approach, due to their ability to operate in real-time even in
embedded computer hardware. Such modifications may have various forms. For
instance, a Cascaded Structure Tensor (CST) is proposed in [17] which took
advantage of contour-based information to extract object proposals; the latter
ones are then classified using a CNN. An alternative lightweight object
detector, called LightRay, is introduced in [18] as a modified version of the
YOLOv4 model for small illicit item detection in complex backgrounds. It
consisted of a fast MobileNetV3 [19] backbone CNN and a feature enhancement
network that includes a Lightweight Feature Pyramid Network (LFPN) [20], to
obtain information of objects at different scales, and a Convolutional Block
Attention Module (CBAM) [21], for refining feature maps through a spatial
attention mechanism.

A different approach was followed in [22], where a novel mechanism called
Foreground and Background Separation (FBS) was proposed for separating illicit
items from complex/cluttered backgrounds. This is achieved by using a feature
extraction DNN combined with Spatial Pyramid Pooling (SPP) and a Path
Aggregation Network, which extracts high-level features. These feature maps
serve as an input to two neural decoders, which reconstruct the background and
the foreground simultaneously. Then, an attention module directs the overall
model’s focus on the foreground objects. In an orthogonal direction, the De-
occlusion Attention Module (DOAM) [23] is a neural module designed to overcome
occlusion in X-ray images; this is important because occlusions are common,
due to the absorption of X-rays by certain materials, such as metals, and the
visual overlap of multiple objects within densely packed parcels. DOAM
consists of two sub-modules, named Edge Guidance (EG) and Material Awareness
(MA), which identify edge and material cues for all visible objects. An
alternative domain-specific module is Lateral Inhibition Module (LIM) [24],
which includes two components: Bidirectional Propagation (BP) and Boundary
Activation (BA). The former one minimizes the impact of neighboring regions,
by isolating irrelevant information and the latter one captures object
boundaries. Both DOAM and LIM have shown promising results in overcoming
object occlusion issues in X-ray scan images.

NMS has also been modified in object detectors for X-ray scan image analysis.
For instance, the framework of[25] is a modified YOLOv4 detector adopting
deformable convolutions [26], the Gradient Harmonizing Mechanism (GHM) loss
[27] and an augmented NMS algorithm combining Soft-NMS [28] with the Distance-
Intersection-over-Union (DIoU) metric. Focusing on real-time performance,
YOLOv5 was modified in [29] using the Stem [30] and CGhost [31] modules,
resulting in a model with reduced number of parameters that still achieves
competitive results in comparison with the baseline method. Finally, the
integrated illicit Object Detection (POD) method [32] for X-ray image analysis
combines a learnable Gabor layer for edge information retrieval, a spatial
attention module for directing focus on low-level features, a Global Context
Feature Extraction (GCFE) module and a Dual Scale Feature Aggregation (DSFA)
module to enhance semantic information from high-level features.

However, to the best of the authors’ knowledge, no object detector devised for
the security domain has attempted to modify one basic building block of most
single-stage detection frameworks: the anchor boxes. This is particularly
important for X-ray screening of luggage or mailed parcels, because better
matching between the anchor boxes and the distribution of object sizes/shapes
in the training dataset leads to better detection performance on test images.
Additionally, despite certain attempts to improve NMS for security
applications, the results remain typically sub-optimal under object
occlusions, which are common in this domain. Thus, this paper proposes a two-
fold improvement of anchor-based, single-stage object detectors for
automatically detecting contraband items in X-ray scan images, contributing
the following two novelties:

  * •

Anchor box optimization by applying Hierarchical Clustering (HC) on the
ground-truth object RoIs of the training set. By clustering the ground-truth
bounding boxes based on their similarity in terms of size, shape, and
position, the resulting clusters can be used to define a natural hierarchy,
with larger clusters representing more general object shapes and smaller
clusters capturing finer details and variations. The resulting hierarchy can
also provide information about the relationships between different object
classes.

  * •

NMS modification to handle occluded object detection and to reduce false
predictions, by computing richer geometrical correlations among candidate RoIs
before final bounding box prediction. This is achieved by incorporating the
Efficient-IoU metric into the Weighted-Cluster NMS method [33].

The remainder of the paper is organized as follows. Section II briefly
presents the specific baseline algorithms which is adopted for implementing
the proposed novelties (YOLOv5, Weighted-Cluster NMS). Section III details the
proposed method, consisting of an anchor box refinement approach and a
modified NMS algorithm. Section IV outlines the experimental evaluation
process, which was conducted on a well-known public dataset, and discusses the
obtained results. Section V concludes the preceding discussion by identifying
the implications of these findings, the limitations of this study and
directions for future research.

##  II Preliminaries

In order to evaluate the proposed two-fold method, YOLOv5 [34] was adopted as
a baseline object detector. The reason behind this choice was solely
practical; in principle, the proposed method can be used to augment any other
variant of the general anchor-based, single-stage object detection framework,
as well.

###  II-A YOLOv5 Architecture

You Only Look Once (YOLO) [16] is a series of fast anchor-based, single-stage
object detectors, where object localization and classification are performed
using a single CNN. This architecture can, however, be divided into a backbone
network, a succeeding neck network and a final prediction head. YOLOv5 [34],
which is an update of YOLOv4 [35], was inspired by EfficientNet [36] and,
thus, can be easily reconfigured for different network complexity profiles.
Out of the common variants (YOLOv5s, YOLOv5m, YOLOv5l, YOLOv5x) the one
employed in this paper is YOLOv5l.

The backbone CNN of YOLOv5 is CSP-Darknet53, a modified version of Darknet53
[37] combined with Cross Stage Partial (CSP). As presented in Fig. 1, the main
convolutional block of CSP-Darknet53 consists of convolutional layers,
residuals and the SiLU activation function, while the final feature maps are
refined using a Spatial Pyramid Pooling-Fast (SPPF) module [38]. The neck
network consists of a Feature Pyramid Network (FPN) [20] and a Path
Aggregation Network (PAN) [39]. These modules repeatedly fuse feature maps
from different scales and depth levels, thus leading to final image
representations, which are simultaneously characterized by accurate spatial
localization details, rich semantics and high invariance regarding object
detection. Finally, the prediction head outputs the candidate detected RoIs
through a set of convolutional operations. Overall YOLOv5 architecture is
presented in Fig. 2.

Figure 1: Main YOLO components.

Figure 2: YOLOv5 overall architecture.

###  II-B Non-Maximum Suppression

Similarly to the majority of object detectors, YOLO generates a large set of
overlapping object proposals, in the form of RoIs in pixel coordinates, along
with the corresponding class labels and confidence scores. Thus, these
candidate RoIs are filtered in a post-processing step based on certain
criteria; this is called Non-Maximum Suppression (NMS). The conventional
Greedy NMS algorithms processes the generated candidate bounding boxes and
their corresponding confidence scores for each input image, sorting RoIs in
descending confidence order. At first, the box with the highest confidence
score is selected and the IoU between itself and all other boxes is
calculated. All significantly overlapping RoIs, with an IoU greater than a
threshold, are removed. This process is repeated until no bounding boxes
remain in the sorted list. The NMS algorithm is presented in Fig. 3.

Figure 3: The Greedy Non-Maximum Suppression algorithm.

##  III Proposed Method

The proposed method is a two-fold improvement of anchor-based, single-stage
object detectors, which is highly suitable for the X-ray security scan image
analysis domain, due to the peculiarities of such images (e.g., heavy
occlusions, cluttered background, etc.).

###  III-A Anchor boxes refinement

Most anchor-based single-stage object detectors utilize reference anchor boxes
of different sizes and aspect ratios, which are placed at various positions
across the input image. The goal of these anchor boxes is to capture the
variation in object shapes and sizes present in the dataset. Typically, they
are predefined (e.g., in the case of YOLOv5, they have been calculated based
on prior knowledge of the sizes, aspect ratios, and distributions of ground-
truth objects in the COCO dataset [40]). In many implementations (e.g.,
YOLOv5) the match between these predefined anchor boxes and the training
dataset is verified before training, by computing the achievable recall rate
if the object detector using these anchors had access to the ground-truth for
all objects in the dataset. If this recall rate is too low, the predefined
anchors are assumed to be unfit and a new set of dataset-specific anchor boxes
is estimated. This is performed via a run of K-Means to group the ground-truth
dataset RoIs into clusters, based on their dimensions in pixel space. The
resulting cluster centers are selected as the new anchor boxes, with
additional optimizations being possible (e.g., in YOLO a genetic algorithm
refines them further). K-Means++ [41] can be utilized instead of classic
K-Means [42].

In the current work, Hierarchical Clustering (HC) is used to obtain anchor
boxes in a dataset-specific manner, according to Algorithm in Fig. 4. The goal
is to generate anchor boxes that both fit the distribution of ground-truth
object sizes/shapes and reflect their arrangement into a natural hierarchy,
aligned with the spatial interrelations between the dataset’s object classes.
Thus, the hierarchy of anchor boxes can provide the detector an explicit
template of the dataset’s semantic structure, as expressed by the spatial
relationships between different object classes. For example, in illicit item
detection, RoIs corresponding to classes such as “knife” and “wrench” will
likely fall under different sub-clusters of a common super-cluster, containing
all small-sized handheld items. Since a selection of anchor boxes that better
match the dataset’s object sizes and shapes is known to lead to better object
detection accuracy [42], it is reasonable to expect further improvements by
obtaining an arrangement of anchors that not only fit the distribution of the
dataset’s object sizes and shapes, but also reflect the natural hierarchy of
the dataset’s semantic classes (at least in terms of RoI shape/size).

Algorithm of Fig. 4 adopts an agglomerative HC method [43] and adapts it to
the anchor box refinement task. Its goal is to hierarchically group the
training dataset’s ground-truth bounding boxes, where each RoI is described as
a 2D feature vector: $[w,h]^{T}$, where $w$/$h$ is the RoI width/height,
respectively. First, the pairwise Euclidean distances between all RoIs in the
entire training dataset are computed and a linkage matrix is constructed using
Ward’s minimum variance [44]. Then, bounding boxes are assigned to clusters
using the maximum cluster criterion and the mean of each cluster is calculated
to obtain a new set of corresponding anchor boxes (one per cluster). The total
number of target clusters is set to 9, according to experimental evaluation.
HC generates a tree-like arrangement of clusters and sub-clusters. The leaves
and the root of the clustering tree are not included in the final set of
formed anchor boxes.

Figure 4: Anchor boxes refinement using Hierarchical Clustering.

###  III-B Modified Non-Maximum Suppression

The default Greedy NMS method suffers in the presence of occlusions and gives
rise to false positives, triggering various improvements that have been
proposed over the years. For instance, Soft-NMS [28] modifies the candidate
RoI scores by a Gaussian decay based on the degree of overlap, instead of
directly setting the score of all overlapping bounding boxes to zero. This
generates more accurate RoIs, even if they are occluded by other objects.
Weighted-NMS [45] utilizes the weighted combination of scores and IoU values
to define the merged coordinates of the predicted bounding boxes; the result
is higher accuracy at the expense of increased time complexity, due to the
number of iterations. To mitigate this, Weighted-Cluster NMS (WC-NMS) [33] has
been developed: WC-NMS groups the detected candidate bounding boxes according
to the IoU values and then selects the final RoIs according to the maximum
score within each group. Implementation-wise this is done with SIMD
parallelism and by exploiting cache locality, thanks to formulating the
process as a series of matrix operations instead of naive iterative loops,
resulting in the fast NMS Algorithm presented in Fig. 5. Thus, suppression is
implemented by calculating the so-called IoU matrix. The latter is a symmetric
matrix M, where $m_{i,j}$ is the IoU between the $i$-th and the $j$-th
candidate RoI. Exploiting the symmetry of M, Algorithm in Fig. 5 retains only
its upper triangular part.

Figure 5: EIoU Weighted-Cluster NMS algorithm.

Considering the importance of an efficient NMS method in the X-ray security
image analysis domain, due to the densely packed nature of typical luggage and
mailed parcels, this paper adopts WC-NMS and improves it, by employing the
Efficient-IoU (E-IoU) [46] as the overlap metric. E-IoU is an improvement of
Complete-IoU (C-IoU) [33] and captures richer geometrical information about
overlapping candidate bounding boxes, taking into account their overlapping
area, their distance and their aspect ratios. Eqs. (1)-(2) define the
overlapping criterion according to the so-called E-IoU matrix.

$\mathbf{X}=\mathbf{M}_{IoU}-\mathbf{R}_{EIoU},$ (1)

where $\textbf{M}_{IoU}$ is the IoU matrix, while $\textbf{R}_{EIoU}$ derives
from the E-IoU loss [46]. Both matrices are calculated using the predicted
candidate RoIs (in top-left, bottom-right pixel coordinates format), as
follows:

$\mathbf{R}_{EIoU}=\frac{\mathbf{D}_{centers}}{(\mathbf{W}^{c})^{2}+(\mathbf{H}^{c})^{2}}+\frac{\mathbf{D}^{w}}{(\mathbf{W}^{c})^{2}}+\frac{\mathbf{D}^{h}}{(\mathbf{H}^{c})^{2}},$
(2)

where , $\textbf{D}_{centers}\in\mathbb{R}^{N\times N}$ is a matrix containing
the pairwise Euclidean distances between all $N$ candidate RoIs, with
$d_{i,j}$ being the Euclidean distance between the centers of the $i$-th and
the $j$-th predicted bounding boxes. $\textbf{W}^{c}\in\mathbb{R}^{N\times N}$
and $\textbf{H}^{c}\in\mathbb{R}^{N\times N}$ contain the width and the
height, respectively, of the smallest enclosing box covering each pair of
candidate ROIs. Thus, $w^{c}_{i,j}$ and $h^{c}_{i,j}$ are the width and height
of the smallest bounding box that can contain both the $i$-th and the $j$-th
RoI. $\textbf{D}^{w}\in\mathbb{R}^{N\times N}$ contains the pairwise Euclidean
distances between the widths, whereas $\textbf{D}^{h}\in\mathbb{R}^{N\times
N}$ contains the pairwise Euclidean distances between the heights of all $N$
candidate RoIs. Similarly to X, all of the matrices mentioned above are
symmetric and only their upper triangular part is important for computations
and the final result is calculated using element-wise division between the
nominator and the denominator, at each term of the sum in Eq. (2). Power
operations in the denominators are element-wise as well. To sum up, the first
term of in Eq. (2) captures information about the distance between candidate
regions, while the remainder about aspect ratio.

Grouping the candidate RoIs based on their E-IoU values means that rich
geometrical information, regarding the spatial interrelations and the
arrangement of the detected bounding boxes in 2D pixel space, are inherently
considered for grouping them more efficiently. In comparison with D-IoU, which
only describes the overlapping area and the distance between candidate boxes,
E-IoU captures information about the overlapping area, distance and aspect
ratios between the compared RoIs. Although C-IoU also takes into account such
geometrical factors, it is computed by estimating the simple difference in
aspect ratio between the compared bounding boxes. Such a naive approach may
not accurately reflect the actual relationship between the RoI shapes [46]. In
order to handle these issues, the second and the third term in Eq. (2) capture
more sufficiently the similarity between the compared bounding boxes, in terms
of their aspect ratio.

##  IV Experimental Evaluation

This Section overviews the experimental setups used for evaluating the
proposed method and discusses the evaluation results. As previously described,
YOLOv5 was selected as the baseline to be improved using the proposed method.

###  IV-A Experimental Dataset

The proposed method was implemented and tested using SIXray [7], a publicly
available X-ray security dataset consisting of 1,059,231 X-ray images from
subway stations. The 6 classes of illicit objects contained in these images
are “gun”, “knife”, “wrench”, “pliers” and “scissors”. Additionally, a
“negative” class includes all images without any illicit item. Three different
dataset subsets are typically utilized in different experimental setups,
namely SIXray10, SIXray100 and SIXray1000, where the number indicates the
ratio of negative against positive samples. SIXray contains ground-truth
whole-image class label annotations manually set by human security inspectors,
while their ground-truth object RoIs/bounding boxes are available only for the
test set. This paper uses the revised object detection annotations for the
training subset provided by [47]. Despite the fact that only images containing
at least one contraband item were utilized, official training-test set split
was adopted.

###  IV-B Evaluation Metrics

The effectiveness of the proposed method is measured using the precision,
recall and mean Average Precision (mAP) metrics. In object detection tasks,
IoU is used to measure the overlap between the predicted and the corresponding
ground-truth RoI. In addition, a threshold value was defined in order to
decide whether the prediction is actually correct. True Positives (TP), False
Positives (FP), and False Negatives (FN) depend on the IoU, the predicted
label and the ground-truth label. These elementary metrics are utilized to
calculate Precision and Recall:

$Precision=\frac{TP}{TP+FP}.$ (3)

$Recall=\frac{TP}{TP+FN}.$ (4)

The Precision-Recall (PR) curve depicts the trade-off between precision and
recall for different discrimination thresholds. Average Precision (AP) is the
area under the PR curve and its range is between 0 to 1. AP is defined as:

$AP=\int_{0}^{1}p(r)\,dr.$ (5)

mAP is calculated as the mean of AP over all classes:

$mAP=\frac{1}{N}\sum_{i}^{N}AP_{i}.$ (6)

###  IV-C Experimental Evaluation

Evaluation of all competing methods in the SIXray dataset was conducted using
the mAP metric at a 0.5 IoU threshold and the average mAP value at a range of
different IoU thresholds. Comparisons were made against the baseline detector
implementation before it was augmented with the proposed method (default
YOLOv5 with anchor boxes obtained from the COCO dataset and Greedy NMS), as
well as with variations using K-Means and K-Means++ clustering for obtaining
the anchor boxes, or employing basic (IoU-based) WC-NMS. Additionally, a
published YOLOv5 result on SIXray is included for completeness.

Table I summarizes the accuracy of the baseline method, which achieved a
precision of 92.1%, a recall of 82.3%, a mAP of 87.6% and an average mAP
across different IoU thresholds (from 0.5 to 0.95) of 72.3%. Table II compares
the proposed method against this baseline and against competing approaches
based on YOLOv5. The first method is a competing one published in [32], using
YOLOv5 with default anchor boxes, conventional Greedy NMS and a different
mini-batch size during training. The next two approaches use Greedy NMS, with
one of them employing K-Means-derived and one employing K-Means++-derived
anchor boxes. As it can be seen, the proposed method outperforms all other
approaches in terms of mAP.

Additionally, Table III presents an ablation study of the proposed method. The
first three variants adopt HC-derived anchor boxes in combination with three
different IoU metrics for WC-NMS, namely IoU, D-IoU and C-IoU, respectively.
Evidently, D-IoU outperforms the other two metrics in terms of mAP. The last
line of Table II demonstrates the results of the full proposed method,
integrating both HC-based anchors and advanced E-IoU-based WC-NMS into YOLOv5,
which outperforms all other (partial) variants in terms of mAP.

TABLE I: Baseline YOLOv5 results across all classes. | Precision | Recall | mAP | mAP50-95
---|---|---|---|---
Overall | 0.921 | 0.823 | 0.876 | 0.723
Class |  |  | AP | AP50-95
Gun | 0.978 | 0.916 | 0.944 | 0.882
Knife | 0.925 | 0.758 | 0.813 | 0.659
Wrench | 0.877 | 0.768 | 0.835 | 0.659
Pliers | 0.919 | 0.845 | 0.916 | 0.728
Scissors | 0.908 | 0.828 | 0.874 | 0.687
TABLE II: Comparative Evaluation. Method | mAP | mAP50-95
---|---|---
Baseline | 87.6 | 72.3
YOLOv5 baseline of [32] | 86.7 | -
K-Means anchors + default NMS | 88 | 73
K-Means++ anchors + default NMS | 88.3 | 72.9
HC anchors + E-IoU WC-NMS (proposed) | 89.9 | 75.7
TABLE III: Ablation study. Method | mAP | mAP50-95
---|---|---
HC anchors + IoU WC-NMS | 89.2 | 75
HC anchors + D-IoU WC-NMS | 89.7 | 75.4
HC anchors + C-IoU WC-NMS | 89.5 | 75.2
HC anchors + E-IoU WC-NMS (full proposed) | 89.9 | 75.7

Table IV presents the complete evaluation of the proposed method across all
classes, using HC-derived anchor boxes and the E-IoU-based WC-NMS. It
outperforms the default YOLOv5-large by 2.3% in terms of mAP. Its mAP is 89.9%
at a 0.5 IoU threshold, while the average mAP across a range of different IoU
thresholds is 75.7%. The proposed method significantly reduces the number of
false predictions and is more accurate in detecting contraband items,
especially in cases of occluded object detection. In Figure 6, the precision-
recall curve of the proposed framework is presented, which highlights the
performance of the model in the desired task. The curve shows a high precision
score for low recall values, indicating that the model is very selective in
its predictions. However, as recall increases, the precision score decreases,
suggesting that the model struggles with correctly classifying some samples.
However, the Area Under the Curve (AUC) value indicates that the model
performs robust predictions. The above findings suggest that the model
exhibits potentials for use in specific applications where high precision is
critical, such as contraband detection. Finally, our model was deployed in the
test subset of SIXray dataset and the predictions are presented in Fig. 7.
Notably, neither HC-derived anchor boxes nor E-IoU-based WC-NMS have been
proposed/investigated before for object detection.

TABLE IV: Accuracy of the proposed method across all classes. Class | Precision | Recall | mAP | mAP 50-95
---|---|---|---|---
Overall | 0.949 | 0.837 | 0.899 | 0.757
|  |  | AP | AP 50-95
Gun | 0.977 | 0.945 | 0.971 | 0.917
Knife | 0.954 | 0.789 | 0.841 | 0.692
Wrench | 0.904 | 0.781 | 0.86 | 0.695
Pliers | 0.955 | 0.833 | 0.924 | 0.75
Scissors | 0.956 | 0.84 | 0.899 | 0.73

Figure 6: Precision-Recall curve of the proposed method.

Figure 7: Predictions on the SIXray test subset.

##  V Conclusions

The large volume and high throughput of passengers or mailed parcels during
rush hours, in airports, subways or post/customs offices, make the automated
detection of contraband items in X-ray images a Big Data analysis task that is
critical for public safety. This paper proposed a novel approach to improve
the performance of single-stage, anchor-based object detectors in the X-ray
domain. It incorporated two complementary improvements: dataset-specific
hierarchical clustering of ground-truth training RoIs, so that the derived
anchor boxes better match the distribution and semantic hierarchy of object
sizes/shapes, and a modification of an efficient NMS algorithm, so as to
better handle occluded objects and to reduce false predictions. According to a
thorough experimental evaluation on a relevant public dataset, the proposed
method outperforms both the baseline and various competing approaches.

Future research will focus on addressing other limitations of existing
methods, such as low generalization ability under domain shifts (e.g., if the
detector has been trained with X-rays from one type of scanner and is deployed
in an airport using a different scanner), as well as integration of the
proposed approach to more recent object detectors (e.g., YOLOv7).

## Acknowledgment

The research leading to these results has received funding from the European
Union’s Horizon Europe research and innovation programme under grant agreement
No 101073876 (Ceasefire). This publication reflects only the authors views.
The European Union is not liable for any use that may be made of the
information contained therein.

## References

  * [1] D. Mery, D. Saavedra, and M. Prasad, “X-ray baggage inspection with computer vision: A survey,” _IEEE Access_ , vol. 8, pp. 145 620–145 633, 2020.
  * [2] S. Akcay and T. Breckon, “Towards automatic threat detection: A survey of advances of deep learning within X-ray security imaging,” _Pattern Recognition_ , vol. 122, p. 108245, 2022.
  * [3] S. Thermos, G. T. Papadopoulos, P. Daras, and G. Potamianos, “Deep affordance-grounded sensorimotor object recognition,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , July 2017.
  * [4] C. Symeonidis, I. Mademlis, N. Nikolaidis, and I. Pitas, “Improving neural Non-Maximum Suppression for object detection by exploiting interest-point detectors,” in _Proceedings of the IEEE International Workshop on Machine Learning for Signal Processing (MLSP)_ , 2019.
  * [5] C. Symeonidis, I. Mademlis, I. Pitas, and N. Nikolaidis, “Neural attention-driven Non-Maximum Suppression for person detection,” _IEEE Transactions on Image Processing_ , 2023, accepted for publication.
  * [6] S. Akçay, M. E. Kundegorski, M. Devereux, and T. P. Breckon, “Transfer learning using convolutional neural networks for object classification within X-ray baggage security imagery,” in _Proceedings of the IEEE International Conference on Image Processing (ICIP)_ , 2016.
  * [7] C. Miao, L. Xie, F. Wan, C. Su, H. Liu, J. Jiao, and Q. Ye, “SIXray: A large-scale security inspection X-ray benchmark for prohibited item discovery in overlapping images,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * [8] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016.
  * [9] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, “Rethinking the inception architecture for computer vision,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016.
  * [10] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, “Densely connected convolutional networks,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017.
  * [11] Y. F. A. Gaus, N. Bhowmik, S. Akcay, and T. Breckon, “Evaluating the transferability and adversarial discrimination of convolutional neural networks for threat object detection and classification within X-ray security imagery,” in _Proceedings of the IEEE International Conference On Machine Learning And Applications (ICMLA)_ , 2019.
  * [12] S. Ren, K. He, R. Girshick, and J. Sun, “Faster R-CNN: Towards real-time object detection with region proposal networks,” _Proceedings of the Advances in Neural Information Processing Systems (NIPS)_ , 2015.
  * [13] K. He, G. Gkioxari, P. Dollár, and R. Girshick, “Mask r-cnn,” in _2017 IEEE International Conference on Computer Vision (ICCV)_ , 2017, pp. 2980–2988.
  * [14] T.-Y. Lin, P. Goyal, R. Girshick, K. He, and P. Dollár, “Focal loss for dense object detection,” in _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_ , 2017.
  * [15] W. Liu, D. Anguelov, D. Erhan, C. Szegedy, S. Reed, C.-Y. Fu, and A. C. Berg, “Ssd: Single shot multibox detector,” in _Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11–14, 2016, Proceedings, Part I 14_. Springer, 2016, pp. 21–37.
  * [16] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi, “You only look once: Unified, real-time object detection,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016.
  * [17] T. Hassan, S. Akcay, M. Bennamoun, S. Khan, and N. Werghi, “Cascaded structure tensor framework for robust identification of heavily occluded baggage items from X-ray scans,” _arXiv preprint arXiv:2004.06780_ , 2020.
  * [18] Y. Ren, H. Zhang, H. Sun, G. Ma, J. Ren, and J. Yang, “LightRay: Lightweight network for prohibited items detection in X-ray images during security inspection,” _Computers and Electrical Engineering_ , vol. 103, p. 108283, 2022.
  * [19] A. Howard, M. Sandler, G. Chu, L.-C. Chen, B. Chen, M. Tan, W. Wang, Y. Zhu, R. Pang, V. Vasudevan _et al._ , “Searching for MobileNetv3,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , 2019.
  * [20] T.-Y. Lin, P. Dollár, R. Girshick, K. He, B. Hariharan, and S. Belongie, “Feature pyramid networks for object detection,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 2117–2125.
  * [21] S. Woo, J. Park, J.-Y. Lee, and I. S. Kweon, “Cbam: Convolutional block attention module,” in _Proceedings of the European conference on computer vision (ECCV)_ , 2018, pp. 3–19.
  * [22] F. Shao, J. Liu, P. Wu, Z. Yang, and Z. Wu, “Exploiting foreground and background separation for prohibited item detection in overlapping X-ray images,” _Pattern Recognition_ , vol. 122, p. 108261, 2022.
  * [23] Y. Wei, R. Tao, Z. Wu, Y. Ma, L. Zhang, and X. Liu, “Occluded prohibited items detection: An X-ray security inspection benchmark and de-occlusion attention module,” in _Proceedings of the ACM International Conference on Multimedia (ACM MM)_ , 2020.
  * [24] R. Tao, Y. Wei, X. Jiang, H. Li, H. Qin, J. Wang, Y. Ma, L. Zhang, and X. Liu, “Towards real-world X-ray security inspection: A high-quality benchmark and lateral inhibition module for prohibited items detection,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , 2021.
  * [25] C. Zhou, H. Xu, B. Yi, W. Yu, and C. Zhao, “X-ray security inspection image detection algorithm based on improved YOLOv4,” in _Proceedings of the IEEE Eurasia Conference on IOT, Communication and Engineering (ECICE)_ , 2021.
  * [26] J. Dai, H. Qi, Y. Xiong, Y. Li, G. Zhang, H. Hu, and Y. Wei, “Deformable convolutional networks,” in _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_ , 2017.
  * [27] B. Li, Y. Liu, and X. Wang, “Gradient harmonized single-stage detector,” in _Proceedings of the AAAI Conference on Artificial Intelligence (AAAI)_ , 2019\.
  * [28] N. Bodla, B. Singh, R. Chellappa, and L. S. Davis, “Soft-NMS: improving object detection with one line of code,” in _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_ , 2017.
  * [29] B. Song, R. Li, X. Pan, X. Liu, and Y. Xu, “Improved YOLOv5 detection algorithm of contraband in x-ray security inspection image,” in _Proceedings of the International Conference on Pattern Recognition and Artificial Intelligence (PRAI)_. IEEE, 2022\.
  * [30] R. J. Wang, X. Li, and C. X. Ling, “Pelee: A real-time object detection system on mobile devices,” _Proceedings of Advances in Neural Information Processing Systems (NIPS)_ , 2018.
  * [31] K. Han, Y. Wang, Q. Tian, J. Guo, C. Xu, and C. Xu, “Ghostnet: More features from cheap operations,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020.
  * [32] C. Ma, L. Zhuo, J. Li, Y. Zhang, and J. Zhang, “Occluded prohibited object detection in X-ray images with global context-aware multi-scale feature aggregation,” _Neurocomputing_ , vol. 519, pp. 1–16, 2023.
  * [33] Z. Zheng, P. Wang, D. Ren, W. Liu, R. Ye, Q. Hu, and W. Zuo, “Enhancing geometric factors in model learning and inference for object detection and instance segmentation,” 2021.
  * [34] G. Jocher, “YOLOv5 by Ultralytics.” [Online]. Available: https://github.com/ultralytics/yolov5
  * [35] A. Bochkovskiy, C.-Y. Wang, and H.-Y. M. Liao, “YOLOv4: Optimal speed and accuracy of object detection,” _arXiv preprint arXiv:2004.10934_ , 2020.
  * [36] M. Tan and Q. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _International conference on machine learning_. PMLR, 2019, pp. 6105–6114.
  * [37] J. Redmon and A. Farhadi, “YOLOv3: An incremental improvement,” _arXiv preprint arXiv:1804.02767_ , 2018.
  * [38] K. He, X. Zhang, S. Ren, and J. Sun, “Spatial pyramid pooling in deep convolutional networks for visual recognition,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 37, no. 9, pp. 1904–1916, 2015\.
  * [39] S. Liu, L. Qi, H. Qin, J. Shi, and J. Jia, “Path aggregation network for instance segmentation,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 8759–8768.
  * [40] T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick, “Microsoft COCO: Common objects in context,” in _Proceedings of the European Conference on Computer Vision (ECCV)_. Springer, 2014.
  * [41] D. Arthur and S. Vassilvitskii, “K-means++: The advantages of careful seeding,” Stanford, Tech. Rep., 2006.
  * [42] X. Luo, Y. Wu, and F. Wang, “Target detection method of UAV aerial imagery based on improved YOLOv5,” _Remote Sensing_ , vol. 14, no. 19, p. 5063, 2022.
  * [43] S. Landau, M. Leese, D. Stahl, and B. S. Everitt, _Cluster analysis_. John Wiley & Sons, 2011.
  * [44] J. H. Ward Jr, “Hierarchical grouping to optimize an objective function,” _Journal of the American statistical association_ , vol. 58, no. 301, pp. 236–244, 1963.
  * [45] H. Zhou, Z. Li, C. Ning, and J. Tang, “CAD: Scale-invariant framework for real-time object detection,” in _Proceedings of the IEEE International Conference on Computer Vision Workshops (ICCVW)_ , 2017.
  * [46] Y.-F. Zhang, W. Ren, Z. Zhang, Z. Jia, L. Wang, and T. Tan, “Focal and efficient IOU loss for accurate bounding box regression,” _Neurocomputing_ , vol. 506, pp. 146–157, 2022.
  * [47] H. D. Nguyen, R. Cai, H. Zhao, A. C. Kot, and B. Wen, “Towards more efficient security inspection via deep learning: A task-driven X-ray image cropping scheme,” _Micromachines_ , vol. 13, no. 4, p. 565, 2022.


handled by Lemma 3.1. In the unweighted case, the size of the independent set
for each type-$d$ subproblem lies in the range from 1 to $B:=O(b^{d})$ (after
removing empty subproblems). In the weighted case, the weight of the
independent set for the $i$-th type-$d$ subproblem lies between $w_{i}$ and
$Bw_{i}$ where $w_{i}$ is the largest rectangle weight in the $i$-th
subproblem. We can approximate the sum to within a constant factor by summing
over $O(B)$ terms, via basic random sampling in the unweighted case or
importance sampling in the weighted case (see Fact A.1), with error
probability at most $1/4$ (which can be lowered by by repeating for
logarithmically many trials and returning the median). Thus, the number of
calls to Lemma 3.1 is $\widetilde{O}(b^{d})$. The total running time is now
$\widetilde{O}(n+b^{O(1)})$, which is $\widetilde{O}(n)$ by setting
$b=n^{\delta}$ for a sufficiently small constant $\delta$. ∎

######  Theorem A.5.

Let $d$ and $c$ be constants. Given a set $S$ of $n$ weighted objects in
$\mathbb{R}^{d}$ of constant description complexity from a $c$-fat collection
${\cal C}$, we can compute an $O(1)$-approximation to the weight of the
maximum-weight independent set for $S$ in $\widetilde{O}(n)$ time w.h.p. by a
static Monte Carlo randomized algorithm.

###### Proof.

This follows by a similar modification to the proof of Theorem 3.6 with
importance sampling. ∎

##  Appendix B MVC: Fast Static Algorithms

###  B.1 Disks in $\mathbb{R}^{2}$

In this subsection, we give a static, near-linear-time,
$(1+O(\varepsilon))$-approximation algorithm $\mathcal{A}$ for MVC for disks
under the promise that $n=O(\mbox{\rm OPT})$. As we are aiming for a
$(1+O(\varepsilon))$-factor approximation, we may tolerate an additive error
of $O(\varepsilon n)$ because of the promise. MVC with $O(\varepsilon n)$
additive error reduces to MIS with $O(\varepsilon n)$ additive error, by
complementing the solution.

PTASs are known for MIS for disks, but allowing $O(\varepsilon n)$ additive
error, there are actually EPTASs that run in near linear time. For example, we
can adapt an approach by Chan [Cha03] based on divide-and-conquer via
separators. Specifically, we can use the following variant of Smith and
Wormald’s geometric separator theorem [SW98]:

######  Lemma B.1 (Smith and Wormald’s separator).

Given $n$ fat objects in a constant dimension $d$, there is an axis-aligned
hypercube $B$, such that the number of objects inside $B$ and the number of
objects outside $B$ are both at most $(1-\beta)n$ for some constant $\beta>0$
(dependent on $d$), and the objects intersecting $\partial B$ can be stabbed
by $O(n^{1-1/d})$ points. Furthermore, $B$ can be constructed in $O(n)$ time.

###### Proof.

(The following is a modification of a proof in [CH24], which is based on Smith
and Wormald’s original proof [SW98].)

For each object, first pick an arbitrary “center” point inside the object. Let
$B_{0}$ be the smallest hypercube that contains at least $\frac{n}{2^{d}+1}$
center points. Let $r$ be the side length of $B_{0}$. Let $h$ be a parameter
to be set later. For each
$t\in\\{\frac{1}{h},\frac{2}{h},\ldots,\frac{h-1}{h}\\}$, let $B_{t}$ be the
hypercube with a scaled copy of $B_{0}$ with the same center and side length
$(1+t)r$. Since $B_{t}$ can be covered by $2^{d}$ quadtree boxes of side
length $<r$, the number of center points inside $B_{t}$ is at most
$\frac{2^{d}n}{2^{d}+1}$.

An object of diameter $\leq r/h$ intersects $\partial B_{t}$ for at most
$O(1)$ choices of $t$. Thus, we can find a value $t$ in $O(n)$ time such that
$\partial B_{t}$ intersects $O(n/h)$ objects of diameter $\leq r/h$. On the
other hand, the objects of diameter $>r/h$ intersecting $\partial B_{t}$ can
all be stabbed by $O(h^{d-1})$ points, because of fatness. Set $h=n^{1/d}$,
and $B_{t}$ satisfies the desired properties with
$\beta=\frac{2^{d}}{2^{d}+1}$.

One remaining issue is how to find $B_{0}$ quickly. Let $C$ be a sufficiently
large constant. Form a $C\times\cdots\times C$ grid, so that the number of
center points between any two consecutive parallel grid hyperplanes is $n/C$.
This takes $O(Cn)$ time by a linear-time selection algorithm. Round each
center point to its nearest grid point. The rounded point set is a multiset
with $O(C^{d})$ distinct elements. Redefine $B_{0}$ as the smallest hypercube
that contains at most $\frac{n}{2^{d}+1}$ rounded center points (multiplicity
included), which can now be computed in $C^{O(1)}$ time. For any box, the
number of center points inside the box changes by at most $O(n/C)$ after
rounding. So, we just need to increase $\beta$ by $O(1/C)$. ∎

The MIS algorithm is simple: we just recursively compute an independent set
for the disks inside $B$ and an independent set for the disks outside $B$, and
take the union of the two sets. If $n$ is below a constant $b$, we solve the
problem exactly by brute force in $2^{O(C)}$ time. (This is analogous to
Lipton and Tarjan’s original EPTAS for independent set for planar graphs
[LT80].)

Since the optimal independent set can have at most $O(\sqrt{n})$ disks
intersecting $\partial B$, the total additive error satisfies a recurrence of
the form

$E(n)\leq\left\\{\begin{array}[]{ll}\max_{n_{1},n_{2}\leq(1-\beta)n:\
n_{1}+n_{2}\leq n}(E(n_{1})+E(n_{2})+O(\sqrt{n}))&\mbox{if $n\geq b$}\\\
0&\mbox{if $n<b$.}\end{array}\right.$

This yields $E(n)=O(n/\sqrt{b})$. We set $b=1/\varepsilon^{2}$. The running
time of the resulting static algorithm $\mathcal{A}$ is
$T(n)=\widetilde{O}(2^{O(1/\varepsilon^{2})}n)$.

###  B.2 Rectangles in $\mathbb{R}^{2}$

In this subsection, we give a static, near-linear-time,
$(\frac{3}{2}+O(\varepsilon))$-approximation algorithm for MVC for rectangles,
by adapting Bar-Yehuda, Hermelin, and Rawitz’s previous, polynomial-time,
$(\frac{3}{2}+O(\varepsilon))$-approximation algorithm [BHR11].

#### Triangle-free case.

We first start with the case when the intersection graph is triangle-free, or
equivalently, when the maximum depth of the rectangles is at most 2. (The
depth of a point $q$ among a set of rectangles refers to the number of
rectangles containing $q$; the maximum depth refers to the maximum over all
$q\in\mathbb{R}^{2}$.)

We will design a static algorithm $\mathcal{A}$ for MVC for rectangles, under
the promise that $n\leq(2+O(\varepsilon))\,\mbox{\rm OPT}$.

We say that a rectangle $s$ dominates another rectangle $s^{\prime}$ if
$\partial s$ intersects $\partial s^{\prime}$ four times, with $s$ having
larger height than $s^{\prime}$. Bar-Yehuda, Hermelin, and Rawitz’s algorithm
proceeds as follows:

  1. 1.

Decompose the input set $S$ into 2 subsets $S_{1}$ and $S_{2}$, such that
there are no dominating pairs within each subset $S_{i}$.

The existence of such a decomposition follows easily from Dilworth’s theorem,
but for an explicit construction, we can just define $S_{1}$ to contain all
rectangles in $S$ that are not dominated by any other rectangle in $S$, and
define $S_{2}$ to be $S\setminus S_{1}$. Correctness is easy to see (since the
maximum depth is 2).

We can compute $S_{1}$ (and thus $S_{2}$) by performing $O(n)$ orthogonal
range queries, after lifting each rectangle to a point in $\mathbb{R}^{4}$. In
fact, computing $S_{1}$ corresponds to the maxima problem in $\mathbb{R}^{4}$,
for which $O(n\log n)$-time algorithms are known [CLP11].

  2. 2.

For each $i\in\\{1,2\\}$, compute a vertex cover $X_{i}$ of $S_{i}$ which
approximates the minimum with additive error at most $\varepsilon n$.

Since we tolerate $\varepsilon n$ additive error, it suffices to approximate
the MIS of $S_{i}$ with $\varepsilon n$ additive error.

Because $S_{i}$ has no dominating pairs, it forms a pseudo-disk arrangement
(each pair’s boundaries intersect at most twice). It is known [BHR11, KLPS86]
that the intersection graph of any set of pseudo-disks that have maximum depth
2 and have no containment pair is in fact planar. (We can easily eliminate
containment pair, by removing rectangles that contain another rectangle; and
we can detect such rectangles again by orthogonal range queries.) So, we can
use a known near-linear-time EPTAS for MIS for planar graphs; for example,
Lipton and Tarjan’s divide-and-conquer algorithm via planar-graph separators
runs in $\widetilde{O}(2^{O(1/\varepsilon^{2})}n)$ time [LT80].

  3. 3.

Return $X$, the smaller of the two sets: $X_{1}\cup S_{2}$ and $X_{2}\cup
S_{1}$.

Let $X^{*}$ be the minimum vertex cover. Then

$\displaystyle|X|$ $\displaystyle\leq$
$\displaystyle\tfrac{1}{2}(|X_{1}|+|S_{2}|+|X_{2}|+|S_{1}|)$
$\displaystyle\leq$ $\displaystyle\tfrac{1}{2}(|X^{*}\cap S_{1}|+|X^{*}\cap
S_{2}|+|S_{1}|+|S_{2}|)+\varepsilon n\ =\ \tfrac{1}{2}(|X^{*}|+n)+\varepsilon
n\ \leq\ (\tfrac{3}{2}+O(\varepsilon))|X^{*}|$

under the assumption that $n\leq(2+O(\varepsilon))|X^{*}|$. The running time
is $T(n)=\widetilde{O}(2^{O(1/\varepsilon^{2})}n)$.

At this point, if we apply Theorem 4.4 with $c=\frac{3}{2}+O(\varepsilon)$,
$\gamma=\varepsilon$, and $\delta=\varepsilon^{3}$ (and
$\tau_{0}(n)=O(\log^{O(1)}n)$ as noted in Section 4.4), we obtain:

######  Corollary B.2.

There is a dynamic data structure for $n$ rectangles in $\mathbb{R}^{2}$ that
maintains a $(\frac{3}{2}+O(\varepsilon))$-approximation of the minimum vertex
cover of the intersection graph, under insertions and deletions, in
$O(2^{O(1/\varepsilon^{2})}\log^{O(1)}n)$ amortized time, under the assumption
that the intersection graph is triangle-free.

#### General case.

We now design our final static algorithm $\mathcal{A}$ for MVC for rectangles
that avoids the triangle-free assumption, again by building on Bar-Yehuda,
Hermelin, and Rawitz’s approach:

  1. 1.

Remove vertex-disjoint triangles $T_{1},\ldots,T_{\ell}$ so that the remaining
set $S^{\prime}:=S\setminus(T_{1}\cup\cdots\cup T_{\ell})$ is triangle-free.

We can implement this step greedily in polynomial time, but a faster approach
is via a plane sweep. Namely, we modify the standard sweep-line algorithm to
computing the maximum depth of $n$ rectangles in $\mathbb{R}^{2}$ (similar to
Klee’s measure problem) [PS85]. The algorithm uses a data structure for a 1D
problem: maintaining the maximum depth of intervals in $\mathbb{R}^{1}$,
subject to insertions and deletions of intervals. Simple modification of
standard search trees achieve $O(\log n)$ time per insertion and deletion. As
we sweep a vertical line $\ell$ from left to right, we maintain the maximum
depth of the $y$-intervals of the rectangles intersected by $\ell$. When
$\ell$ hits the left side of a rectangle, we insert an interval. When $\ell$
hits the right side of a rectangle, we delete an interval. When the maximum
depth becomes 3, we remove the 3 intervals containing the maximum-depth point,
which correspond to a triangle in the intersection graph, and continue the
sweep. The total time is $O(n\log n)$.

  2. 2.

Compute a vertex cover $X^{\prime}$ of $S^{\prime}$ which is a
$(\frac{3}{2}+O(\varepsilon))$-approximation to the minimum. This can be done
by Corollary B.2 in $\widetilde{O}(2^{O(1/\varepsilon^{2})}n)$ time (we
actually don’t need the full power of a dynamic data structure, since we just
want a static algorithm for this step).

  3. 3.

Return $X=X^{\prime}\cup T$ with $T:=T_{1}\cup\cdots\cup T_{\ell}$.

Let $X^{*}$ be the minimum vertex cover. The key observation is that $X^{*}$
must contain at least 2 of the 3 vertices in each triangle $T_{i}$, and so
$|X^{*}\cap T|\geq\frac{2}{3}|T|$. Thus,

$|X|\ =\ |X^{\prime}|+|T|\ \leq\ (\tfrac{3}{2}+O(\varepsilon))|X^{*}\cap
S^{\prime}|+\tfrac{3}{2}|X^{*}\cap T|\ \leq\
(\tfrac{3}{2}+O(\varepsilon))|X^{*}|.$

The running time of our final static algorithm $\mathcal{A}$ is
$T(n)=\widetilde{O}(2^{O(1/\varepsilon^{2})}n)$.


# Bridging disconnected networks of first and second lines of biologic
therapies in rheumatoid arthritis with registry data: Bayesian evidence
synthesis with target trial emulation

Sylwia Bujkiewicz  Biostatistics Research Group, Department of Health
Sciences, University of Leicester, University Road, Leicester, LE1 7RH, UK
<EMAIL_ADDRESS>Janharpreet Singh  Biostatistics Research Group,
Department of Health Sciences, University of Leicester, University Road,
Leicester, LE1 7RH, UK  Lorna Wheaton  Biostatistics Research Group,
Department of Health Sciences, University of Leicester, University Road,
Leicester, LE1 7RH, UK  David Jenkins  Biostatistics Research Group,
Department of Health Sciences, University of Leicester, University Road,
Leicester, LE1 7RH, UK  Centre for Health Informatics, Division of
Informatics, Imaging and Data Science, Faculty of Biology, Medicine and
Health, University of Manchester, Manchester, UK  Reynaldo Martina
Biostatistics Research Group, Department of Health Sciences, University of
Leicester, University Road, Leicester, LE1 7RH, UK  Kimme Hyrich  NIHR
Manchester Biomedical Research Centre, Manchester University NHS Foundation
Trust, Manchester Academic Health Sciences Centre, Manchester, UK  Versus
Arthritis Centre for Epidemiology, Centre for Musculoskeletal Research, The
University of Manchester, Manchester, M13 9PL, UK  Keith R. Abrams
Biostatistics Research Group, Department of Health Sciences, University of
Leicester, University Road, Leicester, LE1 7RH, UK  Department of Statistics,
University of Warwick, , Coventry, CV4 7AL, UK  Centre for Health Economics,
University of York, York, YO10 5DD, UK

###### Abstract

Objective: We aim to utilise real world data in evidence synthesis to optimise
an evidence base for the effectiveness of biologic therapies in rheumatoid
arthritis in order to allow for evidence on first-line therapies to inform
second-line effectiveness estimates.
Study design and setting: We use data from the British Society for
Rheumatology Biologics Register for Rheumatoid Arthritis (BSRBR-RA) to
supplement RCT evidence obtained from the literature, by emulating target
trials of treatment sequences to estimate treatment effects in each line of
therapy. Treatment effects estimates from the target trials inform a bivariate
network meta-analysis (NMA) of first and second-line treatments.
Results: Summary data were obtained from 21 trials of biologic therapies
including 2 for second-line treatment and results from six emulated target
trials of both treatment lines. Bivariate NMA resulted in a decrease in
uncertainty around the effectiveness estimates of the second-line therapies,
when compared to the results of univariate NMA, and allowed for predictions of
treatment effects not evaluated in second-line RCTs.
Conclusion: Bivariate NMA provides effectiveness estimates for all treatments
in first- and second-line, including predicted effects in second-line where
these estimates did not exist in the data. This novel methodology may have
further applications, for example for bridging networks of trials in children
and adults.

_Keywords: real world evidence, bivariate network meta-analysis, target trial
emulation, treatment lines, rheumatoid arthritis, biologic therapies_

##  1 Introduction

The evidence base for health care decision making traditionally consisted of
data from randomised controlled trials (RCTs), considered as a gold standard
in evaluation of health technologies. In recent years, there has been growing
interest in the use of real world data (RWD) from observational studies in
health care evaluation. Routinely collected data, from electronic health
records or patients’ registries, can provide useful information about
effectiveness of treatments, where data from RCTs may be sparse or are not
available at all for some treatment comparisons. Considerable methodological
research has focussed on inclusion of RWD in evidence synthesis with the aim
of overcoming some limitations of RCT data [1, 2, 3]. The focus of such
research has been particularly in circumstances where RCT evidence was sparse
and combining RCT data with RWD aimed to increase the evidence base to improve
the precision of effectiveness estimates [4] and sometimes bridge disconnected
networks.

Whilst research to date has largely focussed on exploitation of RWD to mimic
or replicate RCT data [5, 6, 7], we take a step further to explore use of RWD
in a scenario of data generation not typical for the RCT setting. In this
paper, we explored how RWD can be used to optimise an evidence base by using
evidence on first-line therapies to inform second-line effectiveness estimates
in evidence synthesis. When data from RCTs are available on effectiveness of a
particular treatment, but only in the first line of therapy, a costly trial
needs to be undertaken to also evaluate the effectiveness of the new therapy
used in patients as a second line treatment (or vice versa). We investigated
the added value of registry data, which provides evidence on both first and
second lines in each individual, when amalgamating these data in a network of
RCTs for both lines of therapies. We developed this approach for incorporating
RWD into clinical and HTA decision-making using a case study in rheumatoid
arthritis (RA).

We made use of data from the British Society for Rheumatology Biologics
Register for Rheumatoid Arthritis (BSRBR-RA) to supplement the RCT evidence
available only for either first- or the second-line of therapy. We did so by
emulating target trials using the approach developed by Hern$\acute{a}$n and
Robins [7]. We estimated treatment effects of biologic therapies based on the
data in emulated target trials, which we then used to inform a bivariate
network meta-analysis (NMA) model of first- and second-line treatments. The
estimates from the registry data were used to “bridge” disconnected networks
for the two lines of therapy. The American College of Rheumatology response
criteria (ACR20) was used as an outcome measure.

The remainder of this paper is structured as follows. Data sources and
statistical methods are described in Section 2. The results are presented in
Section 3, which are followed by discussion and conclusion in Section 4.

##  2 Methods

###  2.1 Data sources

####  2.1.1 Summary data from randomised controlled trials

Summary data from a literature review of RCTs of biologic therapies in
patients with rheumatoid arthritis were obtained for the effectiveness of
adalimumab, etanercept, infliximab, golimumab, abatacept and rituximab used as
first-line biologic therapies (in biologic naive patients) and the
effectiveness of golimumab and rituximab used as second-line biologic
therapies in patients who switched from a previous biologic treatment. Data
were obtained from 20 trials including 18 for the first-line treatments and
two trials for second-line treatments. When constructing a network, placebo
arms with methotrexate as concomitant therapy and the arms including a
combination of methotrexate and placebo were treated as the same treatment arm
in the network. Methotrexate, used in many trials as part of the combination
therapy in the biologic arm, was ignored (for some studies methotrexate was
included as concomitant therapy where percentage of patients with addition of
methotrexate varied across studies, similarly as in the BSRBR-RA target
trials).

####  2.1.2 Registry data

We made use of data from the British Society for Rheumatology Biologics
Register for Rheumatoid Arthritis (BSRBR-RA) to supplement randomised trial
evidence. Whilst RCTs included only either first- or second-line therapy,
registry data provided evidence on both lines of therapy for each patient.
BSRBR-RA data consisted of 19410 individuals, 15636 of whom had data recorded
on biologic treatment. The data were used to emulate trials of both lines of
therapy.

###  2.2 Emulation of target trials

We used the BSRBR-RA data to emulate a series of trials of first- and second-
line treatments for a range of biologic therapies using a target trial
approach [7]. In the first instance we specified the key components of the
target trial protocol, which (following recommendation by Hern$\acute{a}$n and
Robins [7]) included: eligibility criteria, treatment strategies, assignment
procedures, the follow-up period, outcome, causal contrast, and statistical
analysis. Note that no RCT included both lines of therapy in sequence, whilst
the proposed protocol of the target trial did include the treatment sequence.
Therefore, we did not aim for the emulated target trials to resemble any
existing RCT (an approach previously used in target trial emulation).

####  2.2.1 Eligibility criteria

Study participants were 18 years of age or older, who had a diagnosis of RA.
Patients who were treated with a biologic disease-modifying antirheumatic drug
(DMARD) prior to the registration with the BSRBR-RA were excluded.

####  2.2.2 Treatment strategies

Patients had to have received at least two lines of therapy, which could be
any of the biologic DMARDs or methotrexate, which is a synthetic DMARD often
used as a combination therapy and/or control treatment in trials of biologic
therapies in RA patients. Data from patients who switched from first line
biologic therapy to no therapy (or to therapies that are neither biologic
DMARDs nor methotrexate) were not included.

####  2.2.3 Assignment procedures

Patients were grouped into treatment arms according to the sequence of
treatment in two lines of therapy. These groups of patients (sequence
treatment arms) were matched to form experimental and control treatment
groups. Matching was conducted based on size of the trial, ensuring well
balanced treatment contrasts, with methotrexate always taken as the control
treatment and rituximab as an experimental treatment. Other biologic therapies
could be used as either experimental treatment or control. The matching
procedure had to ensure unique treatments in experimental and control arms for
each line. The process is schematically described in Figure 1. This procedure
resulted in target trials of two lines of therapy recorded on the same
patients who switched treatment in both treatment arms. For example, patients
in the experimental arm receiving first line adalimumab switched to infliximab
and those in the control arm receiving first line etanercept switched to
methotrexate, thus resulting in the first line comparison of adalimumab versus
etanercept and in the second line comparison of infliximab versus
methotrexate. Since patients were not randomly allocated, we assumed no
unmeasured confounding at baseline conditional on a number of prognostic
factors, measured at baseline or initiation of each treatment, that could
influence the response. The prognostic factors included age, gender, duration
of the disease and a number of clinical measures related to the disease
activity including the number of tender and swollen joints, serology (being
positive for rheumatoid factor), acute-phase reactants (CRP and ESR) and
28-joint count disease activity score (DAS-28).

####  2.2.4 The follow-up period

The minimum follow up time had to ensure that data were collected 24 weeks
after initiation of each line of therapy. Start of the second-line therapy
varied depending when patients needed to switch to second-line treatment,
which was typically due to either lack of response or adverse reactions.

####  2.2.5 Outcome

Patients were assessed according to ACR20 response criteria, which classified
them as responders if they had at least 20% improvement according to ACR
criteria. Due to a large number of missing values on some of the components of
ACR within BSRBR-RA data (the register did not capture patient pain or
physician global score), the definition of response was relaxed allowing
patients to be classed as responders if they had at least 20% improvement in
at least one of the joint count components (tender or swollen joint count) and
at least one of the remaining five components of the ACR measure (physician
global assessment, patient global assessment, pain, HAQ, ESR (or CRP)) [8].

####  2.2.6 Causal contrast and statistical analysis

Baseline characteristics for each group were summarised to ensure that the
covariates were similarly distributed across the treatment arms. The numbers
of responders were then adjusted for covariates using regression adjustment in
each line of therapy. Logistic regression was applied to the data with
baseline characteristics and treatment as covariates. We estimated the per-
protocol effect in all emulated trials.

$\Rightarrow$$\Rightarrow$$\Rightarrow$$\Rightarrow$$\Rightarrow$$\Rightarrow$First
lineSecond line

$\bullet$ IFX, $\bullet$ ADA, $\bullet$ ETA, $\bullet$ RTX, $\bullet$ MTX.

Figure 1: Schematic diagram representing the process of matching sequence
treatment arms. Each row represents patients assigned to unique treatment
sequence depicted by different colours. Gray arrows on the left hand side show
how the sequence treatment groups were matched; to ensure a balance in terms
of the size of each arm and that for both lines the treatment in each arm was
different. ADA – adalimumab, ETA – etanercept, IFX – infliximab, RTX –
rituximab, MTX – methotrexate.

###  2.3 Bivariate network meta-analysis

We used bivariate NMA to model jointly the treatment effects on ACR20 for
first- and second-lines of therapies. A standard approach to any multivariate
meta-analysis is to use a hierarchical model with a multivariate normal
distribution used to describe variability at two levels: within-study (where
the correlation occurs due to the modelled multivariate quantities, such as
treatment effects on multiple outcomes, being measured in the same
individuals) and between-studies (where the correlation is a result of
heterogeneity between the average effects, measured on each outcome in each
study, varying across studies due to, for example, differences in population
or treatment doses). Accounting for the within-study correlation is important
in such analysis [9]. However, modelling jointly non-normal outcomes, such as
Binomial responses, would require transforming data, which can lead to biased
results [10]. Papanikos et al, carried out a simulation study showing that,
when the within-study correlation is weak, a multivariate meta-analysis model
with independent binomial likelihoods is preferable [11]. Exploratory analysis
of the BSRBR-RA data set, estimating the within-study correlation using the
bootstrapping approach [12, 13], showed that the within-study correlation
between the treatment effects for the two lines of therapy transformed onto
the log odds ratio (OR) scale was close to zero. We, therefore, adapted the
approaches to multivariate/bivariate NMA by Achana et al. [14] and Bujkiewicz
et al. [15] by assuming independent binomial likelihoods at the within-study
level, as in Papanikos et al [11], to model the number of responders $r_{ijk}$
to treatment $k$ in line of therapy $j=1,2$ out of the number of participants
$n_{ik}$ in study $i$, with the probability of response denoted by $p_{ijk}$;

$\displaystyle r_{ijk}$ $\displaystyle\sim$ $\displaystyle
Binomial(p_{ijk},n_{ik})$ $\displaystyle logit(p_{ijk})$ $\displaystyle=$
$\displaystyle\theta_{ijk}$ $\displaystyle\theta_{ijk}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{lll}\mu_{ijb}&if&k=b\\\
\mu_{ijb}+\delta_{ij,bk}&if&k>b\end{array}\right.$ (3)

$\left(\begin{array}[]{c}\delta_{i1,bk}\\\
\delta_{i2,bk}\end{array}\right)\sim\rm{MVN}\left(\left(\begin{array}[]{c}d_{1,bk}\\\
d_{2,bk}\end{array}\right),\;\mathbf{T}=\left(\begin{array}[]{cc}\tau_{1}^{2}&\tau_{1}\tau_{2}\rho\\\
\tau_{1}\tau_{2}\rho&\tau_{2}^{2}\end{array}\right)\right)$
$\left(\begin{array}[]{c}d_{1,bk}\\\ d_{2,bk}\\\
\end{array}\right)=\left(\begin{array}[]{c}d_{1,1k}-d_{1,1b}\\\
d_{2,1k}-d_{2,1b}\\\ \end{array}\right)$

In this hierarchical model, $\theta_{ijk}$ denotes log odds of response to
treatment $k$ in line $j$ and study $i$, $\mu_{ijb}$ is the baseline treatment
effect in study $i$ line $j$, and $\delta_{ij,bk}$ are true (random) treatment
effects in each study $i$ comparing treatment $k$ with the baseline treatment
$b$ in line $j$. We assume that the true effects follow a bivariate normal
distribution that is common to the studies investigating the same treatment
contrast $bk$ with mean effect $d_{j,bk}$ and the between-studies covariance
matrix $\mathbf{T}$ (here assumed homogeneous across the treatment contrasts).
The network structure of the data is taken into account by assuming that the
pooled effects $d_{j,bk}$ for each contrast $bk$ and treatment line $j=1,2$
satisfy the consistency assumption (2.3), where $d_{j,1k}$ denote the basic
parameters (average treatment effects of each treatment $k=1,\dots,n_{t}$ in
the network relative to the reference treatment $1$ and treatment line
$j=1,2$, with $n_{t}$ denoting the number of treatments in the network.

Prior distributions are placed on the between-studies heterogeneity parameters
$\tau_{j}\sim U(0,2)$, the between-studies correlation $\rho=r*2-1$, $r\sim
Beta(1.5,1.5)$, the baseline effects $\mu_{ijb}\sim N(0,10^{3})$ and the basic
parameters $d_{j,1k}\sim N(0,10^{3})$. When the between-studies correlation is
zero, the model reduces to two univariate models for the two outcomes modelled
independently; with $\delta_{ij,bk}\sim
N((d_{j,1b}-d_{j,1k}),\tau_{j}^{2}),\;j=1,2$.

To predict treatment effect in the second line when data are only available
for the therapy in first line, additional assumptions of exchangeability need
to be made, where instead of placing a prior distributions on basic
parameters, we add another level of hierarchy to the model as in Bujkiewicz et
al [15]. For each treatment arm $k$; ancillary parameters $\vartheta_{jk}$ for
the two treatment lines $j=1,2$, such that
$d_{j1k}=\vartheta_{jk}-\vartheta_{j1}$, are assumed exchangeable (for the
biologic therapies only) and correlated:

$\left(\begin{array}[]{c}\vartheta_{1k}\\\ \vartheta_{2k}\\\
\end{array}\right)\sim N\left(\left(\begin{array}[]{c}\eta_{1}\\\ \eta_{2}\\\
\end{array}\right),\frac{1}{2}\left(\begin{array}[]{cc}\omega_{1}^{2}&\omega_{1}\omega_{2}\rho_{t}\\\
\omega_{1}\omega_{2}\rho_{t}&\omega_{2}^{2}\\\ \end{array}\right)\right)$

for $k=2,\dots,n_{t}$ (across the biologic therapies only).
Prior distributions are placed on the parameters $\omega_{j}\sim Unif(0,2)$
and $\rho_{t}=r_{t}*2-1$ with $r_{t}\sim Beta(1.5,1.5)$.

##  3 Results

###  3.1 Summary of data and the network structure

Summary data were obtained from 20 RCTs of biologic therapies with 18 trials
for first line treatment (including adalimumab, etanercept, infliximab,
golimumab, abatacept and rituximab) and two for second-line treatment
(including golimumab and rituximab). BSRBR-RA data included 12657 individuals
given first line biologic at the time of registration. Follow up data included
112983 observations, which was on average 8.93 follow-ups per individual. For
a large proportion of the visits, methotrexate was recorded as a concomitant
therapy to a biologic treatment. Target trial emulation using the BSRBR-RA
data led to generation of six trials of biologic therapies in two lines. Table
1 includes information on the number of studies for each treatment contrast,
line of therapy and the type of study.

|  | number |  | first line |  | second line
---|---|---|---|---|---|---
source |  | of studies |  | treatments |  | treatments
first line RCTs | aaaaa | 6 | aaaaa | ADA vs MTX | aaaaa |
| 1 |  | ETA vs MTX |  |
| 4 |  | IFX vs MTX |  |
| 5 |  | GOL vs MTX |  |
| 2 |  | ABT vs MTX |  |
| 1 |  | RTX vs MTX |  |
second line RCTs |  | 1 |  |  |  | GOL vs MTX
| 1 |  |  |  | RTX vs MTX
target trials |  | 1 |  | IFX vs ADA |  | RTX vs IFX
| 1 |  | ADA vs ETA |  | RTX vs IFX
| 1 |  | ETA vs IFX |  | RTX vs ADA
| 1 |  | IFX vs ETA |  | ETA vs MTX
| 1 |  | ETA vs ADA |  | ADA vs MTX
| 1 |  | ADA vs IFX |  | ETA vs MTX
Table 1: Data sources and treatments; ADA – adalimumab, ETA – etanercept, IFX
– infliximab, GOL – golimumab, ABT – abatacept, RTX – rituximab, MTX –
methotrexate.

Figure LABEL:net-alla shows the network structure of RCT data for the first
and second lines of therapy and Figure LABEL:net-allb illustrates the network
structure of target trials emulated from BSRBR-RA data for both lines of
therapy. For the target trials, both treatment lines correspond to the same
trial, in contrast to the RCTs which report only either first or second line
of treatment. To emphasise this in Figure LABEL:net-all, we used the same
colour of the network edges for both treatment lines for the target trials, in
contrast to the RCTs where different colours of edges for different treatment
lines represent different trials. In this paper, we aimed to demonstrate the
value of the registry data in estimating the effect of the biologic therapies
when used as second-line treatments. The network of RCT data for the second-
line therapy was particularly sparse, including only two trials; for golimumab
and rituximab. BSRBR-RA data gave additional information about adalimumab,
etanercept, infliximab as well as rituximab used as second-line therapies. The
network structure for RCT and BSRBR-RA data combined is shown in Figure
LABEL:net-allc.

###  3.2 Results of network meta-analyses

Results of two univariate NMAs of biologic therapies used as second-line
treatments are shown in Table 2, where the lower triangle includes the results
from the NMA using data from RCTs alone and the upper triangle shows results
of NMA combining data from RCTs and BSRBR-RA register. Including the registry
data allowed for estimation of treatment effects for second line biologic
therapies which were not included in the RCT network. There was also some
improvement in the precision of treatment effect estimates for those already
included in the RCT data. For example, comparing rituximab with methotrexate
gave OR=11.34 (95% CrI: 0.37, 59.4) when using RCT data alone, whilst addition
of the registry data resulted in OR=3.81 (0.48, 14.5) thus reducing the
uncertainty by 76% in terms of the width of the credible interval. The
between-studies correlation was weak; $\rho$=-0.17 (95% CrI: -0.86, 0.8).

| MTX | ADA | ETA | IFX | GOL | ABT | RTX
---|---|---|---|---|---|---|---
|  | 4.69 | 4.41 | 11.0 | 4.89 |  | 3.81
MTX |  | (0.6, 17.5) | (0.79, 14.5) | (0.56, 50.9) | (0.32, 21.5) | – | (0.48, 14.5)
|  |  | 2.14 | 3.42 | 2.66 |  | 1.18
ADA | – |  | (0.11, 9.67) | (0.17, 15.5) | (0.05, 12.1) | – | (0.15, 4.41)
|  |  |  | 4.90 | 2.17 |  | 1.60
ETA | – | – |  | (0.11, 23.1) | (0.06, 10.1) | – | (0.08, 7.36)
|  |  |  |  | 2.76 |  | 0.66
IFX | – | – | – |  | (0.02, 10.8) | – | (0.11, 2.19)
| 6.41 |  |  |  |  |  | 3.31
GOL | (0.21, 33.4) | – | – | – |  | – | (0.07, 15.0)
ABT | – | – | – | – | – | (0.10, 00.1) | –
| 11.34 |  |  |  | 15.7 |  |
RTX | (0.37, 59.4) | – | – | – | (0.05, 64.3) | – |
Table 2: Results of a univariate NMA of biologic therapies used as second-line
treatments using data from RCTs alone (lower triangle) and combining data from
RCTs and BSRBR-RA registry (upper triangle).

The results of a bivariate NMA combining data from RCTs and BSRBR-RA of
biologic therapies in both lines of therapy are shown in Table 3, with the
results using the “standard” bivariate NMA model in the upper triangle and the
results from the analysis assuming exchangeability of biologic therapies in
the lower triangle. Combining data from the first and second lines of therapy
through the use of the bivariate NMA led to a further decrease in uncertainty
for most of the individual treatments when compared to the results of the
univariate NMA of second-line therapy alone. For example, ACR20 response to
adalimumab vs methotrexate from the bivariate NMA was OR=4.36 (0.67, 15.5)
compared to OR=4.69 (0.6, 17.5) from the univariate NMA.

The bivariate NMA approach assuming the additional exchangeability of the
absolute effects of the biologic therapies allowed for predictions of
treatment effects that had not been evaluated in trials in a second-line
setting. In this case, it produced effectiveness estimates for abatacept in
the second line of therapy against all other treatments in the network.
Moreover, this additional exchangeability led to a noticeable reduction in
uncertainty around the remaining estimates of effect for other therapies. This
was a result of additional borrowing of information. However, there may have
been some degree of smoothing of the effects across the biologic therapies,
which was difficult to judge due to the large uncertainty. A sensitivity
analysis was carried out using a _t_ -distribution in place of the normal
distribution in model (2.3), which largely produced very similar results but
inflated the uncertainty around the effectiveness estimate for abatacept in
the second line.

| MTX | ADA | ETA | IFX | GOL | ABT | RTX
---|---|---|---|---|---|---|---
MTX |  | 4.36 | 4.03 | 10.2 | 5.07 |  | 3.54
|  | (0.67, 15.5) | (0.86, 12.3) | (0.53, 46.7) | (0.3, 23.0) | – | (0.44, 13.5)
| 3.33 |  | 1.91 | 3.0 | 2.51 |  | 1.06
ADA | (1.37, 6.97) |  | (0.12, 8.22) | (0.18, 13.3) | (0.06, 12.1) | – | (0.16, 3.75)
| 3.36 | 1.12 |  | 4.63 | 2.12 |  | 1.50
ETA | (1.5, 6.58) | (0.44, 2.41) |  | (0.11, 21.7) | (0.06, 10.4) | – | (0.08, 6.87)
| 4.2 | 1.34 | 1.33 |  | 2.89 |  | 0.62
IFX | (1.42, 10.6) | (0.53, 3.26) | (0.48, 3.43) |  | (0.02, 12.1) | – | (0.12, 1.96)
| 3.28 | 1.07 | 1.05 | 0.91 |  |  | 3.18
GOL | (1.16, 7.38) | (0.36, 2.45) | (0.35, 2.36) | (0.24, 2.03) |  | – | (0.06, 14.6)
| 4.08 | 1.31 | 1.29 | 1.07 | 1.42 |  |
ABT | (0.77, 12.3) | (0.26, 3.87) | (0.24, 3.76) | (0.18, 2.92) | (0.26, 4.43) |  | –
| 2.82 | 0.89 | 0.9 | 0.75 | 0.97 | 1.1 |
RTX | (1.09, 5.96) | (0.39, 1.65) | (0.32, 1.8) | (0.3, 1.31) | (0.32, 2.19) | (0.21, 2.96) |
Table 3: Results of a bivariate NMA combining data from RCTs and BSRBR-RA of
biologic in both lines of therapy using the “standard” bivariate NMA model
(upper triangle) and assuming exchangeability of biologic therapies (lower
triangle).

##  4 Discussion and conclusions

In this paper, we provide a conceptual approach for using RWD, such as from
registries or electronic health records, to generate estimates of
effectiveness of treatments in first and second lines of therapy and combining
them with RCT data to enhance the evidence base and provide effectiveness
estimates of therapies in the second line, where data on effectiveness in the
second line are not available from RCTs. In such circumstances, producing
these estimates would require conducting expensive and time-consuming
additional clinical trials. The proposed approach can be used to carry out
feasibility analysis or provide inputs to the trial design or even be used for
evidence based decision-making where evidence is sufficiently robust.

When carrying out this research, we came across a number of limitations. Some
of them were related to data. In particular, the RCT data were relatively
sparse with a star-shaped network for the first line treatments and only two
trials reporting the effectiveness of biologic therapies in the second line.
The data set was simplified by combining the control arms (including
methotrexate as either combination therapy or concomitant therapy with
placebo) into the same control arm denoted as methotrexate. This was done to
strengthen the network structure to better illustrate the methodological
aspect of this work. Most of the biologic arms also included methotrexate.
Considering that for a large proportion of visits in the BSRBR-RA data
methotrexate was recorded as a concomitant therapy to a biologic treatment, an
assumption was made that a large proportion of patients receiving biologic
therapy, across all studies, also received methotrexate. The registry data set
contained a substantial amount of missing data, in particular for some of the
components of the ACR20 response criteria, which was not due to the issues of
quality of the data but owing to the fact that some of the components are not
routinely collected by the register. In order to estimate the response to the
biologic therapies, we chose to relax the definition of the response. In
addition, the register only captures 28 joint counts, which may be different
from some of the trials. Considering these potentially strong assumptions
around the data sources, the results of our analysis should not be used for
clinical interpretation, but only as an illustration of the proposed
methodology.

There were only six target trials generated from the registry data, which
resulted in substantial uncertainty around the between-studies correlation, as
these were the only studies contributing data to estimating the correlation.
The combined network was still limited with lack of data on each contrast and
line across study designs. Target trial data were incorporated in the NMAs at
face value, assuming they were equivalent with RCT data. Extensions of the
analysis could include a power prior approach [16], allowing for down-
weighting of RWD, or hierarchical modelling to differentiate between the two
study designs [4]. Further investigation into data scenarios and model
assumptions needs to be carried out to understand when this framework can be
most efficient.

###  4.1 Conclusions

Registry data can be used to bridge networks of first and second lines of
therapy which are disconnected when using RCT data alone. Bivariate NMA of
combined data from RCTs and RWD can be used to predict effectiveness of a
treatment in second line use when the therapy is only investigated in a RCT as
first line (or vice versa). The approach can be applied to other settings
where RCT data are available for disjoint subsets of population, such as, for
example, children and adults and registries may provide data covering follow-
up period from childhood to adulthood for each individual.

##  5 Contributions

SB and KRA conceived the research idea. SB conceptualised the study. SB, JS,
LW, DJ and KH curated the data. SB and JS undertook formal analyses. KH
provided clinical input. SB and KRA provided supervision for JS, LW and DJ. SB
contributed the original draft of the manuscript and visualisation. All
authors contributed to manuscript revisions.

##  6 Acknowledgements

The early work leading to these results has received support from the
Innovative Medicines Initiative Joint Undertaking under grant agreement no
[115546], resources of which are composed of financial contribution from the
European Union’s Seventh Framework Programme (FP7/2007- 2013) and EFPIA
companies’ in kind contribution. In subsequent research KRA and SB were
supported by the UK Medical Research Council (MRC) Methodology Research Panel
(MRP) [grant no. MR/R025223/1]. SB was also partially supported by the MRC MRP
[MR/L009854/1]. JS and LW were supported by UK National Institute for Health
Research (NIHR) Methods Fellowship [award no. RM-FI-2017-08-027]. KLH is
supported by the NIHR Manchester Biomedical Research Centre.

The authors would like to thank the British Society for Rheumatology for use
of this data. The BSR Biologics Register in Rheumatoid Arthritis (BSRBR-RA) is
a project owned and run by the BSR on behalf of its members. BSR would like to
thanks its funders for their support, currently AbbVie, Amgen, Celltrion, Eli
Lilly, Pfizer, Samsung, Sandoz, Sanofi, and in the past Roche, UCB, SOBI and
Merck. This income finances a wholly separate contract between BSR and The
University of Manchester to host the BSRBR-RA. All decisions concerning
analyses, interpretation and publication are made autonomously of any
industrial contribution. The BSRBR-RA would like to gratefully acknowledge the
support of BSR members and the National Institute for Health Research, through
the Comprehensive Local Research Networks at participating centres.

##  7 Disclosures

SB has served as a paid consultant, providing methodological advice, to NICE,
Roche and RTI Health Solutions, has received payments for educational events
from Roche and has received research funding from European Federation of
Pharmaceutical Industries & Associations (EFPIA) and Johnson & Johnson.
JS,LW,DJ do not have conflict of interest. RM has served as a consultant for
Roche,Takeda, Freeline and Khondrion KH received honoraria from Abbvie
(Speaker fees) and grant income from Pfizer and BMS. KRA has served as a paid
consultant, providing methodological advice, to; Abbvie, Amaris, Allergan,
Astellas, AstraZeneca, Boehringer Ingelheim, Bristol-Meyers Squibb, Creativ-
Ceutical, GSK, ICON/Oxford Outcomes, Ipsen, Janssen, Eli Lilly, Merck, NICE,
Novartis, NovoNordisk, Pfizer, PRMA, Roche and Takeda, and has received
research funding from Association of the British Pharmaceutical Industry
(ABPI), European Federation of Pharmaceutical Industries & Associations
(EFPIA), Pfizer, Sanofi and Swiss Precision Diagnostics. He is a Partner and
Director of Visible Analytics Limited, a healthcare consultancy company.

##  8 References

## References

  * [1] Orestis Efthimiou, Dimitris Mavridis, Thomas PA Debray, Myrto Samara, Mark Belger, George CM Siontis, Stefan Leucht, Georgia Salanti, and GetReal Work Package 4.  Combining randomized and non-randomized evidence in network meta-analysis.  Statistics in medicine, 36(8):1210–1226, 2017.
  * [2] Pablo E Verde and Christian Ohmann.  Combining randomized and non-randomized evidence in clinical research: a review of methods and applications.  Research synthesis methods, 6(1):45–62, 2015.
  * [3] Nicky J Welton, Phillippo David M, Rhiannon Owen, Hayley J Jones, Sofia Dias, Sylwia Bujkiewicz, Anthony E Ades, and Keith R Abrams.  Chte2020 sources and synthesis of evidence; update to evidence synthesis methods.  DSU Report, March 2020.
  * [4] Susanne Schmitz, Roisin Adams, and Cathal Walsh.  Incorporating data from various trial designs into a mixed treatment comparison model.  Statistics in medicine, 32(17):2935–2949, 2013.
  * [5] Victoria L Bartlett, Sanket S Dhruva, Nilay D Shah, Patrick Ryan, and Joseph S Ross.  Feasibility of using real-world data to replicate clinical trial evidence.  JAMA network open, 2(10):e1912869–e1912869, 2019.
  * [6] Reynaldo Martina, David Jenkins, Sylwia Bujkiewicz, Pascale Dequen, and Keith Abrams.  The inclusion of real world evidence in clinical development planning.  Trials, 19(1):1–12, 2018.
  * [7] Miguel A Hernán and James M Robins.  Using big data to emulate a target trial when a randomized trial is not available.  American journal of epidemiology, 183(8):758–764, 2016.
  * [8] David T Felson, Jennifer J Anderson, Maarten Boers, Claire Bombardier, Daniel Furst, Charles Goldsmith, Linda M Katz, Robert Lightfoot Jr, Harold Paulus, Vibeke Strand, et al.  American college of rheumatology preliminary definition of improvement in rheumatoid arthritis.  Arthritis & Rheumatism: Official Journal of the American College of Rheumatology, 38(6):727–735, 1995.
  * [9] Richard D Riley.  Multivariate meta-analysis: the effect of ignoring within-study correlation.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 172(4):789–811, 2009.
  * [10] Taye H Hamza, Hans C van Houwelingen, and Theo Stijnen.  The binomial distribution of meta-analysis was preferred to model within-study variability.  Journal of clinical epidemiology, 61(1):41–51, 2008.
  * [11] Tasos Papanikos, John R Thompson, Keith R Abrams, and Sylwia Bujkiewicz.  A novel approach to bivariate meta-analysis of binary outcomes and its application in the context of surrogate endpoints.  arXiv preprint arXiv:2004.02007, 2020.
  * [12] Michael J Daniels and Michael D Hughes.  Meta-analysis for the evaluation of potential surrogate markers.  Statistics in medicine, 16(17):1965–1982, 1997.
  * [13] Sylwia Bujkiewicz, John R Thompson, Alex J Sutton, Nicola J Cooper, Mark J Harrison, Deborah PM Symmons, and Keith R Abrams.  Multivariate meta-analysis of mixed outcomes: a bayesian approach.  Statistics in Medicine, 32(22):3926–3943, 2013.
  * [14] Felix A Achana, Nicola J Cooper, Sylwia Bujkiewicz, Stephanie J Hubbard, Denise Kendrick, David R Jones, and Alex J Sutton.  Network meta-analysis of multiple outcome measures accounting for borrowing of information across outcomes.  BMC medical research methodology, 14(1):1–16, 2014.
  * [15] Sylwia Bujkiewicz, Dan Jackson, John R Thompson, Rebecca M Turner, Nicolas Städler, Keith R Abrams, and Ian R White.  Bivariate network meta-analysis for surrogate endpoint evaluation.  Statistics in medicine, 38(18):3322–3341, 2019.
  * [16] Joseph G Ibrahim and Ming-Hui Chen.  Power prior distributions for regression models.  Statistical Science, pages 46–60, 2000.