Multiple imputation methods are widely used in practice for missing data. An important consideration for a multiple imputation method is the choice of an imputation model which generates the imputations for each missing value, especially when the missing rate is not low. Mixed effects models are commonly used for modelling longitudinal data which exhibit large between-individual variations. In this case, a good imputation model should generate imputations at the individual level to incorporate the large between-individual variations. In this article, we propose a multiple imputation method for nonlinear mixed effects models with missing responses. We consider an iterative linearization method where the imputations are generated based on a “working” linear mixed effects model. We evaluate the proposed method via simulations and apply the method to a real dataset.