Properties of the full random‐effect modeling approach with missing covariate data

During drug development, a key step is the identification of relevant covariates predicting between‐subject variations in drug response. The full random effects model (FREM) is one of the full‐covariate approaches used to identify relevant covariates in nonlinear mixed effects models. Here we explore the ability of FREM to handle missing (both missing completely at random (MCAR) and missing at random (MAR)) covariate data and compare it to the full fixed‐effects model (FFEM) approach, applied either with complete case analysis or mean imputation. A global health dataset (20 421 children) was used to develop a FREM describing the changes of height for age Z‐score (HAZ) over time. Simulated datasets (n = 1000) were generated with variable rates of missing (MCAR) covariate data (0%‐90%) and different proportions of missing (MAR) data condition on either observed covariates or predicted HAZ. The three methods were used to re‐estimate model and compared in terms of bias and precision which showed that FREM had only minor increases in bias and minor loss of precision at increasing percentages of missing (MCAR) covariate data and performed similarly in the MAR scenarios. Conversely, the FFEM approaches either collapsed at ≥$$ \ge $$ 70% of missing (MCAR) covariate data (FFEM complete case analysis) or had large bias increases and loss of precision (FFEM with mean imputation). Our results suggest that FREM is an appropriate approach to covariate modeling for datasets with missing (MCAR and MAR) covariate data, such as in global health studies.

identification of patients sub-populations with lower drug metabolism).Therefore, during model development the choice of which method to use for covariate identification should be considered carefully.
The methods for covariate identification can generally be classified into (i) covariate selection methods and (ii) full covariate model approaches.Selection methods, among which the stepwise procedures are the most common, use objective selection criteria (eg, a P-value) to identify the covariates that should be included in a model.Despite their recognized advantages, 2,3 these automated selection approaches have several weaknesses, mainly related to multiple testing issues, risk for selection bias, inflated type-I error, 4,5 and sensitivity to covariate correlations. 6,7On the other hand, full covariate model approaches have been suggested as an alternative to stepwise procedure to avoid selection bias and improve model inference. 7Briefly, with the full fixed-effects model (FFEM) approach, a set of non-correlated covariates is identified prior to model development; this full set is implemented to simultaneously estimate the fixed-effect coefficients for each covariate-parameter relation, together with the rest of the model parameters, thus avoiding the multiple testing issue.
Among the full model methods, an additional approach has recently been introduced 8 : the full random effects model (FREM), described in detail by Yngman et al. 9 Similarly to FFEM, also this approach uses a predefined set of covariates in a full model.However, in FREM covariates are considered as responses.This implies that while in FFEM the covariate effects are estimated by fixed-effects coefficients on the same parameter, in FREM they are captured in estimated covariances between individual parameters and covariates.This characteristic makes FREM able to better handle correlated covariates in the analysis.
Despite their specific advantages, both covariate selection methods and full covariate model approaches have a common challenge: handling missing covariate values.How to approach such a situation depends on whether covariate data are missing at random (MAR) or not.If covariate data are Missing Not At Random (MNAR), 10 a common method is to model the mechanisms determining the missing data.On the other hand, if covariate data are MAR (missing completely at random [MCAR] or MAR), other approaches may be used, such as including only subjects with non-missing covariate data (complete case analysis), mean or median imputation of missing data, 11 or more sophisticated methods such as multiple imputation (MI). 10,12In addition to these, FREM may be a new alternative: in this approach, the covariates are treated as responses and linked to the model for the parameters by jointly modeling their random effects and any missingness is handled by the maximum likelihood estimation. 12hile the FREM approach has been introduced, 8,9 and already applied in several published studies, [13][14][15] its ability to handle missing (MCAR and MAR) covariate data has never been fully characterized. 16n this article, the ability of FREM to handle missing covariate data has been investigated.To this aim, the FREM approach is first used to model covariate effects in a global health data set, pooled from 5 separate studies.The performance of FREM is then compared to either a FFEM with complete case analysis (FFEM CC) or a FFEM with mean imputation (FFEM ME).This comparison is made by performing simulation studies with variable degrees of missing covariate information (MCAR) (between 0% and 90%), then re-estimating the parameter coefficients either with one of the two types of FFEM or with FREM, and finally comparing the bias and precision of the three methods.In addition, handling of MAR scenarios conditioned on observed covariates and responses are investigated with 30% missing covariate information with both FREM and FFEM CC.
The article starts with a methods section where the FREM methodology is described in detail (Section 2).A data section follows where the global health data used in this analysis is described, together with the strategy for building the model using FREM (Section 3).After this, we describe the simulation study used to investigate the missing covariate data properties of the FREM (Section 4).The results section (Section 5) first presents the model developed to describe the global health data (Section 5.1) and then analyzes the ability of FREM to handle missing covariate data and compares it to FFEM (Section 5.2).In the final section, the developed model is examined and the results of applying FREM, the implications, and future perspectives are discussed (Section 6).

FULL RANDOM EFFECTS MODEL
The FREM is a method for estimating covariate effects in mixed-effects models.It selects a predefined set of covariates in a full model and it models each covariate as a separate response.A full covariance matrix between random effects for parameters and covariates is estimated together with the other model components.Coefficients for the covariate-parameter relations are obtained by calculating the ratio between covariance (combined parameter and covariate variability) and variance of each covariate.Unlike FFEM, the implementation of covariate effects in FREM does not rely on estimating multiple fixed-effects on the same parameter and it might therefore be an advantage (over FFEM) when estimating highly correlated covariates in FREM.Before introducing the equations used for a FREM, consider first a nonlinear mixed-effects model (NLMEM) for subject i and time point j without any covariates: where f (.) is a nonlinear function with respect to the structural parameters  par and the random effects corresponding to the inter-individual variability (IIV) b par,i ∼  ( 0, D par ) ; t ij is the independent variable (eg, time of observation); h(.) is a nonlinear function describing how the individual observations differ from the model, that is, the residual error function; and  ij is the random effect corresponding to the residual unexplained variability (RUV)  ij ∼  (0, Σ).In addition, D par and Σ are the covariance matrices respectively of the individual random effects (b par,i ) and the residual errors ( ij ).
A multi-response NLMEM is defined as a NLMEM with several types of responses.For example, in drug development modeling, one type may be PK observations that are driving an outcome variable (eg, drug effect response), which is in turn dependent on the drug concentration PK model.A multi-response NLMEM of K responses can be defined as: where each k response can have different parameters, functional forms, and independent variables.The FREM is a multi-response NLMEM and can thus consider several types of responses.Accordingly, all covariates are described by separate functions in which they are treated as responses, without any parameters in common with the structural model.Instead, the covariates are linked to the parameters of the structural model via the covariances of the individual random effects: with the full inter-individual covariance matrix D defined as: where  cov 1∶K and b cov 1∶K ,i ∼  (0, D cov ) are the parameters describing the K covariates;  cov,i is a residual with a variance (Σ cov ) that could be fixed to a value close to machine epsilon (Σ cov → ) if the covariates are assumed to be without measurement error.Equation ( 3) describes a single response; however, since FREM can model several responses, the complete set of responses (see Equation ( 2)) can be linked by the FREM approach.Of note, the K covariates are assumed to be time-invariant (hence the index i for the covariates).In fact, time-varying covariates (eg, food intake) are not handled using this definition of FREM and, if present, they should be included in the base model using fixed-effects.Or alternatively, be specified as a separate response including the variable time or implemented using a multilevel random effects models, where for example, the second level of random effects changes with each covariate observation in time.Therefore, in its simplest form, the model for the k covariate is described as: hence the b cov k ,i is estimated to be close to the difference between the observed covariate value and the mean of the covariate ( cov k ).
An important feature of FREM, especially relevant when assessing its performance in case of missing covariate data, is the flexibility for the covariate distribution.In the case of complete data, the FREM method does work well with regards to the underlying true covariate distribution and, in fact, both continuous and categorical covariates can be parametrized in the same way.This is not surprising since the b cov k ,i is the empirical Bayes estimate estimated with a negligible residual error.An example of this is presented below for the bi-modal SEX covariate.

Distribution of the bi-modal covariate SEX
Assume that the sex of 100 subjects is observed as a covariate and that 70 subjects are males and 30 are females.If males are encoded as 1 and females as 2, the mean of the covariate (or the estimated unbiased mean) is βsex ≈ 1.3.Furthermore, the random effect b sex,i for each subject i is estimated to bsex,i ≈ −0.3 for males and bsex,i ≈ 0.7 for females.Hence the underlying assumption of normality in the random effect b sex,i is not distorting the bi-modal nature of the sex covariate in the estimated covariate distribution.

Linking FREM to a full fixed-effects model
A relevant feature of FREM is that such a model can be transformed into any combination of fixed-effects covariate model, that is, a NLMEM with covariates handled as fixed effects.This implies that any FREM model with K covariates can be expressed as any combination of the 2 K − 1 possible fixed-effects covariate models without re-estimating the parameters.
In the simplest case with linear models, this is done by recalculating the parameters ( par ) and the corresponding random effect b par conditional on the covariates cov according to: where the full inter-individual covariance matrix D is defined as: Note that only the rightmost part in Equation ( 5) is dependent on the actual value of the covariates through cov, which enables efficient calculations of the conditional parameters for different covariate values.Another remark is that the rows and columns of covariates that should be excluded in the fixed-effects model are removed from D cov and D par,cov but their correlations to other included covariates will still influence the  par|cov and b par|cov .

Implicit handling of missing data
In FREM, missing values in covariates are handled as missing values in dependent variables that is, they are simply excluded in the analysis.However, the analysis is still implicitly informed by the missing covariate information through the correlation to other covariates and dependent variables since the variance-covariance matrix of all parameters and covariates is estimated for the whole population.This can be illustrated by the following example: consider a model that includes 2 highly correlated covariates and the value of 1 of the 2 covariates is missing for 1 subject.This missing value can be estimated by using the information provided by the correlation between the non-missing covariate observation and the model parameters.In contrast to modeling complete data, in the case of missing covariate data, the distribution of the missing covariate is influenced by the parametrization of the covariate model.Missing categorical covariate data Referring to the example of the bi-modal distribution of the covariate SEX, the random effect estimated by FREM may now be of for example, bsex,i ≈ 0.4.Meaning that the estimate indicates that it is more likely that the subject is a female than a male even though the estimated imputed value is neither female (ŷ sex,i = 2) nor male (ŷ sex,i = 1) but somewhere in between (ŷ sex,i ≈ 1.7).

Missing continuous covariate data
Similar issues might arise when estimating missing continuous covariates.As an example, a missing value of a log-normally distributed covariate will be imputed based on the assumption of normality if the simple additive form of the covariate model is used (see Equation ( 4)).A way to avoid this is to either log-transform the data of the covariate or change the covariate model to a log-normal model.However, both of these modifications have an impact on how the covariate-parameter model is interpreted, as described in the next section.

Interpretation of fixed-effects covariate-parameter relationships
The basic NLMEM equation (Equation 1) can be further expanded by defining the functional form of the relation between the inter-individual random effects (b par,i ) and the structural parameters  par where g(.) is a (potentially) nonlinear function.Examples of distributions that corresponds to specific functions of g(.) are Since the relationship with the conditional normal distribution mean is linear with respect to the random effect (see Equation ( 5)), the parameter-covariate relationship will be a function of the parameter function g(.) and the covariate model parametrization f cov (.) (or covariate data transformation).This way, covariate effects replace the random effects of the parameter model.For example, a normal covariate model (f cov =  cov + b cov,i +  cov ) with an exponential parameter model (g par =  par e b par,i ) yields an additive parameter-covariate relationship on the log scale: where the covariate contribution is calculated according to Equation (5).
Similarly, a parameter on the normal scale, will yield a parameter-covariate relationship on the normal scale: Most of the parameter-covariate relationships that are available in the full fixed-effects modeling approach can be implemented (by covariate and/or parameter transformations) in FREM but some are only possible to approximate.For example, the log-normal proportional effect is not trivial to derive with the FREM approach, that is,

Simulating with a FREM model
Data can be directly simulated using a FREM, together with simulations of the covariates, since the covariate distributions are explicitly encoded in FREM.However, when simulating covariates, there is no guarantee that the observed covariate distributions are kept.For example, to ensure that the covariate SEX is a dichotomous variable, one can either round the samples into a dichotomous distribution or use transformations of normal samples into binomial samples via the covariate parametrization function f cov (.).A characteristic of FREM when simulating data and covariates is that covariates are simulated using their estimated covariances, both with the response (via the parameter variances) and between covariates.
An alternative way to simulate data is to transform the model into a FFEM.After the transformation, data can be easily simulated using the observed covariates and the FFEM.However, in order to simulate using the observed covariates, a FFEM for each missing covariate pattern is needed, since the covariate coefficients and IIV (see Equation ( 5)) depend on which covariates are missing.Hence, to correctly simulate using K covariates, as many as 2 K − 1 FFEMs might be needed if all combinations of missing covariates are present.On the other hand, if several subjects have the same missing covariate pattern, these values can all be simulated using the same FFEM.This process might decrease the computational burden when simulating large studies.
The latter technique is useful, for example, when creating visual predictive check (VPCs). 17,18Then one can often use the realized design and observed covariates; hence, it is important to acknowledge each subject missing covariate pattern to be able to simulate with the correct IIV for each subject.

Relate covariate effects to explained variability
A useful way to evaluate a model's performance is to quantify how much the model is able to explain the variability in observed features of the data.For the dataset used here, an important feature are the observed, individual, minimum HAZ values (ie, the nadir) across the whole age range.The maximum explainable variability is obtained by comparing the variance of the model-predicted nadir values and the variance of the observed nadir values.Similarly, the ability of each separate covariate to explain variability can be obtained by using a model with only a single covariate included.Since FREM can be reformulated into any combination of fixed-effects, including single covariate models (as described in Section 2.2), one can quantify the ability of each covariate to explain variability.This can be calculated either for the observed variability (eg, for the nadir) or for the maximum explainable variability.

DATA: COHORTS STUDY
The COHORTS collaboration includes birth cohorts from five low and middle-income countries, namely Brazil, Guatemala, India, the Philippines, and South Africa. 19The main characteristics of the cohorts are hereby described by study site: • Brazil: the cohort includes individuals born in the city of Pelotas in 1982.• South Africa: the cohort includes individuals born in public hospital in Johannesburg-Soweto in 1990. 24 the COHORTS study, children were followed for up to 19 years and, among other measurements, longitudinal height for age Z-score (HAZ) observations were collected.The sampling design differed between the sites in number of observations, duration of observations, and observations' times (Figure 1 and Table 1).The information collected from the different sites was harmonized and pooled to increase the power for a joint analysis of the data.Data sharing regulations do not apply to this article as no new data were created or analyzed in this study.
The proportion of missing covariate data (Figure 2) shows that a few covariates are completely missing in some of the cohorts, while other covariates have variable rates of missing data in different cohorts.Study site and SEX were the only covariates that did not have any missing values across the whole COHORTS study.F I G U R E 2 Proportion (%) of complete covariate data in the COHORTS study (all data), shown by covariate and stratified by study site.

Modeling of the COHORTS HAZ data
A NLMEM was developed (see Section 2) to describe the longitudinal HAZ data.The model building steps are summarized below: • The data were randomly split into an analysis dataset including 80% of the subjects (dataset used to build the model), and an hold-out dataset with the remaining 20% of the subjects (dataset used for possible model refinements and predictions).
• Starting form a previously developed model, 25 a structural model was established to describe the trends of the HAZ data.
• Covariate effects where included in the model using the FREM approach, described in Section 2.
• Predictions were performed into the hold-out data.
The covariates that were included in the model are: study site-SITEID; sex (males/females)-SEXN; maternal years of education-MEDUCYRS; maternal age at birth-MAGE; number of previous births-PARITY; social economic status-SESN; the number of people in house, crowding index-CRWDINDX; use of preventive health care-HCUTILIZ; birth weight -BIRTHWT; all children to adult ratio-CHLDADLT; sanitation-SANITATN; Water availability at water source-H2OAVAIL; maternal empowerment quartile-MATEMPQT; maternal marriage status-MMARITN; maternal height-MHTCM; health care access-HCACCESS; maternal age at first birth-MAGE1ST; paternal years of education-FEDUCYRS; gestational age at birth-GAGEBIRTH; child dependency ratio-CHLDDR; normalized total income in household-INCTOT; duration of breast feeding-DURBRST; maternal number of male children-NMCHILD; maternal number of female children-NFCHILD; number of children in the home-NCHLD; number of male adults in the home-NADULTM; number of female adults in the home-NADULTF;-paternal age at child birth-FAGE; mother smoked during pregnancy-SMOKED; child length at birth-BIRTHLEN; number of rooms in house-NROOMS; fathers height-FHTCM.The predictive performance of the models was evaluated by VPCs.To this aim, data were simulated 100 times using the same study design and covariate data from the subjects in the analysis dataset.For the VPCs, simulated data were plotted with the observed data of either the analysis dataset or the hold-out dataset (Figure 3).
A NLMEM estimation tool, NONMEM 7.3 26 with the importance sampling method (IMPMAP) was used for all models evaluated.

SIMULATION SETUP
Simulations were used to assess the performance of the FREM methodology with various rates of missing (MCAR and MAR) covariate data.The FREM approach was compared to the FFEM CC or FFEM ME (only with MCAR) in terms of bias and precision of estimations with the final aim to evaluate their ability to handle missing covariate information.
As already mentioned, missing covariate data are handled differently by the methods we compared: they are treated as missing data in FREM, they are imputed as mean values in FFEM ME, or subjects are disregarded in the FFEM CC.Simulations were performed using only the covariates from the Indian site of the COHORT dataset, since it had rich HAZ observations and large proportion of complete covariate data.Three baseline covariates were considered: birth weight (BW), birth length (BL), and SEX.Only subjects that had complete data for all the three mentioned covariates (n = 6626) were used.
First, the parameter coefficients for the Indian cohort were re-estimated, using maximum likelihood, as a FREM model, using only the three selected baseline covariates as predictors.Then, these estimates were recalculated (Equation 5) to a FFEM model assuming that all covariates were known.This FFEM was then used to simulate the new HAZ data (with either 1000 or 100 subjects sampled with replacement, including the three covariates and sampling design, from the Indian site) with various rates of missing (MCAR) covariate information at random, that is, with a uniform deletion of the covariates according to the missing level.In addition, in the MAR scenario, 30% missing (MCAR) data (in BW and SEX) were considered and the proportions of the 30% missing (MAR) data in BL were dependent on: (1) the observed BL where 90%, 80%, 70%, or 60% of the 30% missing data were missing in subjects with an observed BL above the median.Consequently 10%, 20%, 30%, or 40% of the missing 30% data in BL was missing in subjects with an observed BL below or equal to the median.(2) The likelihood of having missing BL data was dependent on the simulated HAZ, where 80% of the 30% subjects with missing BL had a simulated HAZ below −2 at any observation time while the rest (20%) never had a simulated HAZ below −2.In the MAR scenarios, FREM were compared to FFEM CC.
In brief, the steps of the simulation study are summarized below: 1. Sample subjects (n = 1000, RICH; or n = 100, SPARSE) from the Indian COHORTS data with no missing covariate information.2. Simulate data for the sampled subjects, using the simulation FFEM.A bi-exponential model with drifting baseline, parametrized with six parameters: BASE, BASESL, BP, HLKOFF, HLKON, and PLMAX, and with sex (male or female) (SEX), BL and BW as covariates on all parameters.Simulations describe the HAZ in children (0-15 years) in India using the observed design, for each child with various rates of missing (MCAR) covariate data (0%, 10%, 30%, 70%, and 90%) and 30% missing in the MAR scenarios.3. From the simulated datasets, re-estimate, using maximum likelihood estimation, parameter coefficients using FFEM CC (both MCAR and MAR scenarios) or FFEM ME (only MCAR scenarios) with the various percentages of missing covariate data (0%, 10%, 30%, 70%, and 90%).4. Re-estimate, using maximum likelihood estimation, parameter coefficients using FREM with the various rates of missing covariate data (0%, 10%, 30%, 70%, and 90%) (MCAR scenarios) and 30% in the MAR scenarios.5. Transform the estimated FREM to a FFEM for comparison.
The procedure was repeated 1000 times and 2 metrics estimates were plotted in box plots for comparison of the methods: bias, defined as the difference between estimated and simulated values of the parameter, as well as precision, defined as the interquartile range of the box plots.In the box plots, the horizontal line in a box indicates the median value; the box edges represent the 25th and 75th percentiles; and the whiskers extend to the furthest data point that is no greater than 1.5 times the inter-quartile range.Data beyond the end of the whiskers are plotted as points.

The COHORTS model
The complete COHORTS dataset included a total of 133 471 observations for 20 421 children, with an average of 6.5 observations per child (Table 1).Model development was based on a partial analysis dataset (see Section 3.1), which included 106 568 observations for 16 291 children (Table 1 and Figure 1) as well as all covariates listed in Section 3.1.
Starting from a previously developed model, 25 the structural base model, without covariates, included six parameters: BASE, BASESL, BP, HLKOFF, HLKON, and PLMAX.The base model described changes in HAZ score over time and was a bi-exponential function with a time-varying baseline score (Equation 12).The baseline score changed linearly up to an estimated age (BP i ), after which it was constant.The slope of the linear change was associated with an additive random effect, meaning that the individual change in baseline could be either positive or negative.The model was defined by the following equations: where age i is the age of the ith child at the time of the observation (in years since birth).Child-specific parameters B i (baseline intercept), PL MAX i (the scale for the nadir), and B SL i (the slope for the time-varying baseline) were assumed to be normally distributed, while BP i (the age at which the baseline becomes time-invariant), HL k on,i (the onset half-life for nadir), and HL k off,i (the half-life for the recovery from nadir) were assumed to be log-normally distributed (see Equation ( 9)).
The structural parameter estimates of the final FREM model are presented in Table 2.The final model described the data well both in the analysis data set as well as in the hold-out data set (Figure 3).
All covariates together could explain 22% of the observed variability in HAZ at 2 years of age.When considered in a univariate manner, the covariate that explained most variability was the study site, followed by sanitation and height of the mother (Figure 4).In this univariate setting, the covariate with the strongest effect that is, study site, could explain up to 60% of the explainable variability (Figure 4).

Simulation results-missing data
For this part of the analysis, the normalized covariates' coefficients (Figure 5), the typical population parameters (Figure 6), and the variance estimates (Figure 7) for the RICH dataset design in the MCAR scenarios were considered.These results were stratified by percentage of missing covariate data and were compared among the three methods used for covariate modeling.In addition, in the MAR scenarios, typical population parameters are illustrated in Figure 8 and normalized covariates' coefficients and variance estimates are presented in the supporting information.Accordingly, bias is indicated by how close the values are to the zero, and precision is indicated by the standard deviation of the values in the box plots.Similar illustrations of the models' performance with the SPARSE dataset in the MCAR scenarios are presented in the supporting information.Our results show that, among the three methods considered with MCAR, FREM provides the most stable model with the best performance at each rate of missing (MCAR) covariate data considered.This is shown by the ability of FREM not only to model datasets that have high percentages of missing (MCAR) covariate data but also to provide an unbiased prediction of covariate coefficients regardless of the rate of missing (MCAR) covariate data (up to 90%).By contrast, the FFEM CC analysis, which initially shared with FREM a modest increase in bias (Figure 5), collapsed at 70% of missing (MCAR) covariate data due to the few subjects with complete data (see Supporting information Figures S4, S5 and  S6).On the other hand, when considering FFEM ME, the increase in bias at higher degrees of missing (MCAR) covariate information was clearly more pronounced for the covariate coefficients estimated with FFEM ME than with FREM (Figure 5).For example, in the simulations with 30% of missing (MCAR) data, the bias for the median normalized coefficient of BW on the parameter BASE was −0.0059 for FFEM CC, 1.5 for FFEM ME, and 0.0036 for FREM.The precision of the covariate coefficient estimates decreased for all methods with increasing rates of missing covariate information.This is expected since the missing (MCAR) information affects directly the sample size of the observed covariates.However, also in this case, a difference between methods could be observed: the loss in precision with increasing rates of missing (MCAR) data was lower when using FREM compared to the other two FFEM methods.For example, in the simulations with 30% of missing (MCAR) data, the precision of the estimated coefficient of SEX on the parameter HL k off was 0.27 for FFEM CC, 0.19 for FFEM ME, and 0.05 for FREM.
As expected, the estimations of the structural parameters were almost independent of the rate of missing (MCAR) covariate data in the dataset (Figure 6).For some of the parameters, a slight trend towards increasing bias with increasing   rates of missing (MCAR) data could be observed both in FFEM and FREM.By contrast, the precision was unaffected by the percentage of missing (MCAR) covariate information, with the exception of the FFEM CC analysis, which showed modeling issues already at 30% missing (MCAR) data.The superiority of the FREM approach in handling missing (MCAR) covariate data is even more evident when considering the estimated IIV variance-covariance matrix (D par ) (Figure 7).With FFEM, the bias of most of these estimates increased with increasing rates of missing (MCAR) covariate information.Conversely, with FREM, this bias was much lower and mainly visible in the datasets missing (MCAR) 90% of the covariate data (Figure 7).In addition, with FREM, the precision of the variance-covariance matrix was not affected, while with both FFEM methods the precision decreased with increasing rates of missing (MCAR) covariate information.
Similar trends as those described above, but slightly more pronounced, were also seen when SPARSE data were used (Supporting information, Figures S1, S2 and S3).F I G U R E 8 Estimated population typical parameters with the RICH design, colored by the different percentages of missing (missing at random (MAR)) covariate proportions in birth length (BL) under the two different MAR scenarios; MAR conditioned on BL and MAR conditioned on the response height for age Z-score (HAZ).Stratified by modeling technique: FFEM CC and FREM.In all scenarios, 30% missing data for each covariate is assumed.See Figure 5 for explanation of abbreviations.
In the scenarios with MAR data, the FFEM CC approach shows bias in the structural parameters due to the proportions of missing (MAR) birth length covariates, going from minor bias in the MCAR approach (50/50 MAR scenario) towards increasing bias with increasing proportions of missing high vs low birth lengths (Figure 8).This is not the same with FREM which seems to be quite consistent with a general low bias regardless of the proportions of missing (MAR) birth length.A similar trend but less pronounced is seen for the typical parameters estimated where the proportions of missing (MAR) birth length are related to the predicted HAZ (Figure 8).The MAR proportions impact of the normalized covariates' coefficients (Supporting information, Figures S7 and S8) and the estimated IIV variance-covariance matrix (D par ) (Supporting information, Figures S9 and S10) are mainly that FREM gives estimates with higher precision than FFEM CC in general.When it comes to bias the performance seems very similar between the methods and the estimates are in general unbiased.

DISCUSSION
In this manuscript, we first describe the development of a model that uses a full covariate approach, the FREM, to identify relevant covariates in the COHORTS study, a global dataset with large data variability.The model describes well the changes in HAZ over time in children from five different low and middle-income countries.Next, starting from the same data, we explored the ability of FREM to handle variable rates of missing (MCAR) covariate data, and we compared this feature to the more traditional FFEM CC and FFEM ME methods.Our analyses shows that FREM had the best results in terms of both bias and precision of the estimates of covariates' coefficients and variance.Also, the bias and precision of the estimates of the typical population parameters was similar between the three methods.In general, compared to the two FFEM methods considered, the FREM approach showed the best estimation properties leading to nearly unbiased estimates and minor decreases in precision even when 90% of the covariate data was missing (MCAR), and in both the RICH and SPARSE datasets.We also looked at two different MAR scenarios where we compare the results of FREM and FFEM CC where the proportions of missing birth length is either conditioned on the observed birth length or the predicted response (HAZ).In both scenarios FREM had the best results for both bias and precision for the typical parameters while bias was similar for the estimated covariate coefficients and variance.However precision was improved with FREM in all estimated parameters.The validity of this model as well as its implication and future use are further discussed; the ability of handling missing (MCAR and MAR) covariate data is also explored.The FREM developed in this article to describe changes in HAZ over time in the COHORTS dataset was a bi-exponential function with a time-varying baseline score (Equation 12).The model described the data well and was also able to predict into a hold-out dataset with high predictive performance.This allows for it to be used in similar studies: it can be adopted to design future studies using either optimal design techniques or simulations, assuming that the population in the new study is somewhat similar to the COHORTS population.Our analysis identified the most important covariates in the COHORTS study (Figure 4): this information can also be used to design new studies and give further insight into factors important to growth development in low and middle-income countries.
Variability in the COHORTS data was large and, despite our findings, all covariates together could explain only about 22% of the observed variability in HAZ at 2 years of age.This indicates the need to further investigate additional covariates that can help identifying the most important factors influencing children's physical growth.The low rate of explainable variability also opens up for performing studies where designs are optimized for each individual, similar to what has been proposed in precision medicine.FREM models, and in general NLMEMs, are highly suitable as a prior information when optimizing individual designs and to identify potential interventions.Moreover, the FREM approach can be transformed to investigate designs using any combination of covariates that are included in the model and is therefore suitable for new designs where only part of the FREM-included covariates are available for inclusion.An example of how the model could be adapted to the specific analysis is given when considering the study site.Study site was the covariate that explained most variability in the FREM, however it is also a potentially changing covariate since for example, Brazil in the 80's might be very different from Brazil today.Similarly, different regions in Brazil might have very different outcomes with regard to physical growth.Hence, using the developed model to extrapolate to other study sites might need the adjustment on how a site (or country) has developed since the COHORTS study had started.
FREM appears to be a suitable approach for handling missing (MCAR and MAR) covariate data.The main reason for this is that in FREM covariates are included as responses and linked to the model parameters through their random effects.In the maximum likelihood estimation, the correlation between the covariates and the other model parameters are explicitly acknowledged via the variance-covariate matrix of all the random effects (covariate and non-covariate).Therefore, FREM compensates for the missing information via the correlation between the full block of IIV with both the non-missing covariates and the parameters.Also, compared to other methods, such as MI, 27 or multiple imputation using mixed-effects linear models, 28 the handling of missing data is implicit in FREM, hence no additional computational burden is added when handling missing data.Besides, in the examples showed in this article, despite that the distribution of the missing data was misspecified (as for example for the covariate SEX, which is not normally distributed, see Section 2.1), FREM gives reasonably precise and unbiased estimates even with 90% of missing (MCAR) covariate information.This is clearly not the case for FFEM CC, which collapsed at more than 70% of missing (MCAR) covariate data, or for FFEM ME, which produced increasingly biased and imprecise results at increasing rates of missing (MCAR) data.
In the examples used in the article, the bias was low for all covariates estimated with FREM, which indicates that the distributional assumptions are sufficiently similar to the underlying distribution to prevent the creation of substantial bias in the fixed-effects, regardless of the rate of missing covariate data.However, other examples will likely show larger consequences of the assumed distributions of the covariate model.
In the current analysis, the ability of the FREM approach to handle missing (MCAR and MAR) data was compared to FFEM CC and FFEM ME (only MCAR scenarios).Alternative FFEM-based methods for NLMEM analyses, such as MI and likelihood-based methods, have been used in previous studies. 29,30However, these methods are not particularly common, due to their complex implementation in NLMEM, and were therefore not included in this comparison.Nevertheless, given the current knowledge, one can expect them to perform better than FFEM ME and probably similarly to the FREM approach.However, a drawback of the MI method is that it typically demands preparation of the dataset prior to the analysis, followed by repeated re-estimations of the models using the datasets.This could be very time-consuming since NLMEMs are typically slow to estimate due to their nonlinear nature.In addition, since MI is not a deterministic approach, the results would vary slightly between different analyses of the same dataset.
It is worth pointing out that the present manuscript deals with missing predictor data and not missing response data.This is one of the reasons why the model for the covariates (Equation 4) does not have to be more complicated than the 10970258, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/sim.9979by Uppsala University Library, Wiley Online Library on [22/05/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License one shown.This is in contrast with more complicated situations, for example that of missing response data in a multi-level model, which may require more advanced models for the relevant covariates. 28o conclude, for datasets with missing (MCAR and MAR) covariate data, FREM is an appropriate method for modeling covariates.FREM has been here shown to improve the estimation properties compared to both FFEM CC and FFEM ME.In both the RICH and SPARSE datasets, FREM was able to estimate unbiased and sufficiently precise coefficients, even when 90% of the covariate data were missing (MCAR).In addition, FREM gives precise and unbiased parameters with the investigated MAR scenarios.This suggests that FREM is a valid tool for handling dataset with missing covariate data such as those in global health studies.

F
I G U R E 1 Distribution of HAZ observations by subject age.Observations are stratified by study site and are shown separately for the analysis data set (used to develop the model) and the hold-out dataset (used for model refinement and predictions).HAZ, height for age Z-score.

10970258, 2024, 5 ,F I G U R E 3
Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/sim.9979by Uppsala University Library, Wiley Online Library on [22/05/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Visual predictive checks of HAZ vs age stratified by study site.The red lines are the 2.5th, 50th, and 97.5th percentiles of the observed data; the gray shaded areas represent the simulations and are the 95% confidence intervals around the observed median and the prediction interval.The two upper panels are based on the analysis dataset, used to build the base model and the final FREM, the lower panel shows predictions into the hold-out dataset.FREM, full random effects model; HAZ, height for age Z-score.

4
Fraction of the maximum explainable variability that is explained by each covariate considered in a univariate manner.Results are stratified for the nadir (ie, the observed individual minimum HAZ values), HAZ at 2 years of age, and HAZ at 11 years of age.HAZ, height for age Z-score.

F I G U R E 6 F I G U R E 7
Estimated population typical parameters with the RICH design, colored by the different percentages of missing (MCAR) covariate data.Stratified by modeling technique: FFEM CC, FFEM ME, and FREM.See Figure5for explanation of abbreviations.Note that some boxes, whiskers, and outliers are cropped since the figure is zoomed in to focus on the smaller deviations around zero.Estimated population variance parameters with the RICH design, colored by the different percentages of missing (MCAR) covariate data.Stratified by modeling technique: FFEM CC, FFEM ME, and FREM.See Figure5for explanation of abbreviations.Note that some boxes, whiskers, and outliers are cropped since the figure is zoomed in to focus on the smaller deviations around zero. 20

•
Guatemala: the cohort includes individuals who participated as children in the INCAP Nutrition Supplementation Trial conducted in four villages in El Progreso district between 1969 and 1977.
21• India: the cohort includes individuals born in a defined region of New Delhi between 1969 and 1972. 22 23 COHORTS dataset, stratified by site and analysis/hold-out dataset.
TA B L E 1Abbreviation: HAZ, height for age Z-score.a Number of HAZ observations.b Number of children.c Minimal and maximal children ages (in years) for the HAZ observations.10970258, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/sim.9979by Uppsala University Library, Wiley Online Library on [22/05/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 10970258, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/sim.9979by Uppsala University Library, Wiley Online Library on [22/05/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 10970258, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/sim.9979by Uppsala University Library, Wiley Online Library on [22/05/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Parameter estimates of the structural part of the final FREM model.PL MAX i , scale for the nadir; BP i , age at which the baseline becomes time-invariant; B SL i , slope for the time-varying baseline; HL k off,i , half-life for the recovery from nadir; HL k on,i , onset half-life for nadir; D par covariance matrix of the individual random effects (b par,i ).
TA B L E 2Note: B i , baseline intercept; a Inter-individual variability (on variance scale).bD par = 10970258, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/sim.9979by Uppsala University Library, Wiley Online Library on [22/05/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Estimated normalized covariate coefficients with the RICH design, colored by the different percentages of missing (MCAR) covariate data.Stratified by modeling technique: FFEM CC, FFEM ME, and FREM.MCAR, Missing Completely At Random; FFEM CC, FFEM with complete case analysis; FFEM ME, FFEM with mean imputation; FREM, full random effects model.Note that some boxes, whiskers, and outliers are cropped since the figure is zoomed in to focus on the smaller deviations around zero.