A Bayesian Approach to Multistage Fitting of the Variation of the Skeletal Age Features

Abstract

Accurate assessment of skeletal maturity is important clinically. Skeletal age assessment is usually based on features encoded in ossification centers. Therefore, it is critical to design a mechanism to capture as much as possible characteristics of features. We have observed that given a feature, there exist stages of the skeletal age such that the variation pattern of the feature differs in these stages. Based on this observation, we propose a Bayesian cut fitting to describe features in response to the skeletal age. With our approach, appropriate positions for stage separation are determined automatically by a Bayesian approach, and a model is used to fit the variation of a feature within each stage. Our experimental results show that the proposed method surpasses the traditional fitting using only one line or one curve not only in the efficiency and accuracy of fitting but also in global and local feature characterization.

1. Introduction

Hand X-ray shown in Figure 1 is commonly used for skeletal age assessment in pediatric radiology. A discrepancy between skeletal maturity and the chronical age may indicate the presence of some abnormality in skeletal growth. This abnormality has been found to be related to various diseases such as endocrine disorders [1], metabolic/growth abnormalities [2], malformations and bone dysplasias [3], and gonadal dysgenesis [4]. Therefore, the assessment of skeletal maturity has become more and more important clinically. Clearly the accuracy in assessment is of the first concern.

Figure 1 

Hand X-ray used in skeletal age assessment.

Features encoded in ossification centers form the basis for assessment. If we know the exact characteristics of the features with regard to different stages of ages, we can do the best job on assessment. In reality, one needs a mechanism to capture such characteristics of features. Given data of a feature with respect to skeletal ages, a simple and common approach is to fit a line or a curve, which in turn is used for future prediction of new patients or assisting radiologists to understand the variation rules of the feature.

For instance, Figure 2(a) shows the variation of a ratio feature [5, 6] in vertical axis with regard to the increasing skeletal age along the horizontal axis from newborn to 19 year old boys. (More details on this ratio are provided in Section 3.2.) Here in the figure, a single line is used for fitting the values of the feature. Obviously, a line is not enough to capture the characteristic of the values of the feature. A quadratic curve, shown in Figure 2(c), does not do a good job either. Fitting a more complex curve does not seem to be a feasible approach. This is because sometimes there are available only a small amount of data which could restrict the learning of complex curves, and local properties (with respect to the time) of the feature are often lost when fitting a global complex curve, and thus leading to inaccurate future prediction.

(a)

(b)

(c)

(d)

Figure 2 

Examples of fitting the variation of the ratio feature. The horizontal axis represents the skeletal age and the vertical axis corresponds to the values of the feature.

In this paper, we propose to fit the variation of features of the skeleton age via a multistage fitting approach. With our approach, we divide the skeletal age axis into several stages or phases, and within each stage, a relative simple model (line or curve) is employed for the purpose of fitting. Usually, the variation of a feature does not follow a simple rule when skeletal age increases. Instead, it often shows different variation patterns among different stages of age. As shown in Figures 2(b) and 2(d), multistage fitting not only can capture the entire pattern of feature variation but also carry the local properties regarding the skeletal age. A critical question is then, how does one determine the appropriate positions to separate the stages? The proposed Bayesian cut in this paper provides an answer via a Bayesian approach.

The rest of the paper is organized as follows. In Section 2, we describe our models for fitting, where the Bayesian cut is introduced. In Section 3, we present our experimental results on multi-stage fitting for artificial and real data. We conclude our paper in Section 4.

2. The Proposed Method

In this section, we first describe our proposed method for a simple case and then extend it to a general scenario.

Given a sequence of values 𝑓₁,𝑓₂,…,𝑓_𝑛, which denotes the skeletal age 𝑓 in an ascending order, consider the linear relationship between 𝑓 and one feature 𝑦 found in the hand X-ray (e.g., length of digit). Usually, such a linear relationship varies as the skeletal age increases. That is, one linear form established for one interval of the skeletal age may not hold for the next interval, where a different linear form should be used. The time where two linear forms differ is called a change point. Our model that takes into account linear relationships and change points is stated as follows: 𝑦_𝑖=𝛽₁₁+𝛽₁₂𝑓_𝑖+𝜖_1𝑖,𝑖=1,…,𝑡₁(𝑡₀𝑦=0),_𝑖=𝛽₂₁+𝛽₂₂𝑓_𝑖+𝜖_2𝑖,𝑖=𝑡₁+1,…,𝑡₂,⋮𝑦_𝑖=𝛽_𝑘1+𝛽_𝑘2𝑓_𝑖+𝜖_𝑘𝑖,𝑖=𝑡_𝑘−1+1,…,𝑡_𝑘(𝑡_𝑘=𝑛), (1) where 𝑡₁,…,𝑡_𝑘−1 (correspondingly 𝑓₁,…,𝑓_𝑘−1) indicate the sequential change points, 𝑡_𝑗−𝑡_𝑗−1≥3   (𝑗=1,…,𝑘), and 𝜖_𝑗𝑖 (for all 𝑖) are independent 𝑁(0,𝜎²_𝑗) and 𝜖_𝑗𝑖 (for all 𝑖,𝑗) are independent of each other. In the model, the parameters 𝛽_𝑗1,𝛽_𝑗2,𝜎²_𝑗,𝑡_𝑗 are all unknown, which will be estimated in light of the given data. The interval [𝑡_𝑗−𝑡_𝑗−1] represents the 𝑗th stage or phase, denoted by 𝑝ℎ_𝑗. The main task here is to estimate the times 𝑡_𝑗. Given the estimates of 𝑡_𝑗, the linear forms and the associated parameters can be obtained through the traditional regression technique. We note that the requirement 𝑡_𝑗−𝑡_𝑗−1≥3(𝑗=1,…,𝑘) is needed for estimation of the regression lines. When 𝑘=2, the model will be reduced to the two-phase regression with a single change point in [7].

The above model that uses only one dependent variable 𝑓 can be generalized to include multiple independent variables. This generalization leads to the following model: 𝑦_𝑖=𝛽^𝑇₁𝐟_𝐢+𝜖_1𝑖,𝑖=1,…,𝑡₁(𝑡₀𝑦=0),_𝑖=𝛽^𝑇₂𝐟_𝐢+𝜖_2𝑖,𝑖=𝑡₁+1,…,𝑡₂,⋮𝑦_𝑖=𝛽^𝑇_𝑘𝐟_𝐢+𝜖_𝑘𝑖,𝑖=𝑡_𝑘−1+1,…,𝑡_𝑘(𝑡_𝑘=𝑛), (2) where 𝐟_𝐢 is a 𝑝-dimensional vector of variables, 𝛽_𝑗(𝑗=1,…,𝑘) is a 𝑝-dimensional vector of parameters, 𝑡_𝑗−𝑡_𝑗−1≥𝑝+1, and 𝜖_𝑗𝑖 are as the same as before. We refer 𝑝 as the cardinality of the input vector 𝐟_𝐢, denoted by 𝐶(𝐟_𝐢), and the number of sample points in 𝑝ℎ_𝑗 as the cardinality of [𝑡_𝑗−𝑡_𝑗−1], denoted by 𝐶(𝑝ℎ_𝑗). We note that though linear regression is used for each phase in model (2), this model certainly encompasses other nonlinear cases such as polynomial forms.

We now describe a Bayesian approach to estimate the change points. Denote (𝐟_{𝐭𝐣−1+1},…,𝐟_𝐭𝐣)^𝑇 by 𝐹_𝑗, (𝐹^𝑇₁,…,𝐹^𝑇_𝑘)^𝑇 by 𝐹, (𝑦_{𝑡𝑗−1+1},…,𝑦_𝑡𝑗)^𝑇 by 𝐲_𝐣, (𝐲₁^𝑇,…,𝐲_𝐤^𝑇)^𝑇 by 𝐲, and (𝑡₁,…,𝑡_𝑘−1) by 𝐭. For simplicity, we assume the noninformative or uniform prior for 𝛽_𝑗(𝑗=1,…,𝑘), ln(𝜎_𝑗²) and 𝐭. Noninformative priors are used when information about parameters is completely unknown or when proper priors such as conjugate priors do not apply. (For a vigorous discussion on the choice of priors, see [8].) We can show the following main result (see the Appendix). Given the data 𝐲 and the uniform prior for 𝛽_𝑗(𝑗=1,…,𝑘), ln(𝜎_𝑗²) and 𝐭, where the number 𝑘 is predetermined, the posterior probability that change points occur at 𝐭 is 𝑝(𝐭∣𝐲)=𝐽2^{(𝑛−𝑘𝑝)/2}_𝑗|𝐹_𝑗^𝑇𝐹_𝑗|^−1/2𝑡×Γ(^𝑗−𝑡_𝑗−1−𝑝2)𝑆^{−(𝑡𝑗−𝑡𝑗−1}_𝑗^−𝑝)/2, (3) where 𝐽=(∑_𝐭2^{(𝑛−𝑘𝑝)/2}∏_𝑗|𝐹^𝑇_𝑗𝐹_𝑗|^−1/2Γ((𝑡_𝑗−𝑡_𝑗−1−𝑝)/2)×𝑆^{−(𝑡𝑗−𝑡𝑗−1}_𝑗^−𝑝)/2)⁻¹, and 𝑆_𝑗=(𝐲_𝐣−𝐹_𝑗𝛽_𝑗)^𝑇(𝐲_𝐣−𝐹_𝑗𝛽_𝑗) with 𝛽_𝑗=(𝐹^𝑇_𝑗𝐹_𝑗)⁻¹𝐹^𝑇_𝑗𝐲_𝐣 denoting the least-squares estimator of 𝛽_𝑗. Using this result, we estimate 𝐭 by 𝑡^∗ at which 𝑝(𝐭∣𝐲) has its maximum, that is, 𝑡^∗=argmax_t𝑝(𝐭∣𝐲). We call 𝑡^∗ the Bayesian cut, and the value 2^{(𝑛−𝑘𝑝)/2}∏_𝑗|𝐹^𝑇_𝑗𝐹_𝑗|^−1/2Γ((𝑡_𝑗−𝑡_𝑗−1−𝑝)/2)𝑆^{−(𝑡𝑗−𝑡𝑗−1}_𝑗^−𝑝)/2 the proportional posterior (𝑝𝑝).

3. Experiments

In this section, we perform the Bayesian cut on two data sets: one is synthesized and the other is real. We use the synthesized data for performance evaluation in terms of recovery of changing points. The real data are used to discover the Bayesian cut and describe the feature in a multistage way which has more accurate prediction of the skeletal age compared with fitting by a single line or curve. Both linear and nonlinear regression are used for comparison. For convenience, we call the fitting with a single line or curve the single fitting and the fitting with the Bayesian cut the Bayesian cut fitting.

3.1. Synthesized Data

We consider five cases or models describing the relationship between the dependent and independent variables. These are shown in Table 1 where the input vector 𝐟_𝐢 for models 𝑚₁, 𝑚₂, 𝑚₃, 𝑚₄, and 𝑚₅ is (1,𝑓_𝑖)^𝑇, (1,𝑓_𝑖,𝑓²_𝑖)^𝑇, (1,𝑓_𝑖,𝑓²_𝑖,𝑓³_𝑖)^𝑇, (1,𝑓_𝑖,𝑓²_𝑖,𝑓³_𝑖,𝑓⁴_𝑖)^𝑇, and (1,𝑓_𝑖,𝑓²_𝑖,𝑓³_𝑖,𝑓⁴_𝑖,𝑓⁵_𝑖)^𝑇, respectively. The data are generated according to the setting given in Table 2. Specifically, 𝛽_𝑗𝑖 is randomly chosen from (−5.0, 5.0). 𝜖_𝑗𝑖 is generated from a normal distribution with mean 0 and variance 𝜎²_𝑗 randomly selected from (0,5^{𝐶(𝐟𝐢)−1}). The number of sample points of the 𝑗th phase 𝐶(𝑝ℎ_𝑗) is randomly selected from the set {(𝐶(𝐟_𝐢)+1),…,(𝐶(𝐟_𝐢)+1)+𝑠}, where 𝑠 is predetermined. 𝑓_𝑖 takes the value of 𝑖 for 𝑖=1,2,…,𝑡_𝑘. Note that we use a variable bound for 𝜎²_𝑗 for taking into account the influence of the highest degree of the polynomial. Also, we use the variable number of sample points for each phase by introducing unbalance and scalability factors such that the performance evaluation will be more objective.To present a quantity on the performance of the Bayesian cut, we use the metric absolute deviation (AD), defined as ∑AD=^𝑗|𝑡^∗_𝑗−𝑡_𝑗|𝑘−1,𝑗=1,…,𝑘−1, (4) where 𝑡^∗_𝑗 represents the 𝑗th element of 𝐭^∗ (the Bayesian cut). Intuitively, the smaller AD is, the closer is the Bayesian cut 𝐭^∗ to the true change points 𝐭.

Table 1 

Models for testing the performance of the Bayesian cut.

Table 2 

Experimental setting.

Table 3 shows the AD values. They are obtained by ranging 𝑘 from 2 to 4 and 𝑠 from 1 to 10. For given 𝑘, 𝑠, and a given model, 50 trials are performed to generate data, leading to 50 datasets {(𝐹,𝐲)}. We find the Bayesian cut 𝐭^∗ for each (𝐹,𝐲) and a given model. The final AD score is obtained by averaging the 50 runs.

Table 3 

AD scores for models in Table 1.

Our findings can be summarized as follows. Regardless of linear or nonlinear regression, the Bayesian cut performs well with low AD scores. Introducing the unbalance and scalability factors does not deteriorate the performance of the Bayesian cut significantly. The Bayesian cut scales well when the number of change points increases.

3.2. Real Data

In this part, we apply the Bayesian cut fitting to some real data from our database shown in Table 4. This table describes feature values with regard to the increasing skeletal age that ranges from newborn to 19-year-old boys (shown in column 1) labeled by radiology experts. In order to obtain features independent of the size and the length of digits, two ratio features are used according to the paper [5]. One is 𝐿₁/𝐿₂, the ratio of the length of distal phalanx 𝐿₁ to that of middle phalanx 𝐿₂ of the middle digit, and the other is 𝐿₂/𝐿₃, the ratio of the length of middle phalanx 𝐿₂ to that of proximal phalanx 𝐿₃. See Figure 3 for illustration of 𝐿₁, 𝐿₂, and 𝐿₃. These two features correspond to columns 2 and 3 which are generated in the light of the algorithm in [6]. Columns 4 and 5 represent normalized values of 𝐿₁/𝐿₂ and 𝐿₂/𝐿₃, respectively. This normalization is done according to (𝑥−𝜇)/𝜎, where 𝜇 is the expectation of 𝑥 and 𝜎 is the variance. In our experiments, only normalized values are used. Figure 4 shows some of the Bayesian cut fitting, where features 𝑛(𝐿₁/𝐿₂) and 𝑛(𝐿₂/𝐿₃) are used, models describing the relationship between the feature and the skeletal age are 𝑚₁ and 𝑚₂ from Table 1, and 𝑘 takes values of 2, 3, and 4. In Figure 4, the horizontal axis represents the age and the horizontal axis indicates the feature. For model 𝑚₁, the blue straight line across the entire age range is from the single (line) fitting. For model 𝑚₂, the blue curve across the entire age range is from the single (quadratic) fitting. All red (broken) lines are from the Bayesian cut fitting.

Table 4 

Some features of the skeletal age.

Figure 3 

Illustration Of 𝐿1, 𝐿2 and 𝐿3.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 4 

Illustration of the Bayesian cut fitting applied to the real data on features of the skeletal age.

4. Conlcusion

In this paper, we propose the Bayesian cut fitting to describe features in response to the skeletal age. In the semantic space derived by our approach, the axis of skeletal age is divided into meaningful stages, within each of which the variation pattern of a feature is consistent so that a traditional regression technique can apply to model the relationship between the skeletal age and the feature. Our approach is inspired by the observation that the variation pattern of a feature can differ in different periods of the skeletal age. A critical issue is to determine the times or change points when the variation pattern of a feature changes. This is handled by the Bayesian cut proposed in this paper. Simulations have been used to demonstrate the efficiency of the Bayesian cut fitting in terms of recovery of change points. The experiments on real data show that given a type of relationship (e.g., linear or quadratic) between the skeletal age and a feature, the Bayesian cut fitting surpasses the traditional single fitting when the consistency of the variation pattern (over the entire skeletal age range) of the feature is suspected. One major issue which is not addressed in this paper is the determination of 𝑘, the number of stages. Selection of 𝑘 depends on the given data and the practical need. We leave this as our future research work.

Appendix

A. Derivation of (3)

Proof. According to the Pythagorean theorem, we have the following likelihood 𝑙(𝛽_𝑗,𝜎²_𝑗1∣𝐲)∝(𝜎²_𝑗)^{(𝑡𝑗−𝑡𝑗−1)/2}1×exp{−2𝜎_𝑗²[(𝑦_{𝑡𝑗−1+1}−𝛽^𝑇_𝑗𝐟_{𝐭𝐣−1+1})²+⋯+(𝑦_𝑡𝑗−𝛽^𝑇_𝑗𝐟_𝐭𝐣)²∝1]}(𝜎²_𝑗)^{(𝑡𝑗−𝑡𝑗−1)/2}1×exp{−2𝜎_𝑗²[𝑆_𝑗+(𝛽_𝑗−𝛽_𝑗)^𝑇×𝐹^𝑇_𝑗𝐹_𝑗(𝛽_𝑗−𝛽_𝑗)]}, (A.1) where 𝑆_𝑗=(𝐲_𝐣−𝐹_𝑗𝛽_𝑗)^𝑇(𝐲_𝐣−𝐹_𝑗𝛽_𝑗) and 𝛽_𝑗=(𝐹^𝑇_𝑗𝐹_𝑗)⁻¹𝐹^𝑇_𝑗𝐲_𝐣. Since 𝜖_𝑗𝑖 are independent of each other, the likelihood function of 𝛽₁,…,𝛽_𝑘,𝜎²₁,…,𝜎²_𝑘,𝐭 is then 𝑙(𝛽₁,…,𝛽_𝑘,𝜎²₁,…,𝜎²_𝑘,𝐭∣𝐲)∝_𝑗1(𝜎²_𝑗)^{(𝑡𝑗−𝑡𝑗−1)/2}1×exp{−2𝜎²_𝑗[𝑆_𝑗+(𝛽_𝑗−𝛽_𝑗)^𝑇𝐹^𝑇_𝑗𝐹_𝑗(𝛽_𝑗−𝛽_𝑗)]}. (A.2) Due to the assumption of the uniform prior for 𝛽_𝑗, ln(𝜎²_𝑗) and 𝐭, we have 𝑝(𝛽₁,…,𝛽_𝑘,𝜎²₁,…,𝜎²_𝑘1,𝐭)∝𝜎²₁⋯𝜎²_𝑘. (A.3) Using (A.2) and (A.6), we have ∝𝑝(𝐭∣𝐲)_𝛽1⋯_{𝛽𝑘𝜎21}⋯_{𝜎2𝑘𝑗}1(𝜎²_𝑗)^{(𝑡𝑗−𝑡𝑗−1)/2+1}1×exp{−2𝜎_𝑗²[𝑆_𝑗+(𝛽_𝑗−𝛽_𝑗)^𝑇×𝐹^𝑇_𝑗𝐹_𝑗(𝛽_𝑗−𝛽_𝑗𝛽)]}×𝑑₁⋯𝑑𝛽_𝑘𝑑𝜎²₁⋯𝑑𝜎²_𝑘. (A.4)
Note that _𝛽𝑗1(𝜎²_𝑗)^{(𝑡𝑗−𝑡𝑗−1)/2+1}1×exp{−2𝜎_𝑗²[𝑆_𝑗+(𝛽_𝑗−𝛽_𝑗)^𝑇𝐹^𝑇_𝑗𝐹_𝑗(𝛽_𝑗−𝛽_𝑗)]}𝑑𝛽_𝑗=exp(−𝑆_𝑗/2𝜎²_𝑗)(𝜎²_𝑗)^{(𝑡𝑗−𝑡𝑗−1)/2+1}(2𝜋)^𝑝/2(2𝜎²_𝑗)^𝑝/2|𝐹^𝑇_𝑗𝐹_𝑗|^−1/2. (A.5) This equation exploits the fact exp{(𝐱−𝜇)^𝑇Σ⁻¹(𝐱−𝜇)}𝑑𝐱=(2𝜋)^𝑝/2|Σ|^1/2, (A.6) from the normal density for the 𝑝-dimensional random vector 𝑋1𝑓(𝐱)=(2𝜋)^𝑝/2|Σ|^1/2exp{(𝐱−𝜇)^𝑇Σ⁻¹(𝐱−𝜇)}𝑑𝐱, (A.7) where 𝜇 is the expected value of 𝑋 and Σ is the variance-covariance matrix of 𝑋.
Substituting (A.10) into (A.7), we have 𝑝(𝐭∣𝐲)∝_𝑗|𝐹^𝑇_𝑗𝐹_𝑗|^−1/2_𝜎2𝑗exp(−𝑆_𝑗/2𝜎²_𝑗)(𝜎²_𝑗)^{(𝑡𝑗−𝑡𝑗−1−𝑝)/2+1}𝑑𝜎²_𝑗. (A.8)
In addition, we have 𝑎exp(−2𝑥)𝑥^{−𝑚/2−1}𝑑𝑥=2^𝑚/2𝑚Γ(2)𝑎^−𝑚/2, (A.9) from the probability density function of 𝑋=𝑎𝑈𝑓(𝑥)=2^−𝑚/2𝑚[Γ(2)]⁻¹𝑎^𝑚/2𝑥^{−𝑚/2−1}𝑎exp(−2𝑥), (A.10) where the constant 𝑎>0 and 𝑈⁻¹∼𝜒²_𝑚.
By applying (A.9) to (A.11), we get 𝑝(𝐭∣𝐲)∝𝐽2^{(𝑛−𝑘𝑝)/2}_𝑗|𝐹^𝑇_𝑗𝐹_𝑗|^−1/2𝑡×Γ(^𝑗−𝑡_𝑗−1−𝑝2)𝑆^{−(𝑡𝑗−𝑡𝑗−1}_𝑗^−𝑝)/2, (A.11) where 𝐽=(∑_𝐭2^{(𝑛−𝑘𝑝)/2}∏_𝑗|𝐹^𝑇_𝑗𝐹_𝑗|^−1/2Γ((𝑡_𝑗−𝑡_𝑗−1−𝑝)/2)×𝑆^{−(𝑡𝑗−𝑡𝑗−1}_𝑗^−𝑝)/2)⁻¹. This completes the proof.

Acknowledgments

Dechang Chen was partially supported by the National Science Foundation grant CCF-0729080.