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
where is a discriminability parameter for indicator , is a threshold parameter, and is the cumulative probability that boy with skeletal age is assigned a more mature grade than .
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