We consider Bayesian estimation of a hierarchical linear model (HLM) from small sample sizes. The continuous response Y and covariates C are partially observed and assumed missing at random. With C having linear effects, the HLM may be efficiently estimated by available methods. When C includes cluster-level covariates having interactive or other nonlinear effects given small sample sizes, however, maximum likelihood estimation is suboptimal, and existing Gibbs samplers are based on a Bayesian joint distribution compatible with the HLM, but impute missing values of C by a Metropolis algorithm via a proposal density having a constant variance while the target conditional distribution has a non-constant variance. Therefore, the samplers are not guaranteed to be compatible with the joint distribution and, thus, always produce unbiased estimation of the HLM. We introduce a compatible Gibbs sampler that imputes parameters and missing values directly from the exact conditional distributions. We illustrate analysis of repeated measurements nested within each of 37 patient-physician encounters by our sampler, and compare our estimators with those of existing methods by simulation.