In twin studies, a widespread method of estimating variance components is through structural equation modelling (SEM). For continuous traits with normal distributions, this is a flexible approach in that it is able to accommodate all linear models and allows for testing of equality of means, variances, covariances and variance components across subpopulations. However, with more elaborate models with discrete or categorical observed variables, SEM maximum likelihood (ML) estimation or ML procedures for estimating generalised linear mixed models such as GLAMM (Rabe-Hesketh and Skrondal 2005) soon reach computational boundaries. An alternative method is Bayesian statistical modelling with Markov chain Monte Carlo (MCMC) estimation algorithms (see also Eaves et al. 2005).