Bones: Latent Trait Model for Multiple Ordered Categorical Responses¶
An example from OpenBUGS [43], Roche et al. [81], and Thissen [92] 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 = 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 => 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}(5)
UnivariateDistribution[
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 = [
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