Bones: latent trait model for multiple ordered categorical responses

The concept of skeletal age (SA) arises from the idea that individuals mature at different rates: for any given chronological age (CA), the average SA in a sample of individuals should equal their CA, but with an inter-individual spread which reflects the differential rate of maturation. Roche et al (1975) have developed a model for predicting SA by calibrating 34 indicators (items) of skeletal maturity which may be observed in a radiograph. Each indicator is categorized with respect to its degree of maturity: 19 are binary items (i.e. 0 = immature or 1 = mature); 8 items have 3 grades (i.e. 0 = immature; 1 = partially mature; 2 = fully mature); 1 item has 4 ordered grades and the remaining 6 items have 5 ordered grades of maturity. Roche et al. calculated threshold parameters for the boundarys between grades for each indicator. For the binary items, there is a single threshold representing the CA at which 50% of individuals are mature for the indicator. Three-category items have 2 threshold parameters: the first corresponds to the CA at which 50% of individuals are either partially or fully mature for the indicator; the second is the CA at which 50% of individuals are fully mature. Four and five-category items have 3 and 4 threshold parameters respectively, which are interpreted in a similar manner to those for 3-category items. In addition, Roche et al. calculated a discriminability (slope) parameter for each item which reflects its rate of maturation. Part of this data is shown below. Columns 1--4 represent the threshold parameters (note the use of the missing value code NA to `fill in' the columns for items with fewer than 4 thresholds); column 5 is the discriminability parameter; column 6 gives the number of grades per item.

   
Thissen (1986) (p.71) presents the following graded radiograph data on 13 boys whose chronological ages range from 6 months to 18 years. (Note that for ease of implementation in BUGS we have listed the items in a different order to that used by Thissen):


Some items have missing data (represented by the code NA in the table above). This does not present a problem for BUGS: the missing grades are simply treated as unknown parameters to be estimated along with the other parameters of interest such as the SA for each boy.

Thissen models the above data using the logistic function. For each item j and each grade k, the cumulative probability Q_jk that a boy with skeletal age θ is assigned a more mature grade than k is given by

   logitQ_jk = δ_j(θ - γ_jk)
   
where δ_j is the discriminability parameter and the γ_jk are the threshold parameters for item j. Hence the probability of observing an immature grade (i.e. k =1) for a particular skeletal age θ is p_j,1 = 1 - Q_j,1. The probability of observing a fully mature grade (i.e.k = K_j, where K_j is the number of grades for item j is p_j,Kj = Q_{j,Kj -1}. For items with 3 or more categories, the probability of observing an itermediate grade is p_j,k = Q_j,k-1 - Q_j,k (i.e. the difference between the cumulative probability of being assigned grade k or more, and of being assigned grade k+1 or more).

The BUGS language for this model is given below. Note that the θ_i for each boy i is assigned a vague, independent normal prior theta[i] ~ dnorm(0.0, 0.001). That is, each boy is treated as a separate problem with is no `learning' or `borrowing strength' across individuals, and hence no hierachical structure on the θ_i's.


BUGS language for bones example

   model
   {
      for (i in 1 : nChild) {
         theta[i] ~ dnorm(0.0, 0.001)
         for (j in 1 : nInd) {
   # Cumulative probability of > grade k given theta
            for (k in 1: ncat[j] - 1) {
               logit(Q[i, j, k]) <- delta[j] * (theta[i] - gamma[j, k])
            }
         }

   # Probability of observing grade k given theta
         for (j in 1 : nInd) {
            p[i, j, 1] <- 1 - Q[i, j, 1]
            for (k in 2 : ncat[j] - 1) {
               p[i, j, k] <- Q[i, j, k - 1] - Q[i, j, k]
            }
            p[i, j, ncat[j]] <- Q[i, j, ncat[j] - 1]
            grade[i, j] ~ dcat(p[i, j, 1 : ncat[j]])
         }
      }
   }
   
We note a couple of tricks used in the above code. Firstly, the variable p has been declared as a 3-way rectangular array with the size of the third dimension equal to the maximum number of possible grades (i.e.5) for all items (even though items 1--28 have fewer than 5 categories). The statement

   grade[i, j] ~ dcat(p[i, j, 1 :ngrade[j]])

is then used to select the relevant elements of p[i,j, ] for item j, thus ignoring any `empty' spaces in the array for items with fewer than the maximum number of grades. Secondly, the final section of the above code includes a loop indexed as follows
Results