What’s this about?

Stata now fits nonlinear mixed-effects models, also known as nonlinear multilevel models and nonlinear hierarchical models. These models can be thought of in two ways. You can think of them as nonlinear models containing random effects. Or you can think of them as linear mixed-effects models in which some or all fixed and random effects enter nonlinearly. However you think of them, the overall error distribution is assumed to be Gaussian.

These models are popular because some problems are not, says their science, linear in the parameters. These models are popular in population pharmacokinetics, bioassays, and studies of biological and agricultural growth processes. For example, nonlinear mixed-effects models have been used to model drug absorption in the body, intensity of earthquakes, and growth of plants.

The new estimation command is menl. It implements the popular-in-practice Lindstrom–Bates algorithm, which is based on the linearization of the nonlinear mean function with respect to fixed and random effects. Both maximum likelihood and restricted maximum-likelihood estimation methods are supported.

Let’s see it work

Let’s look at several examples.

Random intercepts

Random coefficients and random-effects covariance structures

Residual covariance structures: Pharmacokinetic (PK) model

Nonlinear three-level model: CES production function

Random intercepts

menl is a serious estimator for serious problems. Nonetheless, for pedagogical reasons, we will begin with a silly example. You will remember it.

Suppose we are interested in modeling brain weight growth of unicorns in the land of Zootopia. The graph below shows the brain weight profiles of 20 unicorns over time.

The variability in brain weight measurements between unicorns increases with time.

Based on previous unicorn studies, a model for unicorn brain growth is believed to be

Parameter ϕ1 represents asymptotic growth—the unicorn’s brain weight as time increases to infinity. Parameter ϕ2 is the time at which the brain weight reaches half of ϕ1, and ϕ3 is a scale parameter that determines the growth rate, the rate at which the brain weight approaches the asymptotic weight ϕ1.

We will initially assume that only the asymptotic growth ϕ1 is unicorn specific by including a random intercept Uj in the equation of ϕ1 for the jth unicorn:

We fit this model using menl as follows:

. menl weight = ({b1}+{U[id]})/(1+exp(-(time-{b2})/{b3}))



Mixed-effects ML nonlinear regression           Number of obs     =        260
Group variable: id                              Number of groups  =         20

						Obs per group:
							      min =         13
							      avg =       13.0
							      max =         13

Linearization log likelihood = -270.54625

You can see that the specification to menl looks much like the equation we wrote earlier. Parameters to be estimated are enclosed in curly braces: {b1}, {b2}, and {b3}. By typing {U[id]}, we specified a random intercept for each unicorn, identified by the group variable id.

The above syntax uses what we call a single-equation or single-stage specification. menl also allows multistage or hierarchical specifications in which parameters of interest can be defined at each level of hierarchy as functions of other model parameters and random effects, such as

. menl weight = {phi1:}/(1+exp(-(time-{phi2:})/{phi3:})),
    define(phi1: {b1}+{U1[id]}) 
    define(phi2: {b2}+{U2[id]})
    define(phi3: {b3}+{U3[id]})



Mixed-effects ML nonlinear regression           Number of obs     =        260
Group variable: id                              Number of groups  =         20

						Obs per group:
							      min =         13
							      avg =       13.0
							      max =         13

Linearization log likelihood = -164.64495

    phi1:  {b1}+{U1[id]}
    phi2:  {b2}+{U2[id]}
    phi3:  {b3}+{U3[id]}

This is the same model except that ϕ2 and ϕ3 are also allowed to vary across unicorns using their own sets of random intercepts.

Random coefficients and random-effects covariance structures

Continuing with our unicorn example, we know that boosting brain weight has been an active research area in Zootopia for the past two decades, and scientists believe that consuming rainbow cupcakes right after birth may help attain higher asymptotic growth. Hence, covariate cupcake, which represents the number of rainbow cupcakes consumed right after birth, is added to the equation for ϕ1j:

We type

. menl weight = {phi1:}/(1+exp(-(time-{b2})/{b3})),
    define(phi1: {b10}+{b11}*cupcake+{U[id]})


Mixed-effects ML nonlinear regression           Number of obs     =        260
Group variable: id                              Number of groups  =         20

						Obs per group:
							      min =         13
							      avg =       13.0
							      max =         13
Linearization log likelihood =  -264.5152

        phi1:  {b10}+{b11}*cupcake+{U[id]}

Consuming rainbow cupcakes after birth indeed leads to higher asymptotic brain growth: /b11 is roughly 0.74 with a 95% CI of [0.38,1.1].

It is bad karma to end a unicorn story without showing how to specify random coefficients or random slopes. Suppose that the equation for ϕ1j is now as follows:

We incorporated a unicorn-specific random slope for variable cupcake. We also assume that the random effects U0j and U1j are correlated. To fit this model, we type

. menl weight = {phi1:}/(1+exp(-(time-{b2})/{b3})),
    define(phi1: {b10}+{U0[id]}+({b11}+{U1[id]})*cupcake)
    covariance(U0 U1, unstructured)


Mixed-effects ML nonlinear regression           Number of obs     =        260
Group variable: id                              Number of groups  =         20

						Obs per group:
                                                              min =         13
                                                              avg =       13.0
                                                              max =         13
Linearization log likelihood = -263.91462

        phi1:  {b10}+{U0[id]}+({b11}+{U1[id]})*cupcake

In addition to the unstructured covariance, menl supports independent, exchangeable, and identity variance–covariance structures for random effects from the same level of hierarchy.

For more complicated NLME models, specifying expressions containing linear combinations may become tedious. menl offers a convenient shorthand specification to handle linear combinations. For example, the above option define() for phi1 may be replaced with

define(phi1: cupcake U0[id] c.cupcake#U1[id])

Residual covariance structures: Pharmacokinetic (PK) model

menl supports flexible variance–covariance structures to model error heteroskedasticity and its within-group dependence. For example, heteroskedasticity can be modeled as a power function of a covariate or even of predicted mean values, and dependence can be modeled using an autoregressive model of any order.

Consider a PK study in which the drug theophylline was administered orally to 12 subjects with the dosage (mg/kg) given on a per weight basis. Serum concentrations (mg/L) were obtained at 11 time points per subject over 25 hours following administration. The graph below shows the resulting concentration-time profiles for the first four subjects.

A PK model used to model the above concentration-time profiles is

Model parameters are the elimination rate constant ke, the absorption rate constant ka, and the clearance Cl. Parameters Cl and ka are allowed to vary across subjects. Because each of the model parameters must be positive to be meaningful, we write

It is common for the PK data to exhibit within-subject heterogeneity that increases with the magnitude of the response variable. In this case, the within-subject error variance is often modeled as a power function of the predicted mean. menl provides the resvariance() option for this.

. menl conc = (dose*{ke:}*{ka:}/({cl:}*({ka:}-{ke:})))*(exp(-{ke:}*time)-exp(-{ka:}*time)),
    define(cl: exp({b0}+{U0[subject]}))
    define(ka: exp({b1}+{U1[subject]}))
    define(ke: exp({b2}))
    resvariance(power _yhat)


Mixed-effects ML nonlinear regression           Number of obs     =        132
Group variable: subject                         Number of groups  =         12

                                                Obs per group:
                                                              min =         11
                                                              avg =       11.0
                                                              max =         11
Linearization log likelihood = -167.67964

          cl:  exp({b0}+{U0[subject]})
          ka:  exp({b1}+{U1[subject]})
          ke:  exp({b2})

In this example, the 95% CI [-0.2,0.8] for power constant delta contains 0, suggesting that within-subject error variance does not depend on the predicted values.

Nonlinear three-level model: CES production function

Production functions are an important component of macroeconomic models. They tell us how economies use inputs such as labor and capital in the production process. A common example is a constant elasticity of substitution (CES) production function.

The CES function defines how production allocates inputs such as capital and labor. It introduces the constant elasticity of substitution parameter that models the allocation, which makes it a very flexible modeling tool. For instance, CES production functions include Cobb–Douglas, Leontief, and perfect substitutes production functions as special cases.

Let’s look at an example.

Suppose that there are five regions, each with seven production firms manufacturing a certain good. Measurements were collected each year over the course of six years. We want to fit the following CES production model,

In the above, U1k and U2jk are random intercepts at the region and firm-within-region levels.

To fit the above nonlinear three-level model using menl, we type

. menl lnoutput = {b0:}-1/{rho}*ln({delta:}*capital^(-{rho})+(1-{delta:})*labor^(-{rho})),
    define(b0: {b0}+{U1[region]})
    define(delta: {delta}+{U2[region>firm]})


Mixed-effects ML nonlinear regression           Number of obs     =        210

Linearization log likelihood = -159.60141

          b0:  {b0}+{U1[region]}
       delta:  {delta}+{U2[region>firm]}

The values of /rho of .54 and /delta of .53 tell us that this economy is allocating capital and labor in roughly equal proportions in the production process.

Regions and firms within regions appear to explain some of the variability in the log-factor productivity parameter and the share parameter, respectively.

We can use nlcom to estimate the elasticity of substitution for the CES production function in this model.

. nlcom (1/(1 + _b[/rho]))


 _nl_1:  1/(1 + _b[/rho])

The elasticity of substitution is estimated to be .65 for these data. If instead we used, for example, a Cobb–Douglas production function, we would have incorrectly constrained the elasticity of substitution to be 1.

Tell me more

Learn more about Stata’s multilevel mixed-effects models features.

Read more about nonlinear mixed-effects models in the Stata Multilevel Mixed-Effects Reference Manual.