Applying Hierarchical Generalized 1 Additive Models to Non-Compartmental 2 Analysis of Pharmacokinetic Data 3

Abstract


Introduction
Non-compartmental analysis (NCA) is a popular strategy for statistical analysis of pharmacokinetic (PK) data.Under specific assumptions, the utility of NCA-based estimation of pharmacokinetic parameters is well-established (Gabrielsson and Weiner, 2012).Though population pharmacokinetic methods, comprising parametric nonlinear multilevel models, have become standard methodology in many applications, important roles for NCA remain.NCA offers simple, direct estimation of key pharmacokinetic parameters in cases where limitations of data, knowledge, expertise, or resources constrain the use of more explicit models.Further, NCA is a primary technique in experimental pharmacokinetics, especially in bioequivalence studies, in cases where selection of a structural model is not desired.
The statistical simplicity of NCA, at least within a typical workflow, carries some relevant practical disadvantages.Numeric methods such as trapezoidal rules for the determination of AUC, the key primary parameter estimated in NCA, are reliant on dense sampling, as they have no mechanism to share information between subjects.The sharing of information between subjects (partial pooling), is a key characteristic of population pharmacokinetics (Gelman, Bois and Jiang, 1996), and more generally, multilevel modelling (Greenland, 2000).The numeric methods produce only a point estimate, and do not indicate their uncertainty in the estimated AUC; regardless of the accuracy of the method, the resulting AUC are drawn from finite data and so are not absolutely known.This uncertainty is important for reporting, to reflect the inherent limitation of results of any individual study, but also to the downstream usage of the parameter estimates, as typical in the two-stage approach, for evaluation of bioequivalence for example.Censored data are also not accommodated in the typical workflow for NCA.Recent work highlights the use of various error and imputation models for estimation of population NCA parameters in the presence of left-censoring (Barnett et al., 2021), but not individual NCA estimation.In veterinary applications for example, arbitrary decisions such as dropping censored observations are generally abundant (Woodward and Whittem, 2019).In population pharmacokinetics censored observations may be rigorously handled via an appropriate likelihood function (Beal, 2001), but this of course is not present in NCA, which may lead to poor performance (Hughes, Upton and Foster, 2017).
The primary parameter estimates in NCA are integrals; the area under the concentration-time curve (AUC), and the area under the moment concentration-time curve (AUMC), from which other parameters are secondarily derived (Veng-Pedersen, 2001).These are classically obtained using trapezoidal methods, which provide point estimation.Though accurate under suitable conditions (Chiou, 1978;Dunne and King, 1989), they determine only one 'best' value, and do not indicate uncertainty.To obtain uncertainty estimation of such integrals, encapsulation of the area estimation within a statistical model, with a suitable likelihood function, is a natural approach.A clear barrier is the specification of the form of the statistical model, given the general principle of non-compartmental modelling to not require a model to be selected by the user.A sensible candidate is a spline-based model, in which the shape of some relationship is a simple smooth function (Tian, Yu and Kim, 2020).
Spline methods can be applied directly to single-level concentration-time data, and examples including cubic splines have been investigated for NCA, with middling results (Purves, 1992).They may also be utilized as elements of a larger statistical model, especially the generalized additive model (GAM; Hastie and Tibshirani, 1986).GAM can be considered a semi-parametric extension of generalized linear models (GLM), in which the relationship between the linear predictor and one or more predictors is expressed by some form of spline (Wood, 2020).In contrast to GLM, which require the investigator to completely specify the form of the relationship between the response variable and all predictors, HGAM allows the form of the effects of one or more of the predictors to be determined from the data during model fitting, allowing for great flexibility.These methods have extensive application in relevant disciplines, such as ecology (Simpson, 2018) and epidemiology (Dominici, 2002), with scope for biomedical applications (Mundo, Tipton and Muldoon, 2022).These methods have a close relationship with other non-parametric techniques such as Gaussian process regression and similar models (Wood, 2020).Current implementations make these models quite accessible.Package 'mgcv' (Wood, 2023) implements a variety of computationally efficient HGAM methods (Wood, 2011) using the penalized likelihood, which are also supported in an explicitly Bayesian form in package 'brms' (Bürkner, 2017), an interface to the Bayesian programming environment Stan (Carpenter et al., 2017).
Various semi-parametric techniques have been widely used in pharmacometrics, generally as components of broader models.Various examples utilize splines to capture the form of covariates (Lai, Shih and Wong, 2006), or to describe time-varying parameters in cases where constant models are insufficient (Li et al., 2002).More recently, Gaussian process models were applied to evaluate potentially non-linear maturation effects based on age (Siivola, Weber and Vehtari, 2021).Semiparametric models for primary analysis of concentration-time functions in pharmacokinetics have been presented sporadically, but appear to have received limited attention.Notably, Park et al., (1997) reported semiparametric modelling using splines, as a convenient and flexible alternative to nonlinear multilevel models, including random effects for between-subject variation.Jullion et al., (2009) described application of Bayesian P-splines, including extensive mathematical details and consideration of priors, with a focus on sparse problems.Most recently, Willemsen et al., (2017) described the application of a specific variant of multilevel spline model to bioequivalence determination with NCA, and nominated several potential advantages, including handling of sparsity and implementation of censoring, which are directly relevant to the HGAM approach evaluated here.Though these evaluations demonstrate substantial utility of a semi-parametric approach to NCA, they are technically sophisticated and demand a relatively high degree of skill to implement.
Broadly, HGAM may facilitate probabilistic, rather than deterministic, estimation in NCA, allowing direct quantification of uncertainty, while aiding the flexibility of analysis.A key goal must be to retain the ease-of-use and data-driven spirit of NCA, relative to other pharmacometric methods which are highly dependent on investigator expertise and nomination of appropriate models.Use of readily available and accessible implementations is thus a key priority.The objective of this manuscript is to explore the application of the hierarchical generalized additive model (HGAM) to NCA, building upon previous applications of semiparametric models for NCA (Park et al., 1997;Jullion et al., 2009;Willemsen et al., 2017).A Bayesian workflow in the open software environment R (R Core Team, 2022) via 'brms' (Bürkner, 2017)

Example 1: Application to Descriptive Pharmacokinetics
The classical theophylline dataset, originally reported by Upton and utilized by Beal et al., (2009) in early NONMEM work, and by Pinheiro and Bates (1995), among many other methodologic applications (Comets, Lavenu and Lavielle, 2017), was obtained from R (R Core Team, 2022).The data comprise observations of the plasma concentration of theophylline from 12 subjects over 24 hours, following oral administration.For the current work these data were selected as an archetypal example of a descriptive PK study with intensive sampling.For this application, the concentrations reported as zero were removed from the dataset; later, censoring will be implemented using HGAM, but in this dataset the censoring limit is not reported.This paper emphasises the practical application and interpretation of hierarchical generalized additive models to non-compartmental analyses, and does not discuss mathematical details.Many resources are available that describe mathematical and statistical principles of HGAM and applications, at various levels of complexity (Pedersen et al., 2019;Van Rij et al., 2019;Wood, 2020).A hierarchical generalized additive model for the theophylline concentration at any point in time after administration was specified as:   (  )~( 0 , ) + (  ( + 1)) +   (  ( + 1)) + (0, ) This model was expressed in R and 'brms' syntax as: brm(Conc_log ~ 1 + s(Time_log1, k = 6) + s(Time_log1, Subject, bs = 'fs', m = 2, k = 6)… In which Conc_log is the log-transformed   , Time_log1 is the pseudo-log-transformed , and Subject is a categorical variable indicating the subject.The model specification and notation here follows a highly recommended introduction to hierarchical GAM in R (Pedersen et al., 2019), based on the model syntax from package 'mgcv' (Wood, 2023).This specification denotes a population ('global') smooth effect of  (), between-subject variation in the effect of  using the factor-smooth interaction term   (), between-subject variation (location shift)  in the intercept  0 , and normally-distributed constant conditional variation , all to predict the concentration   on   scale.The population smooth function in this case, where the basis is not specified, is a thin-plate regression spline (Wood, 2003).Note that the pseudo-logarithmic transformation of  allows  = 0 to be included (Jullion et al., 2009), but is somewhat non-linear especially for small values of , which requires some caution to select an appropriate unit.This corresponds to a group-level model with shared wiggliness, as expressed in Pedersen et al., (2019).Analogous to the more general multilevel models, in which group-level parameters are drawn towards a common average, in this model the group-level functions are drawn towards a common average function, but with group-level variation in shape.
Parameters for this model include the smooth function variability for  and a coefficient for the linear trend of , the intercept  0 and its between-subject variation , and the conditional standard deviation σ.An important consideration is the specification of prior information for these parameters.Some parameters are readily interpretable.A prior distribution for the intercept represents the average   for average  (intercept priors are centered in 'brms'), and  is between-subject variation in this average, representing location-shift of the subjects.For this model these were (0,5) for  0 , intended to be weakly informative;  is implemented within the smooth term.The prior distribution for σ was (0,1), which denotes that the for a sensible population pharmacokinetic model the scale of residual error, on logarithmic scale, should be reasonably small (less than 1).Priors for the smooth components are more difficult because they do not map to very interpretable statements about any real quantity, and each were retained using the package defaults; the linear trend of time was flat (−∞, ∞), and the smooth function variance was (0,2.5,3).
Estimation of the parameters was conducted via MCMC, by the No-U-Turn Sampler as implemented in Stan (Carpenter et al., 2017).Convergence and sufficiency of sampling were assessed using the ^ statistic (Vehtari et al., 2021) and effective sample sizes.Primary parameter estimates are reported in Table 1.The final model was strongly descriptive of the total variation in concentration (Bayes R 2 : 0.946, 90%CrI: 0.933, 0.957; Gelman et al., 2019).The fitted smooth functions, for each subject, are expressed visually (via ggplot2; Wickham, 2016) in Figure 1; overall the apparent fit of the model to the individual concentration data was adequate.There was no clear evidence of any major lack-of-fit for any subject.Residual analysis demonstrated that the constant-variation model was compromised by an apparent effect of time on residual variability (Figure 2), with residual variability declining with time.This pattern suggests relative difficulty in capturing the absorption phase.
A revised model was implemented: in R and 'brms' syntax as: brm(bf(Conc_log ~ 1 + s(Time_log1, k = 6) + s(Time_log1, Subject, bs Regarding the expected value   , the model is equivalent.Instead of constant conditional variation this model included a linear effect of time (on the   scale) on the conditional standard deviation σ.
In this distributional model (Kneib, Silbersdorff and Säfken, 2023) both an intercept and linear coefficient for the effect of time on σ are estimated, on the   scale.Prior distributions were (−1,1) for  ,0 and (0,2) for  , , intended to denote that the expected σ is probably not larger than  1 , and the effect of a 1-unit change in time on the   scale is most plausibly between -2 and 2 units; other priors were as for the constant variation model.
Overall the distributional model generated similar posterior predictions of the time-course of theophylline concentrations (Figure 3), to the constant-variation model (Figure 1).As the residual variation declined with time in this model, the fit is strongest for the later observations and weakest for the earliest observations, which may be reasonable given the greater complexity of the underlying biological processes earlier in time, but the scale of the conditional variation seems excessively small in this case.Residual analysis from the distributional model showed that the predicted conditional variation was not apparently contradicted by the observations (Figure 2).The proportion of described variation remained large (Bayes R 2 : 0.850, 90% CrI: 0.781, 0.904) but is more difficult to directly interpret under the non-constant variation.Though the time-varying model did improve the pattern of residuals, particularly around the early observations, the presence of time-varying error is difficult to explain.It would be surprising that analytic measurement error be time-varying; model lack-of-fit, or error in the exact sampling time, seem like more realistic possibilities.
An additional modelling option is to allow greater flexibility for the form of the subject-level functions.A varying-wiggliness model (Pedersen et al., 2019) requires more parameters than the shared-wiggliness models, but may improve the fit for the individuals by reducing the degree of shrinkage.This model was specified as: The time-course of the individual theophylline concentrations predicted by the model are shown in Figure 4, and are visually of better fit to the individual data than the shared-wiggliness models, particularly considering some of the more poorly-predicted observations (Figure 1).The proportion of variability in concentration explained was comparable (Bayes R 2 : 0.950, 90% CrI: 0.937, 0.960).
Population predictions of the theophylline concentration from each model, conditional on an average subject (i.e.ignoring the subject-specific smooth terms) are expressed in Figure 5.The form of the population smooth functions was similar between the constant-error and varying-wiggliness models, but visibly different for the varying-error model.All were visually reasonable with respect to the aggregate data.Compared with classical graphics for pharmacokinetic data using sample statistics by nominal time, which do not represent a viable model for the expected concentration based on the data (Fossler, 2017), this graphic style using HGAM shows simultaneously the expected concentration for an average subject, its uncertainty, and the variability of the observed subjects.
The key NCA parameter   was obtained via post-processing of the joint posteriors from each of the candidate models, by numeric integration on the interval 0-24h.Posterior distributions, both for the individuals and the population (a hypothetical average subject) were obtained by estimation from each sample.Because the per-sample numerical estimation is somewhat expensive this was fairly slow on a personal computer, taking several hours of computation while utilizing parallel processing by package 'future' (Bengtsson, 2021).For comparison, the   were generated by the trapezoidal method as implemented in the 'ncappc' package (Acharya et al., 2016).
The posterior distributions of population   , from each model, are presented as histograms as implemented with 'ggdist' (Kay, 2022), in Figure 6 along with the trapezoidal   and its confidence interval from percentile bootstrap via package 'boot' (Canty and Ripley, 2022).
Summaries of the population posterior distributions for the   ,   , and   are described in Table 1.The constant and time-varying error models generated similar posterior distributions, and intervals comparable to the bootstrap confidence intervals, while the posterior distributions from the varying-wiggliness model had greater width and skewness.As the subject-level estimates are drawn less strongly to the average in the varying-wiggliness model (less shrinkage), it is reasonable to expect that less information was available to support the population smooth, and its secondary parameters.
Overall the posterior distributions for the individual   converged closely on the point estimates from the trapezoidal method (Figures 7-9), and posterior distributions were narrow, reflecting low, though not negligible, uncertainty.
Noting the weak nature of the priors, especially for the smooth terms, the resulting posterior distributions for the pharmacokinetic parameters may not have strong interpretation in terms of belief, but rather represent mostly likelihood-based uncertainty statements.Compared to the various methods for confidence interval generation from NCA, mostly based on bootstrapping (Jaki, Wolfsegger and Ploner, 2009) which attempt to characterize uncertainty in a population parameter, this approach generates uncertainty statements corresponding to individual parameters.Because these statements are generated from a joint probability distribution for all relevant parameters, a credible region about   , and a dependent credible interval about   , which result from simultaneous uncertainty in multiple primary parameters, are relatively easy (though computationally expensive) to obtain.Similar statements may be available via the penalized likelihood approach in package 'mgcv'.This contrasts with the relatively onerous procedures required for valid simultaneous confidence (i.e.frequentist) statements from population pharmacokinetic (i.e.nonlinear multilevel) models (Kümmel et al., 2018); notwithstanding the arguably poorer utility and interpretability of frequentist uncertainty statements (Morey et al., 2016).Some computational challenges were noted during development of this model.Despite careful selection of the prior distributions, a small step size (adapt_delta) was necessary to avoid divergences (Betancourt, 2020), which are a threat to the validity of the results.In the author's experience this is generally necessary for hierarchical additive models in 'brms', perhaps due to complex posterior geometry.A max_treedepth larger than the default was also required, which reflects sampling inefficiency (Stan Development Team, 2022); this may be resolvable using alternate parameterization, but none were clearly available for this case to the author's knowledge.

Example 2: Implementing Covariates
A popular workflow involving NCA is a 'two-stage' approach.Given some observed concentrationtime data, the first stage is to estimate pharmacokinetic parameters for each subject-occasion by NCA.
These parameter estimates (especially AUC) are then passed as data to a second stage, generally some linear model.Bioequivalence and bioavailability studies are a typical application.Notably, the first stage does not have to be NCA.Compartmental modelling by non-linear regression could be used instead, in which case the transition to a one-stage ('population') model could be implemented using a parametric, non-linear multilevel model as currently popular in pharmacometrics (Mould and Upton, 2012).
The hierarchical GAM suggests the possibility of an easy-to-implement, multilevel, 'one-stage' approach utilizing NCA, analogously to the relationship between nonlinear regression and nonlinear multilevel models.The theophylline model above, though emphasizing the conditional estimates of parameters for each subject, is a multilevel model; the subjects are variants from an estimated population effect and their estimates are drawn to it by partial pooling.Next, we will pursue an analysis with more explicit parameter estimation, including covariates, to illustrate some capabilities of a multilevel approach.
Warner and colleagues described the pharmacokinetics of meloxicam in dairy cattle after intravenous or oral administration (Warner et al., 2020b).The cows were at either mid-lactation or post-partum at the time of the study.Data were obtained from the repository (Warner et al., 2020a).From these data, a HGAM can be defined to estimate the effects of lactation stage  and route of administration  on meloxicam concentration over time.The model was be specified as: (  )~( 0 , ) + (  ( + 1),  * ) +   (  ( + 1)) + (0, ) The specification of this model in R and 'brms' syntax was: The terms have a similar meaning to the theophylline model above using the shared-wiggliness, constant-error structure.The additional  *  term denotes, as in the linear model setting, an interaction of the  and  predictors; as well as the additive effect of each, the combination of both variables has an additional unique effect.In other words, there is a distinct population smooth predicted for each of the 4 combinations of the 2 binary predictors, and the predicted concentration-time relationship for some subject is implemented as the sum of one of these population smooths and a subject-specific smooth, where the subject-specific smooths vary commonly across categories.
Priors for the model were defined as in the theophylline model.An important note on the priors for this model is that they contain no information that reflects an expectation that the bioavailability of the oral route is lower than that of the intravenous route.Information about this feature must be obtained entirely from the data.Of course, no such prior information is available in any frequentist approach to this question (the majority of pharmacokinetic studies), so this is not ultimately a loss, except in the sense of the lost opportunity to provide additional information to the analysis.In the case where strong prior information is desired, fully-parametric pharmacokinetic models are available in which parameters have direct interpretation (Wakefield, 1996;Margossian, Zhang and Gillespie, 2022).
The presence of censoring in the concentration data is a complicating issue (Woodward and Whittem, 2019).The data contain numerous left-censored observations which represent meloxicam concentration below the lower performance limit of the analytic method.Fortunately, it is relatively simple to define these observations as left-censored in the model specification using the cens term (Bürkner, 2017).This implements the censoring in the likelihood function (Stan Development Team, 2022) such that the likelihood contribution from the censored observation is equal for any value below the censoring limit (representing that any value below the limit should be considered equally plausible), comparable to the method highlighted by Beal (2001).
A feature of this dataset is that the absorption phase was poorly observed; the earliest observations have concentrations well above the censoring limit.This initial phase is likely to be quite impactful on the pharmacokinetic parameter estimates.In the absence of a pre-dosing observation, the model is not constrained to any particular prediction at time zero, which may generate problematic predictions under the experimental assumption that drug is absent at time zero.The censoring model presents an opportunity to resolve this.By adding a left-censored observation at time zero, the model will be forced to conform to the a priori expected function shape with a starting point somewhere below the smallest drug concentration that could be detected.This synthetic data behaves essentially as prior information that expresses the knowledge (or the assumption) that drug concentration was below the censoring limit before drug was administered.This seems analogous to catalytic priors (Huang et al., 2022), in which the synthetic data is drawn from predictions by a simpler model.Of course, in cases where time zero observations were made, those observations can be utilized in the same way.The estimation of parameters followed the same procedure as the theophylline model.No particular complications were encountered during sampling except for the requirement for a small step size to avoid divergent transitions.The conditional concentration-time relationship for each subject, highlighting their covariates, is expressed in Figure 10.The shape of the predicted function differed clearly between the IV and PO data, with the synthetic time zero observations apparently successful in generating the a-priori expected form of the absorption phase.Overall the model appeared to have reasonably low uncertainty in the conditional concentration-time function for each subject, though some apparently extreme observations were poorly fitted, reflecting the influence of shrinkage (Figure 10).Visualization of the 'population' smooths, conditional on the subject-level smooths set to zero, was utilized to highlight the covariate effects.These population smooths represent the model's estimate of the underlying 'average' concentration-time relationship, by covariate condition.Those estimates, with credible regions, are expressed in Figure 11.The graphic shows a distinct population smooth function for each combination of the covariates, and a reasonable overall fit to the complete data.As for the theophylline example, this graphic style is a clean way to represent both the observed data and the estimated population effects, without the compromises inherent in traditional presentations (Fossler, 2017).
A clear objective for the analysis is to express the covariate effects on the pharmacokinetic parameters.As for the theophylline model, the posterior distributions for the individual (conditional on subject)   (from 0 to 144h) were visualized as histograms (Figure 12), highlighting the covariates.Using posterior predictions from the population smooths, conditional on subject-level smooths set to zero, the covariate effects were expressed as contrasts in predicted   .The Bayesian approach readily facilitates generation of such contrasts, via working directly with the samples from the posterior distribution.The approximate   was obtained by numeric integration of the predicted concentration-time function, for each posterior sample, for each condition (taking several hours of computation).Comparisons were generated of the   between stage for each of IV and PO administration , and between route for each of mid-lactation and post-partum .Comparisons of   were generated between route, for PO administration.The resulting posterior distributions for   by condition, and for selected   ratios (contrasts) by condition, are expressed in figure 13.
For comparison, a two-stage approach using a conventional NCA as implemented in the 'ncappc' package (Acharya et al., 2016), and a linear model, was applied.Censored observations were dropped.The resulting linear predictions of the geometric mean   by condition, and their ratios, were generated by percentile bootstrap as implemented in package 'boot' (Canty and Ripley, 2022) along with confidence intervals.Notably, these uncertainty statements do not correspond closely to the population predictions from the HGAM model.The HGAM population estimates correspond to the estimate of the   of underlying population function, i.e. the expected   for a hypothetical average subject, whereas the two-stage estimates correspond to the geometric mean of the subjects'   , i.e. the observed subjects'   on average.A comparable statement from the HGAM can be generated easily, after some computation; the posterior distribution of the conditional   for each subject was generated, and the geometric mean and ratios for those   result in posterior distributions for those parameters.These are expressed in figure 14.
Summarized as point estimates and intervals, these are of comparable value to the two-stage bootstrap results (noting of course that these credible intervals are not equivalent to confidence intervals; Morey et al., 2016), and both are much narrower than the population estimate from the HGAM.

Discussion
These explorations suggest potential utility of the HGAM approach for pharmacokinetic data analysis.
Bringing the partial pooling, likelihood-based parameter estimation, and uncertainty quantification that are generally features of population pharmacokinetics to NCA, the method facilitates a more modern, single-stage statistical approach utilizing all the available information.In these examples, the estimation of parameters was apparently accurate in comparison to a conventional workflow.These models were constructed via a Bayesian approach, allowing posterior distributions for functions of the model, such as areas or maxima, to be easily obtained as functions of the samples.As the prior information is overall quite weak, especially regarding the parameters of the smooth functions which lack a clear interpretation (and were simply left at the package defaults), the resulting posteriors may lack a strong interpretation in terms of belief.However, as a pragmatic choice this may be reasonable, if the model is computationally efficient enough.
Examples focussed on the analysis of multiple-subject data, analogously to the typical workflow for nonlinear multilevel modelling.Though it is possible to simply generate a separate GAM for each subject, in principle allowing for some features such as uncertainty quantification and censoring, this is not likely to be an efficient usage of the data in most instances, so was not characterized in detail here.Similarly possible is the implementation of HGAM as the first stage in an otherwise two-stage approach, perhaps for evaluation of covariates, in which case the potential advantages compared to the trapezoidal method include partial pooling, implementation of censoring, and uncertainty quantification.This may be particularly advantageous in sparse-data situations, if NCA is desired rather than a parametric model.Additionally, uncertainty in the pharmacokinetic parameter estimates, expressed for example as the posterior standard deviation, could be transferred to the second stage via a measurement error model, or by multiple imputation.Though joint models would more readily facilitate downstream usage of the parameter estimates from this model in some other estimation, to my knowledge this is currently technically limited by the pre-processing generation of smooth terms in 'brms' which is not easily accessible in custom Stan programs, precluding definitions of smooths which depend on parameters.Future developments of HGAM directly in Stan would facilitate such applications.Of course, spline models such as penalized splines may be directly implemented by the user in Bayesian analysis software, as by Jullion et al., (2009), but this undermines accessibility for non-mathematical users.
Interesting extensions to the model could be considered in future.The thin-plate regression spline (Wood, 2003) as applied for these analyses is capable of modelling interactions of smooths; this could, for example, be used to implement an effect of dose, in cases where the influence of dose variability on the concentration-time relationship is required or of interest.The form of the dose-exposure relationship may present a challenge in this case, suggesting the need for careful selection of priors or spline degrees of freedom.A key application for NCA is the evaluation of bioequivalence (BE).Similarly to the MSITAR approach explored by Willemsen et al., (2017), the HGAM approach may be suitable for BE studies, where formulation or route are implemented as predictors, similarly to the meloxicam example.For more complex BE designs, with crossover for example, extensions of this model for additional levels to implement within-subject varying splines would be required.Such models could be well-suited to BE statements, either conditional on the tested subjects or conditional on the population smooth, potentially acting as an intermediate case, between the two-stage approach using typical NCA and a population pharmacokinetic approach with a fully-parametric model (Dubois et al., 2011).Population BE statements marginalizing the between-subject variation as suggested by (Willemsen et al., 2017), from 'brms' models are a more complex consideration, especially regarding the potential cost of computation and the random-effects definition of the spline terms to conduct this marginalization in post-processing (Wiley and Hedeker, 2022).Estimation of model-predicted AUC under the current workflow is computationally expensive due to multiple calls to the 'brms' model during numeric integration, which could perhaps be improved using a trapezoidal method instead.
A specific limitation of the model warrants specific mention.These splines perform poorly in extrapolation beyond the observed data, simply because the form of the function no longer has any support.Though it is computationally feasible to simply generate  ∞ by numeric integration from the fitted model just as for any other endpoint, in practice the improper integrals, conditional on subject, converged poorly across samples, generating sufficient divergent estimates in many cases that these posterior distributions were of suspect reliability.In typical NCA the extrapolation to infinity is achieved using a simple linear function on the log-scale; this was previously implemented in a spline approach to NCA (Jullion et al., 2009), but its validity for the current HGAM model is not clear.
Future evaluations could consider shape-constrained additive models (SCAM) which present the possibility of enforcing monotone or concave splines for example, including for multiple predictors (Pya and Wood, 2015).Of course, in studies with suitably sensitive analytic methods and observations over a sufficient length of time,   approaches  ∞ , such that resulting bias in simply using   is likely to be small and generally irrelevant, though this is a matter for careful evaluation by the user.Though this study focussed on 'primary' parameters ,   , and   , some of the more complex NCA parameters (Veng- Pedersen, 1984), including those determined from multiple dosing events, would fit naturally in the multilevel approach as they could be jointly estimated, and could be considered in future applications.
In summary, the proposed HGAM approach to PK analysis overcomes many of the disadvantages of more traditional two-stage approaches for descriptive PK studies, but without the demands of model selection imposed by typical non-linear multilevel modelling.Compared to previous demonstrations of semi-parametric modelling for NCA (Park et al., 1997;Jullion et al., 2009;Willemsen et al., 2017), the HGAM approach appears to be at least as flexible, whilst also being supported by an accessible open-source implementation.Priorities for further evaluation of this model include application to sparse datasets, inclusion of dose variation, and extension to within-subject variation via additional grouping levels.1: posterior-predicted theophylline population   ,   , and   .Columns are the percentiles of the posterior distribution (credible limits and posterior median).m2: constant-error, constant-wiggliness model, sd: time-varying-error, constant-wiggliness model, by: varying wiggliness model.

Parameter
Results are from the HGAM model for theophylline concentration over time, with subject-varying wiggliness.The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022). .Interval estimates are credible interval limits (posterior percentiles), and for the trapezoidal method, confidence limits from percentile bootstrap † .
Results are from the HGAM model for meloxicam concentration over time with covariates expressing dose and stage effects.The raw data were obtained from Warner et al., (2020a).Results are from the HGAM model for theophylline concentration over time, with constant residual error.The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022).

598
Figure 2: posterior mean residuals and conditional variation functions for the constant-error ('m2'), time-varying error ('sd'), and subject-varying wiggliness ('by') models, by time on either linear or pseudologarithmic scale.Solid lines in residual scatterplots represent ±1 and ±2 times (at the posterior mean) of the conditional variation.Solid line on the conditional standard deviation line plots represent the posterior mean of the conditional standard deviation on the log-scale, and the grey field its 90% credible interval.
Results are from the HGAM models for theophylline concentration over time.The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022).Results are from the HGAM model for theophylline concentration over time, with time-varying conditional variation (residual error).The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022).Results are from the HGAM model for theophylline concentration over time, with subject-varying wiggliness.The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022).This style of graphic is a possible replacement of the classic graphic style for a single-dose PK experiment, with the sample mean by nominal time as the central points, and the error bars its confidence interval, which is almost entirely uninformative (Fossler, 2017); no data is presented, and the confidence intervals for the sample means do not represent a sensible model.
Results are from the HGAM model for theophylline concentration over time.The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022).Results are from HGAM models for theophylline concentration over time.The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022).Results are from the HGAM model for theophylline concentration over time, with subject-varying wiggliness.The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022).Results are from the HGAM model for theophylline concentration over time, with subject-varying wiggliness.The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022).Results are from the HGAM model for theophylline concentration over time, with subject-varying wiggliness.The raw data were from the classical theophylline dataset reported by Upton, utilized by Beal et al., (2009) and Pinheiro and Bates (1995), and obtained here from R (R Core Team, 2022).
Figure 10: posterior-predicted meloxicam concentrations over time, conditional on subject, for intravenous or oral administration, in mid-lactation or post-partum cows.The line represents the posterior mean prediction, the shaded region the 90% credible intervals (which are narrow and sit close to the line on this scale), the filled circles the observations of meloxicam plasma concentration, and the open circles the left-censored observations (concentration below the analytic limit of quantification).Colors represent the covariates for route and stage effects.Note that the observations for PO at time 0 are synthetic and intended to express a prior belief that the concentration at time zero was too low, in principle, to be measured.
Results are from the HGAM model for meloxicam concentration over time with covariates expressing route and stage effects.The raw data were obtained from Warner et al., (2020a).
Figure 11: posterior-predicted population meloxicam concentrations over time (conditional on an 'average' subject), for intravenous or oral administration, in mid-lactation or post-partum cows, with concentration on logarithmic scale.The line represents the population posterior mean prediction, the shaded regions the 50% and 90% credible intervals, the filled circles the observations of meloxicam plasma concentration, and the open circles the left-censored observations (concentration below the analytic limit of quantification).The grey lines join the observations by subject.Note that the observations for PO at time 0 are synthetic and intended to express a prior belief that the concentration at time zero was too low, in principle, to be measured.
Results are from the HGAM model for meloxicam concentration over time with covariates expressing route and stage effects.The raw data were obtained from Warner et al., (2020a).Results are from the HGAM model for meloxicam concentration over time with covariates expressing route and stage effects.The raw data were obtained from Warner et al., (2020a).
Figure 13: in the upper panel, posterior distributions, as 100 dots, of the population-predicted (conditional on a hypothetical 'average' subject, i.e. ignoring the between-subject variation), dosecorrected (  :dose), for intravenous or oral administration, in mid-lactation or post-partum cows.In the lower panel, the population-predicted   ratio (relative exposure), for intravenous or oral administration, in mid-lactation or post-partum cows; the solid line represented   ratio of 1 (no effect).The points and intervals underneath the distributions represent the posterior mean, 50% and 90% credible intervals (posterior quantiles).
Results are from the HGAM model for meloxicam concentration over time with covariates expressing dose and stage effects.The raw data were obtained from Warner et al., (2020a).
Figure 14: in the upper panel, posterior distributions, as 100 dots, of the predicted mean   , conditional on the observed subjects (i.e.estimated average of the studied subjects), dosecorrected (  :dose), for intravenous or oral administration, in mid-lactation or post-partum cows.In the lower panel, the predicted mean   ratio (relative exposure) conditional on the observed subjects, for intravenous or oral administration, in mid-lactaion or post-partum cows; the solid line represented   ratio of 1 (no effect).The points and intervals underneath the distributions represent the posterior mean, 50% and 90% credible intervals (posterior quantiles).
Results are from the HGAM model for meloxicam concentration over time with covariates expressing dose and stage effects.The raw data were obtained from Warner et al., (2020a).

Figure 1 :
Figure 1: posterior-predicted theophylline concentrations over time, on linear or logarithmic time scale, conditional on subject.The line represents the posterior mean prediction, the shaded regions the 50% and 90% credible intervals, and the filled circles the observations of theophylline plasma concentration.

Figure 3 :
Figure 3: posterior-predicted theophylline concentrations over time, on linear or logarithmic time scale, conditional on subject.The line represents the posterior mean prediction, the shaded regions the 50% and 90% credible intervals, and the filled circles the observations of theophylline plasma concentration.

600 601Figure 4 :
Figure 4: posterior-predicted theophylline concentrations over time, on linear or logarithmic time scale, conditional on subject.The line represents the posterior mean prediction, the shaded regions the 50% and 90% credible intervals, and the filled circles the observations of theophylline plasma concentration.

Figure 5 :
Figure 5: posterior-predicted population theophylline concentrations over time from the HGAM (various models by panel).The line represents the population posterior prediction (population smooth), the shaded regions the 50% and 90% credible intervals, the points the observations of theophylline concentration, and the grey lines join the observations by subject.

Figure 6 :
Figure6: estimation of the theophylline population   under various models.Dot histograms (100 dots) represent posterior distributions of the theophylline   (from the global smoother), the solid points the posterior mean, and the bars the 50A% and 90% credible intervals.The solid point is the sample mean   via the trapezoidal method, and the error bar its 90% confidence interval from percentile bootstrap.

Figure 7 :
Figure 7: posterior distributions of the theophylline   , conditional on subject.The posterior distribution is represented as the histogram.The point represents the posterior mean   and the bar the 90% credible interval.The dashed line is the corresponding point estimate from the trapezoidal method (conventional NCA).

605 606Figure 8 :
Figure 8: posterior distributions of the theophylline   , conditional on subject.The posterior distribution is represented as the histogram.The point represents the posterior mean   and the bar the 90% credible interval.The dashed line is the corresponding point estimate from the trapezoidal method (conventional NCA).

Figure 9 :
Figure 9: posterior distributions of the theophylline   , conditional on subject.The posterior distribution is represented as the histogram.The point represents the posterior mean   and the bar the 90% credible interval.The dashed line is the corresponding point estimate from the trapezoidal method (conventional NCA).

Figure 12 :
Figure 12: posterior distributions of the meloxicam   (dose-corrected to 1mg/kg), conditional on subject.The posterior distribution is represented as the histogram.The point represents the posterior mean   and the bar the 90% credible interval.The dashed line is the corresponding point estimate from the trapezoidal method (conventional NCA), in which censored observations were dropped.Colors represent the covariates for route and stage effects.
is illustrated for common PK studies including simple descriptive, and parallel designs with covariates, via two practical examples utilizing published data.Strategies for statistical inference and visualization of results are explored.

Table 2 :
posterior-predicted meloxicam   by STAGE and ROUTE conditions, for intravenous (IV) or oral (PO) administration, in mid-lactation (ML) or post-partum (PP) cows.Point estimates are the posterior median of the population  A , posterior median of the geometric mean   across subjects B , or maximum likelihood estimate of the geometric mean   C .Interval estimates are credible interval limits, and for the trapezoidal method, confidence limits from percentile bootstrap * .Results are from the HGAM model for meloxicam concentration over time with covariates expressing dose and stage effects.The raw data were obtained fromWarner et al., (2020a).

Table 3 :
posterior-predicted meloxicam   ratio (contrasts) by STAGE and ROUTE conditions, for intravenous (IV) or oral (PO) administration, in mid-lactation (ML) or post-partum (PP) cows.Point estimates are the posterior median of ratios of population   A , posterior median of the geometric mean   across subjects B , or maximum likelihood estimate of the geometric mean   C