You might have a mannequin that’s nonlinear within the parameters. Maybe it’s a mannequin of tree progress and due to this fact asymptotes to a most worth. Maybe it’s a mannequin of serum concentrations of a drug that rise quickly to a peak focus after which decay exponentially. Straightforward sufficient, use nonlinear regression ([R] nl) to suit your mannequin. However … what when you have repeated measures for every tree or repeated blood serum ranges for every affected person? You would possibly wish to account for the correlation inside tree or affected person. You would possibly even imagine that every tree has its personal asymptotic progress. You want nonlinear mixed-effects fashions—additionally referred to as nonlinear hierarchical fashions or nonlinear multilevel fashions.

The menl command, launched in Stata 15, matches NLME fashions. Classical nonlinear fashions assume there may be one remark per topic and that topics are impartial. You possibly can consider NLME fashions as an extension of nonlinear fashions to the case the place a number of measurements could also be taken over a topic and these within-subject observations are normally correlated. You might also consider NLME fashions as a generalization of linear mixed-effects fashions the place some or all random results enter the mannequin in a nonlinear trend. No matter the way you consider them, NLME fashions are used to explain a response variable as a (nonlinear) perform of covariates, accounting for the correlation amongst observations on the identical topic.

On this weblog, I’ll stroll you thru a number of examples of becoming nonlinear mixed-effects fashions utilizing menl. I’ll first match a nonlinear mannequin to the information on the expansion of bushes, ignoring the correlation amongst measurements on the identical tree. I’ll then show alternative ways of accounting for this correlation and the way to incorporate random results into the mannequin parameters to offer the parameters tree-specific interpretation.

We’ve knowledge (Draper and Smith 1998) on trunk circumference of 5 completely different orange bushes the place trunk circumference (in mm) was measured on seven completely different events, over roughly a four-year interval of progress. We wish to mannequin the expansion of orange bushes. Let’s plot the information first.

. webuse orange
. twoway scatter circumf age, join(L) ylabel(#6)

There may be some variation within the progress curves, however the identical general form is noticed for all bushes. Particularly, the expansion tends to stage off towards the top. Pinheiro and Bates (2000) instructed the next nonlinear progress mannequin for these knowledge:

[
{mathsf{circumf}}_{ij}=frac{phi_{1}}{1+expleft{ -left( {mathsf{age}}_{ij}-phi_{2}right)/phi_{3}right}} + epsilon_{ij}, quad j=1,dots,5; i=1,dots,7
]

Parameters of NLME fashions usually have scientifically significant interpretations and kind the idea of analysis questions. On this mannequin, (phi_{1}) is the common asymptotic trunk circumference of the bushes as (mathsf{age}_{ij} to infty), (phi_{2}) is the common age at which the tree attains half the asymptotic trunk circumference (phi_{1}) (half-life), and (phi_{3}) is a scale parameter.

For now, let’s ignore the repercussions (much less correct estimates of uncertainty in parameter estimates, probably much less highly effective assessments of hypotheses, and fewer correct interval estimates) of not accounting for the correlation and potential heterogeneity among the many observations. As an alternative, we are going to deal with all observations as being i.i.d.. We are able to use nl to suit the above mannequin, however right here we are going to use menl. (As of 15.1, menl now not requires that you just embrace random results within the mannequin specification.)

. menl circumf = {phi1}/(1+exp(-(age-{phi2})/{phi3})), stddeviations

Acquiring beginning values:

NLS algorithm:

Iteration 0:   residual SS =  17480.234
Iteration 1:   residual SS =  17480.234

Computing normal errors:

Combined-effects ML nonlinear regression           Variety of obs     =         35
Log Chance = -158.39871
------------------------------------------------------------------------------
     circumf |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       /phi1 |   192.6876   20.24411     9.52   0.000     153.0099    232.3653
       /phi2 |   728.7564   107.2984     6.79   0.000     518.4555    939.0573
       /phi3 |   353.5337   81.47184     4.34   0.000     193.8518    513.2156
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
                sd(Residual) |   22.34805   2.671102      17.68079    28.24734
------------------------------------------------------------------------------

The stddeviations possibility instructs menl to report the usual deviation of the error time period as a substitute of variance. As a result of one progress curve is used for all bushes, the person variations as depicted within the above graph are included within the residuals, thus inflating the residual normal deviation.

We’ve a number of observations from the identical tree which might be prone to be correlated. One approach to account for the dependence between observations throughout the similar tree is by together with a random impact, (u_j), shared by all of the observations throughout the (j)th tree. Thus, our mannequin turns into

$$
mathsf{circumf}_{ij}=frac{phi_{1}}{1+expleft{ -left( mathsf{age}_{ij}-phi_{2}proper)/phi_{3}proper}} + u_j + epsilon_{ij}, quad j=1,dots,5; i=1,dots,7
$$

This mannequin could also be match by typing

. menl circumf = {phi1}/(1+exp(-(age-{phi2})/{phi3}))+{U[tree]}, stddev

Acquiring beginning values by EM:

Alternating PNLS/LME algorithm:

Iteration 1:    linearization log probability = -147.631786
Iteration 2:    linearization log probability = -147.631786

Computing normal errors:

Combined-effects ML nonlinear regression           Variety of obs     =         35
Group variable: tree                            Variety of teams  =          5

                                                Obs per group:
                                                              min =          7
                                                              avg =        7.0
                                                              max =          7
Linearization log probability = -147.63179
------------------------------------------------------------------------------
     circumf |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       /phi1 |   192.2526   17.06127    11.27   0.000     158.8131    225.6921
       /phi2 |   729.3642   68.05493    10.72   0.000      595.979    862.7494
       /phi3 |    352.405   58.25042     6.05   0.000     238.2363    466.5738
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
tree: Identification               |
                       sd(U) |   17.65093   6.065958      8.999985    34.61732
-----------------------------+------------------------------------------------
                sd(Residual) |    13.7099    1.76994      10.64497    17.65728
------------------------------------------------------------------------------

The fixed-effects estimates are related in each fashions, however their normal errors are smaller within the above mannequin. The residual normal deviation can also be smaller as a result of a few of the particular person variations are actually being accounted for by the random impact. Let (y_{sj}) and (y_{tj}) be two observations on the (j)th tree. The inclusion of the random impact (u_j) implies that ({rm Cov}(y_{sj},y_{tj}) > 0) (in actual fact, it is the same as ({rm Var}(u_j))). Subsequently, (u_j) induces a constructive correlation amongst any two observations throughout the similar tree.

There may be one other approach to account for the correlation amongst observations inside a gaggle. You possibly can mannequin the residual covariance construction explicitly utilizing menl‘s rescovariance() possibility. For instance, the above random-intercept nonlinear mannequin is equal to a nonlinear marginal mannequin with an exchangeable covariance matrix. Technically, the 2 fashions are equal solely when the correlation of the within-group observations is constructive.

. menl circumf = {phi1}/(1+exp(-(age-{phi2})/{phi3})),
    rescovariance(exchangeable, group(tree)) stddev

Acquiring beginning values:

Alternating GNLS/ML algorithm:

Iteration 1:    log probability = -147.632441
Iteration 2:    log probability = -147.631786
Iteration 3:    log probability = -147.631786
Iteration 4:    log probability = -147.631786

Computing normal errors:

Combined-effects ML nonlinear regression           Variety of obs     =         35
Group variable: tree                            Variety of teams  =          5

                                                Obs per group:
                                                              min =          7
                                                              avg =        7.0
                                                              max =          7
Log Chance = -147.63179
------------------------------------------------------------------------------
     circumf |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       /phi1 |   192.2526   17.06127    11.27   0.000     158.8131    225.6921
       /phi2 |   729.3642   68.05493    10.72   0.000      595.979    862.7494
       /phi3 |    352.405   58.25042     6.05   0.000     238.2363    466.5738
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
Residual: Exchangeable       |
                          sd |   22.34987    4.87771      14.57155    34.28026
                        corr |   .6237137   .1741451      .1707327    .8590478
------------------------------------------------------------------------------

The fixed-effects portion of the output is an identical within the two fashions. These two fashions additionally produce the identical estimates of the marginal covariance matrix. You possibly can confirm this by working estat wcorrelation after every mannequin. Though the 2 fashions are equal, by together with a random impact on bushes, we cannot solely account for the dependence of observations inside bushes, but in addition estimate between-tree variability and predict tree-specific results after estimation.

The earlier two fashions indicate that the expansion curves for all bushes have the identical form and differ solely by a vertical shift that equals (u_j). This can be too sturdy an assumption in apply. If we glance again on the graph, we are going to discover there may be an growing variability within the trunk circumferences of bushes as they strategy their limiting peak. Subsequently, it could be extra affordable to permit (phi_1) to fluctuate between bushes. This can even give our parameters of curiosity a tree-specific interpretation. That’s, we assume

$$
mathsf{circumf}_{ij}=frac{phi_{1}}{1+expleft{ -left( mathsf{age}_{ij}-phi_{2}proper)/phi_{3}proper}} + epsilon_{ij}
$$

with

$$
phi_1 = phi_{1j} = beta_1 + u_j
$$

Chances are you’ll interpret (beta_1) as the common asymptotic peak and (u_j) as a random deviation from the common that’s particular to the (j)th tree. The above mannequin could also be match as follows:

. menl circumf = ({b1}+{U1[tree]})/(1+exp(-(age-{phi2})/{phi3}))

Acquiring beginning values by EM:

Alternating PNLS/LME algorithm:

Iteration 1:    linearization log probability = -131.584579

Computing normal errors:

Combined-effects ML nonlinear regression           Variety of obs     =         35
Group variable: tree                            Variety of teams  =          5

                                                Obs per group:
                                                              min =          7
                                                              avg =        7.0
                                                              max =          7
Linearization log probability = -131.58458
------------------------------------------------------------------------------
     circumf |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /b1 |    191.049   16.15403    11.83   0.000     159.3877    222.7103
       /phi2 |    722.556   35.15082    20.56   0.000     653.6616    791.4503
       /phi3 |   344.1624   27.14739    12.68   0.000     290.9545    397.3703
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
tree: Identification               |
                     var(U1) |   991.1514   639.4636      279.8776    3510.038
-----------------------------+------------------------------------------------
               var(Residual) |   61.56371   15.89568      37.11466    102.1184
------------------------------------------------------------------------------

If desired, we will moreover mannequin the within-tree error covariance construction utilizing, for instance, an exchangeable construction:

. menl circumf = ({b1}+{U1[tree]})/(1+exp(-(age-{phi2})/{phi3})),
    rescovariance(exchangeable)

Acquiring beginning values by EM:

Alternating PNLS/LME algorithm:

Iteration 1:    linearization log probability = -131.468559
Iteration 2:    linearization log probability = -131.470388
Iteration 3:    linearization log probability = -131.470791
Iteration 4:    linearization log probability = -131.470813
Iteration 5:    linearization log probability = -131.470813

Computing normal errors:

Combined-effects ML nonlinear regression           Variety of obs     =         35
Group variable: tree                            Variety of teams  =          5

                                                Obs per group:
                                                              min =          7
                                                              avg =        7.0
                                                              max =          7
Linearization log probability = -131.47081
------------------------------------------------------------------------------
     circumf |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /b1 |   191.2005   15.59015    12.26   0.000     160.6444    221.7566
       /phi2 |   721.5232   35.66132    20.23   0.000     651.6283    791.4182
       /phi3 |   344.3675   27.20839    12.66   0.000     291.0401     397.695
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
tree: Identification               |
                     var(U1) |   921.3895    582.735      266.7465    3182.641
-----------------------------+------------------------------------------------
Residual: Exchangeable       |
                         var |   54.85736   14.16704      33.06817    91.00381
                         cov |  -9.142893   2.378124     -13.80393   -4.481856
------------------------------------------------------------------------------

We didn’t specify possibility group() within the above as a result of, within the presence of random results, rescovariance() mechanically determines the lowest-level group based mostly on the required random results.

The NLME fashions we used up to now are all linear within the random impact. In NLME fashions, random results can enter the mannequin nonlinearly, identical to the fastened results, they usually usually do. As an example, along with (phi_1), we will let different parameters fluctuate between bushes and have their very own random results:

. menl circumf = {phi1:}/(1+exp(-(age-{phi2:})/{phi3:})),
    outline(phi1:{b1}+{U1[tree]})
    outline(phi2:{b2}+{U2[tree]})
    outline(phi3:{b3}+{U3[tree]})
(output omitted)

The above mannequin is just too sophisticated for our knowledge containing solely 5 bushes and is included purely for demonstration functions. However from this instance, you may see that the outline() possibility is beneficial when you could have a sophisticated nonlinear expression, and also you wish to break it down into smaller items. Parameters specified inside outline() might also be predicted for every tree after estimation. For instance, after we match our mannequin, we might want to predict the asymptotic peak, (phi_{1j}), for every tree (j). Under, we’re requesting to assemble a variable named phi1 that incorporates the expected values for the expression {phi1:}:

. predict (phi1 = {phi1:})
(output omitted)

If we had extra bushes, we might have additionally specified a covariance construction for these random results as follows:

. menl circumf = {phi1:}/(1+exp(-(age-{phi2:})/{phi3:})),
    outline(phi1:{b1}+{U1[tree]})
    outline(phi2:{b2}+{U2[tree]})
    outline(phi3:{b3}+{U3[tree]})
    covariance(U1 U2 U3, unstructured)
(output omitted)

I demonstrated only some fashions on this article, however menl can do far more. For instance, menl can incorporate larger ranges of nesting similar to three-level fashions (for instance, repeated observations nested inside bushes and bushes nested inside planting zones). See [ME] menl for extra examples.

References

Draper, N., and H. Smith. 1998. Utilized Regression Evaluation. third ed. New York: Wiley.

Pinheiro, J. C., and D. M. Bates. 2000. Combined-Results Fashions in S and S-PLUS. New York: Springer.