LAUSR

We consider Bayesian estimation of a hierarchical linear model (HLM) from partially observed data, assumed to be missing at random, and small sample sizes. A vector of continuous covariates C$$ C $$ includes cluster‐level partially observed covariates with interaction effects. Due to small sample sizes from 37 patient–physician encounters repeatedly measured at four time points, maximum‐likelihood estimation is suboptimal. Existing Gibbs samplers impute missing values of C$$ C $$ by a Metropolis algorithm using proposal densities that have constant variances while the target posterior distributions have nonconstant variances. Therefore, these samplers may not ensure compatibility with the HLM and, as a result, may not guarantee unbiased estimation of the HLM. We introduce a compatible Gibbs sampler that imputes parameters and missing values directly from the exact posterior distributions. We apply our Gibbs sampler to the longitudinal patient–physician encounter data and compare our estimators with those from existing methods by simulation.