Bones: Latent Trait Model for Multiple Ordered Categorical Responses¶

An example from OpenBUGS [44], Roche et al. [82], and Thissen [93] 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 Distributed
@everywhere using Mamba

## Data
bones = Dict{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 => permutedims(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,
    (ncat, delta, theta, gamma, nChild, nInd) ->
      begin
        p = Array{Float64}(undef, 5)
        UnivariateDistribution[(
          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])) for i in 1:nChild, j in 1:nInd
        ]
      end,
    false
  ),

  theta = Stochastic(1,
    () -> Normal(0, 100)
  )

)


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


## Sampling Scheme
scheme = [MISS(:grade),
          AMWG(:theta, 0.1)]
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.20640874 0.0023834028 0.0039448110 2737.81276
 theta[2]  1.37861692 0.25824308 0.0029819342 0.0058663965 1937.82503
 theta[3]  2.35227822 0.27998526 0.0032329913 0.0067161153 1737.93707
 theta[4]  2.90165730 0.29713320 0.0034309987 0.0078730921 1424.33353
 theta[5]  5.54427283 0.50242324 0.0058014839 0.0169090038  882.88350
 theta[6]  6.70804782 0.57206890 0.0066056827 0.0221532973  666.83738
 theta[7]  6.49138381 0.60154625 0.0069460578 0.0219158412  753.39330
 theta[8]  8.93701249 0.73636136 0.0085027686 0.0336199950  479.71875
 theta[9]  9.03585289 0.65172497 0.0075254717 0.0233182299  781.15561
theta[10] 11.93125529 0.69360918 0.0080091090 0.0282955741  600.88678
theta[11] 11.53686992 0.92271657 0.0106546132 0.0493587234  349.46912
theta[12] 15.81482824 0.54261736 0.0062656056 0.0210976666  661.48275
theta[13] 16.93028146 0.72458739 0.0083668145 0.0323302069  502.30161

Quantiles:
              2.5%       25.0%       50.0%       75.0%       97.5%
 theta[1] -0.11215555  0.19557824  0.33881555  0.45840506  0.7174563
 theta[2]  0.91705346  1.19969433  1.36116575  1.53751273  1.9466119
 theta[3]  1.78287586  2.17136780  2.35623189  2.53035766  2.9211580
 theta[4]  2.32940825  2.69621746  2.89121336  3.10758151  3.4945343
 theta[5]  4.59142954  5.20543314  5.53392246  5.86525435  6.5867245
 theta[6]  5.56649066  6.30983797  6.70666338  7.09569168  7.8229872
 theta[7]  5.38663728  6.07628064  6.46533033  6.88636840  7.7051374
 theta[8]  7.47304526  8.43125608  8.96072241  9.45344704 10.2856733
 theta[9]  7.80477915  8.60559136  9.01498109  9.46962522 10.3024722
theta[10] 10.64129157 11.48379528 11.89611699 12.37737647 13.3873043
theta[11]  9.83558611 10.88717498 11.49029895 12.15757004 13.4263451
theta[12] 14.79250437 15.45889470 15.79840132 16.15824313 16.9593310
theta[13] 15.61843069 16.42289725 16.90719268 17.41900248 18.3895761