Bones: Latent Trait Model for Multiple Ordered Categorical Responses

An example from OpenBUGS [38], Roche et al. [65], and Thissen [72] concerning skeletal age in 13 boys predicted from 34 radiograph indicators of skeletal maturity.

Model

Skeletal ages are modelled as

\operatorname{logit}(Q_{i,j,k}) &= \delta_j (\theta_i - \gamma_{j,k}) \quad\quad i=1,\ldots,13; j=1,\ldots,34; k=1,\ldots,4 \\
\theta_i &\sim \text{Normal}(0, 100),

where \delta_j is a discriminability parameter for indicator j, \gamma_{j,k} is a threshold parameter, and Q_{i,j,k} is the cumulative probability that boy i with skeletal age \theta_i is assigned a more mature grade than k.

Analysis Program

using Mamba

## Data
bones = (Symbol => Any)[
  :gamma => reshape(
    [ 0.7425,     NaN,     NaN,    NaN,  10.2670,     NaN,     NaN,     NaN,
     10.5215,     NaN,     NaN,    NaN,   9.3877,     NaN,     NaN,     NaN,
      0.2593,     NaN,     NaN,    NaN,  -0.5998,     NaN,     NaN,     NaN,
     10.5891,     NaN,     NaN,    NaN,   6.6701,     NaN,     NaN,     NaN,
      8.8921,     NaN,     NaN,    NaN,  12.4275,     NaN,     NaN,     NaN,
     12.4788,     NaN,     NaN,    NaN,  13.7778,     NaN,     NaN,     NaN,
      5.8374,     NaN,     NaN,    NaN,   6.9485,     NaN,     NaN,     NaN,
     13.7184,     NaN,     NaN,    NaN,  14.3476,     NaN,     NaN,     NaN,
      4.8066,     NaN,     NaN,    NaN,   9.1037,     NaN,     NaN,     NaN,
     10.7483,     NaN,     NaN,    NaN,   0.3887,  1.0153,     NaN,     NaN,
      3.2573,  7.0421,     NaN,    NaN,  11.6273, 14.4242,     NaN,     NaN,
     15.8842, 17.4685,     NaN,    NaN,  14.8926, 16.7409,     NaN,     NaN,
     15.5487, 16.8720,     NaN,    NaN,  15.4091, 17.0061,     NaN,     NaN,
      3.9216,  5.2099,     NaN,    NaN,  15.4750, 16.9406, 17.4944,     NaN,
      0.4927,  1.3556,  2.3016,  3.2535,  1.3059,  1.8793,  2.4970,  3.2306,
      1.5012,  1.8902,  2.3689,  2.9495,  0.8021,  2.3873,  3.9525,  5.3198,
      5.0022,  6.3704,  8.2832, 10.4988,  4.0168,  5.1537,  7.1053, 10.3038],
     4, 34)',
   :delta =>
     [2.9541, 0.6603, 0.7965, 1.0495, 5.7874, 3.8376, 0.6324,
      0.8272, 0.6968, 0.8747, 0.8136, 0.8246, 0.6711, 0.978,
      1.1528, 1.6923, 1.0331, 0.5381, 1.0688, 8.1123, 0.9974,
      1.2656, 1.1802, 1.368,  1.5435, 1.5006, 1.6766, 1.4297,
      3.385,  3.3085, 3.4007, 2.0906, 1.0954, 1.5329],
   :ncat =>
     [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
      3, 3, 3, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5, 5, 5],
   :grade => reshape(
     [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,2,1,1,2,1,1,
      2,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,3,1,1,2,1,1,
      2,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,4,3,3,3,1,1,
      2,1,1,1,2,2,1,1,1,1,1,1,NaN,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,4,5,4,3,1,1,
      2,1,1,1,2,2,1,1,2,1,1,1,1,1,1,1,2,1,1,3,2,1,1,1,1,1,3,1,5,5,5,4,2,3,
      2,1,1,1,2,2,1,2,1,1,1,1,1,2,1,1,2,NaN,1,3,2,1,1,1,1,1,3,1,5,5,5,5,3,3,
      2,1,1,1,2,2,1,1,1,NaN,NaN,1,1,1,1,1,2,NaN,1,3,3,1,1,1,1,1,3,1,5,5,5,5,3,3,
      2,1,2,2,2,2,2,2,1,NaN,NaN,1,2,2,1,1,2,2,1,3,2,1,1,1,1,1,3,1,5,5,5,5,3,4,
      2,1,1,2,2,2,2,2,2,1,1,1,2,1,1,1,2,1,1,3,3,1,1,1,1,1,3,1,5,5,5,5,4,4,
      2,1,2,2,2,2,2,2,2,1,1,1,2,2,2,1,2,NaN,2,3,3,1,1,1,1,1,3,1,5,5,5,5,5,5,
      2,1,NaN,2,2,2,NaN,2,2,1,NaN,NaN,2,2,NaN,NaN,2,1,2,3,3,NaN,1,NaN,1,1,3,1,5,5,5,5,5,5,
      2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,1,NaN,2,1,3,2,5,5,5,5,5,5,
      2,2,2,2,2,2,2,2,2,2,NaN,2,2,2,2,2,2,2,2,3,3,3,NaN,2,NaN,2,3,4,5,5,5,5,5,5],
     34, 13)',
  :nChild => 13,
  :nInd => 34
]


## Model Specification

model = Model(

  grade = Stochastic(2,
    @modelexpr(ncat, delta, theta, gamma, nChild, nInd,
      begin
        p = Array(Float64, 5)
        Distribution[
          begin
            n = ncat[j]
            p[1] = 1.0
            for k in 1:n-1
              Q = invlogit(delta[j] * (theta[i] - gamma[j,k]))
              p[k] -= Q
              p[k+1] = Q
            end
            Categorical(p[1:n])
          end
          for i in 1:nChild, j in 1:nInd
        ]
      end
    ),
    false
  ),

  theta = Stochastic(1,
    :(Normal(0, 100))
  )

)


## Initial Values
inits = [
  [:grade => bones[:grade], :theta => [0.5,1,2,3,5,6,7,8,9,12,13,16,18]],
  [:grade => bones[:grade], :theta => [1,2,3,4,5,6,7,8,9,10,11,12,13]]
]


## Sampling Scheme
scheme = [MISS([:grade]),
          AMWG([:theta], fill(0.1, bones[:nChild]))]
setsamplers!(model, scheme)


## MCMC Simulations
sim = mcmc(model, bones, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)

Results

Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750

Empirical Posterior Estimates:
              Mean        SD      Naive SE       MCSE        ESS
 theta[1]  0.32603385 0.2064087 0.0023834028 0.005562488 1376.94938
 theta[2]  1.37861692 0.2582431 0.0029819342 0.007697049 1125.66422
 theta[3]  2.35227822 0.2799853 0.0032329913 0.008869417  996.50667
 theta[4]  2.90165730 0.2971332 0.0034309987 0.009730398  932.48357
 theta[5]  5.54427283 0.5024232 0.0058014839 0.022915189  480.72039
 theta[6]  6.70804782 0.5720689 0.0066056827 0.029251931  382.46138
 theta[7]  6.49138381 0.6015462 0.0069460578 0.030356237  392.68306
 theta[8]  8.93701249 0.7363614 0.0085027686 0.042741328  296.81508
 theta[9]  9.03585289 0.6517250 0.0075254717 0.031204552  436.20719
theta[10] 11.93125529 0.6936092 0.0080091090 0.039048344  315.51820
theta[11] 11.53686992 0.9227166 0.0106546132 0.065584299  197.94151
theta[12] 15.81482824 0.5426174 0.0062656056 0.028535483  361.59041
theta[13] 16.93028146 0.7245874 0.0083668145 0.043348628  279.40285

Quantiles:
             2.5%       25.0%      50.0%       75.0%        97.5%
 theta[1] -0.1121555  0.1955782  0.33881555  0.45840506  0.717456283
 theta[2]  0.9170535  1.1996943  1.36116575  1.53751273  1.946611947
 theta[3]  1.7828759  2.1713678  2.35623189  2.53035766  2.921158047
 theta[4]  2.3294082  2.6962175  2.89121336  3.10758151  3.494534287
 theta[5]  4.5914295  5.2054331  5.53392246  5.86525435  6.586724493
 theta[6]  5.5664907  6.3098380  6.70666338  7.09569168  7.822987248
 theta[7]  5.3866373  6.0762806  6.46533033  6.88636840  7.705137414
 theta[8]  7.4730453  8.4312561  8.96072241  9.45344704 10.285673300
 theta[9]  7.8047792  8.6055914  9.01498109  9.46962522 10.302472161
theta[10] 10.6412916 11.4837953 11.89611699 12.37737647 13.387304288
theta[11]  9.8355861 10.8871750 11.49029895 12.15757004 13.426345115
theta[12] 14.7925044 15.4588947 15.79840132 16.15824313 16.959330990
theta[13] 15.6184307 16.4228972 16.90719268 17.41900248 18.389576101