Existing command menl has new features
for fitting nonlinear mixed-effects models (NLMEMs) that may include
lag, lead (forward), and difference operators. One important class of such
models is the class of pharmacokinetic (PK) models and, specifically,
multiple-dose PK models. menl's new features can
also be used to fit other models, such as certain growth
models and time-series nonlinear multilevel models.

Let's see it work

Growth models

Growth models based on difference equations are often used to model animal
growth. For instance,
Andrews (1982)
used the so-called logistic-by-weight model to model the growth of small
lizards. Wright et al. (2013) applied this model to study the length of
brown anole lizards. The authors also extended the model to account for the
between-lizard variation by including random effects.

Inspired by the above studies, we created our own fictional data on lengths
of small lizards. In what follows, we will use menl to fit a
random-effects logistic-by-weight model to these data.

First, let's take a quick look at the data. We list the observations for the
first two lizards.

. list if id<=2, sepby(id)

id time length

1.

1 103 22.83874

2.

1 208 29.38543

3.

1 310 35.20067

4.

1 411 39.46608

5.

1 512 40.801

6.

2 5 28.97618

7.

2 212 38.83975

8.

2 313 41.70709

Lizards were captured and their length measured. Then, they were
marked and released. After that, the lizards were recaptured at different
times, and their length was recorded again. Variable time contains the
number of days since the beginning of the study, and length records the
length (in mm) of lizards.

The interval version of the random-effects logistic-by-weight model can be
defined as

where LijLij is the length of lizard jj at time ii and
DijDij is the difference in days between measurements i−1i−1 and
ii for lizard jj. Model parameters are lizard-specific growth rate,
rjrj, and asymptotic length, aa. Here
ϵij∼N(0,σ2ϵ)ϵij∼N(0,σϵ2) and rj=r0+ujrj=r0+uj and
uj∼N(0,σ2r)uj∼N(0,σr2), where r0r0 is the mean growth rate,
σ2rσr2 is the between-lizards variance, and σ2ϵσϵ2 is
sampling error variance.

Let's fit the above model using menl. Our menl specification is
as follows.

We specify the above nonlinear model using menl's standard nonlinear
specification with model parameters enclosed in curly braces, {}. New
in the above specification are the uses of the lag operator L. (with
variable length) and the difference operator D. (with variable
time).

The define() option models the lizard-specific growth rate, rjrj,
as a random effect, ujuj, plus a constant, r0r0. tsorder()
specifies the time variable that will be used by the L. and D.
operators to create the lagged and differenced versions of the variables.
Variable id in the term U[id] identifies the grouping structure
of our dataset. initial() specifies the initial value for the
asymptotic length.

The estimated asymptotic length is 42.5 mm with a standard error of 0.16,
and the mean growth rate is 0.0097 with a standard error of 0.0002. The
between-lizards standard deviation for the growth rate is 0.0007 with a 95% CI
of [0.0005, 0.0011], which suggests that there is variability in lizards' growth
rates compared with their mean value.

Let's visualize the predicted growth rates for lizards. We can use
predict to compute the predictions not only for the response but
also for any parameter or function defined in define(). So we use
predict to compute lizard-specific and mean growth rates.

The growth rates of lizards vary around the mean value of 0.0097, with the
largest value of 0.0113 for lizard 40 and the smallest value of 0.0083 for
lizard 22.

Multiple-dose PK models

PK models concern drug absorption, distribution, metabolism, and
excretion. Single-dose and multiple-dose models are of special interest.
Each comes in two flavors distinguished by how long it takes the drug to enter
the body: orally or intravenously. menl can fit all those models.

Let's use menl to fit a multiple-dose model with intravenous
administration. Consider a study of the neonatal PKs of phenobarbital
(Grasela and Donn 1985). Each infant received one or more intravenous doses, dose (mg/kg).
One to six blood serum phenobarbital concentration measurements,
conc (mg/L), were obtained from each infant, subject.
Other variables of interest are birthweight and a dichotomized
Apgar score (fapgar), a measure of the physical condition.

A one-compartment open model with intravenous administration and first-order
elimination was used to model the PKs of this phenobarbital study.

Model parameters are the clearance, CljClj (L/h), and volume of distribution,
VjVj (L), for each infant jj. E(conci−1,j)E(conci−1,j) is the expected
concentration at measurement time i−1i−1 for infant jj.

In addition, CljClj and VjVj are parameterized as follows,

CLj=β1weightj×exp(u1j)CLj=β1weightj×exp⁡(u1j)

Vj=β2weightj{1+β3fapgarj×exp(u2j)}Vj=β2weightj{1+β3fapgarj×exp⁡(u2j)}

where u1ju1j and u2ju2j are two independent Gaussian random effects
with zero means and the corresponding variances σ21σ12 and
σ22σ22.

Option define() was used to incorporate the parameterizations for
CljClj and VjVj. For instance, in option define(Cl:...),
{cl:weight} models β1β1 and {U1[subject]} models
u1ju1j.

The difference operator D.time models (timeij−ij−timei−1,ji−1,j).

The expression L.{conc:} models the predicted value for
E(conci−1,j)E(conci−1,j). It represents the predicted value for the
concentration corresponding to the previous time value.

Option tsorder(time) specified the variable, time, that will
be used to determine the time order for the time-series operators.

Option tsinit() specifies the initial condition for the
concentration at time 0.

Option tsmissing specifies that observations containing system
missing values (.) in the outcome conc be retained in the
computation. In this dataset, these observations represent intermittent
measurements at which the dose was administered but the concentration was
not measured.

Heavier babies have a higher clearance and volume of distribution. There is a
positive association between the volume of distribution and the Apgar score:
healthier babies have a better ability to eliminate the drug.

After estimation, let's use predict to predict concentrations.

. predict yhat
(option yhat assumed)

The default prediction is the predicted mean function, including the estimates
of subject-specific random effects. So variable yhat contains
observation-specific mean concentration for each infant.

Let's plot the predicted concentration against the observed one for, say,
the first four infants.

The model predicted the missing intermittent concentrations and provided the
entire trajectory of the concentration for the observed measurement times.

Other features

In addition, menl provides two new within-group error structures:
Toeplitz and banded. They can be specified within options
rescovariance() and rescorrelation() as toeplitz# and banded#, respectively, where # defines
the order of the error structures.

menl also provides new option lrtest, which reports a
likelihood-ratio test comparing the nonlinear mixed-effects model with the
model fit by ordinary nonlinear regression.