Population pharmacokinetic-pharmacodynamic analysis involves nonlinear hierarchical modelling where the mean response in a population and the variability in response from different sources are studied. It generally consists of two model hierarchies: a model for residual error and a model for heterogeneity termed between subject variance (BSV). The overall variability in a parameter within a population termed population parameter variance (PPV) consists of within subject variance (WSV) and BSV. Both these variances can further be split into random and predictable components. The predictable component of BSV (termed BSVP) is explained by covariates, individual characteristics e.g. weight. As BSVP increases, the remaining unpredictable (or random) between subject variability (BSVR) decreases since BSV = BSVP + BSVR, and BSV is a constant in any given data set. Since BSV and BSVR are estimated from the base and full covariate models, respectively, then BSVP = BSV-BSVR. The aim of this study was to explore the hypothesis, that a significant covariate may not always decrease BSVR. The specific aims were: (1) to explore circumstances where BSVR may not be reduced when adding a significantly correlated covariate and (2) to explore whether specific models for covariates may eliminate this anomaly when assessing BSVR. Simulations were performed using MATLAB (2011a) and estimation using NONMEM (ver 7.2) with FOCE and INTERACTION. A 1-compartment intravenous bolus PK model was used for simulation following a single unit dose (d = 1). The BSV of clearance [BSV(CL)] was described according to a log-normal distribution model with mean zero and variance ω². An additive random unexplained variability was assumed. Initially, we show through a simple simulation that BSVR can increase when a significantly correlated covariate is added to the model. We follow this with five simulation scenarios, A to E, that have various levels of correlation between the continuous covariate (Z) and CL ranging from 0 to 100 %. Each simulated scenario was replicated 100 times and estimated by a base model (i.e. without covariate addition) and six covariate models (M1-M6) which included non-nested (M1), nested (M2), and two types of interaction models for each of M1 and M2; non-nested interaction (M3, M5), nested interaction (M4, M6). Initially, through a motivating example we show that BSVR may not reduce even when there is 50 % correlation between the covariate Z and CL. It was found that with 0 % correlation M1, the non-nested covariate model (NNCM) resulted in negative BSVP (inflated BSVR) whereas M2, the nested covariate model (NCM), resulted in a calculated BSVP of zero. NNCM (M1) shows negative BSVP (BSVR > BSV) with correlation as high as 50 % and this model needs a minimum of 75 % correlation to show a positive BSVP. NCM (M2) shows positive but downwardly biased BSVP with 25, 50 and 75 % correlations. However, inclusion of a covariate-eta interaction term for both types of covariate models resulted in greater BSVP for 25, 50 and 75 % correlation scenarios compared to NNCM and NCM respectively. For 100 % correlation, it was found that covariate-eta interaction models show the same BSVP as the models without the interaction term, i.e. under perfect positive correlation all models perform similarly and correctly. It was found that a significantly correlated covariate may not reduce BSVR and in fact it may inflate the BSVR due to statistical misspecification of the covariate model. Incorporating statistical models that account for the covariate-eta interaction may be useful diagnostically in identifying the variability explained by covariates.