Consider the conventional multilevel model Y=C[gamma]+Zu+e where [gamma] represents fixed effects and (u,e) are multivariate normal random effects. The continuous outcomes Y and covariates C are fully observed with a subset Z of C. The parameters are [theta]=([gamma],var(u),var(e)). Dempster, Rubin and Tsutakawa (1981) framed the estimation as a missing data problem, where (Y,u) are the complete data and the random effects u are conceived as missing data. Viewed in this way, the Expectation-Maximization (EM) algorithm has proven to be a natural and popular approach to estimation. However, when C is partially observed or subject to measurement error, it is natural to formulate a multilevel model for C that includes random effects, [nu]. In this article, we extend this thinking to allow estimation of the joint distribution of data Y=(Y,C)=(U[subscript o],U[subscript m]) and random effects b=(u,v) from observed data Y[subscript o]=(Y[subscript o],C[subscript o]) and to generate multiple imputations of missing data (Y[subscript m],b) based on the estimated distribution under the assumption that the data Y are missing at random. This approach contributes to the literature on multiple imputation in three ways: (a) it allows random effects [nu] to be conceived as latent covariates, thus addressing measurement errors of C; (b) it allows non-linearities, including random coefficients, interaction effects, and other polynomial effects involving partially observed covariates; (c) it imputes (Y[subscript m],b) using two-step importance sampling. In these cases, the joint distribution of Y is not analytically tractable even if the analytic multilevel model of interest to the analyst follows a multivariate normal distribution. We prove that our method of maximizing the likelihood and imputing missing data ensures compatibility of the non-normal joint distribution with the analytic normal theory multilevel model via provisionally known random effects. We present and evaluate a sufficient condition under which the produced imputations are compatible with the analytic model. [This paper will be published in the "Journal of Educational and Behavioral Statistics."]