nlmeVPC: Visual Model Diagnosis for the Nonlinear Mixed Effect Model

A nonlinear mixed effects model is useful when the data are repeatedly measured within the same unit or correlated between units. Such models are widely used in medicine, disease mechanics, pharmacology, ecology, social science, psychology, etc. After fitting the nonlinear mixed effect model, model diagnostics are essential for verifying that the results are reliable. The visual predictive check (VPC) has recently been highlighted as a visual diagnostic tool for pharmacometric models. This method can also be applied to general nonlinear mixed effects models. However, functions for VPCs in existing R packages are specialized for pharmacometric model diagnosis, and are not suitable for general nonlinear mixed effect models. In this paper, we propose nlmeVPC , an R package for the visual diagnosis of various nonlinear mixed effect models. The nlmeVPC package allows for more diverse model diagnostics, including visual diagnostic tools that extend the concept of VPCs along with the capabilities of existing R packages.


Introduction
After fitting a model, diagnosing the fitted model is essential for verifying that the results are reliable (Nguyen et al. 2017).For linear models, the residuals are usually used to determine the goodness of fit of the fitted model.However, due to random effects and nonlinearity, residuals are less useful for diagnosing fit in nonlinear mixed effects models.Therefore, various diagnostic tools for this type of model have been developed.The nonlinear mixed effect model is useful when the data are repeatedly measured within the same unit, or the relationship between the dependent and independent variables is nonlinear.It is widely used in various fields, including medicine, disease mechanics, pharmacology, ecology, pharmacometrics, social science, and psychology (J.C. Pinheiro and Bates 2000;Davidian 2017).Recently, among the various diagnostic tools applicable to nonlinear mixed models, simulationbased diagnostic methods have been developed in the field of pharmacology (M.Karlsson and Savic 2007).The visual predictive check (VPC; M. O. Karlsson and Holford (2008)) is a critical diagnostic tool that visually tests for the adequacy of the fitted model.It allows for the diagnosis of fixed and random effects in mixed models (M.O. Karlsson and Holford 2008) in the original data space.Currently, it is a widely used method for diagnosing population pharmacometric models: Heus et al. (2022) used the VPC method to evaluate their population pharmacokinetic models of vancomycin, Mudde et al. (2022) checked their final population PK model for each regimen of antituberculosis drug using the VPC method, and Otto et al. (2021) compares the predictive performance of parent-metabolite models of (S)-ketamine with the VPC method.This VPC method can be used for general nonlinear mixed effect models, including hierarchical models.
The psn (Lindbom, Ribbing, and Jonsson 2004) and xpose4 (R. J. Keizer, Karlsson, and Hooker 2013) packages provide various diagnostic methods for pharmacometric models in R (R Core Team 2023), including the VPC plot.These packages were developed only for the pharmacometricians, and the users needed to use the NONMEM software (Bauer 2011) for generating inputs of functions in these two packages.However, NONMEM is licensed software, and it is mainly designed towards the analysis of population pharmacometric models.Therefore, it is not easy for nonpharmacometrician to use NONMEM and to draw the VPC plot through psn and xpose4.Recently, the vpc (R. Keizer 2021) package and nlmixr2 (Fidler et al. 2019) package have been developed to draw VPC plots in R without results from NONMEM.However, vpc package provides the function for drawing the original VPC plot only.nlmixr2 was developed initially for fitting general dynamic models.To check the fitted model, nlmixr2 provides a function to use graphical diagnostics in xpose4 with the nlmixr2 model object.Also, nlmixr2 uses the VPC plot in the vpc package with the nlmixr2 model object.Therefore, both packages only provide a function to draw the basic VPC plot, and the other newly developed simulation-based methods, including extensions of VPC, are not provided.
We have developed a new R package, nlmeVPC, to provide a suite for various visual checking methods for the nonlinear mixed models.This package includes various state-of-the-art model diagnostic methods based on the visual comparison between the original and simulated data.Most methods compare the statistics calculated from the observed data to the statistics from the simulated data.Percentiles, for example, the 10 th , 50 th , and 90 th percentiles, are widely used to summarize relevant statistics from the observed and simulated data.We compare the similarities between the statistics from the observed data and those from the simulated data in two different spaces: the data space and the model space (Hadley Wickham, Cook, and Hofmann 2015).The original data comprise The R Journal Vol.15/1, March 2023 ISSN 2073-4859 the data space.Usually, the data space is represented by the independent and dependent variables, such as time and blood concentration, in the pharmacokinetic data.On the other hand, the model space is composed of quantities obtained from the fitted model, for example, residuals and summary values from the fitted model.From this viewpoint, we categorize the well-known visual diagnostic tools into two categories.One method compares the observed and simulated data in the original data space, and the other in the model space.In the method with the data space, we developed functions for the original VPC (VPCgraph), the additive quantile regression VPC (aqrVPC), and the bootstrap VPC (bootVPC).In addition, we proposed a new VPC method: the average shifted VPC (asVPC).In the method with the model space, the coverage plot (coverageplot) and the quantified VPC (quantVPC) are included.For a more detailed diagnosis, we developed a coverage detailed plot (coverageDetailplot).In nlmeVPC, the ggplot2 package (H.Wickham 2016) is used to create all plots.

Nonlinear mixed effect model
The general nonlinear mixed effect model is defined as follow: where y ij is the dependent variable of the j th observation of the i th individual, x ij is the independent variable, f is the nonlinear function of x ij , θ is the population parameter, η i represents the variability of the individual i, and ϵ ij represents the random error.From the data, θ , Ω, and Σ are estimated.For the simulated data, the fitted model

Validation in the data space
The VPC is the most popular model validation method in the pharmacometrics area.It was developed to diagnose population pharmacokinetic/pharmacodynamic models visually.The main idea of the VPC is to compare the distribution of observations and the distribution of predicted values, where the distribution of predicted values is obtained from simulated data drawn from the fitted model.If the fitted model explains the observed data well, these two distributions should be similar.Both distributions can be represented in the original data space that consists X axis as the independent variable and Y axis as the dependent variable.It allows us to compare the observed data with the fitted model in the original data space.In nlmeVPC, we include an original VPC plot, an additive quantile regression VPC (Jamsen et al. 2018), and a bootstrap VPC (Post et al. 2008).We also proposed a new approach to draw the VPC: the average shifted VPC.

Visual Predictive Check
The visual predictive check (VPC; M. O. Karlsson and Holford (2008)) is based on the principle that if the fitted model adequately describes the observed data, the distribution of the simulated data from the fitted model should be similar to the distribution of the observed data.There are several ways to compare the similarities between the distributions.In the VPC approach, profiles of quantiles are used.Two profiles are mainly used to compare the distributions of observations and predictions.One profile is from the upper bound of the prediction intervals, and the other is from the lower bound.These prediction intervals are calculated from the simulated data.90% prediction intervals are usually used.For small and sparse samples, 80% prediction interval is also used.The lower and upper bounds of 80% prediction interval are the 10 th and 90 th percentiles of the simulated data.Figure 1(A) shows the "scatter" type of the VPC plot.Dots indicate the observed data.Two dashed blue lines represent profiles of the 10 th and 90 th percentiles of the simulated data, and the solid blue line represents the 50 th percentile.If the fitted model represents the observed data well, most observed data should lie between profiles of 10 th and 90 th percentiles.
Figure 1(B) is the "percentile" type of the VPC plot.In this plot, profiles of percentiles from the observed data are compared to profiles of percentiles from the simulated data.Two dashed red lines represent profiles of the 10 th and 90 th percentiles of the observed data, and the solid red line represents profiles of the 50 th percentile of the observed data.If the fitted model represents the observed data well, two profiles in each percentile -one from the original data and the other from the simulated data -are similar.Light blue areas represent the 95% confidence areas of the 10 th and 90 th percentiles, and pink areas represent the 95% confidence areas of the 50 th percentile.These confidence areas were calculated from the simulated data.After calculating percentiles in each simulated data, we find 95% confidence intervals for each percentile and use this to draw the areas.In this plot, it is necessary to verify that the profiles of the original data are in confidence areas of each profile from the simulated data in each percentile.If each percentile line of the observed data is in the corresponding confidence area, this can be evidence that the fitted model represents the observed data quite well.Otherwise, the fitted model needs to be improved.The "CI" type of the VPC plot is the most widely used type in pharmacometrics.
The percentiles of the dependent variable are calculated in each bin.To estimate the percentiles accurately, enough data points need to lie within each bin.No binning means that the number of bins is equal to the number of different independent variable values.In this case, the VPC plot shows all details of the relationship between the independent and dependent variables.However, the resulting lines become too irregular to show meaningful trends or patterns.
As the number of bins decreases, the lines become smoother and more regular, however this can come at the loss of information if too much smoothing is used.Therefore, the selection of the best number of bins is crucial.The way of determining cutoffs for bin also plays an important role.Lavielle and Bleakley (2011) proposed a procedure for finding the optimal bin configuration automatically, including the optimal number of bins and the optimal bin cutoffs.VPCgraph provides the automatic binning with optK and makeCOVbin; here, optK finds the optimal number of bins, and makeCOVbin finds the optimal cutoffs of bins using Lavielle and Bleakley's method.

Additive quantile regression VPC
To overcome the difficulties of making bins as well as determining the number of bins, Jamsen et al. (2018) used additive quantile regression to calculate the quantiles of the observed and simulated data.This regression method makes it possible to estimate quantiles without discrete binning, which is especially useful when the data are insufficient, irregular, or inappropriate to configure the bins.To fit the additive quantile regression, we used the rqss function in the quantreg (Koenker 2023) package and developed the aqrVPC function to draw the VPC type plot with additive quantile regression.Figure 2 shows the additive quantile regression VPC plot.The solid and dashed lines represent the 10 th , 50 th , and 90 th additive quantile regression lines of the observed data, and the pink and light blue areas represent the confidence areas of the additive quantile regression lines of the simulated data.Lines and areas in the additive quantile regression VPC plot are much smoother than those in the original VPC plot.

Bootstrap VPC
The bootstrap VPC (Post et al. 2008) compares the distribution of the simulated data to the distribution of the bootstrap samples drawn from the observed data.This plot reflects the uncertainty of the observed data and allows for more objective comparisons with the predicted median.
Figure 3 shows the bootstrap VPC plot using bootVPC.The solid and dashed blue lines represent the 10 th , 50 th , and 90 th percentiles of the simulated data.The solid red line represents the 50 th percentile line, and the pink areas represent the 95% confidence areas of the 50 th percentile line, calculated from the bootstrap samples of the observed data.If the solid blue line and the solid red line are similar, the solid blue line is in the pink area, and the pink area is located between two dashed blue lines, then this is evidence that the fitted model fit the observed data well.

Average shifted VPC
Even though binning mitigates the problem with highly irregular data, the VPC plot still has a precision problem with sparse data.In this paper, we propose a new approach to draw the adapted VPC plot from the average shifted histograms (Scott 1985).A histogram is a widely used method for displaying the density of a single continuous variable.However, histograms can look quite different based on different choice of bin width and anchor.This requires computing an optimal bin width.To overcome the problem with the choice of bin width and to obtain smoother estimates, Scott (1985) proposed the averaged shifted histogram (ASH).The idea behind ASH is to make K different histograms with different anchors and to combine them via a weighted average.For each histogram, the starting point is shifted by d/K, where d is the width of the bin.We extend this idea to the VPC graph and propose the average shifted visual predictive check (asVPC) plot.The binning of the VPC varies from the traditional binning in a histogram.The width of the bins is determined such that each bins contains the same, fixed number of observations, and as such the width of each bin is different.To apply the ASH idea to the VPC, we divide our data into K * B bins along with the independent variable, where each bin has a different width but the same number of observations.Here, B is the number of the original bins in the VPC plot, and K is the number of the histograms averaged in the asVPC plot.
To draw the VPC plots, the following information is needed from each bin: (A) the median of the independent variable and percentiles of the dependent variable of the observed data.
(B) the median of the independent variable and percentiles of the dependent variable of the simulated data.
(C) 95% confidence interval of percentiles (usually 10 th , 50 th , and 90 th percentiles) of the dependent variable, calculated from the simulated data.
These values are also needed to produce confidence areas and percentile lines in the asVPC plot.Furthermore, additional steps are needed for asVPC.To find the median of the independent variable and percentiles of the dependent variable using the ASH algorithm, the following procedures are needed: 1. Divide the independent variable into N = K * B bins.

For
Calculate the median of the independent variable and the weighted percentiles of the dependent variable in the combined bin 3. Collect the medians of the independent variable and the weighted percentiles of the dependent variable from 2, and connect them to the lines.
We can implement (A) and (B) by applying these procedures separately to the observed and simulated data.Additionally, (C) can be implemented using these procedures for each simulated dataset.First, we find the weighted percentiles, combine the results from each simulated dataset, and then calculate the 95% confidence intervals of each percentile.Using these three quantities (A), (B), and (C), we can draw the VPC type plot with the ASH approach, producing the asVPC plot.

Determining the weights
In the asVPC plot, the observations in each bin are combined using weights.Typically, the data near the center of the integrated bin have higher weights, and the data far from the center have smaller weights.This idea is used in the ASH algorithm as well as the density estimation literature.We suggest two different ways to apply weights for the asVPC calculation.
• Bin-related weight: For each bin in the combined bin, K-1 bins in both directions have weights and the other bins have no weight.Use 1 K as weights for 2K-1 bins.This approach is the same as the original ASH approach.
• Distance-related weight: Use the reciprocal of the distance from the center of the independent variable in each combined bin so that the points near the center have higher weights and the points far from the center have lower weights.
Figure 4 shows the results from the asVPC function using bin-related weights and distance-related weights.The solid and dashed lines represent the average shifted quantile lines of the observed data, and the pink and light blue areas represent the confidence areas of the simulated data.The lines in the asVPC plot are smoother than those in the original VPC plot, and the confidence areas in the asVPC plot are thinner than those in the original VPC plot.

Validation in the model space
To validate the fitted model in the model space, we need to choose appropriate statistics to visualize and describe them in the model space.In nlmeVPC, we use the same statistics, that is quantiles, as the VPC plot in the data space, and compare them in two ways -numerically and visually.predictive check and the coverage plot are two commonly used methods in pharmacometrics.However, there is a limitation to detecting the illness of the fitted model.These methods combine the results of all ranges of the independent variable, and it is helpful to see the overall fitness.However, it is not easy to detect the detailed discrepancy between the fitted model and the observed data.To overcome this limitation, we developed a new plot, the coverage detailed plot, to compensate for the shortness of the coverage plot.

Numerical predictive check
The VPC plot visually compares the observed data to the simulated percentiles in the data space.On the other hand, the numerical predictive check (NPC; Wang and Zhang (2012); M. Karlsson and Savic (2007)) numerically compares the observed data to the simulated data.For a given level of prediction (for example, 90%), the predicted interval is calculated using the simulated data, and the number of the observed data points outside of the prediction interval is counted, both below the lower bound and above the upper bound.The expected number of points below the lower bound of the predicted interval (for example, 5% of observations) is also calculated and compared to the observed number.If these two numbers are similar, the suggested model is suitable (Maharaj et al. 2019).
NumericalCheck provides information summarizing the results of the NPC for various prediction levels.PI is the prediction level × 100 and n is the number of observations."Expected points below PI" and "Expected points above PI" are respectively the expected numbers below and above the PI; "points below PI" and "points above PI" respectively represent the numbers of points below and above the PI; "95%CIBelowFrom(%)" and "95%CIBelowTo(%)" represent the 95% confidence interval of "points below PI(%)"; and "95%CIAboveFrom(%)" and "95%CIAboveTo(%)" represent the 95% confidence interval of "points above PI(%)".If "points below PI(%)" is in the 95% confidence intervals of "points below PI(%)" and is similar to "Expected points below PI", and if "points above PI(%)" is in the 95% confidence intervals of "points above PI(%)" and is similar to "Expected points above PI", this is the evidence that the fitted model explains the data well. NumericalCheck(origdata,simdata,pred.level=c(0,0.2,0.4,0.6,0.8,0.9),N_xbin=8)$NPC

Coverage plot
The result of the NPC is a table with many values, which, while useful, can be difficult to parse visually.The coverage plot (M.O. Karlsson and Holford 2008) was developed to help visually check the fitted model with the NPC result.In each level of the predicted interval, the ratios between the expected number of data points (Exp) outside the prediction interval and the observed number of data points (Obs) outside the prediction interval are calculated.These ratios are calculated separately for the upper and lower bound of the prediction interval.For example, when the prediction level is 90, a 90% prediction interval is used, and 10% of the observations are expected to locate outside this prediction interval.To be more precise, 5% of observations are expected to be above the upper limit, and 5% of observations are expected to be below the lower limit.
The R Journal Vol.15/1, March 2023 ISSN 2073-4859 The coverage plot with the NPC result can diagnose a model using multiple prediction intervals.The X-axis represents the prediction level, and the Y-axis represents Obs/Exp.The closer Obs/Exp is to 1, the more appropriate the model is.Furthermore, the confidence intervals of Obs/Exp values are obtained using simulated data and then expressed together in the plot for a more objective comparison.This plot can provide more information than the VPC plot, which interprets only a couple of quantiles -usually the 10%, 50%, and 90% percentiles.
The xpose4 package (R.J. Keizer, Karlsson, and Hooker 2013) provides a coverage plot.However, to draw the coverage plot using xpose4, PsN (Lindbom, Ribbing, and Jonsson 2004) software is needed to calculate the NPC result.Therefore, we developed coverageplot to draw the coverage plot using the results from NumericalCheck.
Figure 5(A) shows the coverage plot.The X-axis shows the level of the prediction interval.The Y-axis show the ratio between the observed number of data and the expected number of data of the lower and upper parts in each level of the prediction interval.The white line is the reference line, and the gray area represents the confidence areas of the ratios in each prediction level.If the solid lines are near the white line or in the gray area, we can conclude that the suggested model is suitable.

Coverage detailed plot
Unlike the VPC plot, which represents the data space, the information in the observed data space does not come together in the coverage plot, which makes it difficult to determine whether the model is overestimated, underestimated, or adequate in the specific region of the data space.To overcome this limitation of the coverage plot, we propose a new method called the coverage detailed plot.The percentages of observations above the prediction interval are calculated in each bin of the independent variable.Additionally, the percentages of observations below the prediction interval are calculated.The white dots in the plot represent the expected percentages.If the percentages of upper and lower observations are near the white dots, we can conclude that the suggested model is suitable for the specific prediction interval.
Figure 5(B) is the result of coverageDetailplot when the prediction level is 80%.The white dots represent the expected percentages of the lower and upper the prediction intervals, 10%, and 90%, respectively.The upper and lower percentages of observation in each time bin are shown in darker gray.The left bin(before 0.045 hours) shows all light gray in the coverage detailed plot, and it is quite different patterns from the expected one.However, it is mainly due to the characteristics of this example data.All observations in this bin are 0. It makes the lower and upper bound of the prediction interval all 0, and the lower and upper percentages become 0. Post et al. (2008) proposed the quantified VPC (QVPC).It expanded the existing VPC, including information about missing data.The QVPC plot visually represents actual and unavailable observations around the predicted medians through the form of percent, regardless of the observed data's density or shape.Here, "unavailable observations" refer to all kinds of missing values and unobserved data that occur for various reasons, including deviations and unexpected omissions.

Quantified VPC
If the model is appropriate, observations at each time bin are allocated randomly around the predicted median of the model.In the QVPC plot, white dots represent the model's predicted median in each bin.If the borderlines above and below are close to the white dots, we can conclude that the fitted model describes the observed data well.
Figure 6 shows the result of quantVPC.The darker gray areas represent the percentages below the median.The lighter gray areas represent the percentage above.The brightest gray areas represent the percent unavailable in each time bin.The white dots represent the ideal location where the above and the below percentages meet.In this example, there is no missing value.

The nlmeVPC package: structure and functionality
Table ?? shows the list of functions and Table ?? shows the list of arguments used in the functions of nlmeVPC.The following codes are for Figure 1 to Figure 6.The coverage plot and the coverage detailed plot for the 80% prediction interval.In the coverage plot, the X-axis is the level of the prediction interval.The Y-axis is the ratio between the number of observed data and the number of expected data of the lower and upper parts in each level of the prediction interval.The white line is the reference line, and the gray area represents the confidence area of the ratios.If the solid lines are near the white line, we can conclude that the suggested model is suitable.In the coverage detailed plot, the white dots represent the expected percentages of lower and upper prediction intervals of, 10%, and 90%, respectively.The upper and lower percentages of observation in each time bin are darker gray.
Figure 8 shows the results of aqrVPC, and Figure 9 shows the results of asVPC for Model 1 and Model 2. The additive quantile regression VPCs and the average shifted VPCs show similar patterns to the original VPCs.Model 1 shows all quantile lines in the confidence area, and the 10 th and 90 th percentiles in Model 2 are mostly outside the confidence area.
Figure 10 shows the results of bootVPC for Model 1 and Model 2. The solid and dashed blue lines show the 10 th , 50 th , and 90 th percentiles of the simulated data.The solid red line represents the 50 th percentile line of the observed data, and the pink areas represent the 95% confidence areas of the 50 th percentile line, calculated from the bootstrap samples of the observed data.The solid blue line is in the pink area, and the two solid red and blue lines are almost identical in both models.The dashed blue lines in Model 1 cover most of the observed data.However, the dashed blue lines in Model 2 do not Figure 11 shows the coverageplot results for Model 1 and Model 2. The lines are in the gray area and close to the white line in Model 1.However, the lines in Model 2 are not in the gray area, especially when the PI value is large.Figures 12 and 13 show the results of coverageDetailplot for Model 1 and Model 2 when PIs are 50% and 80%.The upper and lower percentages in both figures are close to the white points in Model 1.On the other hand, the upper percentages of the most time bins are far from the white points in Model 2, especially the time bin (3.54,5.28]when PI = 50%.When PI = 80%, most upper and lower percentages are far from the white points. Figure 14 shows the results of quantVPC for Model 1 and Model 2. In Model 2, the right tail area (after 11 hours) looks quite different from the expected pattern.The above percentages are much larger than the below percentages.

Figure 1 (
Figure 1(C) is the "CI" type of the VPC plot.The solid red line represents the 50 th percentile of the observed data, and dashed red lines represent the 10 th and 90 th percentiles of the observed data.

TheFigure 1 :Figure 2 :Figure 3 :
Figure 1:The visual predictive check plot.The solid red line represents the 50 th percentile of the observed data, and dashed red lines represent the 10 th and 90 th percentiles of the observed data.The solid blue line represents the 50 th percentile of the simularted data, and dashed blue lines represent the 10 th and 90 th percentiles of the simulated data.Light blue and pink areas represent the 95% confidence areas of the 10 th , 50 th and 90 th percentile lines.

Figure 4 :
Figure 4: The average shifted VPC plot.Dots indicate the observed data.The solid line represents the 50th quantiles of the observed data, and dashed lines represent the 10 th and 90 th percentiles of the observed data.Light blue and pink areas represent the 95% confidence areas of the 10 th , 50 th , and 90 th percentiles.

Figure 5 :
Figure 5:The coverage plot and the coverage detailed plot for the 80% prediction interval.In the coverage plot, the X-axis is the level of the prediction interval.The Y-axis is the ratio between the number of observed data and the number of expected data of the lower and upper parts in each level of the prediction interval.The white line is the reference line, and the gray area represents the confidence area of the ratios.If the solid lines are near the white line, we can conclude that the suggested model is suitable.In the coverage detailed plot, the white dots represent the expected percentages of lower and upper prediction intervals of, 10%, and 90%, respectively.The upper and lower percentages of observation in each time bin are darker gray.

Figure 6 :
Figure 6:The quantified VPC plot.The darker gray areas represent the percentages below the median, the lighter gray areas represent the percentage above, and the brightest gray areas represent the percent unavailable in each time bin.The white dots represent the ideal location where the above and the below percentages meet.