Enhancement of Non-Linear Deep Learning Model by Adjusting Confounding Variables for Bone Age Estimation in Pediatric Hand X-rays

Dataset and Preprocessing

A large dataset for skeletal age prediction consisting of hand X-ray images with bone age and sex information has been released by the Radiological Society of North America (RSNA) [26, 29]. The RSNA dataset is composed of 12,611 training images, 1,425 validation images, and 200 testing images. To train the bone age prediction model, only the RSNA dataset was used. For this dataset, the Greulich and Pyle [30] (GP) method, which interprets bone age by comparing hand X-rays with representative images, was used to measure the ground truth bone age. For external validation, 227 hand X-rays of clinically normal participants with sex information from the South Korean referral hospital, Asan Medical Center (AMC), were used. For the AMC dataset, two expert radiologists labeled the bone ages using the GP method. The mean of the bone age label and chronological age was used for evaluation. The bone ages of the RSNA test dataset were similarly distributed to those in the train dataset; the study populations were sampled according to the uniform distribution of chronological age in the AMC dataset. These two datasets showed significantly different distributions according to two-sample Kolmogorov–Smirnov test for goodness of fit (P-value < 0.001). The bone ages of participants ranged from 11 to 219 months in the RSNA dataset, while they ranged from 24 to 228 in the AMC dataset. The age histograms of the RSNA and AMC datasets are depicted in Fig. 1.

When preprocessing the images, they were first adjusted from the top 1% pixel value to avoid L/R markers and other high intensity artifacts acting as confounding factors [31]. To preserve the original image shape, the images were padded into a square. And they were reshaped into 512 × 512 pixels, considering the sophisticated nature of medical images and resource limitations, such as training time and GPU memory. The pixel values were normalized from 0 to 1 using min–max scaling; Additionally, contrast limited adaptive histogram equalization was used to emphasize the bony part of hand X-rays.

U-Net [32] was trained for hand-image segmentation using a semi-automated data labeling strategy. 200 randomly sampled hand X-ray images were binarized and small clusters were removed. The first U-Net segmentation network was trained using these masks and the intersection over union (IoU) for the baseline 200 images was 0.982. Subsequently, additional masks were generated from this network using the training dataset. Next, well-generated masks were sampled and the U-Net segmentation network was trained with these masks. The final IoU for the entire RSNA dataset was 0.975. Finally, a hand mask was generated for all training, validation, and testing datasets. Hand images only were extracted and used for bone age prediction. The hand mask generation and hand extraction processes are summarized in Fig. 2.

Network Architectures

U-Net based multi-task learning (MTL) architecture was adopted for bone age estimation. The network simultaneously conducted bone age estimation and image reconstruction tasks; therefore, the model could learn good representations when reconstructing the original image from the augmented image. An auxiliary regression layer was added to the bottleneck layer of the U-Net architecture to estimate bone age. Furthermore, because bone age estimation was the main target for the entire architecture, additional encoder layers were added to accelerate the target auxiliary regression task. The schematic figure of U-Net MTL and auxiliary-accelerated MTL (AA-MTL) architecture is depicted in Fig. 3.

Handling Methods of Confounding Variables

In the typical statistical model (i.e., generalized linear models), the relationship between the observation Y_i and independent variable X_ij is formulated as a linear combination of independent variables and parameters as:

given a random sample (Y_i, X_i1, X_i2, …, X_ij-1, X_ij), where W₁, W₂, …, W_j-1, W_j are the parameters for the variables X_i1, X_i2, …, X_ij-1, X_ij; and β_i is the random bias.

However, considering the non-linear properties of deep learning, the relationship between the input vector v_ij and output vector u_ij cannot be defined as in Eq. (1). The relationship between vectors v_ij and u_ij in a non-linear deep learning model $f(\bullet )$ can be formulated as:

Therefore, three different handling methods of the important confounding variable (sex) in non-linear deep learning models were compared. The model that did not use sex information was regarded as the baseline. The models that used sex information as the input vector and predicted sex information as the output vector were trained. When sex was set as an input to a deep learning model, sex was concatenated just before the fully connected layer for bone age estimation; when sex was set as an output, an auxiliary classifier was added to predict it.

Furthermore, the backpropagation from one loss function was hypothesized to affect another loss function, and vice versa. Therefore, sex classification and bone age estimation at the diverse encoder level were investigated when sex was set as the prediction output. The hierarchical model structures of the following were compared. Bone age estimation after bottleneck block, followed by sex classification after second encoder block (the two-convolutional layer followed by the max pooling layer or input image was considered as an encoder block, see Fig. 3); sex classification after bottleneck block, followed by bone age estimation after second encoder block; and sex classification and bone age estimation simultaneously after bottleneck block.

Average Class Activation Maps

Saliency maps using gradient-weighted class activation map (Grad-CAM) [33] were acquired to assess the explainability of each model. The Grad-CAM images were stratified by pubertal stage because the growth patterns of children change with the pubertal stage. Therefore, the hand X-ray images were grouped into pre-puberty (0–6 years), puberty (6–14 for boys, 6–12 for girls), and post-puberty (14–18 for boys, 12–18 for girls) [34–36].

Next, the image registration technique was applied to the pubertal stage-stratified Grad-CAMs to show each model focus on. First, the representative image of each pubertal stage was selected. The center, width, and height from the hand mask of each image was acquired. Subsequently, an affine transformation using the center, width, and height of the hand mask was applied to the Grad-CAM and averaged. Finally, the averaged Grad-CAM was superimposed on the representative image of each pubertal stage.

Statistical Analyses

Parametric and nonparametric statistics were used on the RSNA and AMC datasets, respectively, as the RSNA dataset was similarly distributed to the train dataset, while the AMC dataset was not. The mean average error (MAE) and Pearson’s correlation coefficient (r) or the Kendall rank correlation coefficient (τ) between the model-estimated bone age and ground truth age were evaluated and compared. A paired t-test and Wilcoxon signed-rank test were conducted on the RSNA and AMC datasets, respectively, to compare the statistical differences between each model. Correlation plots and Bland–Altman plots [37, 38] were constructed to show the correlation between the estimated bone age and ground truth age. The statistics software R version 4.2.0 (R Foundation for Statistical Computing, Vienna, Austria) was used for all statistical analyses. Two-sided P-values were used to determine statistical significance, which was set at an alpha value of 0.05.

Implementation Details

U-Net [32] based encoder-decoder architecture was adopted. Because deep learning models can learn good representations from self-reconstruction of the original image [39], the U-Net decoder learnt to reconstruct the original image from the augmented image. The models without decoders were trained to isolate the effect of this self-reconstruction. All models were implemented in Python version 3.6.9 with PyTorch version 1.9.0. Adam optimizer was chosen with a learning rate of 5e-3, and the ReduceLROnPlateau learning rate scheduler was used. The models were trained with augmentation techniques plausible for medical images, such as shifting, zooming, rotation, blur, gaussian noise addition, and sharpening. Additionally, strong augmentation, such as cutout [40], was used to the extent that it did not compromise clinical viability. All networks were trained using a single GPU (Titan RTX) and a batch size of eight. For bone age estimation, regression loss was calculated as L1 Loss, as was image reconstruction loss. L1 loss was calculated as:

where y is ground truth age and $\widehat{y}$ is bone age regression output when applied as regression loss. When applied as reconstruction loss, y was regarded as the original image and $\widehat{y}$ as the model reconstruction output. When sex was set as the prediction output, the classification loss was calculated using cross entropy loss:

where $\widehat{y}$ is the softmax probability of the model classification output. As the target task was bone age estimation and the other tasks were performed for confounding adjustment, the loss weights for the image reconstruction and sex classification were set to 0.1 and the sum of loss weights was set to 1 (i.e., bone age estimation weight + image reconstruction weight, if performed + sex classification weight, if performed = 1). The models were trained from the scratch without pretrained weights.

Comparison of Handling Methods of Confounding Variables

A comparison of the confounding handling methods in U-Net backbone, U-Net MTL, and AA-MTL is summarized in Table 1. The absolute errors between the model-predicted bone age and ground truth age label were compared for each confounding adjustment method and the best performing method for each architecture. The correlation analyses between the model-predicted bone age and ground truth age labels are depicted in Supplementary Table 1.

In the U-Net backbone architecture, the confounding adjustment by input method showed the best performances in the RSNA dataset with an MAE of 5.740 months (P-value < 0.001, no adjustment and adjustment by output). For the AMC dataset bone age, confounding adjustment by input showed statistically better performance than confounding adjustment by output (P-value = 0.023). For the AMC dataset chronological age, confounding adjustment by output outperformed no adjustment (P-value < 0.001) and adjustment by input (P-value = 0.002).

The confounding adjustment by input method showed the best performances in bone age estimation of the RSNA dataset when compared with the other methods (P-value < 0.001 in both methods) in the U-Net MTL architecture. For the AMC dataset bone age, a significant difference between the confounding adjustment by input and output (P-value = 0.027) was observed, whereas no statistical difference between input and no adjustment (P-value = 0.205) was observed. For the chronological age, there were no significant differences between the best performing method, no adjustment, and adjustment by input (P-value = 0.542) and output (P-value = 0.329).

Finally, in the proposed AA-MTL architecture, confounding adjustment by input outperformed other methods in the RSNA dataset with an MAE of 5.434 months (P-value < 0.001, no adjustment and adjustment by output). For the AMC dataset bone age, adjustment by output significantly outperformed no adjustment (P-value = 0.009). However, no significant difference between output and input was observed (P-value = 0.814). For the AMC dataset chronological age, adjustment by output outperformed no adjustment (P-value = 0.006), but not adjustment by input (P-value = 0.051).

Comparison of Conducting Multi-task at Diverse Encoder Level

The performance of each MTL network depended on the other tasks and model size [1, 41]. Considering the multi-task nature of this study, confounding adjustments and target task were conducted at the diverse encoder levels. Three different multi-task combinations of confounding adjustments and task were performed: sequential (sex classification at second encoder level, bone age estimation at bottleneck level), simultaneous (sex classification and bone age estimation at bottleneck level), and inverse-sequential (bone age estimation at second encoder level, sex classification at bottleneck level). The U-Net backbone architecture performed only sex classification and bone age estimation. The U-Net MTL and AA-MTL architectures performed sex classification, bone age estimation, and source image reconstruction from augmented images.

For the RSNA dataset bone age, the simultaneous model outperformed the inverse-sequential model (P-value = 0.007); however, it did not outperform the sequential model (P-value = 0.328) with double-task U-Net backbone architecture. In the AMC dataset, the simultaneous model outperformed the inverse-sequential model (P-value = 0.008) but not the sequential model (P-value = 0.075) with bone age labels; whereas the sequential model outperformed the inverse-sequential (P-value < 0.001) and simultaneous (P-value = 0.021) models with chronological age labels.

In the U-Net MTL architecture, the simultaneous model showed the best performance with an MAE of 7.601 months for the RSNA dataset bone age; however, no statistical difference was observed compared to the sequential (P-value = 0.080) or inverse-sequential (P-value = 0.281) models. In the AMC dataset, the simultaneous model outperformed both the inverse-sequential (P-value < 0.001) and sequential (P-value = 0.006) models with bone age labels, and the inverse-sequential model (P-value = 0.003) but not the sequential model (P-value = 0.148) with chronological age labels.

For the RSNA dataset bone age, the simultaneous model showed the best performance but did not have statistically significant results compared to the inverse-sequential (P-value = 0.115) or sequential (P-value = 0.249) models with AA-MTL architecture. For the AMC dataset with bone age labels, the sequential model of AA-MTL architecture showed the best performance with an MAE of 8.190 m. The AA-MTL sequential model significantly outperformed the inverse-sequential model (P-value = 0.001), but no statistical difference between it and the simultaneous model (P-value = 0.990) was observed. The sequential model of AA-MTL architecture showed the best performance for the AMC dataset with chronological age labels: it outperformed the inverse-sequential (P-value = 0.002) and simultaneous (P-value = 0.001) models. Table 2 summarizes the comparison results of conducting multi-task at diverse encoder level.

Comparison of Multiple Tasks in Performance and Generalizability

A comparison of the single-task model, which only performed bone age estimation; double-task model, which performed bone age estimation and source image reconstruction; and triple-task model, which performed bone age estimation, source image reconstruction, and sex classification is summarized in Table 3. In the U-Net MTL architecture, the double-task model with sex information showed the best performance in both datasets. In the AA-MTL architecture, the double-task model with sex information showed the best performance in the RSNA dataset. However, in the AMC dataset, the triple-task model showed the best performance with both bone age and chronological age labels.

Correlation plot and Bland–Altman plot of each multi-task model in the AA-MTL architecture are depicted in Fig. 4. All three models showed decent performances in estimating bone ages, with R² values of 0.969, 0.975, and 0.946 in the single-, double-, and triple-task models in the RSNA dataset, respectively. The Bland–Altman plot of all three models showed good agreement between the ground truth and model-predicted bone ages.

In the AMC dataset, all models showed decent performances in estimating bone ages, with R² values of 0.963, 0.962, and 0.969 in the single-, double-, and triple-task models, respectively. The Bland–Altman plots of all three models showed good agreements between the ground truth and model-predicted bone ages, with mean differences of 1.382 (95% confidence interval [CI], –13.543–16.306), 1.827 (95% CI, –11.423–15.077), and 2.622 (95% CI, –17.146–22.391) for the single-, double-, and triple-task models in the RSNA dataset, respectively; additionally, they showed good agreements with mean differences of –0.198 (95% CI, –22.480–22.085), –0.759 (95% CI, –23.247–21.729), and 0.625 (95% CI, –20.298–21.548) in the AMC dataset, respectively.

Average Class Activation Map

The qualitative results were evaluated using the average Grad-CAMs. Grad-CAMs acquired from each image were coordinated according to the representative image of each puberty group. Subsequently, the coordinated Grad-CAM images were averaged and superimposed on the representative images. The average Grad-CAM images acquired from the single-, double-, and triple-task models based on the AA-MTL architecture are depicted in Fig. 5.

Below is the link to the electronic supplementary material.

Conflicts of Interest

The authors report no conflicts of interest.