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 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 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.