Nonlinear mixed-effects models

Houssein Assaad, PhD

Senior Statistician and Software Developer

StataCorp LLC

SUGM, Colombus OH

19 Jul 2018

Outline

  • What are nonlinear mixed-effects models ? Why and when you should use them ?
  • Exploring the syntax: A decay model
  • Covariance Structures: A Pharmacokinetics example
  • Multilevel specification: A growth model
  • Summary
  •  Nonlinear hierarchical models
  • Goal: To describe a continuous response that evolves nonlinearly over time (or other condition) as a function of covariates, while taking into account the correlation among observations  within the same subject.
  • Nonlinear models are often based on a model for the mechanism producing the response. As a consequence, the model parameters in a nonlinear model generally have a natural physical interpretation.
  • Nonlinear mixed-effects models (NLME), also known as:
  •  Nonlinear multilevel models
  • Popular in agricultural, environmental, and biomedical applications. For example, NLME models have been used in drug absorption in the body, intensity of earthquakes, and growth of plants.

...First, an existential question:

Why would you want to use NLME models ?

Why nonlinear ?

Why mixed-effects ?

ht=β11+exp{(tβ2)/β3}h_t = \frac{\beta_1}{1+\exp{\left\{-\left(t - \beta_2\right)/\beta_3\right\}}}
h_t = \frac{\beta_1}{1+\exp{\left\{-\left(t - \beta_2\right)/\beta_3\right\}}}

Let hth_t be the expected height of a tree at time tt which is modelled using a logistic growth model

  • β1\beta_1 is the asymptotic height as tt\to\infty

Why nonlinear ?

  • β3\beta_3 is the time it takes for a tree to grow from 50% to about 73% of  β1 \beta_1

     

  • β2\beta_2 is the time at which the tree reaches half of its asymptotic height
ht2.3+16.6t44.4t2+56.8t331.5t4+6.3t5h_t \approx -2.3 + 16.6t -44.4t^2 + 56.8t^3-31.5t^4+6.3t^5
h_t \approx -2.3 + 16.6t -44.4t^2 + 56.8t^3-31.5t^4+6.3t^5
ht=31+exp{(t1)/.2}h_t = \frac{3}{1+\exp{\left\{-\left(t - 1\right)/.2\right\}}}
h_t = \frac{3}{1+\exp{\left\{-\left(t - 1\right)/.2\right\}}}
  • The logistic curve

is approximated very well

in the interval [0.4,1.6][0.4, 1.6] by the 5th5^{th} degree polynomial

  • Polynomial approximation is unreliable outside the original interval.
  • Suppose you have a dataset of tree heights measured over time

Why nonlinear ?

  • Interpretability : Coefficients in polynomial model do not have any physical interpretation (e.g. what is the interpretation of the coefficient of t4t^4 ?)
  • Parsimony : Linear polynomial model uses twice as many parameters as the nonlinear model (logistic model)
  • Validity beyond the observed range of the data
  +---------------------+
  |   id      X       Y |
  |---------------------|
  |    1     14     5.6 |
  |    1     21     9.6 |
  |    1     42    15.5 |
  |---------------------|
  |    2     15    12.6 |
  |    2     21    16.4 |
  |    2     42    16.9 |
  |    2     51    18.9 |
  +---------------------+
Y42\texttt{Y}_{42}
\texttt{Y}_{42}
Y31\texttt{Y}_{31}
\texttt{Y}_{31}
Y12\texttt{Y}_{12}
\texttt{Y}_{12}
Yij=β11+exp(Xijβ2)/β3+ϵij\texttt{Y}_{ij} = \frac{\beta_1}{1+\exp{\left(\texttt{X}_{ij} - \beta_2\right)/\beta_3}} + \epsilon_{ij}
\texttt{Y}_{ij} = \frac{\beta_1}{1+\exp{\left(\texttt{X}_{ij} - \beta_2\right)/\beta_3}} + \epsilon_{ij}
  • Data consists of j=1,,Nj=1,\ldots,N subjects, each measured njn_j times
Y11\texttt{Y}_{11}
\texttt{Y}_{11}

Why mixed-effects ?

  • β1\beta_1, β2\beta_2, and β3\beta_3 may be subject-specific, e.g.
β1=β1j=b1+Uj\beta_1 = \beta_{1j} = b_1 + U_j
\beta_1 = \beta_{1j} = b_1 + U_j

fixed effect

random effect

++
+

mixed effect

where

ϵijN(0,σ2)\epsilon_{ij}\sim N\left(0, \sigma^2\right)
\epsilon_{ij}\sim N\left(0, \sigma^2\right)
UjN(0,σu2)U_{j}\sim N\left(0, \sigma^2_u\right)
U_{j}\sim N\left(0, \sigma^2_u\right)

where

Why mixed effects ?

  • Observations within the same subject are correlated
  • Mixed-effects models allow us to account for the correlation of YsjY_{sj} and YtjY_{tj} by including in our model  a random effect UjU_j specific to jthj^{th} subject
  • The inclusion of random effects allow the parameters of interest (e.g. asymptotic growth ) to have subject-specific interpretations. Very important in pharmacokinetics, for example, where the goal is to derive dosage regimen for a specific patient.
cov(Ysj,Ytj)0\mathtt{cov}\left(Y_{sj}, Y_{tj}\right) \neq 0
\mathtt{cov}\left(Y_{sj}, Y_{tj}\right) \neq 0

2 observations on jthj^{th} subject

The Stata 15 menl command

Exploring the syntax

Fun example: Unicorn health study in Zootopia

. webuse unicorn
(Brain shrinkage of unicorns in the land of Zootopia)

. twoway connected weight time if id < 10, connect(L) ///
	title(Brain Shrinkage of unicorn)
     +--------------------------+
     | id       time     weight |
     |--------------------------|
  1. |  1          0   8.266746 |
  2. |  1   .1666667   7.163768 |
  3. |  1   .3333333   6.725025 |
  4. |  1         .5    6.38325 |
  5. |  1   .6666667    6.00481 |
     |--------------------------|
776. | 60   1.333333    3.75486 |
777. | 60        1.5   3.581552 |
778. | 60   1.666667   3.701847 |
779. | 60   1.833333   3.836078 |
780. | 60          2   3.729806 |
     +--------------------------+
. list id time weight if (_n <=5 | _n >775), sep(5)
  • Based on previous studies, a model for unicorn brain shrinkage is believed to be
weightij=β1+(β2β1)exp(β3timeij)+ϵij\texttt{weight}_{ij} = \beta_1 + \left(\beta_2 - \beta_1\right)\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}
\texttt{weight}_{ij} = \beta_1 + \left(\beta_2 - \beta_1\right)\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}

where ϵijN(0,σ2)\epsilon_{ij}\sim N\left(0,\sigma^2\right)

  • β1\beta_1 is average asymptotic brain weight of unicorns as timeij\texttt{time}_{ij} \to \infty
  • β2\beta_2 is average brain weight at birth (at timeij=0\texttt{time}_{ij} = 0 )
  • β3\beta_3 is a scale parameter that determines the rate at which the average brain weight approaches the asymptotic weight β1\beta_1 (decay rate)
weightij=β1+(β2β1)exp(β3timeij)+ϵij\texttt{weight}_{ij} = \beta_1 + \left(\beta_2 - \beta_1\right)\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}
\texttt{weight}_{ij} = \beta_1 + \left(\beta_2 - \beta_1\right)\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}
  • We assume the asymptotic weight β1\beta_1 is unicorn-specific (varies across unicorns) by including a random effect UjU_j as follows
β1=β1j=b1+Uj,\beta_1 = \beta_{1j} = b_1 + U_j,
\beta_1 = \beta_{1j} = b_1 + U_j,
  • Thus our model may be written as
weightij=b1+Uj+{β2(b1+Uj)}exp(β3timeij)+ϵij\texttt{weight}_{ij} = b_1 + U_j + \left\{\beta_2 - \left(b_1 + U_j\right)\right\}\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}
\texttt{weight}_{ij} = b_1 + U_j + \left\{\beta_2 - \left(b_1 + U_j\right)\right\}\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}
UjN(0,σu2)U_j \sim N\left(0, \sigma^2_u\right)
U_j \sim N\left(0, \sigma^2_u\right)
. menl weight = {b1}+{U[id]} + ({beta2} - ({b1}+{U[id]}) )*exp(-{beta3}*time)
  • Here is how you can fit it with menl

( Output is on next slide)

β1j\beta_{1j}
\beta_{1j}
β1j\beta_{1j}
\beta_{1j}

gsem-like syntax

. menl weight = {b1}+{U[id]} + ({beta2} - ({b1}+{U[id]}) )*exp(-{beta3}*time)

Obtaining starting values by EM:

Alternating PNLS/LME algorithm:

Iteration 1:    linearization log likelihood = -56.9757597

Computing standard errors:

Mixed-effects ML nonlinear regression           Number of obs     =        780
Group variable: id                              Number of groups  =         60

                                                Obs per group:
                                                              min =         13
                                                              avg =       13.0
                                                              max =         13
Linearization log likelihood =  -56.97576
------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /b1 |   4.707954   .1414511    33.28   0.000     4.430715    4.985193
      /beta2 |   8.089432   .0260845   310.12   0.000     8.038307    8.140556
      /beta3 |    4.13201   .0697547    59.24   0.000     3.995293    4.268726
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                      var(U) |   1.189809   .2180058      .8308328    1.703887
-----------------------------+------------------------------------------------
               var(Residual) |   .0439199   .0023148      .0396095    .0486995
------------------------------------------------------------------------------

weightij=b1+Uj+{β2(b1+Uj)}exp(β3timeij)+ϵij\texttt{weight}_{ij} = b_1 + U_j + \left\{\beta_2 - \left(b_1 + U_j\right)\right\}\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}
\texttt{weight}_{ij} = b_1 + U_j + \left\{\beta_2 - \left(b_1 + U_j\right)\right\}\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}
. menl weight = {b1}+{U[id]} + ({beta2} - ({b1}+{U[id]}) )*exp(-{beta3}*time)
  • We could have used the two-stage formulation of our model
weightij=β1+(β2β1)exp(β3timeij)+ϵij\texttt{weight}_{ij} = \beta_1 + \left(\beta_2 - \beta_1\right)\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}
\texttt{weight}_{ij} = \beta_1 + \left(\beta_2 - \beta_1\right)\exp\left(-\beta_3\texttt{time}_{ij}\right) + \epsilon_{ij}
β1=β1j=b1+Uj\beta_1 = \beta_{1j} = b_1 + U_j
\beta_1 = \beta_{1j} = b_1 + U_j
. menl weight = {beta1:} + ({beta2} - {beta1:})*exp({beta3}*time), ///
      define(beta1: {b1}+{U[id]})
  • Particularly useful for more complicated models
  • Parameters specified within define() may be easily predicted after estimation

Incorporating covariates

  • cupcake\texttt{cupcake} : number of rainbow cupcakes consumed right after birth
  • female\texttt{female} :     femalej\texttt{female}_j is 11 if the jjth unicorn is a female and 00 otherwise
     +---------------------------------------------+
     | id   cupcake       time   female     weight |
     |---------------------------------------------|
  1. |  1         2          0   female   8.266746 |
  2. |  1         2   .1666667   female   7.163768 |
  3. |  1         2   .3333333   female   6.725025 |
  4. |  1         2         .5   female    6.38325 |
  5. |  1         2   .6666667   female    6.00481 |
     |---------------------------------------------|
776. | 60         4   1.333333     male    3.75486 |
777. | 60         4        1.5     male   3.581552 |
778. | 60         4   1.666667     male   3.701847 |
779. | 60         4   1.833333     male   3.836078 |
780. | 60         4          2     male   3.729806 |
     +---------------------------------------------+

. list  if (_n <=5 | _n >775), sep(5)
  • Reducing brain weight loss has been an active research area in Zootopia, and scientists believe that consuming rainbow cupcakes right after birth may help slow down brain shrinkage.
β3j=b30+b31cupcakej\beta_{3j} = b_{30} + b_{31}\mathtt{cupcake}_j
\beta_{3j} = b_{30} + b_{31}\mathtt{cupcake}_j
  • Also, we would like to investigate whether the asymptote, β1j\beta_{1j}, is gender specific, so we include the factor variable female\mathtt{female} in the equation for β1j\beta_{1j}.
β1j=b10+b11femalej+U1j\beta_{1j} = b_{10} + b_{11}\mathtt{female}_j + U_{1j}
\beta_{1j} = b_{10} + b_{11}\mathtt{female}_j + U_{1j}
define(beta3: {b30} + {b31}*cupcake)
define(beta1: {b10} + {b11}*1.female + {U[id]})

decay rate

. menl weight = {beta1:} + ({beta2:} - {beta1:})*exp(-{beta3:}*time), ///
        define(beta1: {b10}+{b11}*1.female+{U0[id]})                  ///
        define(beta2: {b2})                                           ///
        define(beta3: {b30}+{b31}*cupcake)
  • Based on our results, consuming rainbow cupcakes after birth indeed slows down brain shrinkage: /b31\mathtt{/b31} is roughly 0.2-0.2 with a 95% CI\mathtt{CI} of [.27,.13][-.27, -.13]


Mixed-effects ML nonlinear regression           Number of obs     =        780
Group variable: id                              Number of groups  =         60

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

        beta1:  {b10}+{b11}*1.female+{U0[id]}
        beta2:  {b2}
        beta3:  {b30}+{b31}*cupcake

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        /b10 |   4.072752   .1627414    25.03   0.000     3.753785     4.39172
        /b11 |   1.264407   .2299723     5.50   0.000     .8136694    1.715144
         /b2 |   8.088102   .0255465   316.60   0.000     8.038032    8.138172
        /b30 |   4.706926   .1325714    35.50   0.000      4.44709    4.966761
        /b31 |  -.2007309   .0356814    -5.63   0.000    -.2706651   -.1307966
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                     var(U0) |   .7840578   .1438924      .5471839    1.123473
-----------------------------+------------------------------------------------
               var(Residual) |   .0420763   .0022176      .0379468    .0466551
------------------------------------------------------------------------------

Useful shortcut for your nonlinear expression

define(beta3: {b30} + {b31}*cupcake + {b32}*juice + {U2[id]})
define(beta3: cupcake juice U2[id])
define(beta3: {b30} + {b31}*cupcake)
define(beta3: cupcake, xb)
  • What if you have many covariates ?
  • Model parameters and random effects are bound in braces if specified directly in the expression: {b0},{U0[id]}, etc.\texttt{\{b0\}}, \texttt{\{U0[id]\}} \text{, etc}.
  • Linear combinations of variables can be included using the specification

 

  • Random effects can be specified within a linear combination, in which case they should be included without curly braces, for example, {lc_u: x1 x2 U[id]}\texttt{\{lc\_u: x1 x2 U[id]\}}

e.g., 

{lc: x1 x2, noconstant}\texttt{\{lc: x1 x2, noconstant\}}
\texttt{\{lc: x1 x2, noconstant\}}
{eqname:varlist [ , xb noconstant ]}\texttt{\{eqname:varlist [ , xb noconstant ]\}}
\texttt{\{eqname:varlist [ , xb noconstant ]\}}
{beta: mpg weight i.rep78}\texttt{\{beta: mpg weight i.rep78\}}
\texttt{\{beta: mpg weight i.rep78\}}

and

  • Random-effects names must start with a capital letter as in

    {U0[id]}\texttt{\{U0[id]\}}

Pharmacokinetics (PK) study of Theophylline

  • Pharmacokinetics is  the study of “what the body does with a drug
  • PK seeks to understand how drugs enter the body, distribute throughout the body, and leave the body, by summarizing concentration-time measurements, while accounting for patient-specific characteristics.

     

  • Maintain drug concentrations within the therapeutic window
  • Consider data on the antiasthmatic agent theophylline (Boeckmann, Sheiner, and Beal [1994] 2011).
  • The drug was administrated orally to 1212 subjects with the dosage (mg/kg) given on a per weight basis.
  • Serum concentrations (in mg/L) were obtained at 1111 time points per subject over 2525 hours following administration.
. webuse theoph

. twoway connected conc time, connect(L)

  • Clearance (Cl\mathtt{Cl}) measures the rate at which a drug is cleared from the plasma. It is defined as the volume of plasma cleared of drug per unit time.
  • Elimination rate constant (kek_e) describes the rate at which a drug is removed from the body. It is defined as the fraction of drug in the body eliminated per unit time.
  • Absorption rate constant (kak_a) describes the rate at which a drug is absorbed by the body.
concij=dosejkekajClj(kajke)×\texttt{conc}_{ij} = \frac{\texttt{dose}_j k_{e} k_{a_j}}{\mathtt{Cl}_j\left(k_{a_j} - k_{e}\right)}\times
\texttt{conc}_{ij} = \frac{\texttt{dose}_j k_{e} k_{a_j}}{\mathtt{Cl}_j\left(k_{a_j} - k_{e}\right)}\times
{exp(ketimeij)exp(kajtimeij)}+ϵij\left\{\exp\left(-k_{e}\texttt{time}_{ij}\right)-\exp\left(-k_{a_j}\texttt{time}_{ij}\right)\right\} + \epsilon_{ij}
\left\{\exp\left(-k_{e}\texttt{time}_{ij}\right)-\exp\left(-k_{a_j}\texttt{time}_{ij}\right)\right\} + \epsilon_{ij}
concij=dosejkekajClj(kajke)×\texttt{conc}_{ij} = \frac{\texttt{dose}_j k_{e} k_{a_j}}{\mathtt{Cl}_j\left(k_{a_j} - k_{e}\right)}\times
\texttt{conc}_{ij} = \frac{\texttt{dose}_j k_{e} k_{a_j}}{\mathtt{Cl}_j\left(k_{a_j} - k_{e}\right)}\times
{exp(ketimeij)exp(kajtimeij)}+ϵij\left\{\exp\left(-k_{e}\texttt{time}_{ij}\right)-\exp\left(-k_{a_j}\texttt{time}_{ij}\right)\right\} + \epsilon_{ij}
\left\{\exp\left(-k_{e}\texttt{time}_{ij}\right)-\exp\left(-k_{a_j}\texttt{time}_{ij}\right)\right\} + \epsilon_{ij}

Because each of the model parameters must be positive to be meaningful, we write

Clj=exp(β0+U0j)\mathtt{Cl}_j = \exp\left(\beta_0 + U_{0j}\right)
\mathtt{Cl}_j = \exp\left(\beta_0 + U_{0j}\right)
kaj=exp(β1+U1j)k_{a_j} = \exp\left(\beta_1 + U_{1j}\right)
k_{a_j} = \exp\left(\beta_1 + U_{1j}\right)
ke=exp(β2)k_{e} = \exp\left(\beta_2\right)
k_{e} = \exp\left(\beta_2\right)

(Output is on next slide)

. 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})) 
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 = -177.02233

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

------------------------------------------------------------------------------
        conc |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /b0 |  -3.227212   .0600089   -53.78   0.000    -3.344827   -3.109597
         /b1 |   .4656357   .1985941     2.34   0.019     .0763984     .854873
         /b2 |  -2.454679   .0524896   -46.77   0.000    -2.557556   -2.351801
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
subject: Independent         |
                     var(U0) |    .027864   .0121371      .0118652    .0654355
                     var(U1) |   .4143467   .1945672       .165067    1.040082
-----------------------------+------------------------------------------------
               var(Residual) |   .5030242   .0688883      .3846079    .6578995
------------------------------------------------------------------------------
covariance(U1 U2 U3, unstructured)
  • By default, menl assumes that all random effects are independent, thus, in the previous model, we have assumed that cov(U0,U1)=0\text{cov}\left(U_{0}, U_{1}\right) = 0, but we could have considered other covariance structures for the random effects
Σ=[σ12σ12σ13σ12σ22σ23σ13σ23σ32]\Sigma = \begin{bmatrix}\sigma^2_{1} & \sigma_{12} & \sigma_{13} \\ \sigma_{12} & \sigma^2_{2} & \sigma_{23} \\ \sigma_{13} & \sigma_{23} & \sigma^2_{3}\end{bmatrix}
\Sigma = \begin{bmatrix}\sigma^2_{1} & \sigma_{12} & \sigma_{13} \\ \sigma_{12} & \sigma^2_{2} & \sigma_{23} \\ \sigma_{13} & \sigma_{23} & \sigma^2_{3}\end{bmatrix}
U1U_1
U_1
U2U_2
U_2
U3U_3
U_3
U1U_1
U_1
U2U_2
U_2
U3U_3
U_3
  • For example, when there are 3 random effects U1U_1, U2U_2, and U3U_3 in the model, their covariance structure may be summarized as
cov(U1,U3)\text{cov}\left(U_{1}, U_{3}\right)
\text{cov}\left(U_{1}, U_{3}\right)
var(U1)\text{var}\left( U_{1}\right)
\text{var}\left( U_{1}\right)
Σ=[σ11σ12σ12σ11σ12σ11]\Sigma = \begin{bmatrix}\sigma_{11} & \sigma_{12} & \sigma_{12} \\ & \sigma_{11} & \sigma_{12} \\ & & \sigma_{11}\end{bmatrix}
\Sigma = \begin{bmatrix}\sigma_{11} & \sigma_{12} & \sigma_{12} \\ & \sigma_{11} & \sigma_{12} \\ & & \sigma_{11}\end{bmatrix}
covariance(U*, exchangeable)
covariance(U*, identity)
Σ=[σ1100σ110σ11]\Sigma = \begin{bmatrix}\sigma_{11} & 0 & 0 \\ & \sigma_{11} & 0 \\ & & \sigma_{11}\end{bmatrix}
\Sigma = \begin{bmatrix}\sigma_{11} & 0 & 0 \\ & \sigma_{11} & 0 \\ & & \sigma_{11}\end{bmatrix}
Σ=[σ11000σ22σ230σ23σ33]\Sigma = \begin{bmatrix}\sigma_{11} & 0 & 0 \\ 0 & \sigma_{22} & \sigma_{23} \\ 0 & \sigma_{23} & \sigma_{33}\end{bmatrix}
\Sigma = \begin{bmatrix}\sigma_{11} & 0 & 0 \\ 0 & \sigma_{22} & \sigma_{23} \\ 0 & \sigma_{23} & \sigma_{33}\end{bmatrix}
covariance(U2 U3, unstructured)
  • Other examples of random effects covariance structures are:
U1U_1
U_1
U2U_2
U_2
U3U_3
U_3
. 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})) covariance(U*, unstructured)
  • Going back to our PK model, we may wish to assume that U0jU_{0j} and U1jU_{1j} have a general covariance matrix
U1U_1
U_1
U2U_2
U_2
U3U_3
U_3
. 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}))                                  ///
    covariance(U*, unstructured)                            

Linearization log likelihood = -177.02222

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

--------------------------------------------------------------------------------
          conc |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
           /b0 |  -3.227224   .0600072   -53.78   0.000    -3.344836   -3.109612
           /b1 |   .4656428   .1985911     2.34   0.019     .0764113    .8548742
           /b2 |  -2.454696   .0524895   -46.77   0.000    -2.557574   -2.351819
--------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
subject: Unstructured        |
                     var(U0) |   .0278616   .0121451      .0118566    .0654714
                     var(U1) |   .4143347    .194565      .1650595    1.040069
                  cov(U0,U1) |  -.0002056   .0332183     -.0653122    .0649011
-----------------------------+------------------------------------------------
               var(Residual) |   .5030238   .0688885      .3846073    .6578995
------------------------------------------------------------------------------
  • Combine resvariance() and rescorrelation() to build flexible covariance structures.
  • Use option resvariance() to model heteroskedasticity as a linear, power, or exponential function of other covariates or of predicted values.
  • Use option rescorrelation() to model the dependence of the within-group errors as, e.g., AR or MA processes.
  • menl provides flexible modeling of within-group error or residual covariance structures.

So far, we have assumed that Var(ϵij)=σ2\text{Var}\left(\epsilon_{ij}\right) = \sigma^2 and cov(ϵsj,ϵtj)=0\text{cov}\left(\epsilon_{sj}, \epsilon_{tj}\right)=0

  • Use option rescorrelation() to model the dependence of the within-group errors as, e.g., auto-regressive (AR) or moving-average (MA) processes.

Variance increases with mean

Constant variance

Var(ϵij)=σ2(conc^ij)2δ\text{Var}\left(\epsilon_{ij}\right) = \sigma^2\left(\hat{\texttt{conc}}_{ij}\right)^{2\delta}
\text{Var}\left(\epsilon_{ij}\right) = \sigma^2\left(\hat{\texttt{conc}}_{ij}\right)^{2\delta}
Var(ϵij)=σe2\text{Var}\left(\epsilon_{ij}\right) = \sigma^2_e
\text{Var}\left(\epsilon_{ij}\right) = \sigma^2_e
Var(ϵij)=σ2{(conc^ij)δ+c}2\text{Var}\left(\epsilon_{ij}\right) = \sigma^2\left\{\left(\hat{\texttt{conc}}_{ij}\right)^{\delta} + c\right\}^2
\text{Var}\left(\epsilon_{ij}\right) = \sigma^2\left\{\left(\hat{\texttt{conc}}_{ij}\right)^{\delta} + c\right\}^2
resvariance(power _yhat)
  • For example, heteroskedasticity, often present in PK data, is modeled using the power function plus a constant.

 

Adding a constant avoids the variance of zero at timeij\texttt{time}_{ij} = 0, because the concentration is zero at that time.

 

  • We will also use restricted maximum likelihood estimation (REML) via option reml instead of the default maximum likelihood (ML) estimation to account for a moderate number of subjects.
. 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) reml


Mixed-effects REML nonlinear regression         

Linear. log restricted-likelihood = -172.44384

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

------------------------------------------------------------------------------
        conc |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /b0 |  -3.227295   .0619113   -52.13   0.000    -3.348639   -3.105951
         /b1 |   .4354519   .2072387     2.10   0.036     .0292716    .8416322
         /b2 |  -2.453743   .0517991   -47.37   0.000    -2.555267   -2.352218
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
subject: Independent         |
                     var(U0) |   .0316416    .014531      .0128634    .0778326
                     var(U1) |   .4500585   .2228206      .1705476    1.187661
-----------------------------+------------------------------------------------
Residual variance:           |
  Power _yhat                |
                      sigma2 |   .1015759    .086535      .0191263    .5394491
                       delta |   .3106636   .2466547     -.1727707    .7940979
                       _cons |   .7150935   .3745256      .2561837    1.996063
------------------------------------------------------------------------------
. predict (Cl = {cl:})

. list subject Cl if subject[_n]!=subject[_n-1], sep(0) noobs

  • We can predict the clearance parameter Clj\mathtt{Cl}_j for each subject
  • Interpretation: For subject 1, 0.0270.027 liters of serum concentration are cleared of the theophylline drug per hour per kg body weight.
  • In other words,  if subject 1 weighs 75 kg, 75×0.027275\times0.027 \approx 2 liters of serum concentration are cleared of theophylline every hour
  +--------------------+
  | subject         Cl |
  |--------------------|
  |       1   .0274995 |
  |       2   .0402322 |
  |       3   .0406919 |
  |       4   .0371865 |
  |       5   .0425807 |
  |       6   .0461769 |
  |       7   .0466375 |
  |       8   .0442402 |
  |       9   .0324682 |
  |      10   .0356245 |
  |      11   .0510018 |
  |      12     .03785 |
  +--------------------+

Summary II

yij=func(covariates,βj)+ϵijy_{ij} = \mathtt{func}\left(\texttt{covariates}, {\beta}_j\right) + \epsilon_{ij}
y_{ij} = \mathtt{func}\left(\texttt{covariates}, {\beta}_j\right) + \epsilon_{ij}

Stage I:

βj=func(covariates,b,Uj)\beta_{j} = \mathtt{func}\left(\texttt{covariates}, \mathbf b, \mathbf U_j\right)
\beta_{j} = \mathtt{func}\left(\texttt{covariates}, \mathbf b, \mathbf U_j\right)

Stage II:

menl y = <menlexpr>, define() [define() ...]
covariance(U0 U1 [U2 ...], <covtype>)
resvariance() rescorrelation()

with

Uj=(U0j,U1j,,Uqj)MVN(0,Σ)\mathbf U_j = \left(U_{0j},U_{1j},\ldots,U_{qj}\right) \sim \mathtt{MVN}\left(\mathbf 0, \Sigma\right)
\mathbf U_j = \left(U_{0j},U_{1j},\ldots,U_{qj}\right) \sim \mathtt{MVN}\left(\mathbf 0, \Sigma\right)
(ϵ1j,,ϵnjj)MVN(0,σ2Λj)\left(\epsilon_{1j}, \ldots ,\epsilon_{n_jj}\right) \sim \mathtt{MVN}\left(\mathbf 0, \sigma^2 \Lambda_j\right)
\left(\epsilon_{1j}, \ldots ,\epsilon_{n_jj}\right) \sim \mathtt{MVN}\left(\mathbf 0, \sigma^2 \Lambda_j\right)

and

2-stage specification

Single-stage specification

yij=func(covariates,b,Uj)y_{ij} = \mathtt{func}\left(\texttt{covariates}, \mathbf b, \mathbf U_j\right)
y_{ij} = \mathtt{func}\left(\texttt{covariates}, \mathbf b, \mathbf U_j\right)
+ϵij+ \epsilon_{ij}
+ \epsilon_{ij}
menl y = <menlexpr>
rescovariance()

What about multilevel NLME models ?

  • Consider fictional data on tree circumference of 30 trees planted in 5 different planting zones.  The tree circumferences were recorded over the period of 1990 to 2003.

Zone 1

Zone 5

Tree 27

Tree 6

Tree 30

Tree 1

Tree 2

19901990

19911991

20032003

Zone 2

Level 3

Level 2

Level 1

  • We wish to fit the 3-parameter logistic growth model:
  • circumijt\texttt{circum}_{ijt} and ageijt\texttt{age}_{ijt} are circumference (cm) and age (in years) of tree jj within planting zone ii in year tt
  • Parameters: Asymptotic circumference (β1\beta_1), time at which half of β1\beta_1 is attained (β2\beta_2), and scale parameter (β3\beta_3)
  • For simplicity we assume that only β2\beta_2 may be affected by zones and trees-within-zones:
β2=β2ij=b2+U1i+U2ij\beta_2 = \beta_{2ij} = b_2 + U_{1i} + U_{2ij}
\beta_2 = \beta_{2ij} = b_2 + U_{1i} + U_{2ij}
circumijt=β11+exp{(ageijtβ2)/β3}+ϵijt\texttt{circum}_{ijt} = \frac{\beta_1}{1+\exp{\left\{-\left(\texttt{age}_{ijt} - \beta_2\right)/\beta_3\right\}}} + \epsilon_{ijt}
\texttt{circum}_{ijt} = \frac{\beta_1}{1+\exp{\left\{-\left(\texttt{age}_{ijt} - \beta_2\right)/\beta_3\right\}}} + \epsilon_{ijt}
. menl circum = {b1}/(1 + exp(-(age - {beta2:})/{b3})),  ///
        define(beta2: {b2} + {U1[zone]} + {U2[zone > tree]}) stddeviations




circumijt=β11+exp{(ageijtβ2)/β3}+ϵijt\texttt{circum}_{ijt} = \frac{\beta_1}{1+\exp{\left\{-\left(\texttt{age}_{ijt} - \beta_2\right)/\beta_3\right\}}} + \epsilon_{ijt}
\texttt{circum}_{ijt} = \frac{\beta_1}{1+\exp{\left\{-\left(\texttt{age}_{ijt} - \beta_2\right)/\beta_3\right\}}} + \epsilon_{ijt}
β2=β2ij=b2+U1i+U2ij\beta_2 = \beta_{2ij} = b_2 + U_{1i} + U_{2ij}
\beta_2 = \beta_{2ij} = b_2 + U_{1i} + U_{2ij}
Mixed-effects ML nonlinear regression           Number of obs     =        420

-------------------------------------------------------------
                |     No. of       Observations per Group
        Path    |     Groups    Minimum    Average    Maximum
----------------+--------------------------------------------
           zone |          5         56       84.0        112
      zone>tree |         30         14       14.0         14
-------------------------------------------------------------

Linearization log likelihood =  313.18701

       beta2:  {b2}+{U1[zone]}+{U2[zone>tree]}

------------------------------------------------------------------------------
     circumf |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /b1 |   30.01228   .0228861  1311.38   0.000     29.96742    30.05713
         /b2 |   6.459742   .0757833    85.24   0.000     6.311209    6.608274
         /b3 |   2.902458   .0046722   621.21   0.000     2.893301    2.911615
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
zone: Identity               |
                      sd(U1) |   .1599475    .057271      .0792858    .3226707
-----------------------------+------------------------------------------------
zone>tree: Identity          |
                      sd(U2) |   .1244983   .0178957      .0939312    .1650126
-----------------------------+------------------------------------------------
                sd(Residual) |   .0966751   .0034615      .0901233    .1037034
------------------------------------------------------------------------------
  • The time to reach half of the asymptotic circumference  (β2\beta_2) for jjth tree within planting zone ii:
. predict (beta2_i = {beta2:}), relevel(zone)

. list zone beta2_i if zone[_n]!=zone[_n-1], sep(0) noobs

  +----------------------------+
  |            zone    beta2_i |
  |----------------------------|
  |     New England   6.357335 |
  |    Mid Atlantic   6.273229 |
  | E North Central   6.723731 |
  | W North Central   6.481954 |
  |  South Atlantic    6.46246 |
  +----------------------------+

β^2=β^2ij=b^2+U1i+U2ij\hat{\beta}_2 = \hat{\beta}_{2ij} = \hat{b}_2 + U_{1i} + U_{2ij}
\hat{\beta}_2 = \hat{\beta}_{2ij} = \hat{b}_2 + U_{1i} + U_{2ij}
β^2=β^2i=b^2+U1i+0=b^2+U1i\hat{\beta}_2 = \hat{\beta}_{2i} = \hat{b}_2 + U_{1i} + 0 = \hat{b}_2 + U_{1i}
\hat{\beta}_2 = \hat{\beta}_{2i} = \hat{b}_2 + U_{1i} + 0 = \hat{b}_2 + U_{1i}
  • We can predict the time to reach half of the asymptotic circumference  (β2\beta_2) for each planting zone ii:

Summary

  • menl supports ML and REML estimation and provides flexible random-effects and residual covariance structures.

 

  • menl supports single-stage and two-stage specifications.

     

  • You can predict random effects and their standard errors, group-specific nonlinear parameters, and more after estimation; see [ME] menl postestimation.

 

  • NLME models are known to be sensitive to initial values. menl   provides default, but for some models you may need to specify your own. Use option initial().

 

  • NLME models are known to have difficulty converging or converging to a local maximum. Trying different initial values may help.

 

  • It also supports specifying models with no random effects (nonlinear marginal models)

ALWAYS BE YOURSELF

UNLESS YOU CAN BE A UNICORN

THEN ALWAYS BE A UNICORN

 

 

Thank You!

Stata Conference (Colombus, ohio)

By H Assaad

Stata Conference (Colombus, ohio)

  • 3,840