Appendix

NONMEM 7 Technical Guide

Fifteen appendices (A-O) covering matrix algebra, Fischer score and exact second-derivative matrices, Cholesky parameterization of Sigma, OMEGA prior degrees of freedom, nonparametric analysis, TNPRI, t-distribution sampling, parameter transformations, super-ID etas, distributed-delay output, and the ddexpand utility.

Appendix A: Matrix Algebra Tools

We wish to determine the derivative of $-\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega}))$ with respect to $\boldsymbol{\mu}$ and $\boldsymbol{\Omega}^{-1}$.

Differentiating with respect to $\boldsymbol{\mu}$ is easily done as follows:

\[ \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\mu} } = -\frac{1}{2} \frac{ \partial \log(|\boldsymbol{\Omega}^{-1}|) }{ \partial \boldsymbol{\mu} } + \frac{1}{2} \frac{ \partial (\boldsymbol{\phi} - \boldsymbol{\mu})' \boldsymbol{\Omega}^{-1} (\boldsymbol{\phi} - \boldsymbol{\mu}) }{ \partial \boldsymbol{\mu} } \tag{1}\]

For a particular parameter,

\[ \begin{aligned} \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \mu_j } &= \frac{1}{2} \frac{ \partial (\boldsymbol{\phi} - \boldsymbol{\mu})' \boldsymbol{\Omega}^{-1}(\boldsymbol{\phi} - \boldsymbol{\mu}) }{ \partial \mu_j } \\ &= -\frac{1}{2} \mathbf{i}'_j \boldsymbol{\Omega}^{-1}(\boldsymbol{\phi} - \boldsymbol{\mu}) - \frac{1}{2} (\boldsymbol{\phi} - \boldsymbol{\mu})' \boldsymbol{\Omega}^{-1} \mathbf{i}_j \end{aligned} \tag{2}\]

where $\mathbf{i}_j$ is a vector of 0’s except for the $j$th element, which is 1. But the scalar terms are equal:

\[ \left((\boldsymbol{\phi} - \boldsymbol{\mu})' \boldsymbol{\Omega}^{-1} \mathbf{i}_j\right)' = \mathbf{i}'_j \boldsymbol{\Omega}^{-1}(\boldsymbol{\phi} - \boldsymbol{\mu}) \tag{3}\]

\[ \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \mu_j } = -\mathbf{i}'_j \boldsymbol{\Omega}^{-1}(\boldsymbol{\phi} - \boldsymbol{\mu}) \tag{4}\]

and for the entire vector $\boldsymbol{\mu}$,

\[ \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\mu} } = -\boldsymbol{\Omega}^{-1}(\boldsymbol{\phi} - \boldsymbol{\mu}) \tag{5}\]

To differentiate $-\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega}))$ with respect to $\boldsymbol{\Omega}^{-1}$ is more difficult:

\[ \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\Omega}^{-1} } = -\frac{1}{2} \frac{ \partial \log(|\boldsymbol{\Omega}^{-1}|) }{ \partial \boldsymbol{\Omega}^{-1} } + \frac{1}{2} \frac{ \partial (\boldsymbol{\phi} - \boldsymbol{\mu})' \boldsymbol{\Omega}^{-1}(\boldsymbol{\phi} - \boldsymbol{\mu}) }{ \partial \boldsymbol{\Omega}^{-1} } \tag{6}\]

To do so, we must develop some partial derivative relationships in linear algebra. Consider any non-singular square matrix $\mathbf{Z}$ which is related to its inverse by

\[ \mathbf{Z}\mathbf{Z}^{-1} = \mathbf{I} \tag{7}\]

Therefore, the partial derivative with respect to some variable $x$ yields

\[ \begin{aligned} \frac{\partial (\mathbf{Z}\mathbf{Z}^{-1})}{\partial x} &= \frac{\partial \mathbf{Z}}{\partial x} \mathbf{Z}^{-1} + \mathbf{Z} \frac{\partial \mathbf{Z}^{-1}}{\partial x} \\ &= \frac{\partial (\mathbf{I})}{\partial x} \\ &= \mathbf{0} \end{aligned} \tag{8}\]

It follows that

\[ \frac{\partial \mathbf{Z}^{-1}}{\partial x} = -\mathbf{Z}^{-1} \frac{\partial \mathbf{Z}}{\partial x} \mathbf{Z}^{-1} \tag{9}\]

Suppose

\[ x = z_{jk} \tag{10}\]

then

\[ \frac{\partial \mathbf{Z}^{-1}}{\partial z_{jk}} = -\mathbf{Z}^{-1} \frac{\partial \mathbf{Z}}{\partial z_{jk}} \mathbf{Z}^{-1} = -\mathbf{Z}^{-1} \mathbf{I}_{jk} \mathbf{Z}^{-1} = \left\{-z^{-1}_{mj} z^{-1}_{kp},\ \text{for all}\ m = 1\ \text{to}\ n,\ p = 1\ \text{to}\ n\right\} \tag{11}\]

where $\mathbf{I}_{jk}$ is a matrix that is 0 everywhere except for element $j,k$ which has a value of 1. Or

\[ \frac{\partial z^{-1}_{mp}}{\partial z_{jk}} = -z^{-1}_{mj} z^{-1}_{kp} \tag{12}\]

Similarly,

\[ \frac{\partial \mathbf{Z}^{-1}}{\partial z_{jk}} = \frac{\partial \mathbf{Z}^{-1}}{\partial z_{kj}} \frac{\partial z_{kj}}{\partial z_{jk}} = -\mathbf{Z}^{-1} \frac{\partial \mathbf{Z}}{\partial z_{kj}} \mathbf{Z}^{-1} = -\mathbf{Z}^{-1} \mathbf{I}_{kj} \mathbf{Z}^{-1} = \left\{-z^{-1}_{jm} z^{-1}_{pk},\ \text{for all}\ m = 1\ \text{to}\ n,\ p = 1\ \text{to}\ n\right\} \tag{13}\]

Suppose it is more convenient to differentiate a particular function with respect to the inverse of $\mathbf{Z}$, then

\[ \begin{aligned} \frac{\partial f(\mathbf{Z})}{\partial z_{jk}} &= \sum_{m=1}^{n}\sum_{p=1}^{n} \frac{\partial f(\mathbf{Z})}{\partial z^{-1}_{mp}} \frac{\partial z^{-1}_{mp}}{\partial z_{jk}} \\ &= \sum_{m=1}^{n}\sum_{p=1}^{n} -z^{-1}_{mj} \frac{\partial f(\mathbf{Z})}{\partial z^{-1}_{mp}} z^{-1}_{kp} \\ &= \sum_{m=1}^{n}\sum_{p=1}^{n} -z^{-1}_{jm} \frac{\partial f(\mathbf{Z})}{\partial z^{-1}_{mp}} z^{-1}_{pk} \end{aligned} \tag{14}\]

or in matrix notation,

\[ \frac{\partial f(\mathbf{Z})}{\partial \mathbf{Z}} = -\mathbf{Z}^{-1} \frac{\partial f(\mathbf{Z})}{\partial \mathbf{Z}^{-1}} \mathbf{Z}^{-1} \tag{15}\]

Furthermore, according to linear algebra, for any matrix $\mathbf{Z}$

\[ |\mathbf{Z}| = \sum_{j=1}^{n} z_{jk} Z_{jk}\ \text{for any}\ k = 1\ \text{to}\ n \tag{16}\]

where $z_{jk}$ is an element of matrix $\mathbf{Z}$, and $Z_{jk}$ is its cofactor. The cofactor $Z_{jk}$ is the determinant of the sub-matrix of $\mathbf{Z}$ that does not include row $j$ and column $k$. Therefore, $Z_{jk}$ does not contain the element $z_{jk}$. Again, according to linear algebra,

\[ z^{-1}_{jk} = \frac{Z_{kj}}{|\mathbf{Z}|} \tag{17}\]

where $z^{-1}_{jk}$ is the $j,k$ element of $\mathbf{Z}^{-1}$. Therefore

\[ \begin{aligned} \frac{\partial |\mathbf{Z}|}{\partial z_{lk}} &= \sum_{j=1}^{n} \frac{\partial z_{jk}}{\partial z_{lk}} Z_{jk} \\ &= Z_{lk} \\ &= |\mathbf{Z}| z^{-1}_{kl} \end{aligned} \tag{18}\]

\[ \frac{\partial |\mathbf{Z}|}{\partial \mathbf{Z}} = |\mathbf{Z}| \mathbf{Z}'^{-1} \tag{19}\]

More generally, for $\mathbf{Z}$ raised to any power $p$,

\[ \begin{aligned} \frac{\partial |\mathbf{Z}|^p}{\partial \mathbf{Z}} &= \frac{\partial |\mathbf{Z}|^p}{\partial |\mathbf{Z}|} \frac{\partial |\mathbf{Z}|}{\partial \mathbf{Z}} \\ &= p |\mathbf{Z}|^{p-1} \frac{\partial |\mathbf{Z}|}{\partial \mathbf{Z}} \\ &= p |\mathbf{Z}|^p \mathbf{Z}'^{-1} \end{aligned} \tag{20}\]

It follows that for any variable $x$ which influences the elements of the matrix $\mathbf{Z}$,

\[ \begin{aligned} \frac{\partial |\mathbf{Z}|^p}{\partial x} &= \sum_{j=1}^{n}\sum_{k=1}^{n} \frac{\partial |\mathbf{Z}|^p}{\partial z_{jk}} \frac{\partial z_{jk}}{\partial x} \\ &= p |\mathbf{Z}|^p \sum_{j=1}^{n}\sum_{k=1}^{n} z^{-1}_{kj} \frac{\partial z_{jk}}{\partial x} \\ &= p |\mathbf{Z}|^p \mathrm{tr}\left(\mathbf{Z}^{-1} \frac{\partial \mathbf{Z}}{\partial x}\right) \end{aligned} \tag{21}\]

So also,

\[ \begin{aligned} \frac{\partial \log(|\mathbf{Z}|)}{\partial z_{lk}} &= \frac{1}{|\mathbf{Z}|} \frac{\partial |\mathbf{Z}|}{\partial z_{lk}} \\ &= z^{-1}_{kl} \end{aligned} \tag{22}\]

so that

\[ \frac{\partial \log(|\mathbf{Z}|)}{\partial \mathbf{Z}} = \mathbf{Z}'^{-1} \tag{23}\]

or, for any variable $x$,

\[ \begin{aligned} \frac{\partial \log(|\mathbf{Z}|)}{\partial x} &= \sum_{j=1}^{n}\sum_{k=1}^{n} \frac{\partial \log(|\mathbf{Z}|)}{\partial z_{jk}} \frac{\partial z_{jk}}{\partial x} \\ &= \sum_{j=1}^{n}\sum_{k=1}^{n} z^{-1}_{kj} \frac{\partial z_{jk}}{\partial x} \\ &= \mathrm{tr}\left(\mathbf{Z}^{-1} \frac{\partial \mathbf{Z}}{\partial x}\right) \end{aligned} \tag{24}\]

For any $n \times m$ matrix $\mathbf{Y}$, any $n \times m$ matrix $\mathbf{X}$, and any $n \times n$ matrix $\mathbf{Z}$,

\[ \begin{aligned} \mathrm{tr}(\mathbf{X}'\mathbf{Z}\mathbf{Y}) &= \sum_{i=1}^{m}\sum_{j=1}^{n}\sum_{k=1}^{n} x'_{ij} z_{jk} y_{ki} \\ &= \sum_{j=1}^{n}\sum_{k=1}^{n}\sum_{i=1}^{m} z_{jk} y_{ki} x'_{ij} \\ &= \mathrm{tr}(\mathbf{Z}\mathbf{Y}\mathbf{X}') \\ &= \sum_{k=1}^{n}\sum_{i=1}^{m}\sum_{j=1}^{n} y_{ki} x'_{ij} z_{jk} \\ &= \mathrm{tr}(\mathbf{Y}\mathbf{X}'\mathbf{Z}) \end{aligned} \tag{25}\]

If $m = 1$, then

\[ \begin{aligned} \mathrm{tr}(\mathbf{x}'\mathbf{Z}\mathbf{y}) &= \mathbf{x}'\mathbf{Z}\mathbf{y} \\ &= \mathrm{tr}(\mathbf{Z}\mathbf{y}\mathbf{x}') \\ &= \mathrm{tr}(\mathbf{y}\mathbf{x}'\mathbf{Z}) \end{aligned} \tag{26}\]

Derivatives to trace functions can be derived as follows:

\[ \begin{aligned} \frac{ \partial \mathrm{tr}(\mathbf{X}'\mathbf{Z}\mathbf{Y}) }{ \partial x_{mp} } &= \mathrm{tr}(\mathbf{I}'_{pm}\mathbf{Z}\mathbf{Y}) \\ &= \mathrm{tr}(\mathbf{I}'_p \mathbf{I}_m \mathbf{Z}\mathbf{Y}) \\ &= \mathrm{tr}(\mathbf{I}_m \mathbf{Z}\mathbf{Y}\mathbf{I}'_p) \end{aligned} \tag{27}\]

Where $\mathbf{I}_j$ is a column vector of zeros for all elements except for element $j$, for which it has a value of 1. Or,

\[ \frac{ \partial \mathrm{tr}(\mathbf{X}'\mathbf{Z}\mathbf{Y}) }{ \partial \mathbf{X} } = \mathbf{Z}\mathbf{Y} \tag{28}\]

Similarly (taking some shortcuts in element/matrix nomenclature),

\[ \begin{aligned} \frac{ \partial \mathrm{tr}(\mathbf{X}'\mathbf{Z}\mathbf{Y}) }{ \partial \mathbf{Z} } &= \mathrm{tr}(\mathbf{X}'\mathbf{I}_{mp}\mathbf{Y}) \\ &= \mathrm{tr}(\mathbf{X}'\mathbf{I}_m \mathbf{I}'_p \mathbf{Y}) \\ &= \mathrm{tr}(\mathbf{I}'_p \mathbf{Y}\mathbf{X}'\mathbf{I}_m) \\ &= \mathbf{X}\mathbf{Y}' \end{aligned} \tag{29}\]

and

\[ \frac{ \partial \mathrm{tr}(\mathbf{X}'\mathbf{Z}\mathbf{Y}) }{ \partial \mathbf{Y} } = \mathbf{X}\mathbf{Z}' \tag{30}\]

Using the above relationships for the following,

\[ \frac{ \partial (\boldsymbol{\phi} - \boldsymbol{\mu})' \boldsymbol{\Omega}^{-1}(\boldsymbol{\phi} - \boldsymbol{\mu}) }{ \partial \boldsymbol{\Omega}^{-1} } = (\boldsymbol{\phi} - \boldsymbol{\mu})(\boldsymbol{\phi} - \boldsymbol{\mu})' \tag{31}\]

It follows that

\[ \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\Omega}^{-1} } = -\frac{1}{2} \boldsymbol{\Omega} + \frac{1}{2} (\boldsymbol{\phi} - \boldsymbol{\mu})(\boldsymbol{\phi} - \boldsymbol{\mu})' \tag{32}\]

and incidentally,

\[ \begin{aligned} \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\Omega} } &= -\boldsymbol{\Omega}^{-1}\left(\frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\Omega}^{-1} }\right) \boldsymbol{\Omega}^{-1} \\ &= \frac{1}{2} \boldsymbol{\Omega}^{-1}(\boldsymbol{\Omega} - (\boldsymbol{\phi} - \boldsymbol{\mu})(\boldsymbol{\phi} - \boldsymbol{\mu})') \boldsymbol{\Omega}^{-1} \end{aligned} \tag{33}\]

The above tools also allow us to evaluate the following

\[ \frac{\partial L_{Ni}}{\partial \boldsymbol{\Omega}^{-1}} = -\frac{1}{2} \frac{ \partial \log(|\boldsymbol{\Omega}^{-1}|) }{ \partial \boldsymbol{\Omega}^{-1} } + \frac{1}{2} \frac{ \partial (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu})' \boldsymbol{\Omega}^{-1}(\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu}) }{ \partial \boldsymbol{\Omega}^{-1} } + \frac{1}{2} \frac{ \partial \log(|\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1}|) }{ \partial \boldsymbol{\Omega}^{-1} } \tag{34}\]

Since

\[ \begin{aligned} \frac{ \partial \log(|\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1}|) }{ \partial \omega^{-1}_{jk} } &= \frac{ 1 }{ \det(\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1}) } \frac{ \partial |\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1}| }{ \partial \omega^{-1}_{jk} } \\ &= \frac{ 1 }{ |\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1}| } \frac{ \partial |\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1}| }{ \partial (\omega^{-1}_{jk} + \hat{s}^{-1}_{jki}) } \frac{ \partial (\omega^{-1}_{jk} + \hat{s}^{-1}_{jki}) }{ \partial \omega^{-1}_{jk} } \\ &= \left(\omega^{-1}_{jk} + \hat{s}^{-1}_{jki}\right)^{-1} \end{aligned} \tag{35}\]

\[ \frac{ \partial \log(|\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1}|) }{ \partial \boldsymbol{\Omega}^{-1} } = (\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})^{-1} \tag{36}\]

It follows that

\[ \frac{\partial L_{Ni}}{\partial \boldsymbol{\Omega}^{-1}} = -\frac{1}{2} \boldsymbol{\Omega} + \frac{1}{2} (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu})(\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu})' + \frac{1}{2} (\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})^{-1} \tag{37}\]

and

\[ \frac{\partial L_{Ni}}{\partial \boldsymbol{\Omega}} = -\boldsymbol{\Omega}^{-1}\left(\frac{\partial L_{Ni}}{\partial \boldsymbol{\Omega}^{-1}}\right) \boldsymbol{\Omega}^{-1} \tag{38}\]

\[ = \frac{1}{2} \boldsymbol{\Omega}^{-1}\left[\boldsymbol{\Omega} - (\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu})(\hat{\boldsymbol{\phi}}_i - \boldsymbol{\mu})' - (\boldsymbol{\Omega}^{-1} + \hat{\mathbf{S}}_i^{-1})^{-1}\right] \boldsymbol{\Omega}^{-1} \tag{39}\]

Because $\boldsymbol{\Omega}$ is symmetrical, the independent parameters which must be varied to minimize the objective function consist of only half of the matrix. Let us define the lower triangular matrix $\mathbf{A}$ containing independent parameters which relate to the elements of $\boldsymbol{\Omega}$ such that

\[ \begin{aligned} \omega_{jk} &= a_{jk} \\ \omega_{kj} &= a_{jk},\ \text{for}\ j = 1\ \text{to}\ n,\ k = 1\ \text{to}\ j \end{aligned} \tag{40}\]

\[ \boldsymbol{\Omega} = \mathbf{A} + \mathbf{A}' - \mathrm{diag}(\mathbf{A}) \tag{41}\]

It follows that

\[ \begin{aligned} \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial a_{jk} } & = \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \omega_{jk} } + \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \omega_{kj} } \\ &\quad - \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \omega_{jk} } \delta(j, k) \\ & = 2 \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \omega_{jk} } - \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \omega_{jk} } \delta(j, k) \end{aligned} \tag{42}\]

or in matrix notation,

\[ \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \mathbf{A} } = \mathrm{Lower}\left[2 \frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\Omega} } - \mathrm{diag}\left(\frac{ \partial -\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\Omega} }\right)\right] \tag{43}\]

and equation 30 results. At the minimum,

\[ \begin{aligned} \frac{\partial L}{\partial \mathbf{A}} &= \frac{\partial L}{\partial \boldsymbol{\Omega}} \\ &= \mathbf{0} \end{aligned} \tag{44}\]

Another trace relationship that appears in probability densities is:

\[ \mathrm{tr}(\boldsymbol{\Sigma}^{-1}\mathbf{R}) = \mathrm{tr}(\boldsymbol{\Lambda}'^{-1}\boldsymbol{\Lambda}^{-1}\mathbf{R}) \tag{45}\]

Where $\boldsymbol{\Lambda}$ is the lower triangular cholesky matrix to a symmetrical matrix $\boldsymbol{\Sigma}$, and $\mathbf{R}$ is also a symmetrical matrix. Derivatives with respect to the cholesky elements is often desired, so,

\[ \begin{aligned} \frac{ \partial \mathrm{tr}(\boldsymbol{\Lambda}'^{-1}\boldsymbol{\Lambda}^{-1}\mathbf{R}) }{ \partial \boldsymbol{\Lambda} } & = -\boldsymbol{\Lambda}'^{-1} \frac{ \partial \mathrm{tr}(\boldsymbol{\Lambda}'^{-1}\boldsymbol{\Lambda}^{-1}\mathbf{R}) }{ \partial \boldsymbol{\Lambda}^{-1} } \boldsymbol{\Lambda}'^{-1} = -\boldsymbol{\Lambda}'^{-1}\left(\boldsymbol{\Lambda}^{-1}\mathbf{R} + \boldsymbol{\Lambda}^{-1}\mathbf{R}'\right) \boldsymbol{\Lambda}'^{-1} \\ & = -2 \boldsymbol{\Lambda}'^{-1} \boldsymbol{\Lambda}^{-1} \mathbf{R} \boldsymbol{\Lambda}'^{-1} = -2 \boldsymbol{\Sigma}^{-1} \mathbf{R} \boldsymbol{\Lambda}'^{-1} \end{aligned} \tag{46}\]

Appendix B: Positive Definite Properties

A matrix $\mathbf{A}$ is defined as positive definite if

\[ \mathbf{x}'\mathbf{A}\mathbf{x} > 0 \tag{47}\]

for any vector $\mathbf{x} \neq \mathbf{0}$. Consider a matrix constructed as follows

\[ \mathbf{A} = \mathbf{Y}\mathbf{Z}\mathbf{Y}' \tag{48}\]

where $\mathbf{Y}$ is any non-zero $n \times m$ matrix, $\mathbf{Z}$ is an $m \times m$ positive definite matrix, so that $\mathbf{A}$ is an $n \times n$ matrix. Then, with any non-zero $n \times 1$ vector $\mathbf{x}$,

\[ \begin{aligned} \mathbf{x}'\mathbf{A}\mathbf{x} &= \mathbf{x}'\mathbf{Y}\mathbf{Z}\mathbf{Y}'\mathbf{x} \\ &= \mathbf{v}'\mathbf{Z}\mathbf{v} \end{aligned} \tag{49}\]

where

\[ \mathbf{v} = \mathbf{Y}'\mathbf{x} \tag{50}\]

is a non-zero $m \times 1$ vector. It follows that

\[ \mathbf{x}'\mathbf{A}\mathbf{x} = \mathbf{v}'\mathbf{Z}\mathbf{v} > 0 \tag{51}\]

and $\mathbf{A}$ is therefore positive definite. The sum of two positive definite matrices $\mathbf{A}$ and $\mathbf{B}$ is also positive definite. Let

\[ \begin{aligned} \mathbf{C} &= \mathbf{A} + \mathbf{B} \\ \mathbf{x}'\mathbf{A}\mathbf{x} &= a > 0 \\ \mathbf{x}'\mathbf{B}\mathbf{x} &= b > 0 \end{aligned} \tag{52}\]

then

\[ \begin{aligned} \mathbf{x}'\mathbf{A}\mathbf{x} + \mathbf{x}'\mathbf{B}\mathbf{x} &= \mathbf{x}'(\mathbf{A} + \mathbf{B})\mathbf{x} \\ &= \mathbf{x}'\mathbf{C}\mathbf{x} \\ &= a + b > 0 \end{aligned} \tag{53}\]

Thus, matrices of the form

\[ \sum_{i=1}^{n} \mathbf{Y}_i \mathbf{Z}_i \mathbf{Y}'_i \tag{54}\]

where $\mathbf{Z}_i$ is positive definite, are also positive definite.

The second derivative of any objective function $L$ evaluated at its minimum is positive definite. This can be shown by considering that the derivative of $L$ evaluated at its minimum $L_0$ is equal to 0 (otherwise, it would not be at a minimum/maximum):

\[ \left(\frac{\partial L}{\partial \boldsymbol{\theta}}\right)_0 = \mathbf{0} \tag{55}\]

And any small perturbation $\Delta\boldsymbol{\theta}$ from the minimum results in a change in $L$, called $\Delta L$, that is positive, otherwise it would be at its maximum, not its minimum. By Taylor series expansion we have:

\[ \begin{aligned} 0 < \Delta L &= \left(\frac{\partial L}{\partial \boldsymbol{\theta}}\right)_0 \Delta\boldsymbol{\theta} + \frac{1}{2} \Delta\boldsymbol{\theta}'\left(\frac{\partial^2 L}{\partial \boldsymbol{\theta}^2}\right)_0 \Delta\boldsymbol{\theta} \\ &= \frac{1}{2} \Delta\boldsymbol{\theta}'\left(\frac{\partial^2 L}{\partial \boldsymbol{\theta}^2}\right)_0 \Delta\boldsymbol{\theta} \end{aligned} \tag{56}\]

and therefore the second derivative of the objective function evaluated at its minimum is positive definite.

Appendix C: The Fischer Score Matrix for Error Assessment in EM Problems

We wish to evaluate the expected value, over all possible data $\mathbf{y}$ and over an infinite number of subjects $m$, of the second derivative of the objective function with respect to the population parameters $\mathbf{q} = \boldsymbol{\mu},\boldsymbol{\Omega}$. The inverse of that matrix is then the asymptotic error matrix to the parameters. We do this as follows. Noting that

\[ p(\mathbf{y}_i \mid \mathbf{q}) = \int p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q}) d\boldsymbol{\phi} \tag{57}\]

where

\[ p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q}) = l(\mathbf{y}_i \mid \boldsymbol{\phi}) h(\boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega}) \tag{58}\]

and

\[ \begin{aligned} L &= \sum_{i=1}^{m} -\log(p_i(\mathbf{y}_i \mid \mathbf{q})) \\ &= -\log(p(\mathbf{y} \mid \mathbf{q})) \end{aligned} \tag{59}\]

then

\[ E_y\left[\frac{ \partial^2 -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_j \partial q_k }\middle| \mathbf{q}\right] = \int_y \frac{ \partial^2 -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_j \partial q_k } p(\mathbf{y} \mid \mathbf{q}) d\mathbf{y} = \tag{60}\]

\[ \int_y \left\{\frac{-1}{p(\mathbf{y} \mid \mathbf{q})} \frac{ \partial^2 p(\mathbf{y} \mid \mathbf{q}) }{ \partial q_j \partial q_k } + \frac{ \partial \log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial \log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_k }\right\} p(\mathbf{y} \mid \mathbf{q}) d\mathbf{y} \tag{61}\]

But

\[ \begin{aligned} \int_y \frac{-1}{p(\mathbf{y} \mid \mathbf{q})}\frac{ \partial^2 p(\mathbf{y} \mid \mathbf{q}) }{ \partial q_j \partial q_k } p(\mathbf{y} \mid \mathbf{q}) d\mathbf{y} &= -\int_y \frac{ \partial^2 p(\mathbf{y} \mid \mathbf{q}) }{ \partial q_j \partial q_k } d\mathbf{y} \\ &= -\frac{\partial^2}{\partial q_j \partial q_k}\int_y p(\mathbf{y} \mid \mathbf{q}) d\mathbf{y} \\ &= -\frac{\partial^2 1}{\partial q_j \partial q_k} \\ &= 0 \end{aligned} \tag{62}\]

\[ E_y\left[\frac{ \partial^2 -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_j \partial q_k }\middle| \mathbf{q}\right] = \int_y \frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_k } p(\mathbf{y} \mid \mathbf{q}) d\mathbf{y} = \tag{63}\]

\[ E_y\left[\frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_k }\middle| \mathbf{q}\right] \tag{64}\]

We note that

\[ \int_y \frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_k } p(\mathbf{y} \mid \mathbf{q}) d\mathbf{y} = \tag{65}\]

\[ \int_y \left[\sum_{i_1=1}^{m} \frac{ \partial -\log(p_{i_1}(\mathbf{y}_{i_1} \mid \mathbf{q})) }{ \partial q_j }\right]\left[\sum_{i_2=1}^{m} \frac{ \partial -\log(p_{i_2}(\mathbf{y}_{i_2} \mid \mathbf{q})) }{ \partial q_k }\right] \prod_{i_3=1}^{m} p_{i_3}(\mathbf{y}_{i_3} \mid \mathbf{q}) d\mathbf{y} = \tag{66}\]

\[ \sum_{i_1=1}^{m}\sum_{i_2=1}^{m} \int_y \frac{ \partial -\log(p_{i_1}(\mathbf{y}_{i_1} \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial -\log(p_{i_2}(\mathbf{y}_{i_2} \mid \mathbf{q})) }{ \partial q_k } \prod_{i_3=1}^{m} p_{i_3}(\mathbf{y}_{i_3} \mid \mathbf{q}) d\mathbf{y} \tag{67}\]

Since

\[ \int_{y_i} p_i(\mathbf{y}_i \mid \mathbf{q}) d\mathbf{y}_i = 1 \tag{68}\]

then

\[ E_y\left[\frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_k }\middle| \mathbf{q}\right] = \sum_{i_1=1}^{m}\sum_{i_2=1}^{m} \int_y \frac{ \partial -\log(p_{i_1}(\mathbf{y}_{i_1} \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial -\log(p_{i_2}(\mathbf{y}_{i_2} \mid \mathbf{q})) }{ \partial q_k } \prod_{i_3=1}^{m} p_{i_3}(\mathbf{y}_{i_3} \mid \mathbf{q}) d\mathbf{y} = \tag{69}\]

\[ \sum_{\substack{i_1=1 \\ i_2 \neq i_1}}^{m}\sum_{i_2=1}^{m} \int_{y_{i_1}}\int_{y_{i_2}} \frac{ \partial -\log(p_{i_1}(\mathbf{y}_{i_1} \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial -\log(p_{i_2}(\mathbf{y}_{i_2} \mid \mathbf{q})) }{ \partial q_k } p_{i_1}(\mathbf{y}_{i_1} \mid \mathbf{q}) p_{i_2}(\mathbf{y}_{i_2} \mid \mathbf{q}) d\mathbf{y}_{i_1} d\mathbf{y}_{i_2} + \tag{70}\]

\[ \sum_{i=1}^{m} \int_{y_i} \frac{ \partial -\log(p_i(\mathbf{y}_i \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial -\log(p_i(\mathbf{y}_i \mid \mathbf{q})) }{ \partial q_k } p_i(\mathbf{y}_i \mid \mathbf{q}) d\mathbf{y}_i \tag{71}\]

But

\[ \sum_{\substack{i_1=1 \\ i_2 \neq i_1}}^{m}\sum_{i_2=1}^{m}\left[\int_{y_{i_1}} \frac{ \partial -\log(p_{i_1}(\mathbf{y}_{i_1} \mid \mathbf{q})) }{ \partial q_j } p_{i_1}(\mathbf{y}_{i_1} \mid \mathbf{q}) d\mathbf{y}_{i_1}\right]\left[\int_{y_{i_2}} \frac{ \partial -\log(p_{i_2}(\mathbf{y}_{i_2} \mid \mathbf{q})) }{ \partial q_k } p_{i_2}(\mathbf{y}_{i_2} \mid \mathbf{q}) d\mathbf{y}_{i_2}\right] = \tag{73}\]

\[ \sum_{\substack{i_1=1 \\ i_2 \neq i_1}}^{m}\sum_{i_2=1}^{m}\left[\int_{y_{i_1}} \frac{ \partial p_{i_1}(\mathbf{y}_{i_1} \mid \mathbf{q}) }{ \partial q_j } d\mathbf{y}_{i_1}\right]\left[\int_{y_{i_2}} \frac{ \partial p_{i_2}(\mathbf{y}_{i_2} \mid \mathbf{q}) }{ \partial q_k } d\mathbf{y}_{i_2}\right] = \tag{74}\]

\[ \sum_{\substack{i_1=1 \\ i_2 \neq i_1}}^{m}\sum_{i_2=1}^{m} \frac{\partial 1}{\partial q_j}\frac{\partial 1}{\partial q_k} = 0 \tag{75}\]

Furthermore,

\[ \begin{aligned} \frac{ \partial -\log(p_i(\mathbf{y}_i \mid \mathbf{q})) }{ \partial q_j } & = \frac{1}{p_i(\mathbf{y}_i \mid \mathbf{q})}\frac{ \partial -p_i(\mathbf{y}_i \mid \mathbf{q}) }{ \partial q_j } = -\int \frac{ [\partial p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q}) / \partial q_j] d\boldsymbol{\phi} }{ p_i(\mathbf{y}_i \mid \mathbf{q}) } \\ &= \int \frac{ [\partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) / \partial q_j] p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q}) }{ p_i(\mathbf{y}_i \mid \mathbf{q}) } d\boldsymbol{\phi} \\ &= \int \left[\frac{ \partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) }{ \partial q_j }\right] z(\boldsymbol{\phi} \mid \mathbf{y}_i, \mathbf{q}) d\boldsymbol{\phi} \\ &= E_\phi\left[\frac{ \partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) }{ \partial q_j }\middle| \mathbf{y}_i, \mathbf{q}\right] \end{aligned} \tag{76}\]

So that

\[ E_y\left[\sum_{i=1}^{m} E_\phi\left[\frac{ \partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) }{ \partial q_j }\middle| \mathbf{y}_i, \mathbf{q}\right] E_\phi\left[\frac{ \partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) }{ \partial q_k }\middle| \mathbf{y}_i, \mathbf{q}\right]\middle| \mathbf{q}\right] = \tag{78}\]

\[ \sum_{i=1}^{m} E_y\left[E_\phi\left[\frac{ \partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) }{ \partial q_j }\middle| \mathbf{y}_i, \mathbf{q}\right] E_\phi\left[\frac{ \partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) }{ \partial q_k }\middle| \mathbf{y}_i, \mathbf{q}\right]\middle| \mathbf{q}\right] \tag{79}\]

\[ E_y\left[\frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_j }\frac{ \partial -\log(p(\mathbf{y} \mid \mathbf{q})) }{ \partial q_k }\middle| \mathbf{q}\right] = \sum_{i=1}^{m} E_y\left[E_\phi\left[\frac{ \partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) }{ \partial q_j }\middle| \mathbf{y}_i, \mathbf{q}\right] E_\phi\left[\frac{ \partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) }{ \partial q_k }\middle| \mathbf{y}_i, \mathbf{q}\right]\middle| \mathbf{q}\right] \tag{80}\]

We define

\[ \mathbf{g}_i = E_\phi\left[\frac{ \partial -\log(p_i(\mathbf{y}_i, \boldsymbol{\phi} \mid \mathbf{q})) }{ \partial \mathbf{q} }\middle| \mathbf{y}_i, \mathbf{q}\right] \tag{81}\]

and is the contribution of data from individual $i$ to the total gradient $\mathbf{g}$, where

\[ \begin{aligned} \mathbf{g} &= \sum_{i=1}^{m} \mathbf{g}_i \\ &= \mathbf{0} \end{aligned} \tag{82}\]

at the minimum.

The gradient components are evaluated by methods of differentiation of matrix algebra.

\[ \mathbf{g}_{i\boldsymbol{\theta}} = -\frac{ \partial \boldsymbol{\mu}'_i }{ \partial \boldsymbol{\theta} }\boldsymbol{\Omega}^{-1}(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \tag{83}\]

Now, for the lower triangular part of the covariance matrix,

\[ \mathbf{g}_{i\mathrm{Lower}(\boldsymbol{\Omega})} = \mathrm{Lower}\left[\boldsymbol{\Omega}^{-1}(\boldsymbol{\Omega} - \bar{\boldsymbol{\Omega}}_i)\boldsymbol{\Omega}^{-1} - \frac{1}{2}\mathrm{diag}\left(\boldsymbol{\Omega}^{-1}(\boldsymbol{\Omega} - \bar{\boldsymbol{\Omega}}_i)\boldsymbol{\Omega}^{-1}\right)\right] \tag{84}\]

To summarize,

\[ \mathrm{Var}(\mathbf{q}\mathbf{q}') = \left(\sum_{i=1}^{m} \mathbf{g}_i \mathbf{g}'_i\right)^{-1} \tag{85}\]

where $\sum_{i=1}^{m} \mathbf{g}_i \mathbf{g}'_i$ is known as the Fisher score matrix.

One caveat is in order. Note that the structure of the expected value second derivative can be written in the form $\mathbf{G}\mathbf{G}'$ where $\mathbf{G}$ is a $k \times m$ matrix having $m$ column vectors of $\mathbf{g}_i$ of size $k \times 1$, and $k$ is the total number of population parameters. If the number of subjects $m$ is less than $k$, then the $k \times k$ Fischer score matrix has only a rank of $m$, and the matrix is not invertible. Thus, this manner of constructing the expected value second derivative only holds true as the number of subjects $m$ as well as the number of data points per subject approaches infinity. Put another way, increasing the number of data points towards infinity for each subject, while having only a limited number of subjects, especially $k < m$, will not lead the Fischer score matrix to approach the expected value second derivative. For such conditions, the exact second derivative should be evaluated as given in appendix D.

Appendix D: The Exact Second Derivative Matrix for Error Assessment

The second derivative matrix for the population parameters and the population variance parameters, when the population parameter density is normally distributed, is determined as in the previous appendix, but without eliminating the terms that are canceled when taking the expected value over all $\mathbf{y}$. For the second derivative for a normal population parameter, we note the following:

\[ \begin{aligned} \frac{ \partial^2 O }{ \partial \mathbf{q} \partial \mathbf{q}' } &= \frac{\partial}{\partial \mathbf{q}}\left(\frac{\partial O}{\partial \mathbf{q}'}\right) = \frac{\partial^2}{\partial \mathbf{q} \partial \mathbf{q}'}\sum_{i=1}^{m} -\log(p_i) \int z_i(\mathbf{q}) d\boldsymbol{\phi} \\ &= \sum_{i=1}^{m}\int \frac{ \partial^2 -\log(p_i) }{ \partial \mathbf{q} \partial \mathbf{q}' } z_i(\mathbf{q}) d\boldsymbol{\phi} - \sum_{i=1}^{m}\int \frac{\partial -\log(p_i)}{\partial \mathbf{q}}\frac{\partial -\log(p_i)}{\partial \mathbf{q}'} z_i(\mathbf{q}) d\boldsymbol{\phi} \\ &\quad + \sum_{i=1}^{m}\left[\int \frac{\partial -\log(p_i)}{\partial \mathbf{q}} z_i(\mathbf{q}) d\boldsymbol{\phi}\right]\left[\int \frac{\partial -\log(p_i)}{\partial \mathbf{q}'} z_i(\mathbf{q}) d\boldsymbol{\phi}\right] \end{aligned} \tag{86}\]

where we let $\mathbf{q}$ represent the vector of all population parameters and variances, and

\[ \begin{aligned} p_i &= p(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\mu}_i(\boldsymbol{\theta}), \boldsymbol{\Omega}) \\ &= l(\mathbf{y}_i \mid \boldsymbol{\phi}) h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i(\boldsymbol{\theta}), \boldsymbol{\Omega}) \end{aligned} \tag{87}\]

\[ \mathbf{g}_i = \int \frac{\partial -\log(p_i)}{\partial \mathbf{q}} z_i(\mathbf{q}) d\boldsymbol{\phi} \tag{88}\]

We already know the first derivatives contributed by each individual $i$. For mu modeled thetas,

\[ \frac{\partial -\log(p_i)}{\partial \theta_{\mu_{j_1}}} = -\frac{ \partial \boldsymbol{\mu}'_i }{ \partial \theta_{\mu_{j_1}} }\boldsymbol{\Omega}^{-1}(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \tag{89}\]

Also,

\[ \frac{\partial -\log(p_i)}{\partial \omega_{j_1 j_2}} = -c(j_1, j_2)\mathbf{I}'_{j_1}\boldsymbol{\Omega}^{-1}\left[(\boldsymbol{\Omega} - \bar{\boldsymbol{\Omega}}_i) - (\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)'\right]\boldsymbol{\Omega}^{-1}\mathbf{I}_{j_2} \tag{90}\]

where $\mathbf{I}_{j_1}$ is a vector of 0’s for every element except for element $j_1$, which is valued at 1. Also,

\[ \begin{aligned} c(j_1, j_2) &= 1\ \text{for}\ j_1 \neq j_2 \\ &= 1/2\ \text{for}\ j_1 = j_2 \end{aligned} \tag{91}\]

Finally, for non-mu modeled thetas, a finite difference formula is employed:

\[ \frac{\partial -\log(p_i)}{\partial \theta_j} \approx \frac{ -[\log(p_i(\theta_j + \Delta\theta_j)) - \log(p_i(\theta_j - \Delta\theta_j))] }{ 2\Delta\theta_j } \tag{92}\]

We now derive the second derivatives

\[ \frac{\partial}{\partial \theta_{\mu_{j_2}}}\left(\frac{\partial -\log(p_i)}{\partial \theta_{\mu_{j_1}}}\right) = \frac{ \partial \boldsymbol{\mu}'_i }{ \partial \theta_{\mu_{j_1}} }\boldsymbol{\Omega}^{-1}\frac{ \partial \boldsymbol{\mu}_i }{ \partial \theta_{\mu_{j_2}} } - \frac{ \partial^2 \boldsymbol{\mu}'_i }{ \partial \theta_{\mu_{j_2}}\partial \theta_{\mu_{j_1}} }\boldsymbol{\Omega}^{-1}(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \tag{93}\]

\[ \begin{aligned} \frac{\partial^2}{\partial \omega_{j_3 j_4}}\left(\frac{\partial -\log(p_i)}{\partial \omega_{j_1 j_2}}\right) & = c(j_1, j_2)c(j_3, j_4)\mathbf{I}'_{j_1}\boldsymbol{\Omega}^{-1}(\mathbf{I}_{j_3}\mathbf{I}'_{j_4} + \mathbf{I}_{j_4}\mathbf{I}'_{j_3})\boldsymbol{\Omega}^{-1}\mathbf{I}_{j_2} \\ &\quad + 2c(j_1, j_2)c(j_3, j_4)\mathbf{I}'_{j_1}\boldsymbol{\Omega}^{-1}(\mathbf{I}_{j_3}\mathbf{I}'_{j_4} + \mathbf{I}_{j_4}\mathbf{I}'_{j_3})\boldsymbol{\Omega}^{-1}\left[(\boldsymbol{\Omega} - \bar{\boldsymbol{\Omega}}_i) - (\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)'\right]\boldsymbol{\Omega}^{-1}\mathbf{I}_{j_2} \end{aligned} \tag{94}\]

\[ \frac{\partial}{\partial \omega_{j_2 j_3}}\left(\frac{\partial -\log(p_i)}{\partial \theta_{\mu_{j_1}}}\right) = c(j_2, j_3)\frac{ \partial \boldsymbol{\mu}'_i }{ \partial \theta_{\mu_{j_1}} }\boldsymbol{\Omega}^{-1}(\mathbf{I}_{j_2}\mathbf{I}'_{j_3} + \mathbf{I}_{j_3}\mathbf{I}'_{j_2})\boldsymbol{\Omega}^{-1}(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \tag{95}\]

Also:

\[ \frac{ \partial^2 -\log(p_i) }{ \partial \theta_{j_1}\partial \theta_{j_2} } = \frac{\partial}{\partial \theta_{j_1}}\left(\frac{\partial -\log(p_i)}{\partial \theta_{j_2}}\right) \approx \frac{ -\log(p_i(\theta_{j_1}+\Delta\theta_{j_1},\theta_{j_2}+\Delta\theta_{j_2})) + \log(p_i(\theta_{j_1}+\Delta\theta_{j_1},\theta_{j_2} -\Delta\theta_{j_2})) }{ 4\Delta\theta_{j_1}\Delta\theta_{j_2} } + \frac{ -\log(p_i(\theta_{j_1}-\Delta\theta_{j_1},\theta_{j_2}+\Delta\theta_{j_2})) + \log(p_i(\theta_{j_1}-\Delta\theta_{j_1},\theta_{j_2} -\Delta\theta_{j_2})) }{ 4\Delta\theta_{j_1}\Delta\theta_{j_2} } \tag{96}\]

And, since $\boldsymbol{\phi}$ show up only in the $h()$ portion of the joint density, and $\boldsymbol{\varepsilon}$ show up only in the $l()$ portion of the joint density,

\[ \frac{\partial}{\partial \sigma_{j_1}}\left(\frac{\partial -\log(p_i)}{\partial \omega_{j_2}}\right) = 0 \tag{97}\]

Defining the individual subject central moments:

\[ \phi^{(1)}_{ir_1} = \int -(\phi_{r_1} - \mu_{ir_1}) z(\boldsymbol{\phi} \mid \mathbf{y}_i, \mathbf{q}) d\boldsymbol{\phi} \tag{98}\]

\[ \phi^{(2)}_{ir_1 r_2} = \int -(\phi_{r_1} - \mu_{ir_1})(\phi_{r_2} - \mu_{ir_2}) z(\boldsymbol{\phi} \mid \mathbf{y}_i, \mathbf{q}) d\boldsymbol{\phi} \tag{99}\]

\[ \phi^{(3)}_{ir_1 r_2 r_3} = \int -(\phi_{r_1} - \mu_{ir_1})(\phi_{r_2} - \mu_{ir_2})(\phi_{r_3} - \mu_{ir_3}) z(\boldsymbol{\phi} \mid \mathbf{y}_i, \mathbf{q}) d\boldsymbol{\phi} \tag{100}\]

\[ \phi^{(4)}_{ir_1 r_2 r_3 r_4} = \int -(\phi_{r_1} - \mu_{ir_1})(\phi_{r_2} - \mu_{ir_2})(\phi_{r_3} - \mu_{ir_3})(\phi_{r_4} - \mu_{ir_4}) z(\boldsymbol{\phi} \mid \mathbf{y}_i, \mathbf{q}) d\boldsymbol{\phi} \tag{101}\]

\[ \sigma^{(2)}_{ir_1 r_2} = \int s_{ir_1} s_{ir_2} z(\boldsymbol{\phi} \mid \mathbf{y}_i, \mathbf{q}) d\boldsymbol{\phi} \tag{102}\]

\[ \sigma^{(2)}_{ir_1 r_2} = \int -(\phi_{r_1} - \mu_{ir_1}) s_{ir_2} z(\boldsymbol{\phi} \mid \mathbf{y}_i, \mathbf{q}) d\boldsymbol{\phi} \tag{103}\]

\[ \sigma^{(3)}_{ir_1 r_2 r_3} = \int -s_{ir_1}(\phi_{r_2} - \mu_{ir_2})(\phi_{r_3} - \mu_{ir_3}) z(\boldsymbol{\phi} \mid \mathbf{y}_i, \mathbf{q}) d\boldsymbol{\phi} \tag{104}\]

where

\[ s_{ir_1} = \frac{\partial -\log(p_i)}{\partial \theta_{\mu_{r_1}}} \tag{105}\]

and the overall central moments:

\[ \phi^{(1)}_{r_1} = \frac{1}{m}\sum_{i=1}^{m} \phi^{(1)}_{ir_1} \tag{106}\]

\[ \phi^{(2)}_{r_1 r_2} = \frac{1}{m}\sum_{i=1}^{m} \phi^{(2)}_{ir_1 r_2} \tag{107}\]

\[ \phi^{(3)}_{r_1 r_2 r_3} = \frac{1}{m}\sum_{i=1}^{m} \phi^{(3)}_{ir_1 r_2 r_3} \tag{108}\]

\[ \phi^{(4)}_{r_1 r_2 r_3 r_4} = \frac{1}{m}\sum_{i=1}^{m} \phi^{(4)}_{ir_1 r_2 r_3 r_4} \tag{109}\]

We now have the following:

\[ \frac{ \partial^2 O }{ \partial \theta_{\mu_{j_2}}\partial \theta_{\mu_{j_1}} } = \sum_{i=1}^{m}\sum_{r_1}^{n}\sum_{r_2=1}^{n}\sum_{r_3=1}^{n}\sum_{r_4=1}^{n} \frac{ \partial \boldsymbol{\mu}_{ir_1} }{ \partial \theta_{\mu_{j_1}} }\omega^{-1}_{r_1 r_2}\left(\omega_{r_2 r_3} - \phi^{(2)}_{ir_2 r_3}\right)\omega^{-1}_{r_3 r_4}\frac{ \partial \boldsymbol{\mu}_{ir_4} }{ \partial \theta_{\mu_{j_2}} } - \sum_{i=1}^{m}\sum_{r_1}^{n}\sum_{r_2=1}^{n}\frac{ \partial^2 \boldsymbol{\mu}_{ir_1} }{ \partial \theta_{\mu_{j_2}}\partial \theta_{\mu_{j_1}} }\omega^{-1}_{r_1 r_2}\phi^{(1)}_{ir_2} + \sum_{i=1}^{m} g_{i\theta_{\mu_{j_1}}} g_{i\theta_{\mu_{j_2}}} \tag{110}\]

\[ \begin{aligned} \frac{ \partial^2 O }{ \partial \omega_{j_3 j_4}\partial \omega_{j_1 j_2} } &= mc(j_1, j_2)c(j_3, j_4)\left(\omega^{-1}_{j_1 j_3}\omega^{-1}_{j_4 j_2} + \omega^{-1}_{j_1 j_4}\omega^{-1}_{j_3 j_2}\right) \\ &\quad + 2mc(j_1, j_2)c(j_3, j_4)\sum_{r_1=1}^{n}\sum_{r_2=1}^{n}\omega^{-1}_{j_1 j_3}\omega^{-1}_{j_4 r_1}\left(\omega_{r_1 r_2} - \phi^{(2)}_{r_1 r_2}\right)\omega^{-1}_{r_2 j_2} \\ &\quad + 2mc(j_1, j_2)c(j_3, j_4)\sum_{r_1=1}^{n}\sum_{r_2=1}^{n}\omega^{-1}_{j_1 j_4}\omega^{-1}_{j_3 r_1}\left(\omega_{r_1 r_2} - \phi^{(2)}_{r_1 r_2}\right)\omega^{-1}_{r_2 j_2} \\ &\quad - mc(j_1, j_2)c(j_3, j_4)\sum_{r_1=1}^{n}\sum_{r_2=1}^{n}\sum_{r_3=1}^{n}\sum_{r_4=1}^{n}\phi^{(4)}_{r_1 r_2 r_3 r_4}\omega^{-1}_{j_1 r_1}\omega^{-1}_{j_2 r_2}\omega^{-1}_{j_3 r_3}\omega^{-1}_{j_4 r_4} \\ &\quad + mc(j_1, j_2)c(j_3, j_4)\sum_{r_1=1}^{n}\sum_{r_2=1}^{n}\sum_{r_3=1}^{n}\sum_{r_4=1}^{n}\omega^{-1}_{j_1 r_1}\phi^{(2)}_{r_1 r_2}\omega^{-1}_{r_2 j_2}\omega^{-1}_{j_3 j_4} \\ &\quad + mc(j_1, j_2)c(j_3, j_4)\sum_{r_1=1}^{n}\sum_{r_2=1}^{n}\sum_{r_3=1}^{n}\sum_{r_4=1}^{n}\omega^{-1}_{j_1 j_2}\omega^{-1}_{j_3 r_1}\phi^{(2)}_{r_1 r_2}\omega^{-1}_{r_2 j_4} \\ &\quad - mc(j_1, j_2)c(j_3, j_4)\omega^{-1}_{j_1 j_2}\omega^{-1}_{j_3 j_4} + \sum_{i=1}^{m} g_{i\omega_{j_1 j_2}}g_{i\omega_{j_3 j_4}} \end{aligned} \tag{111}\]

\[ \begin{aligned} \frac{ \partial^2 O }{ \partial \omega_{j_2 j_3}\partial \theta_{\mu_{j_1}} } & = c(j_2, j_3)\sum_{r_1=1}^{n}\sum_{r_2=1}^{n}\sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_{ir_1} }{ \partial \theta_{\mu_{j_1}} }\phi^{(2)}_{ir_2}\left(\omega^{-1}_{r_1 j_2}\omega^{-1}_{j_3 r_2} + \omega^{-1}_{r_1 j_3}\omega^{-1}_{j_2 r_2}\right) \\ &\quad - c(j_2, j_3)\sum_{i=1}^{m}\sum_{r_4=1}^{n} \frac{ \partial \boldsymbol{\mu}_{ir_4} }{ \partial \theta_{\mu_{j_1}} }\omega^{-1}_{r_4 r_1}\phi^{(3)}_{ir_1 r_2 r_3}\omega^{-1}_{j_2 r_2}\omega^{-1}_{j_3 r_3} - a(j_2, j_3)g^c_{j_1} + \sum_{i=1}^{m} g_{i\theta_{\mu_{j_1}}}g_{i\omega_{j_2 j_3}} \end{aligned} \tag{112}\]

where

\[ g_{j_1} = \sum_{i=1}^{m} g_{i\theta_{\mu_{j_1}}} \tag{113}\]

etcetera. And finally:

For the non-mu parameters:

\[ \frac{ \partial^2 O }{ \partial \theta_{j_1}\partial \theta_{j_2} } = \sum_{i=1}^{m}\int \frac{ \partial^2 -\log(p_i) }{ \partial \theta_{j_1}\partial \theta_{j_2} } z_i(\mathbf{q}) d\boldsymbol{\phi} - \sum_{i=1}^{m}\sigma^{(2)}_{ij_1 j_2} + \sum_{i=1}^{m} g_{i\theta_{j_1}}g_{i\theta_{j_2}} \tag{114}\]

\[ \frac{ \partial^2 O }{ \partial \theta_{\mu_{j_1}}\partial \theta_{j_2} } = \sum_{i=1}^{m}\sum_{r_1}^{n} \frac{ \partial \boldsymbol{\mu}_{ir_1} }{ \partial \theta_{\mu_{j_1}} }\omega^{-1}_{j_2 r_1}\sigma^{(2)}_{ij_2} + \sum_{i=1}^{m} g_{i\theta_{\mu_{j_1}}}g_{i\theta_{j_2}} \tag{115}\]

\[ \frac{ \partial^2 O }{ \partial \omega_{j_2 j_3}\partial \theta_{j_1} } = c(j_2, j_3)\sum_{i=1}^{m}\sum_{r_1}^{n}\sum_{r_2}^{n} \sigma^{(3)}_{ij_1 r_1 r_2}\omega^{-1}_{j_2 r_1}\omega^{-1}_{r_2 j_3} + \sum_{i=1}^{m} g_{i\theta_{j_1}}g_{i\omega_{j_2 j_3}} \tag{116}\]

At the minimum of the objective function,

\[ g_{j_1} = 0 \tag{117}\]

\[ g_{j_1} = 0 \tag{118}\]

\[ \omega_{r_1 r_2} = \phi^{(2)}_{r_1 r_2} \tag{119}\]

we then obtain following simplifications:

\[ \begin{aligned} \frac{ \partial^2 O }{ \partial \theta_{\mu_{j_2}}\partial \theta_{\mu_{j_1}} } &= \sum_{i=1}^{m} g_{i\theta_{\mu_{j_1}}}g_{i\theta_{\mu_{j_2}}} - \sum_{i=1}^{m}\sum_{r_1}^{n}\sum_{r_2=1}^{n}\sum_{r_3=1}^{n}\sum_{r_4=1}^{n} \frac{ \partial \boldsymbol{\mu}_{ir_1} }{ \partial \theta_{\mu_{j_1}} }\omega^{-1}_{r_1 r_2}\left(\omega_{r_2 r_3} - \phi^{(2)}_{ir_2 r_3}\right)\omega^{-1}_{r_3 r_4}\frac{ \partial \boldsymbol{\mu}_{ir_4} }{ \partial \theta_{\mu_{j_2}} } \\ &\quad - \sum_{i=1}^{m}\sum_{r_1}^{n}\sum_{r_2=1}^{n}\frac{ \partial^2 \boldsymbol{\mu}_{ir_1} }{ \partial \theta_{\mu_{j_2}}\partial \theta_{\mu_{j_1}} }\omega^{-1}_{r_1 r_2}\phi^{(1)}_{ir_2} \end{aligned} \tag{120}\]

\[ \begin{aligned} \frac{ \partial^2 O }{ \partial \omega_{j_3 j_4}\partial \omega_{j_1 j_2} } &= mc(j_1, j_2)c(j_3, j_4)\left(\omega^{-1}_{j_1 j_2}\omega^{-1}_{j_3 j_4} + \omega^{-1}_{j_1 j_3}\omega^{-1}_{j_4 j_2} + \omega^{-1}_{j_1 j_4}\omega^{-1}_{j_3 j_2}\right) \\ &\quad - mc(j_1, j_2)c(j_3, j_4)\sum_{r_1=1}^{n}\sum_{r_2=1}^{n}\sum_{r_3=1}^{n}\sum_{r_4=1}^{n}\phi^{(4)}_{r_1 r_2 r_3 r_4}\omega^{-1}_{j_1 r_1}\omega^{-1}_{j_2 r_2}\omega^{-1}_{j_3 r_3}\omega^{-1}_{j_4 r_4} + \sum_{i=1}^{m} g_{i\omega_{j_1 j_2}}g_{i\omega_{j_3 j_4}} \end{aligned} \tag{121}\]

\[ \begin{aligned} \frac{ \partial^2 O }{ \partial \omega_{j_2 j_3}\partial \theta_{\mu_{j_1}} } & = c(j_2, j_3)\sum_{r_1=1}^{n}\sum_{r_2=1}^{n}\sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_{ir_1} }{ \partial \theta_{\mu_{j_1}} }\phi^{(2)}_{ir_2}\left(\omega^{-1}_{r_1 j_2}\omega^{-1}_{j_3 r_2} + \omega^{-1}_{r_1 j_3}\omega^{-1}_{j_2 r_2}\right) \\ &\quad - c(j_2, j_3)\sum_{i=1}^{m}\sum_{r_4=1}^{n} \frac{ \partial \boldsymbol{\mu}_{ir_4} }{ \partial \theta_{\mu_{j_1}} }\omega^{-1}_{r_4 r_1}\phi^{(3)}_{ir_1 r_2 r_3}\omega^{-1}_{j_2 r_2}\omega^{-1}_{j_3 r_3} + \sum_{i=1}^{m} g_{i\theta_{\mu_{j_1}}}g_{i\omega_{j_2 j_3}} \end{aligned} \tag{122}\]

Let us see that, if we take the expected value over all $\mathbf{y}$, and as the number of subset of subjects $m_i$ sharing a particular set of covariates (and therefore sharing the same $\mu_{ir}$), approaches infinity, we should obtain the results in Appendix C. First, note that:

\[ \int_y \phi^{(3)}_{ir_1 r_2 r_3} p(\mathbf{y} \mid \mathbf{q}) d\mathbf{y} = \frac{1}{m_i}\sum_{j=1}^{m_i}\int_y -(\phi_{r_1} - \mu_{ir_1})(\phi_{r_2} - \mu_{ir_2})(\phi_{r_3} - \mu_{ir_3}) z(\boldsymbol{\phi} \mid \mathbf{y}_j, \mathbf{q}) d\boldsymbol{\phi} p(\mathbf{y} \mid \mathbf{q}) d\mathbf{y} = \tag{126}\]

\[ \frac{1}{m_i}\sum_{j=1}^{m_i}\int_{y_j}\int -(\phi_{r_1} - \mu_{ir_1})(\phi_{r_2} - \mu_{ir_2})(\phi_{r_3} - \mu_{ir_3}) l(\mathbf{y}_j \mid \boldsymbol{\phi}) h(\boldsymbol{\phi} \mid \mathbf{q}) d\boldsymbol{\phi} d\mathbf{y}_j = \tag{127}\]

\[ \frac{1}{m_i}\sum_{j=1}^{m_i}\int -(\phi_{r_1} - \mu_{ir_1})(\phi_{r_2} - \mu_{ir_2})(\phi_{r_3} - \mu_{ir_3}) h(\boldsymbol{\phi} \mid \mathbf{q})\int_{y_j} l(\mathbf{y}_j \mid \boldsymbol{\phi}) d\mathbf{y}_j = \tag{128}\]

\[ \frac{1}{m_i}\sum_{j=1}^{m_i}\int -(\phi_{r_1} - \mu_{ir_1})(\phi_{r_2} - \mu_{ir_2})(\phi_{r_3} - \mu_{ir_3}) h(\boldsymbol{\phi} \mid \mathbf{q}) d\boldsymbol{\phi} = \hat{\phi}^{(3)}_{ir_1 r_2 r_3} \tag{129}\]

Since $h()$ is a multi-variate normal density, so therefore the $\hat{\phi}^{(3)}_{ir_1 r_2 r_3}$ is the skewness of a normally distributed variable:

\[ \hat{\phi}^{(3)}_{ir_1 r_2 r_3} = 0 \tag{130}\]

Similarly,

\[ \hat{\phi}^{(1)}_{ir_1} = 0 \tag{131}\]

and

\[ \begin{aligned} \int_y \phi^{(4)}_{ir_1 r_2 r_3 r_4} p(\mathbf{y} \mid \mathbf{q}) d\mathbf{y} &= \hat{\phi}^{(4)}_{ir_1 r_2 r_3 r_4} \\ &= \omega_{r_1 r_2}\omega_{r_3 r_4} + \omega_{r_1 r_3}\omega_{r_2 r_4} + \omega_{r_1 r_4}\omega_{r_2 r_3} \end{aligned} \tag{132}\]

is derived from the kurtosis of a normally distributed random variable (Triantafyllopoulos 2003). We may now make the final simplification:

\[ E_y\left[\frac{ \partial^2 O }{ \partial \theta_{\mu_{j_2}}\partial \theta_{\mu_{j_1}} }\right] = \sum_{i=1}^{m} g_{i\theta_{\mu_{j_1}}}g_{i\theta_{\mu_{j_2}}} \tag{133}\]

\[ E_y\left[\frac{ \partial^2 O }{ \partial \omega_{j_3 j_4}\partial \omega_{j_1 j_2} }\right] = \sum_{i=1}^{m} g_{i\omega_{j_1 j_2}}g_{i\omega_{j_3 j_4}} \tag{134}\]

\[ E_y\left[\frac{ \partial^2 O }{ \partial \omega_{j_2 j_3}\partial \theta_{\mu_{j_1}} }\right] = \sum_{i=1}^{m} g_{i\theta_{\mu_{j_1}}}g_{i\omega_{j_2 j_3}} \tag{135}\]

These results are expected based on the general proof for the expected information matrix given in Appendix C for any population parameter density.

For population mixture parameters, the second derivatives are simply:

\[ \frac{\partial^2 O}{\partial \alpha \partial \alpha} = \sum_{i=1}^{m} g_{i\alpha} g_{i\alpha} \tag{136}\]

\[ \frac{\partial^2 O}{\partial \alpha \partial \theta} = \sum_{i=1}^{m} g_{i\alpha} g_{i\theta} \tag{137}\]

\[ \frac{ \partial^2 O }{ \partial \alpha \partial \boldsymbol{\omega} } = \sum_{i=1}^{m} g_{i\alpha} g_{i\boldsymbol{\omega}} \tag{138}\]

\[ \frac{ \partial^2 O }{ \partial \alpha \partial \boldsymbol{\sigma} } = \sum_{i=1}^{m} g_{i\alpha} g_{i\boldsymbol{\sigma}} \tag{139}\]

For Three hierarchical stage analysis, the second derivative matrices would have added to them:

\[ \frac{\partial^2 OP}{\partial(\boldsymbol{\theta})^2} = \boldsymbol{\Theta}^{-1} \tag{140}\]

for all theta parameters (mu or non-mu modeled). For inter-subject variance components:

\[ \begin{aligned} \frac{ \partial^2 OP }{ \partial \omega_{j_3 j_4}\partial \omega_{j_1 j_2} } & = z_Wc(j_1, j_2)c(j_3, j_4)\mathbf{I}'_{j_1}\boldsymbol{\Omega}^{-1}(\mathbf{I}_{j_3}\mathbf{I}'_{j_4} + \mathbf{I}_{j_4}\mathbf{I}'_{j_3})\boldsymbol{\Omega}^{-1}\mathbf{I}_{j_2} \\ &\quad + 2c(j_1, j_2)c(j_3, j_4)\mathbf{I}'_{j_1}\boldsymbol{\Omega}^{-1}(\mathbf{I}_{j_3}\mathbf{I}'_{j_4} + \mathbf{I}_{j_4}\mathbf{I}'_{j_3})\boldsymbol{\Omega}^{-1}(\boldsymbol{\Omega} - z_W\boldsymbol{\Omega})\boldsymbol{\Omega}^{-1}\mathbf{I}_{j_2} \\ & = c(j_1, j_2)c(j_3, j_4)\mathbf{I}'_{j_1}\boldsymbol{\Omega}^{-1}(\mathbf{I}_{j_3}\mathbf{I}'_{j_4} + \mathbf{I}_{j_4}\mathbf{I}'_{j_3})(2\boldsymbol{\Omega}^{-1} - z_W\boldsymbol{\Omega}^{-1})\mathbf{I}_{j_2} \end{aligned} \tag{141}\]

Similarly for Sigma parameters, we add the following:

\[ \begin{aligned} \frac{ \partial^2 OP }{ \partial \sigma_{j_3 j_4}\partial \sigma_{j_1 j_2} } & = z_{\Sigma}c(j_1, j_2)c(j_3, j_4)\mathbf{I}'_{j_1}\boldsymbol{\Sigma}^{-1}(\mathbf{I}_{j_3}\mathbf{I}'_{j_4} + \mathbf{I}_{j_4}\mathbf{I}'_{j_3})\boldsymbol{\Sigma}^{-1}\mathbf{I}_{j_2} \\ &\quad + 2c(j_1, j_2)c(j_3, j_4)\mathbf{I}'_{j_1}\boldsymbol{\Sigma}^{-1}(\mathbf{I}_{j_3}\mathbf{I}'_{j_4} + \mathbf{I}_{j_4}\mathbf{I}'_{j_3})\boldsymbol{\Sigma}^{-1}(\boldsymbol{\Sigma} - z_{\Sigma}\boldsymbol{\Sigma})\boldsymbol{\Sigma}^{-1}\mathbf{I}_{j_2} \\ & = c(j_1, j_2)c(j_3, j_4)\mathbf{I}'_{j_1}\boldsymbol{\Sigma}^{-1}(\mathbf{I}_{j_3}\mathbf{I}'_{j_4} + \mathbf{I}_{j_4}\mathbf{I}'_{j_3})(2\boldsymbol{\Sigma}^{-1} - z_{\Sigma}\boldsymbol{\Sigma}^{-1})\mathbf{I}_{j_2} \end{aligned} \tag{142}\]

If we designate the exact second derivative as described in this appendix as matrix $\mathbf{R}$, and the Fischer score matrix (appendix C) as $\mathbf{S}$, then we can construct a variance matrix in a manner similar to NONMEM:

\[ \mathrm{Var}(\mathbf{q}\mathbf{q}') = \mathbf{R}^{-1}\mathbf{S}\mathbf{R}^{-1} \tag{143}\]

The $\mathbf{R}$ matrix is not always numerically positive definite. For Monte Carlo assessed information matrices, the NONMEM program passes the matrix through a positive definiteness filter that makes small adjustments to the eigenvalues of $\mathbf{R}$, if necessary.

Appendix E: Adjustment of Error Matrix for Constraints and Non-Positive Definiteness

The user supplied subroutine CONSTRAINT allows the user to impose constraints on the population parameters. Therefore the error matrix must be adjusted to account for these constraints, and it is done as follows.

Let the constraint matrix $\mathbf{W}$ be defined as a matrix with elements:

\[ w_{ij} = \frac{\partial q_j}{\partial q_i} \tag{144}\]

where $q_j$ is the $j$th population parameter (or Omega variance), which could have a dependence on some other population parameter $q_i$. If no dependence exists for parameter $j$, then $w_{ij}=0$ for $i \neq j$, and $w_{ii}=1$. If no dependencies are defined for any parameter, then $\mathbf{W}=\mathbf{I}$, and no correction occurs. The Error matrix is corrected as follows:

For the R matrix type covariance:

\[ \mathrm{Var}(\mathbf{q}\mathbf{q}') = \mathbf{W}'(\mathbf{W}\mathbf{R}\mathbf{W}')^{-1}\mathbf{W} \tag{145}\]

The logic behind this equation is as follows. Parameters that are dependent on other parameters are considered secondary parameters, in contrast to the primary parameters that are independent of other parameters. The $(\mathbf{W}\mathbf{R}\mathbf{W}')^{-1}$ term creates the error matrix with rows and columns pertaining to the secondary parameters zeroed out, while the errors of the primary parameters are adjusted to account for the constraint on the model. This matrix is then multiplied on either side by $\mathbf{W}$ and $\mathbf{W}'$, to fill the zeroed secondary parameter rows and columns of $(\mathbf{W}\mathbf{R}\mathbf{W}')^{-1}$ with errors from the primary parameters, in accordance with their dependencies on the primary parameters. The resulting error matrix therefore contains errors to the primary as well as secondary parameters, and this matrix is placed in the var table by the poperr command, and the varc table by the poperr_corr command.

Similarly for the S matrix:

\[ \mathrm{Var}(\mathbf{q}\mathbf{q}') = \mathbf{W}'(\mathbf{W}\mathbf{S}\mathbf{W}')^{-1}\mathbf{W} \tag{146}\]

For the RSR matrix:

\[ \begin{aligned} \mathrm{Var}(\mathbf{q}\mathbf{q}') &= \mathbf{W}'\left[(\mathbf{W}\mathbf{R}\mathbf{W}')^{-1}\mathbf{W}\right]\mathbf{S}\left[\mathbf{W}'(\mathbf{W}\mathbf{R}\mathbf{W}')^{-1}\right]\mathbf{W} \\ &= \mathbf{W}'(\mathbf{W}\mathbf{R}\mathbf{W}')^{-1}(\mathbf{W}\mathbf{S}\mathbf{W}')(\mathbf{W}\mathbf{R}\mathbf{W}')^{-1}\mathbf{W} \end{aligned} \tag{147}\]

If a particular parameter is constrained to a fixed value, then $\mathbf{W}$ will be singular. The matrix $\mathbf{W}\mathbf{R}\mathbf{W}'$ is therefore inverted by the Jacobi method of extracting eignvalues and eigenvector matrices. That is, for any symmetric matrix $\mathbf{A}$, the Jacobian process decomposes the matrix to:

\[ \mathbf{A} = \mathbf{E}\boldsymbol{\Lambda}\mathbf{E}' \tag{148}\]

where $\boldsymbol{\Lambda}$ is the diagonal matrix of eigenvalues and $\mathbf{E}$ is a matrix of eigenvector columns, which has the property:

\[ \mathbf{E}' = \mathbf{E}^{-1} \tag{149}\]

The generalized inverse of $\mathbf{A}$ is then obtained as:

\[ \mathbf{A}^{-} = \mathbf{E}\boldsymbol{\Lambda}^{-}\mathbf{E}' \tag{150}\]

where $\boldsymbol{\Lambda}^{-}$ has diagonal elements of

\[ \begin{aligned} \lambda^{-}_i &= \frac{1}{\lambda_i}\ \text{for}\ |\lambda_i| > 0 \\ &= 0\ \text{for}\ |\lambda_i| = 0 \end{aligned} \tag{151}\]

In addition, the $\mathbf{R}$ matrix itself can at times be not positive definite (has negative eigenvalues), because of the imprecision of evaluating this matrix using random sampling, in the manner decribed in appendix D. It has been found in practice that using the absolute value of the eigenvalues to evaluate the inverse for Monte Carlo constructed information matrices effectively yields satisfactory error matrices. This is because negative eigenvalues are usually close to 0, and arise in the least important portions of the matrix.

Appendix F: Obtaining Analytical Derivatives of Likelihood with Respect to Cholesky of Sigma Parameters

The following is used to provide more rapid analysis for the importance sampling, direct sampling, and SAEM methods.

We concern ourselves with the derivatives of the likelihood:

\[ -\log(l(\mathbf{y}_i \mid \boldsymbol{\phi})) = \frac{1}{2}(\mathbf{y}_i - \mathbf{f}_i)'\mathbf{V}_i^{-1}(\mathbf{y}_i - \mathbf{f}_i) + \frac{1}{2}\log(|\mathbf{V}_i|) \tag{152}\]

With respect to the Sigma parameters. Consider that the Sigma parameters are involved in the construction of the residual variance matrix as follows:

\[ \begin{aligned} \mathbf{V}_i &= \mathbf{H}_i\boldsymbol{\Sigma}\mathbf{H}'_i \\ &= \mathbf{H}_i\boldsymbol{\Lambda}\boldsymbol{\Lambda}'\mathbf{H}'_i \end{aligned} \tag{153}\]

Where $\boldsymbol{\Lambda}$ is the cholesky of $\boldsymbol{\Sigma}$. It is the elements of $\boldsymbol{\Lambda}$ that are actually varied to optimize the objective function, therefore we wish to determine the derivative with respect to the elements of $\boldsymbol{\Lambda}$.

First derivatives:

\[ \frac{ \partial -\log(l(\mathbf{y}_i \mid \boldsymbol{\phi})) }{ \partial \lambda_{j_1 k_1} } = \frac{1}{2}(\mathbf{y}_i - \mathbf{f}_i)'\frac{ \partial \mathbf{V}_i^{-1} }{ \partial \lambda_{j_1 k_1} }(\mathbf{y}_i - \mathbf{f}_i) + \frac{1}{2}\frac{ \partial \log(|\mathbf{V}_i|) }{ \partial \lambda_{j_1 k_1} } = \tag{154}\]

\[ -\frac{1}{2}(\mathbf{y}_i - \mathbf{f}_i)'\mathbf{V}_i^{-1}\frac{\partial \mathbf{V}_i}{\partial \lambda_{j_1 k_1}}\mathbf{V}_i^{-1}(\mathbf{y}_i - \mathbf{f}_i) + \frac{1}{2}\mathrm{tr}\left(\mathbf{V}_i^{-1}\frac{\partial \mathbf{V}_i}{\partial \lambda_{j_1 k_1}}\right) \tag{155}\]

\[ \begin{aligned} \frac{\partial \mathbf{V}_i}{\partial \lambda_{j_1 k_1}} &= \mathbf{H}_i\frac{ \partial \boldsymbol{\Lambda}\boldsymbol{\Lambda}' }{ \partial \lambda_{j_1 k_1} }\mathbf{H}'_i \\ &= \mathbf{H}_i\frac{ \partial \boldsymbol{\Sigma} }{ \partial \lambda_{j_1 k_1} }\mathbf{H}'_i \end{aligned} \tag{156}\]

\[ \begin{aligned} \frac{ \partial \boldsymbol{\Sigma} }{ \partial \lambda_{j_1 k_1} } &= \frac{ \partial \boldsymbol{\Lambda}\boldsymbol{\Lambda}' }{ \partial \lambda_{j_1 k_1} } = \frac{ \partial \boldsymbol{\Lambda} }{ \partial \lambda_{j_1 k_1} }\boldsymbol{\Lambda}' + \boldsymbol{\Lambda}\frac{ \partial \boldsymbol{\Lambda}' }{ \partial \lambda_{j_1 k_1} } = \mathbf{I}_{j_1}\mathbf{I}'_{k_1}\boldsymbol{\Lambda}' + \boldsymbol{\Lambda}\mathbf{I}_{k_1}\mathbf{I}'_{j_1} = \\ &\quad\{\delta(j, j_1)\lambda_{k k_1} + \delta(k, j_1)\lambda_{j k_1}\ \text{for}\ j = 1\ \text{to}\ n,\ k = 1\ \text{to}\ n\} \end{aligned} \tag{157}\]

Second derivatives:

\[ \frac{ \partial^2 -\log(l(\mathbf{y}_i \mid \boldsymbol{\phi})) }{ \partial \lambda_{j_2 k_2}\partial \lambda_{j_1 k_1} } = \tag{158}\]

\[ \begin{aligned} &(\mathbf{y}_i - \mathbf{f}_i)'\mathbf{V}_i^{-1}\frac{\partial \mathbf{V}_i}{\partial \lambda_{j_2 k_2}}\mathbf{V}_i^{-1}\frac{\partial \mathbf{V}_i}{\partial \lambda_{j_1 k_1}}\mathbf{V}_i^{-1}(\mathbf{y}_i - \mathbf{f}_i) - \frac{1}{2}(\mathbf{y}_i - \mathbf{f}_i)'\mathbf{V}_i^{-1}\frac{ \partial^2 \mathbf{V}_i }{ \partial \lambda_{j_2 k_2}\partial \lambda_{j_1 k_1} }\mathbf{V}_i^{-1}(\mathbf{y}_i - \mathbf{f}_i) \\ &\quad + \frac{1}{2}\mathrm{tr}\left(\mathbf{V}_i^{-1}\frac{ \partial^2 \mathbf{V}_i }{ \partial \lambda_{j_2 k_2}\partial \lambda_{j_1 k_1} }\right) - \frac{1}{2}\mathrm{tr}\left(\mathbf{V}_i^{-1}\frac{\partial \mathbf{V}_i}{\partial \lambda_{j_2 k_2}}\mathbf{V}_i^{-1}\frac{\partial \mathbf{V}_i}{\partial \lambda_{j_1 k_1}}\right) \end{aligned} \tag{159}\]

\[ \begin{aligned} \frac{ \partial^2 \mathbf{V}_i }{ \partial \lambda_{j_2 k_2}\partial \lambda_{j_1 k_1} } &= \mathbf{H}_i\frac{ \partial^2 \boldsymbol{\Sigma} }{ \partial \lambda_{j_2 k_2}\partial \lambda_{j_1 k_1} }\mathbf{H}'_i \\ &= \mathbf{H}_i\frac{ \partial^2 \boldsymbol{\Lambda}\boldsymbol{\Lambda}' }{ \partial \lambda_{j_2 k_2}\partial \lambda_{j_1 k_1} }\mathbf{H}'_i \end{aligned} \tag{160}\]

\[ \begin{aligned} \frac{ \partial^2 \boldsymbol{\Sigma} }{ \partial \lambda_{j_2 k_2}\partial \lambda_{j_1 k_1} } & = \frac{ \partial^2 \boldsymbol{\Lambda}\boldsymbol{\Lambda}' }{ \partial \lambda_{j_2 k_2}\partial \lambda_{j_1 k_1} } = \frac{ \partial \boldsymbol{\Lambda} }{ \partial \lambda_{j_1 k_1} }\frac{ \partial \boldsymbol{\Lambda}' }{ \partial \lambda_{j_2 k_2} } + \frac{ \partial \boldsymbol{\Lambda} }{ \partial \lambda_{j_2 k_2} }\frac{ \partial \boldsymbol{\Lambda}' }{ \partial \lambda_{j_1 k_1} } = \mathbf{I}_{j_1}\mathbf{I}'_{k_1}\mathbf{I}_{k_2}\mathbf{I}'_{j_2} + \mathbf{I}_{j_2}\mathbf{I}'_{k_2}\mathbf{I}_{k_1}\mathbf{I}'_{j_1} = \\ &\quad\{\delta(j, j_1)\delta(k_1, k_2)\delta(j_2, k) + \delta(j, j_2)\delta(k_2, k_1)\delta(j_1, k)\ \text{for}\ j = 1\ \text{to}\ n,\ k = 1\ \text{to}\ n\} \end{aligned} \tag{161}\]

Since the Sigma parameters are only in the data likelihood portion of the conditional density, then

\[ \begin{aligned} \frac{ \partial -\log(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\not\mu}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\theta}_{\not\mu} } &= \frac{ \partial -\log(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \lambda_{j_1 k_1}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})) }{ \partial \lambda_{j_1 k_1} } \\ &= \frac{ \partial -\log(l(\mathbf{y}_i, \boldsymbol{\phi} \mid \lambda_{j_1 k_1}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})) }{ \partial \lambda_{j_1 k_1} } \end{aligned} \tag{162}\]

which is used in equation 23. The Sigma-like theta parameters cannot be processed in this way, because the user defines sigma-like parameters in $\mathbf{H}$, with unpredictable functional relationships to that theta.

Appendix G: Degrees of Freedom Assessment for `OMEGA` Priors

The heuristic justification for Mats Karlsson’s formula:

\[ N = 2\left(\frac{\boldsymbol{\Omega}}{E_{\boldsymbol{\Omega}}}\right)^{2} \tag{163}\]

Where $N$ is the effective number of subjects of the previous nlme (maximum likelihood) analysis, $\boldsymbol{\Omega}$ is an omega diagonal element, and $E_{\boldsymbol{\Omega}}$ is its standard error (the error in the estimate of $\boldsymbol{\Omega}$), is as follows.

For a normally distributed random variable $x$, with mean 0, and variance $\boldsymbol{\Omega}$, the following holds:

\[ \begin{aligned} \int x p(x \mid 0, \boldsymbol{\Omega}) dx &= \bar{x} \\ &= 0 \end{aligned} \tag{164}\]

Define the random variable $y$:

\[ y = (x - \bar{x})^2 \tag{165}\]

\[ \begin{aligned} \int (x - \bar{x})^2 p(x \mid 0, \boldsymbol{\Omega}) dx &= \int y p(x \mid 0, \boldsymbol{\Omega}) dx \\ &= \bar{y} \\ &= \mathrm{Var}(x) \\ &= \boldsymbol{\Omega} \end{aligned} \tag{166}\]

Finally, the fourth central moment is (Triantafyllopoulos 2003):

\[ \begin{aligned} \int (x - \bar{x})^4 p(x \mid 0, \boldsymbol{\Omega}) dx &= \int x^4 p(x \mid 0, \boldsymbol{\Omega}) dx \\ &= 3\boldsymbol{\Omega}^2 \end{aligned} \tag{167}\]

Then,

\[ \begin{aligned} \mathrm{Var}(\boldsymbol{\Omega}) = \mathrm{STD}^2(\boldsymbol{\Omega}) &= \int (y - \bar{y})^2 p(x \mid 0, \boldsymbol{\Omega}) dx = \int y^2 p(x \mid 0, \boldsymbol{\Omega}) dx - \bar{y}^2 \\ &= \int (x - \bar{x})^4 p(x \mid 0, \boldsymbol{\Omega}) dx - \bar{y}^2 = 3\boldsymbol{\Omega}^2 - \boldsymbol{\Omega}^2 = 2\boldsymbol{\Omega}^2 \end{aligned} \tag{168}\]

For $N$ normal random deviates, the variance of the estimate of its average variance is

\[ \begin{aligned} \mathrm{Var}(\hat{\boldsymbol{\Omega}}_N) &= \mathrm{Var}(\boldsymbol{\Omega}) / N \\ &= \mathrm{SE}^2(\hat{\boldsymbol{\Omega}}_N) \\ &= E_{\boldsymbol{\Omega}}^2 \end{aligned} \tag{169}\]

Thus,

\[ E_{\boldsymbol{\Omega}}^2 = \frac{2}{N}\boldsymbol{\Omega}^2 \tag{170}\]

That is, the standard error of the variance is related by the above equation, as long as the $N$ items that contribute to its assessment are normally distributed. This is the best error in the inter-subject variance that can be expected in a set of parameters from subjects with rich data for each. In a population analysis, however, some subjects with few data points will not have

much information for their parameter. However, population analysis yields empirical standard errors of Omegas $E_{\boldsymbol{\Omega}}$, that properly reflect the total information available for the Omega. Thus, given $E_{\boldsymbol{\Omega}}$, the “effective N”, or degrees of freedom can be evaluated as:

\[ N = 2\left(\frac{\boldsymbol{\Omega}}{E_{\boldsymbol{\Omega}}}\right)^{2} \tag{171}\]

A more accurate justification on ascertaining the effective degrees of freedom from OMEGAS of a previous maximum likelihood/nlme analysis is as follows. The nlme process estimates variances via the Fisher Information matrix, which for omegas, for a normally distributed set of etas, is:

\[ \frac{\partial}{\partial \omega_{j_3 j_4}}\left(\frac{\partial(-LL)}{\partial \omega_{j_1 j_2}}\right) = \frac{\partial}{\partial \omega_{j_3 j_4}}\left(\frac{ \partial\left(\frac{1}{2}\mathrm{tr}(\boldsymbol{\Omega}^{-1}m\mathbf{R}) - \frac{m}{2}\log(|\boldsymbol{\Omega}^{-1}|)\right) }{ \partial \omega_{j_1 j_2} }\right) = \tag{172}\]

\[ = \frac{m}{2}\mathbf{I}'_{j_1}\boldsymbol{\Omega}^{-1}\mathbf{I}_{j_3 j_4}\boldsymbol{\Omega}^{-1}\mathbf{I}_{j_2} = \frac{m}{2}\omega^{-1}_{j_1 j_3}\omega^{-1}_{j_4 j_2} \tag{173}\]

where $m$ is the number of subjects and

\[ \mathbf{R} = \frac{1}{m}\sum_{i=1}^{m}\boldsymbol{\eta}_i\boldsymbol{\eta}'_i \tag{174}\]

The above equation is based on matrix algebra principles outlined in appendix A, and we have seen a version of this structure in appendix F. Here we will ignore the symmetry of $\boldsymbol{\Omega}$ for the moment, and suppose we have $p^2$ independent omega elements (where $p$ is the dimension of $\boldsymbol{\Omega}$), and seek to package their $p^2 \times p^2$ second partial derivative elements into an extended $p^2 \times p^2$ Kronecker product matrix, so that we can conveniently invert it later:

\[ \boldsymbol{\Omega}^{-1} \otimes \boldsymbol{\Omega}^{-1} = \begin{pmatrix} \omega^{-1}_{11}\omega^{-1}_{11} & \omega^{-1}_{11}\omega^{-1}_{12} & \cdots & \omega^{-1}_{11}\omega^{-1}_{1p} & \cdots & \cdots & \omega^{-1}_{1p}\omega^{-1}_{11} & \omega^{-1}_{1p}\omega^{-1}_{12} & \cdots & \omega^{-1}_{1p}\omega^{-1}_{1p} \\ \omega^{-1}_{11}\omega^{-1}_{21} & \omega^{-1}_{11}\omega^{-1}_{22} & \cdots & \omega^{-1}_{11}\omega^{-1}_{2p} & \cdots & \cdots & \omega^{-1}_{1p}\omega^{-1}_{21} & \omega^{-1}_{1p}\omega^{-1}_{22} & \cdots & \omega^{-1}_{1p}\omega^{-1}_{2p} \\ \vdots & \vdots & \ddots & \vdots & & & \vdots & \vdots & \ddots & \vdots \\ \omega^{-1}_{11}\omega^{-1}_{p1} & \omega^{-1}_{11}\omega^{-1}_{p2} & \cdots & \omega^{-1}_{11}\omega^{-1}_{pp} & \cdots & \cdots & \omega^{-1}_{1p}\omega^{-1}_{p1} & \omega^{-1}_{1p}\omega^{-1}_{p2} & \cdots & \omega^{-1}_{1p}\omega^{-1}_{pp} \\ \vdots & \vdots & & \vdots & \ddots & & \vdots & \vdots & & \vdots \\ \vdots & \vdots & & \vdots & & \ddots & \vdots & \vdots & & \vdots \\ \omega^{-1}_{p1}\omega^{-1}_{11} & \omega^{-1}_{p1}\omega^{-1}_{12} & \cdots & \omega^{-1}_{p1}\omega^{-1}_{1p} & \cdots & \cdots & \omega^{-1}_{pp}\omega^{-1}_{11} & \omega^{-1}_{pp}\omega^{-1}_{12} & \cdots & \omega^{-1}_{pp}\omega^{-1}_{1p} \\ \omega^{-1}_{p1}\omega^{-1}_{21} & \omega^{-1}_{p1}\omega^{-1}_{22} & \cdots & \omega^{-1}_{p1}\omega^{-1}_{2p} & \cdots & \cdots & \omega^{-1}_{pp}\omega^{-1}_{21} & \omega^{-1}_{pp}\omega^{-1}_{22} & \cdots & \omega^{-1}_{pp}\omega^{-1}_{2p} \\ \vdots & \vdots & \ddots & \vdots & & & \vdots & \vdots & \ddots & \vdots \\ \omega^{-1}_{p1}\omega^{-1}_{p1} & \omega^{-1}_{p1}\omega^{-1}_{p2} & \cdots & \omega^{-1}_{p1}\omega^{-1}_{pp} & \cdots & \cdots & \omega^{-1}_{pp}\omega^{-1}_{p1} & \omega^{-1}_{pp}\omega^{-1}_{p2} & \cdots & \omega^{-1}_{pp}\omega^{-1}_{pp} \end{pmatrix} \tag{175}\]

Based on Kronecker product principles, this is readily inverted to form the extended variance matrix:

\[ \begin{aligned} EV(\boldsymbol{\Omega}) &= \left(\frac{m}{2}\boldsymbol{\Omega}^{-1} \otimes \boldsymbol{\Omega}^{-1}\right)^{-1} \\ &= \frac{2}{m}(\boldsymbol{\Omega} \otimes \boldsymbol{\Omega}) \end{aligned} \tag{176}\]

Therefore, each “extended variance” element consists of:

\[ EV(\omega_{j_1 j_2}, \omega_{j_3 j_4}) = \frac{2}{m}\omega_{j_1 j_3}\omega_{j_4 j_2} \tag{177}\]

We wish now have the symmetricized variance version of it as follows, using the functional transformation for symmetricizing:

\[ \begin{aligned} V(\omega_{j_1 j_2}, \omega_{j_3 j_4}) &= V((\omega_{j_1 j_2} + \omega_{j_2 j_1})/2, (\omega_{j_3 j_4} + \omega_{j_4 j_3})/2) \\ &= \frac{2}{m}\left(EV(\omega_{j_1 j_2}, \omega_{j_3 j_4}) + EV(\omega_{j_2 j_1}, \omega_{j_3 j_4}) + EV(\omega_{j_1 j_2}, \omega_{j_4 j_3}) + EV(\omega_{j_2 j_1}, \omega_{j_4 j_3})\right)/4 \\ &= \frac{1}{2m}\left(\omega_{j_1 j_3}\omega_{j_4 j_2} + \omega_{j_2 j_3}\omega_{j_4 j_1} + \omega_{j_1 j_4}\omega_{j_3 j_2} + \omega_{j_2 j_4}\omega_{j_3 j_1}\right) \\ &= \frac{1}{m}\left(\omega_{j_1 j_3}\omega_{j_4 j_2} + \omega_{j_2 j_3}\omega_{j_4 j_1}\right) \end{aligned} \tag{178}\]

where the simplification of the last step is due to Omega actually already being symmetrical.

Some examples:

\[ V(\omega_{ii}, \omega_{ii}) = 2\omega_{ii}\omega_{ii}/m \tag{179}\]

\[ V(\omega_{ik}, \omega_{ik}) = (\omega_{ii}\omega_{kk} + \omega_{ik}\omega_{ik})/m \tag{180}\]

\[ V(\omega_{ij}, \omega_{kl}) = (\omega_{ik}\omega_{lj} + \omega_{jk}\omega_{li})/m \tag{181}\]

And the argument can be made, as before, that the effective degrees of freedom can be extracted from the Omega estimates of an nlme estimation, and their empirical, Fisher information derived variance or SE as:

\[ N = \frac{2(\omega_{11})^2}{V(\omega_{11}, \omega_{11})} \tag{182}\]

When the previous analysis is MCMC Bayesian, the justification can be derived from the known mean and variance formulas for an inverse Wishart distribution of Omega. According to Wikipedia (with source references to (Haff 1979; Press 1979)), the mean of an inverse Wishart distributed matrix of dimension $p$ is:

\[ \mathbf{Q} = \boldsymbol{\Psi}/(\nu - p - 1) \tag{183}\]

Where $\boldsymbol{\Psi}$ is the effective scale and $\nu$ is the effective degrees of freedom of the OMEGA previous analysis. NONMEM multiplies inputted $OMEGAP values $\mathbf{R}$ with the inputted $OMPEGAPD degrees of freedom to serve as the scale, so we can express this scale as

\[ \boldsymbol{\Psi} = \nu\mathbf{R} \tag{184}\]

Thus $\mathbf{R}$ represents the scale per degree of freedom. Then the mean Omega would be expressed as:

\[ \mathbf{Q} = \nu\mathbf{R}/(\nu - p - 1) \tag{185}\]

The variance for a diagonal element is:

\[ \begin{aligned} w_{ii}^2 &= 2(\psi_{ii})^2/(\nu - p - 1)^2/(\nu - p - 3) \\ &= 2(\nu q_{ii})^2/(\nu - p - 1)^2/(\nu - p - 3) \\ &= 2 q_{ii}^2/(\nu - p - 3) \end{aligned} \tag{186}\]

Therefore effective degrees of freedom from a previous analysis could be derived from its diagonal mean ($q$) and variance elements ($w$) of OMEGA:

\[ \nu = 2(q_{ii})^2/w_{ii} + p + 3 \tag{187}\]

and this $\nu$ would be inputted as the $OMEGAPD for the new analysis. The effective scale (per degree of freedom) inputted as $OMEGAP for the new analysis would be:

\[ \mathbf{R} = (\nu - p - 1)\mathbf{Q}/\nu \tag{188}\]

where again, $\mathbf{Q}$ is the Omega matrix of the previous analysis.

It is interesting that because Omega is inverse Wishart distributed, this means that the Omega inverse is Wishart distributed, and its mean is as follows (see (Haff 1979; Press 1979), and Wikipedia entry for Wishart distribution), with respect to the normalized scale $\mathbf{R}$:

\[ \begin{aligned} \mathbf{P} &= \frac{1}{n}\sum_{i=1}^{n} \boldsymbol{\Omega}_i^{-1} \\ &= \mathbf{R}^{-1} \end{aligned} \tag{189}\]

where $\boldsymbol{\Omega}_i^{-1}$ is the inverse of the omega matrix obtained from a sample (for a total of $n$ Bayesian samples) from the stationary distribution of the MCMC Bayesian analysis (non-negative iterations in the .ext file), and $\mathbf{P}$ is the average of these Omega inverse elements. The variance of a diagonal element of $\mathbf{P}$ is:

\[ z_{ii} = 2(p_{ii})^2/\nu \]

Therefore, the effective degrees of freedom to be inputted into $OMEGAPD would be

\[ \nu = 2(p_{ii})^2/z_{ii} \tag{190}\]

and the $OMEGAP input would be the inverse of the average omega inverse:

\[ \mathbf{R} = (\mathbf{P})^{-1} \tag{191}\]

Thus, the inverse of the average omega inverse (harmonic mean of omega) from a Bayesian analysis will be similar to the mode Omega estimate from an nlme (FOCE/EM) analysis. The utility priorget offers convenient transformations of estimates from the previous analysis into properly packaged prior information in the form of ready to use control stream templates.

Appendix H: Technical Note on NonParametric Analysis

Perform a standard FOCE analysis, to produce vectors of $\boldsymbol{\eta}_i$, $i=1$ to $N$, at the mode a posteriori (MAP estimates, or empirical Bayes estimates (EBE)) for each subject $i$, evaluated at the final population parameters $(\boldsymbol{\theta}, \boldsymbol{\Omega}, \boldsymbol{\Sigma})$, where $N$ is the number of subjects, and $N_s$ is the number of support points. These best fit etas for each subject serve as the anchors, or grid points, for the non-parametric analysis, to be evaluated by subroutine NP.

In the subroutine NP, using data of subject $i$, and grid point $k$ (which may have come from EBE of subject $k$ of the FOCE analysis for $k \leq N$, or random creation of extra support points, for $N < k \leq N_s$), the data likelihood is evaluated:

\[ l(\mathbf{y}_i, \boldsymbol{\eta}_k, \boldsymbol{\Sigma}) \tag{192}\]

by subroutine OBJ3, and an initial prior (population density) is evaluated as

\[ \pi^0(\boldsymbol{\eta}_k) = \exp\left(-\frac{1}{2}\boldsymbol{\eta}'_k\boldsymbol{\Omega}\boldsymbol{\eta}_k\right) \tag{193}\]

$\pi(\boldsymbol{\eta}_k)$ is stored in vnonpara(1), and file system of subroutine DAT8.

Like other optimization methods, non-parametric estimation is reiterated, until the objective function no longer changes. At any given iteration of the non-parametric optimization, the following is assessed.

For each subject $i$ with $\mathbf{y}_i$, the set of $\boldsymbol{\eta}$ that yields the largest posterior density is evaluated,

\[ l(\mathbf{y}_i, \boldsymbol{\eta}_{m_i}, \boldsymbol{\Sigma})\pi(\boldsymbol{\eta}_{m_i}) \geq l(\mathbf{y}_i, \boldsymbol{\eta}_k, \boldsymbol{\Sigma})\pi(\boldsymbol{\eta}_k)\ \text{for all}\ k = 1, N_s \tag{194}\]

where $m_i$ produces the largest value of the posterior density for subject $i$. $l(\mathbf{y}_i, \boldsymbol{\eta}_k, \boldsymbol{\Sigma})$ is stored in VNONPARA(3). The $m_i$ is stored in IC(I), and $l(\mathbf{y}_i, \boldsymbol{\eta}_{m_i}, \boldsymbol{\Sigma})\pi(\boldsymbol{\eta}_{m_i})$ is stored in X79(I). The final $m_i$ is stored in VNONPARA(5).

The objective function is evaluated for a given iteration as

\[ O = -2\sum_{i=1}^{N}\log\left(\sum_{j=1}^{N_s} l(\mathbf{y}_i, \boldsymbol{\eta}_j, \boldsymbol{\Sigma})\pi(\boldsymbol{\eta}_j)\right) \tag{195}\]

where summing the probability density over all discrete positions $\boldsymbol{\eta}_j$, $j = 1$ to $N_s$, is the non-parametric or discrete density equivalent to integrating over all $\boldsymbol{\eta}$ for a continuous density function, to obtain a marginal density for each subject $i$. These marginal densities are in turn multiplied among all subjects $i$ to $N$, to obtain the joint marginal density. Since the objective function is typically $-2\log$(joint marginal density), it is more convenient to sum the log of the marginal densities among all subjects $i$, $i = 1$ to $N$, as shown in the above equation. $O$ is stored in OBJNP, and

\[ \sum_{j=1}^{N_s} l(\mathbf{y}_i, \boldsymbol{\eta}_j, \boldsymbol{\Sigma})\pi(\boldsymbol{\eta}_j) \tag{196}\]

is temporarily stored in U(I).

Normalized posterior densities are also evaluated:

\[ p_i(\boldsymbol{\eta}_k) = \frac{ l(\mathbf{y}_i, \boldsymbol{\eta}_k, \boldsymbol{\Sigma})\pi(\boldsymbol{\eta}_k) }{ \sum_{j=1}^{N_s} l(\mathbf{y}_i, \boldsymbol{\eta}_j, \boldsymbol{\Sigma})\pi(\boldsymbol{\eta}_j) } \tag{197}\]

The term

\[ \sum_{i=1}^{N}\sum_{j=1}^{N_s} \frac{ l(\mathbf{y}_i, \boldsymbol{\eta}_k, \boldsymbol{\Sigma}) }{ l(\mathbf{y}_i, \boldsymbol{\eta}_j, \boldsymbol{\Sigma})\pi(\boldsymbol{\eta}_j) } \]

is stored in VNONPARA(2).

The posterior densities are normalized in the sense that

\[ \sum_{k=1}^{N_s} p_i(\boldsymbol{\eta}_k) = 1 \tag{198}\]

as required for a proper probability density of $\boldsymbol{\eta}$. The final $p_i(\boldsymbol{\eta}_k)$ are stored in row subject $i$, column IPROB(K), of the .npi file. These are averaged among all subjects at a given anchor point $\boldsymbol{\eta}_k$ to obtain a posterior, or empirical, assessed “weight” at that anchor:

\[ p(\boldsymbol{\eta}_k) = \frac{1}{N}\sum_{i=1}^{N} p_i(\boldsymbol{\eta}_k) \tag{199}\]

The final values are reported as subject 0, IPROB(K), in the .npi file. If the following test is satisfied:

\[ \delta_k = \left|\frac{ p(\boldsymbol{\eta}_k) - \pi(\boldsymbol{\eta}_k) }{ \pi(\boldsymbol{\eta}_k) }\right| \geq \varepsilon_2 \tag{200}\]

for some small optimization criterion $\varepsilon_2$, then $p(\boldsymbol{\eta}_k)$ serves as the new prior density for the next iteration (this is the default expectation-maximization update method):

\[ \pi(\boldsymbol{\eta}_k) = p(\boldsymbol{\eta}_k) \tag{201}\]

If the test fails, it means that $p(\boldsymbol{\eta}_k)$ is no longer changing sufficiently with respect to its previous value $\pi(\boldsymbol{\eta}_k)$, and the updates no longer need to be performed for future iterations, for that $\boldsymbol{\eta}_k$.

When for all $k$ the following is satisfied:

\[ \delta_k \leq \varepsilon_1\ \text{for}\ k = 1\ \text{to}\ N \tag{202}\]

for some small $\varepsilon_1$, then the non-parametric optimization is ended.

The following final information is stored:

$\pi(\boldsymbol{\eta}_k)$ is stored in VNONPARA(1), which is retrieved for each $k$ via sequential calls to subroutine DAT8.

\[ \bar{\boldsymbol{\eta}} = \sum_{k=1}^{N_s} \boldsymbol{\eta}_k\pi(\boldsymbol{\eta}_k) \tag{203}\]

is stored in EXNPETA(), EXETA().

\[ \bar{\boldsymbol{\Omega}} = \sum_{k=1}^{N_s} (\boldsymbol{\eta}_k - \bar{\boldsymbol{\eta}})(\boldsymbol{\eta}_k - \bar{\boldsymbol{\eta}})'\pi(\boldsymbol{\eta}_k) \tag{204}\]

is stored in COVNPETA(), COVETA().

The expected values EXNPETA(j), and expected covariances COVNPETA(j,k) are reported as ETA(j) and ETC(j,k), respectively, in the .npe file.

If INPETA$\neq 0$, then

$m_i$ is stored in VNONPARA(2) to VNONPARA(1+neta), where neta is the eta vector length, retrieved for each subject $i$ via sequential calls to subroutine DAT8, and also placed in the .npd file, labeled as ETM(j), pertaining to eta(j). That is, the grid point eta vector that best fits subject $i$ is stored in record $i$ of the DAT8 storage system (the entire DAT8 storage contains $N_s$ records, where $N_s \geq N$).

If INPETA=0, then cumulative distributions are stored as follows:

\[ C(\boldsymbol{\eta}_{k(j)}) = \sum_{i : \eta_{i(j)} \leq \eta_{k(j)}}^{N_s} \pi(\boldsymbol{\eta}_i) \tag{205}\]

is stored in VNONPARA(1+j), retrieved for each support point $k$ sequentially from subroutine DAT8, and also placed in the .npd file, labeled as CUM(j). Here, $\eta_{i(j)}$ is the $j$th element of the eta vector belonging to support point $i$. That is, $C(\boldsymbol{\eta}_{k(j)})$ is the sum of densities $\pi(\boldsymbol{\eta}_i)$ for which the $j$th element of eta is less than or equal to the $j$th element of $\boldsymbol{\eta}_k$ (which is $\eta_{k(j)}$). In turn, $\eta_{k(j)}$ is stored in VETA(j), retrieved for each support point $k$ from subroutine DAT4, and reported as eta(j) in the .npd file. The $\pi(\boldsymbol{\eta}_k)$ for support point $k$ is reported in the .npd file in the last column, labeled as PROBABILITY.

See references (Savic and Karlsson 2009; Baverel et al. 2011) for information of supplementary support points and bootstrapping.

Appendix I: Note on `TNPRI`

The statistical basis of the frequentist method (TNPRI which stands for Total Normal PRIor) for priors is that of sampling about some mean, with some measure of dispersion, but not requiring a rigid rule of a particular distribution, other than that it has some semblance of normal distribution where this makes sense.

Regarding THETAS:

Let THETA be the theta estimate, and SE be the standard error of the theta estimate.

When no boundaries are given in the $THETA record, then $SIML will generate random sample thetas that are normally distributed with mean THETA and variance SE$^2$, and with suitable correction for correlation between theta(1) and theta(2), etc., in accordance with the variance-covariance matrix of the estimates.

When boundary is set, an intermediate variable normal deviate $v$ will be generated with mean $\log(\mathrm{THETA})$, and variance $(\mathrm{SE}/\mathrm{THETA})^2$. Again, the random deviate $v$ for each theta(1), theta(2), etc., is corrected to account for correlation (covariance of estimates) between theta(1) and theta(2), etc.

This log-normal deviate $v$ is then transformed to a final theta sample $w$ as follows:

Lower bound only:

w = exp(v) + LB

(range of $v$(-inf,+inf) transposes to range of $w$=(LB,+inf))

Upper bound only:

w = UB - exp(-v)

(range of $v$(-inf,+inf) transposes to range of $w$=(-inf,UB))

Lower and upper bound:

w = LB + (UB - LB) * exp(v) / (1 + exp(v))

(range of $v$(-inf,+inf) transposes to range of $w$=(LB,UB))

When SE/mean is small, then the resulting distribution of the thetas is nearly normally distributed with mean THETA and standard error SE, as reflected from the original estimates and standard errors. When SE/mean is large, this creates considerable non-normal distribution in the samples.

Regarding OMEGAS and Sigmas:

In the case of the OMEGAS, using the mean OMEGA and its standard error from a previous problem, the TNPRI method transposes this into its cholesky form and its equivalent standard error, so that each element in the cholesky matrix has the appropriate “mean and standard error”. Transforming the mean omega elements and their standard errors into the equivalent mean and standard errors for the cholesky elements is not trivial, but it can be done using matrix algebra and the principal of propagation of errors. The principle of propagation of errors itself is accurate only as an asymptotic rule, that is, if the error is sufficiently small, then it provides reasonable results.

Also, while having omega elements be normally distributed is not reasonable, the underlying cholesky elements can be modeled as normally distributed. Further, the off-diagonal elements are allowed to be positive or negative, and these elements are sampled as a strict normal distribution. The diagonal elements of the cholesky matrix can also be normally sampled, but with the proviso that the sample be positive. Thus, as with thetas, if the standard error is small relative to the mean value, very few if any samples will be negative. If the standard error is large, then a great many samples will end up as negative. To avoid negative values altogether, the additional transformation that is done for the cholesky diagonal is to generate a random sample as v-log(cholesky estimate)+(se of cholesky estimate)/(cholesky estimate)r, where r is created as a -N(0,1) deviate, and then exponentiate v, and this is the cholesky diagonal sample. When se/mean is small, exp(log(mean)+SE/meanr) is very nearly normally distributed with the appropriate mean and standard error.

Next, the cholesky matrix is multiplied by its transpose to create the random OMEGA sample, with the result of being positive definite, and having the appropriate mean and dispersion.

Appendix J: T Distribution Sample Generation

The tdist6 and tdist7 examples described in intro7.pdf use the fact that two normal random deviates can be converted into a T distributed normal random deviate. The derivation is as follows.

From Numerical Recipes (Press et al. 1992), using the Box-Muller algorithms for creating a random normal deviate pair is as follows. Two uniform random deviate pairs $u$ and $v$ are generated, and modified:

\[ U = 2u - 1 \tag{206}\]

\[ V = 2v - 1 \tag{207}\]

\[ W = U^2 + V^2 \tag{208}\]

Values of W>=1.0 are rejected.

Normal deviates $x$ and $y$ are generated as follows:

\[ x = U\sqrt{-2\log(W)/W} \tag{209}\]

\[ y = V\sqrt{-2\log(W)/W} \tag{210}\]

According to (Shaw 2006), a t-distribution sample of $n$ degrees of freedom can be generated as

\[ z = \frac{U}{\sqrt{W}}\sqrt{n\left(u^{-2/n} - 1\right)} \tag{211}\]

This suggests that rather than starting with uniform deviates $u$ and $v$, one can use normal deviates $x$ and $y$ to generate the t deviate:

\[ x^2 + y^2 = -2\log(W) \tag{212}\]

\[ W = \exp\left(-\frac{1}{2}(x^2 + y^2)\right) \tag{213}\]

\[ \begin{aligned} \frac{x}{\sqrt{x^2 + y^2}} &= \frac{U\sqrt{-2\log(W)/W}}{\sqrt{-2\log(W)}} \\ &= \frac{U}{\sqrt{W}} \end{aligned} \tag{214}\]

So,

\[ z = \frac{x}{\sqrt{x^2 + y^2}}\sqrt{n\left(\exp((x^2 + y^2)/2) - 1.0\right)} \tag{215}\]

This transformation is used in the tdist6 and tdist7 examples.

From the alternative method of the Box-Muller method ((Press et al. 1992)) of using trigonometric functions, we note that uniform deviates $u$ and $v$ can be generated from the normal deviates:

\[ u = \exp\left(-\frac{1}{2}(x^2 + y^2)\right) \tag{216}\]

\[ \begin{aligned} v &= \frac{1}{2\pi}\tan^{-1}\left(\frac{y}{x}\right) \\ &= \frac{1}{2\pi}\cos^{-1}\left(\frac{x}{\sqrt{x^2 + y^2}}\right) \end{aligned} \tag{217}\]

Therefore,

\[ t = \cos(2\pi v)\sqrt{n\left(u^{-2/n} - 1.0\right)} \tag{218}\]

which is a means of generating a t sample using a uniform deviate pair, without having to reject uniform random samples, and is the method used for t sample generating function TDEV2 in ..\source\GENERAL.f90, for general degrees of freedom $n$.

Appendix K: Transformation of Parameters during Classical NONMEM Estimation

During classical NONMEM estimation, the thetas, omegas, and sigmas are first transformed into “unconstrained” parameters, domain of which is from $-\infty$ to $+\infty$, regardless of the boundaries on the original parameters. This allows movement in the unconstrained domain by the estimation process.

In the following discussion, x=parameter in the user domain, l=lower bound in user domain, w=upper bound in u=unconstrained parameter, u=unconstrained parameter, and s=scaling.

Scaling is evaluated using the initial values of $x$ at the start of the estimation, and remains unchanged throughout the estimation. The $u$ parameters are varied during the estimation, so new values of $x$ are back-calculated from the fixed $s$ and varying $u$.

Thetas:

No lower or upper bound:

\[ u = \mathrm{sign}(x_{\mathrm{initial}}) \cdot 0.1 \tag{219}\]

\[ x = |x_{\mathrm{initial}}/u_{\mathrm{initial}}| \tag{220}\]

Lower bound only:

\[ u = 0.1 \tag{221}\]

\[ s = (x_{\mathrm{initial}} - l)/\exp(u_{\mathrm{initial}}) \tag{222}\]

Upper bound only:

\[ u = 0.1 \tag{223}\]

\[ s = (w - x_{\mathrm{initial}}) \cdot \exp(u_{\mathrm{initial}}) \tag{224}\]

Upper and lower bounds:

\[ u = 0.1 \tag{225}\]

\[ s = \log(x_{\mathrm{initial}} - l) - \log(w - x_{\mathrm{initial}}) - u_{\mathrm{initial}} \tag{226}\]

The Omegas and Sigmas are transformed as follows

\[ \mathbf{C} = \mathrm{Cholesky}(\mathbf{X}) \tag{227}\]

Where

\[ \mathbf{X} = \mathbf{C}\mathbf{C}^T \tag{228}\]

And $\mathbf{C}$ is lower triangular.

Then:

Diagonal elements:

\[ u_{ii} = 0.1 \tag{229}\]

\[ s_{ii} = c_{ii\mathrm{initial}}/\exp(u_{ii\mathrm{initial}}) \tag{230}\]

Off-diagonal elements:

\[ u = \mathrm{sign}(c_{ij\mathrm{initial}}) \cdot 0.1 \tag{231}\]

\[ s_{ij} = |c_{ij\mathrm{initial}}/u_{\mathrm{initial}}| \tag{232}\]

Appendix L: Estimating Super ID Level Etas

Two cycles of MAP or Monte Carlo estimation are evaluated for each total objective function call. In the first cycle, all etas are varied and estimated, in the usual manner, without accounting for $LEVEL information. Then, summary statistics for etas are obtained by evaluating mean etas among individuals belonging to a particular group of a particular level ($LEVEL information is used). For the second cycle, these mean etas calculated from the first cycle are then placed in the higher (super-individual) eta positions and fixed there, while the individual level etas alone are allowed to be newly estimated.

Let’s consider the FOCE method. For the first cycle, all MAP etas are estimated in the usual manner using the following objective function for FOCE, for all etas. For a given individual with K levels (including the individual subject level 1):

\[ L_{\mathrm{indv}} = -\log(l(\mathbf{y}_i \mid \boldsymbol{\eta}_{k(i)}, \boldsymbol{\theta})) + \sum_{k=1}^{K}\left(\frac{1}{2}\log(\det(\boldsymbol{\Omega}_k)) + \frac{1}{2}\boldsymbol{\eta}'_k(\boldsymbol{\theta})\boldsymbol{\Omega}_k^{-1}\boldsymbol{\eta}_k(\boldsymbol{\theta})\right) + \frac{1}{2}\log(\det(\boldsymbol{\Omega}_1^* + \mathbf{S}_i^*(\mathbf{y}_i \mid \boldsymbol{\eta}_{k(i)}, \boldsymbol{\theta}))) \tag{233}\]

Where B stands for block:

\[ \hat{\boldsymbol{\eta}}_B = \begin{pmatrix} \hat{\boldsymbol{\eta}}_1 \\ \hat{\boldsymbol{\eta}}_2 \\ \vdots \\ \hat{\boldsymbol{\eta}}_K \end{pmatrix} \tag{234}\]

\[ \boldsymbol{\Omega}_B^{-1} = \begin{pmatrix} \boldsymbol{\Omega}_1^{-1} & 0 & 0 & 0 \\ 0 & \boldsymbol{\Omega}_2^{-1} & 0 & 0 \\ 0 & 0 & \boldsymbol{\Omega}_3^{-1} & 0 \\ 0 & 0 & 0 & \boldsymbol{\Omega}_K^{-1} \end{pmatrix} \tag{235}\]

\[ \mathbf{S}_B^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) = \begin{pmatrix} \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) & \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) & \cdots & \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) \\ \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) & \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) & \cdots & \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) & \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) & \cdots & \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) \end{pmatrix} \tag{236}\]

since

\[ \begin{aligned} \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) &= \mathbf{S}^*_{ll'}(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) \\ &= \left(\frac{ \partial f_i(\boldsymbol{\eta}_B, \boldsymbol{\theta}) }{ \partial \boldsymbol{\eta}_B }\right)_{\hat{\boldsymbol{\eta}}_B,\boldsymbol{\theta}} \boldsymbol{V}^{-1}\left(\frac{ \partial f_i(\boldsymbol{\eta}_B, \boldsymbol{\theta}) }{ \partial \boldsymbol{\eta}_B }\right)'_{\hat{\boldsymbol{\eta}}_B,\boldsymbol{\theta}} \end{aligned} \tag{237}\]

because each element of $\boldsymbol{\eta}_B$ has the same position in the prediction function $f$ of its corresponding element in every level, $k=1$ to $K$. This is effectively optimizing without the $LEVEL information.

To prepare for the second cycle, the mean eta vector is evaluated, starting with level K, and going down to level 1. Let

\[ \bar{\boldsymbol{\eta}}_{k,i} = \sum_{j \in k_{(l)(k,i)}} \hat{\boldsymbol{\eta}}_j \tag{238}\]

be the total eta vector for subject $i$ summed across all levels. Let also:

\[ \bar{\boldsymbol{\eta}}_{k(i)} = \frac{ \sum_{l \in l(k,i)} \hat{\boldsymbol{\eta}}_l }{ \sum_{l \in l(k,i)} 1 } \tag{239}\]

Then

\[ \hat{\boldsymbol{\eta}}_{(B),i} = \bar{\boldsymbol{\eta}}_{k,i} - \bar{\boldsymbol{\eta}}_{k+1(i)} \tag{240}\]

where $k_{(l)(k,i)}$ is the $j$th group of level $k$ to which subject $i$ belongs, and $(k+1)(k,i)$ is the $k+1$ th group of level $k+1$ to which subject $i$ belongs. The subject subscript $i$ may boomitted for simpler notation:

\[ \hat{\boldsymbol{\eta}}_{k(l)} = \bar{\boldsymbol{\eta}}_{k,i} - \bar{\boldsymbol{\eta}}_{k+1(i)} \tag{241}\]

These are then used to evaluate the OFV of the second cycle (and evaluating eta hat(2) to the individual level):

\[ \begin{aligned} L_{\mathrm{indv},2} &= -\log(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta})) \\ &\quad + \frac{1}{2}\hat{\boldsymbol{\eta}}'_1\boldsymbol{\Omega}_1^{-1}\hat{\boldsymbol{\eta}}_1 + \frac{1}{2}\hat{\boldsymbol{\eta}}'_2\boldsymbol{\Omega}_2^{-1}\hat{\boldsymbol{\eta}}_2 + \frac{1}{2}\log(\det(\boldsymbol{\Omega}_1^* + \mathbf{S}_i^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}))) + \sum_{k=1}^{K}\left(\frac{1}{2}w_k\log(\det(\boldsymbol{\Omega}_k))\right) \\ &\quad + \frac{1}{2}w_k\hat{\boldsymbol{\eta}}'_k\boldsymbol{\Omega}_k^{-1}\hat{\boldsymbol{\eta}}_k \end{aligned} \tag{242}\]

Where subscript (2) refers to cycle 2 optimization, and $w_k$ is a weight parameter. If LEVWT=0 (default), then the weight is such that when the contribution of this term from all subjects is summed, the total OFV will have a weight of 1 for each group in level $k$. If LEVWT=1, then the summed weighting among all subjects, in the total OFV, will be in proportion to their number of LEVCENTER values, $w_k$ is always equal to 1. When the weight is not modeled, or the default expectation-maximization update method is used (individual subject), $w_k$ is always equal to 1. The weight for the higher level) may differ slightly from $w_k$ to account for degrees of freedom corrections, if LEVCENTER=1 (see loss of degrees of freedom table and w and v table at end of this section).

Experience has shown that optimizing the cycle 2 objective function in FOCE (LEVOBJTYPE=1) results in a less accurate assessment of the higher level (k>1) omegas. Instead, the OFV that is optimized and used for advancement of parameters is the average of OFVs to cycle 1 and 2 (LEVOBJTYPE=0, default), providing less biased results:

\[ \begin{aligned} (L_{\mathrm{indv},1} + L_{\mathrm{indv},2})/2 &= -\log(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta})) + \log(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_{B(2)}, \boldsymbol{\theta}))/2 \\ &\quad + \frac{1}{4}\hat{\boldsymbol{\eta}}'_1\boldsymbol{\Omega}_1^{-1}\hat{\boldsymbol{\eta}}_1 + \frac{1}{4}\hat{\boldsymbol{\eta}}'_2\boldsymbol{\Omega}_2^{-1}\hat{\boldsymbol{\eta}}_2 + \frac{1}{4}\log(\det(\boldsymbol{\Omega}_1^* + \mathbf{S}_i^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}))) \\ &\quad + \sum_{k=1}^{K}\left(\frac{1}{2}w_k\log(\det(\boldsymbol{\Omega}_k)) + \frac{1}{2}w_k\hat{\boldsymbol{\eta}}'_k\boldsymbol{\Omega}_k^{-1}\hat{\boldsymbol{\eta}}_k\right) \end{aligned} \tag{243}\]

where (1) portion of subscript refers to cycle 1 optimization.

To show that using the average objective function of cycle 1 and cycle 2 does not lead to a biased estimate, that is, the two cycles are consistent with each other, consider a two-level problem. Then,

\[ \begin{aligned} L_{\mathrm{indv}, w=1} &= -\log(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta})) + \frac{1}{2}\hat{\boldsymbol{\eta}}'_1\boldsymbol{\Omega}_1^{-1}\hat{\boldsymbol{\eta}}_1 + \frac{1}{2}\hat{\boldsymbol{\eta}}'_2\boldsymbol{\Omega}_2^{-1}\hat{\boldsymbol{\eta}}_2 \\ &\quad + \frac{1}{2}\log(\det(\boldsymbol{\Omega}_1)) + \frac{1}{2}\log(\det(\boldsymbol{\Omega}_2)) + \frac{1}{2}\log\left(\det\begin{pmatrix} \boldsymbol{\Omega}_1^{-1} + Z & Z \\ Z & \boldsymbol{\Omega}_2^{-1} + Z \end{pmatrix}\right) \end{aligned} \tag{244}\]

\[ Z = \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) \]

Then,

\[ \begin{aligned} &\log(\det(\boldsymbol{\Omega}_1)) + \log(\det(\boldsymbol{\Omega}_2)) + \log\left(\det\begin{pmatrix} \boldsymbol{\Omega}_1^{-1} + Z & Z \\ Z & \boldsymbol{\Omega}_2^{-1} + Z \end{pmatrix}\right) = \\ &\log\left(\det\begin{pmatrix} \boldsymbol{\Omega}_1 & 0 \\ 0 & \boldsymbol{\Omega}_2 \end{pmatrix}\begin{pmatrix} \boldsymbol{\Omega}_1^{-1} + Z & Z \\ Z & \boldsymbol{\Omega}_2^{-1} + Z \end{pmatrix}\right) = \log\left(\det\begin{pmatrix} \mathbf{I} + \boldsymbol{\Omega}_1 Z & \boldsymbol{\Omega}_1 Z \\ \boldsymbol{\Omega}_2 Z & \mathbf{I} + \boldsymbol{\Omega}_2 Z \end{pmatrix}\right) = \\ &\log(\det((\mathbf{I} + \boldsymbol{\Omega}_1 Z)(\mathbf{I} + \boldsymbol{\Omega}_2 Z) - \boldsymbol{\Omega}_1 Z\boldsymbol{\Omega}_2 Z)) = \log(\det(\mathbf{I} + (\boldsymbol{\Omega}_1 + \boldsymbol{\Omega}_2)Z)) = \\ &\log(\det(\mathbf{I} + \boldsymbol{\Omega} Z)) = \log(\det(\boldsymbol{\Omega}^{-1} + Z)) + \log(\det(\boldsymbol{\Omega})) \end{aligned} \tag{245}\]

Where we define

\[ \boldsymbol{\Omega} = \boldsymbol{\Omega}_1 + \boldsymbol{\Omega}_2 \tag{246}\]

The above derivation of equation 245 is based on the fact that for block matrices of equal dimensions,

\[ \left|\begin{matrix} A & B \\ C & D \end{matrix}\right| = |AD - BC| \]

if $CD = DC$ (ref. (Silvester 2000)).

We can similarly simplify the eta terms if we note that

\[ \begin{aligned} \hat{\boldsymbol{\eta}}_k & = \boldsymbol{\Omega}_k\frac{\partial f_i}{\partial \boldsymbol{\eta}_k}\boldsymbol{\Sigma}^{-1}(\mathbf{y}_i - f_i) \\ &= \boldsymbol{\Omega}_k \mathbf{b}_k \\ &= \boldsymbol{\Omega}_k \mathbf{b} \end{aligned} \tag{247}\]

since

\[ \mathbf{b}_k = \mathbf{b} \tag{248}\]

because as mentioned before, each element of $\hat{\boldsymbol{\eta}}_k$ has the same position in the prediction function $f$ of its corresponding element in every level, $k=1$ to $K$.

Then

\[ \boldsymbol{\Omega}_k^{-1}\hat{\boldsymbol{\eta}}_k = \mathbf{b} \tag{249}\]

\[ \begin{aligned} \hat{\boldsymbol{\eta}}'_1\boldsymbol{\Omega}_1^{-1}\hat{\boldsymbol{\eta}}_1 + \hat{\boldsymbol{\eta}}'_2\boldsymbol{\Omega}_2^{-1}\hat{\boldsymbol{\eta}}_2 &= (\hat{\boldsymbol{\eta}}'_1 + \hat{\boldsymbol{\eta}}'_2)\mathbf{b} \\ &= \hat{\boldsymbol{\eta}}'\mathbf{b} \end{aligned} \tag{250}\]

Where we define

\[ \hat{\boldsymbol{\eta}} = \hat{\boldsymbol{\eta}}_1 + \hat{\boldsymbol{\eta}}_2 \tag{251}\]

Note also,

\[ \begin{aligned} \partial^2\Omega &= (\boldsymbol{\eta}'_1\boldsymbol{\Omega}_1^{-1} + \boldsymbol{\eta}'_2\boldsymbol{\Omega}_2^{-1})(\hat{\boldsymbol{\eta}}_1 + \hat{\boldsymbol{\eta}}_2) \\ &= b'(\boldsymbol{\Omega}_1 + \boldsymbol{\Omega}_2)b \\ &= b'\boldsymbol{\Omega} b \\ &= b'(\boldsymbol{\Omega}_1 + \boldsymbol{\Omega}_2)b \end{aligned} \]

So, in general,

\[ \begin{aligned} L_{\mathrm{(1)}} &= -\log(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta})) + \sum_{k=1}^{K}\left(\frac{1}{2}\log(\det(\boldsymbol{\Omega}_k)) + \frac{1}{2}\hat{\boldsymbol{\eta}}'_k\boldsymbol{\Omega}_k^{-1}\hat{\boldsymbol{\eta}}_k\right) + \frac{1}{2}\log(\det(\boldsymbol{\Omega}_1^* + \mathbf{S}_i^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}))) \end{aligned} \tag{252}\]

\[ = -\log(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}, \boldsymbol{\theta})) + \frac{1}{2}\log(\det(\boldsymbol{\Omega})) + \frac{1}{2}\hat{\boldsymbol{\eta}}'\mathbf{b} + \frac{1}{2}\log(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}, \boldsymbol{\theta}))) \tag{253}\]

So, rephrasing equation 243 to the following equivalent:

\[ \begin{aligned} (L_{\mathrm{indv},1} + L_{\mathrm{indv},2})/2 &= -\log(l(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta})) + \frac{1}{2}\sum_{k=1}^{K}\left(\frac{1}{2}w_k\log(\det(\boldsymbol{\Omega}_k))\right) \\ &\quad + \frac{1}{4}\hat{\boldsymbol{\eta}}'_1\boldsymbol{\Omega}_1^{-1}\hat{\boldsymbol{\eta}}_1 + \frac{1}{4}\hat{\boldsymbol{\eta}}'_2\boldsymbol{\Omega}_2^{-1}\hat{\boldsymbol{\eta}}_2 + \sum_{k=1}^{K}\left(\frac{1}{2}w_k\hat{\boldsymbol{\eta}}'_k\boldsymbol{\Omega}_k^{-1}\hat{\boldsymbol{\eta}}_k\right) \\ &\quad + \frac{1}{4}\log(\det(\boldsymbol{\Omega}_1^* + \mathbf{S}_i^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}))) + \frac{1}{4}\log(\det(\boldsymbol{\Omega}^{-1} + \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}, \boldsymbol{\theta}))) \end{aligned} \tag{254}\]

(Note that $\mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}_B, \boldsymbol{\theta}) = \mathbf{S}^*(\mathbf{y}_i \mid \hat{\boldsymbol{\eta}}, \boldsymbol{\theta})$.)

Thus, the LN(1) portion contributes to assessing the total $\boldsymbol{\Omega}$ (and the total $\boldsymbol{\Omega}$, adding stability to the analysis without introducing bias/inconsistency (since the component level etas must sum to the total eta, and the component level omegas must sum to the total omega anyways)) to the $E$-STEP. The SCOV step is assessed using LN(2) only. When covariate modeling, the main reduction in the OFV will be in reducing $\hat{\boldsymbol{\eta}}_{2}$ and $\boldsymbol{\Omega}_2$, which of course will reduce the total $\boldsymbol{\Omega}$ to the same degree. This means both the LN(1) and LN(2) portions will be equally reduced, and there will be no loss in statistical power to detect a covariate effect when using equation 254.

using (LN(1)+LN(2))/2. In fact, the same reduction in LN(1) will also occur, so it may be actually simpler to use the analysis that is suggested by equation 253, that is, a model that only models the total omegas, and without need to use $LEVEL, as long as one knows that one is modeling, for example, inter-subject variance that can be expected in a set of parameters from subjects with rich data for each. In a population analysis, however, some subjects with few data points will not have properly packaged prior information in the form of ready to use control stream templates.

Extensive numerical verification tables follow (pp. 99–105) demonstrating that the total Omegas from the several Superid LEVEL implementations (LN(1), LN(2), (LN(1)+LN(2))/2, etc.) produce consistent total Omega estimates across FOCE, ITS, IMP, BAYES, and NUTS estimation methods with both LEVCENTER=0 and LEVCENTER=1 settings, and Superid2 (two super levels), Superid3 (one super level), and Superid5 (one super level with a very small Omega) configurations. The tables show that estimating via (LN(1)+LN(2))/2 or by LN(2) are approximately equivalent, as numerically confirmed, and that total Omega is reduced similarly in covariate search whether using $LEVEL or not.

The degrees of freedom corrections are as follows, where $L$=loss of degrees of freedom:

LEVCENTER	Analysis type	Level $k$=1 to K-1	Level K
0	EM/`BAYES`	L=0	L=1 by default, can be changed by user with `$PK *LDF=x` Set “LDF=0” for closer consistency with `FOCE`/Laplace
0	`FOCE`/Laplace	L=0	0
1	EM/`BAYES`	L=Number of groups at $k$=1	L=1 by default, can be changed by user with `$PK *LDF=x` Set “LDF=0” for closer consistency with `FOCE`/Laplace
1	`FOCE`/Laplace	L=Number of groups at $k$=1	0

Coefficients $w$ (on $\boldsymbol{\eta}'\boldsymbol{\Omega}$megalnv$\boldsymbol{\eta}$ term) and $v$ (on $\log(\det(\boldsymbol{\Omega}))$ term) (calculated for each subject):

	LEVWT=0	LEVWT=1
$w_k$	$1/N_g$	$N_g/N$
$v_k$	$(N_g-L)/N_g(1+N_g)$	$(N_g-L)/N_g(N_g/N)$

$N$=number of subjects in entire data set $N_g$=number of groups in level $k$ ($=N$ for level 1) $N_{\mathrm{kg}}$=number of individuals in group $j$ of level $k$ to which subject belongs ($=1$ for level 1) $L$=loss of degrees of freedom based on previous table.

Appendix M: Simple Distributed Delayed Output Problems and How They May Be Solved with Repetition Parameters

In guide VIII, the Repetition Parameters I example shows an example of how a convolution can be performed in NONMEM, using RPTO and RPTI variables, particularly the repeatI example. The following shows the basis for such a problem. Consider a problem in which a state variable N represents the number of cells (progenitor white blood cells, for example) in a particular stage of development, and the output rate is related to the input rate $k(t)$ by some distributed delay $l(t)$ (which represents a probability density) of the input, for example,

\[ \begin{aligned} \frac{dN(t)}{dt} &= k(t) - k(t) \ast L(t) \\ &= k(t) - \int_{0}^{t} l(t-s)k(s)ds \end{aligned} \tag{255}\]

That is, the output is delayed by an amount $t-s$, with possible values of $s$ can vary from 0 to infinity, distributed according to $l(s)$. If we allow the ODE solver to do all the work, it requires 2 independent variables of integration, $s$ and $t$ to solve for $N(t)$:

\[ \begin{aligned} N(t) &= \int_{0}^{t} dN(u) \\ &= \int_{0}^{t} k(u)du - \int_{0}^{t}\int_{0}^{u} l(u-s)k(s)dsdu \end{aligned} \tag{256}\]

which is difficult to perform. We can simplify the problem by carrying further the derivation as follows:

\[ N(t) = \int_{0}^{t} dN(u) = \int_{0}^{t} k(u)du - \int_{0}^{t}\int_{0}^{u} l(u-s)k(s)dsdu = \tag{257}\]

\[ [K(t) - K(-\infty)] - \int_{-\infty}^{t}\int_{-\infty}^{u} l(u-s)k(s)dsdu = \tag{258}\]

(where $K$=integral of $k$)

\[ [K(t) - K(-\infty)] - \left[\int_{-\infty}^{t} l(u-s)k(s)dsdu\right]L(u) = \int_{-\infty}^{t} L(u)\left(\int_{-\infty}^{u} \frac{dk(z-s)}{ds}dz\right)du = \int_{-\infty}^{t} L(u)dkdu = \tag{259}\]

(where $L$=integral of $l$, and using integration of parts)

\[ [K(t) - K(-\infty)] - \left[[K(-\infty) - K(-\infty)]L(u)\right]_{-\infty}^{t} - \int_{-\infty}^{t} [K(u-s) - K(-\infty)]L'(u)du = \int_{-\infty}^{t} L(u)dkdu = \tag{260}\]

\[ [K(t) - K(-\infty)] - [(K(-\infty) - K(-\infty))L(t) + \int_{-\infty}^{t} L(u)k(u)du] - \int_{-\infty}^{t} [K(u-s) - K(-\infty)]l(u)du = \tag{261}\]

\[ \int_{-\infty}^{t} [\bar{k} \cdot s(u)] du - \int_{-\infty}^{t} L(u)k(u)du - \int_{-\infty}^{t} [K(u-s) - K(-\infty)]l(u)du = \int_{-\infty}^{t} L(u)k(u)du \tag{262}\]

\[ \int_{0}^{t} (1 - L(s))k(s)ds \tag{263}\]

This last step is true because

\[ \begin{aligned} \int_{0}^{t} k(s)ds &= \int_{-\infty}^{t} k(s)ds \\ &= \int_{-\infty}^{t} k(s)ds \\ &= K(t) - K(-\infty) \end{aligned} \]

letting $z=t-s$, $dz=-ds$, and noting that

\[ \begin{aligned} L(0) &= 0;\ L(\infty) \\ &= 1 \end{aligned} \tag{264}\]

since $l$ is a probability density, and $L$ is its cumulative probability function.

Therefore, the state variable $N(t)$ can be solved as:

\[ N(t) = \int_{0}^{t} (1 - L(s))k(s)ds \tag{265}\]

Thus, by using the CDF $L$ instead of $l$, we remove one variable of integration. CDF functions to common distributions are generally available, including in the NONMEM library. We wish to re-arrange the variable of integration as follows, and we consider change in input rate after $t_0$:

\[ z=t-s \tag{266}\]

\[ k'(\tau) = k(\tau) - \bar{k}_0 \tag{267}\]

\[ k(t) = \bar{k}_0\ \tau < 0 \tag{268}\]

\[ \begin{aligned} N(t) &= \int_{0}^{t} (1 - L(s))k(s)ds + N_0 = -\int_{0}^{t} (1 - L(t-z))k(z)dz + N_0 = \\ &= \int_{0}^{t} (1 - L(t-z))k'(z)dz + N_0 \end{aligned} \tag{269}\]

\[ \int_{0}^{t} (1 - L(t-z))dz \cdot \bar{k}_0 + N_0 \tag{270}\]

We could then submit the above integration to the ODE solver. If we let $s$ be the variable of integration T, and let $t$ be the time TI of the present data record, then we could submit the problem as:

\[ \frac{N(t)}{dt} = k(t)(1 - L(TI - t))k(t) \tag{271}\]

In the repeat1 problem:

; This example is discussed in Repetition_I example (repetI.mak)
; in the online help directory.
; It demonstrates the use of RPTO variable

$PROB repeat1
$INPT ID TIME AMT DV EVID
$DATA repeat1.dat NULL=.
$SUB ADVAN6 TOL=5
$MODEL COMP=(DEPOT DEFOSDE) COMP=(CENTRAL DEFOBS)

$PK
IF (ICALL.EQ.0) RPTO=1  ;enables use of repetition feature
$PK
KA=THETA(1)*EXP(ETA(1))

K=THETA(2)*EXP(ETA(2))
C=THETA(3)
V=THETA(4)*EXP(ETA(4))
Z2=V
IF (RPTI.EQ.0) TI=TIME
IF (NEWIND.EQ.2) RPTO=1

$DES
DADT(1)=-KA*A(1)
D=EXP(-K*(TI-T)*C)
DADT(2)=KA*A(1)-D*D

$ERROR
Y=F+EPS(2)

$THETA 2 1 1 2
$OMEGA 1 2 (1 FIX) 1
$SIGMA 1
$TABLE TIME TI

Note that

\[ D=1-L(TI-T) \]

and that this is 1-WEIBULLCDF(1,1), while

\[ KA*A(1)=k(t) \]

In another problem, shown in intro7.pdf in the distributed delay section, the following problem is considered, based on (Koch and Schropp 2015), with a Michaelis-Menten elimination of drug from central compartment, while drug in central compartment stimulates cell production as an smax/sc50 manner, and with the output rate being the production rate Weibull distribution delayed. This can readily be solved by the RPTO/RPTI method, as shown (smax_weibull_rep.ctl):

$PROBLEM Weibull1 Distributed Delay on Smax Drug Responsive Input function
$ABBR DECLARE INTEGER I,RPTO=MODEL(COMPONENT 1: DEPOT)=<A1>
$INPUT ID AMT TIME DV EVID MDV
$DATA smax_weibull.dat IGNORE=C
$SUBROUTINES ADVAN13 TOL=6 ATOL=6 ; OTHER=ddmain.f90
$MODEL NCOMPARTMENTS=2

$PK
IF (ICALL.EQ.0) RPTO=1  ;enables use of repetition feature
$PK
CALL=-2
MAXIT=200000000
SMAX=THETA(1)
SCS=THETA(2)
KINS=THETA(3)
KA=THETA(4)
BB=THETA(5)

VMAX = 2
V50 = 3
V = 1
IF (RPTI.EQ.0) TI=TIME
IF (NEWIND.EQ.2) RPTO=1

; initial conditions
A_0(1)=0
A_0(2)=5.0

$DES
C = A(1)/V
kin = kins*(1.0+Smax*c/(SC50+c))
DD=D
D=(1-t-T).D=EXP(-(T-T-T/BB)**AA)  ; DD=1-WEIBULLCDF(T-T)
DADT(1) = -(Vmax*A(1))/(V50+A(1))
DADT(2) = DD*(KIN-KIN8) + A(1)*0(1)

$ERROR
X1=A(2)
IPRED=A(2)
Y=IPRED+IPRED*EPS(1)+EPS(1)*EPS(1)

$THETA
0.5
10.0
0.8
7.0
6.5

$OMEGA (0.0 fixed)x4
$SIGMA (0.0 fixed)
$SIML (122345) ONLYSIM SUBP=1
$TABLE TIME MDV EVID NOPRINT NOAPPEND FILE=smax_weibull_rep.tab
$NONHEADER FORMAT=S10E23.15

The RPTO/RPTI setup causes the integration to repeatedly integrate from T=0 to T=TIME of the present record for each record, thus performing a convolution of KIN-KIN8, with the integrated result A(2) being the total number of cells at time $t$, a result of an output rate that is equal to production rate Weibull distribution delayed. What is important to consider here is that this convolution process can be implemented only if the derivative of the state variable as a function of time can be represented in a k-k*l format, where k is the input or production rate. Any other form of distributed delay will not work, therefore, this convolution method using the repetition variables RPTO/RPTI will not work for general convolution problems. For example, suppose we have the following:

\[ DADT(1) = -KT \cdot (1 - L(TI-T)) \cdot k(t) / TDD \cdot KA \cdot A(1) \]

The above cannot be put into a k-k*l format, and therefore cannot be treated with a (1-L)*k type of solution, which is the only type that the RPTO/RPTI method can solve properly. For these more general distributed problems, we need the distributed delay methods described in intro7.pdf (intro7_8_x_y/ADG_x_y.ctl), which provide some type of double variable of integration process, by creating an extended partitioning via transit compartments, or by discrete time partitions using the DDE solver.

Appendix N: Distributed Delay Assessments by the ddexpand Utility, Using a Discrete Delay Solver

The ddexpand utility can utilize discrete delay differential equation solvers to approximate the convolution integral, in the following manner. Consider a state variable $A(x)$ that the user wishes to have convolved with a distribution. The user identifies its initial value with:

\[ \mathrm{APG\_x\_y} = \{\text{expression equivalent to}\ x(0)\} \tag{272}\]

(not AP_x_y), which is for discrete delay), where $x$ is the variable name, and $y$ is an integer index equal to or greater than 1 referring to the $y$th convolution description (the user may have more than one convolution description). If $x$ is an integer, then $x$ represents $A(x)$, a state variable of the series of differential equations in the user’s model.

In the model, the user may refer to

\[ \mathrm{ADG\_x\_y} \]

To represent the convolution integral

\[ \int_{0}^{t} l(s)x(t-s)ds \]

and it is up to ddexpand to create an approximation to this integral and define ADG_x_y accordingly.

Given the further identification by the user of the following (as comments in the control stream file):

;NNy={number of equations to create for convolution approximation}
;PRCy={imprecision value}
;DISTRIBy={distribname}CDF
...
$DES
APG_x_y=x0
...
{equations that may use ADG_x_y, perhaps such as}
DADT(1)=kin-ADG_1_1

where CDF stands for cumulative distribution function. Then these values are used to represent the following: $m$={NNy}

\[ \alpha = 1 - \{\mathrm{PRC}\ y\} \tag{273}\]

\[ L(s) = \{\mathrm{distribname}\}\mathrm{CDF} \tag{274}\]

\[ \begin{aligned} l(s) &= \{\mathrm{distribname}\} \\ &= \frac{dL(s)}{ds} \end{aligned} \tag{275}\]

\[ \tau_e = L^{-1}(\alpha) \tag{276}\]

\[ \begin{aligned} k &= m/\tau_e \\ \Delta\tau &= \frac{1}{k} \end{aligned} \tag{277}\]

\[ \tau_i = i\Delta\tau,\ i = 1\ \text{to}\ m \tag{278}\]

\[ L_i = L(\tau_i) = L(i\Delta\tau),\ i = 1\ \text{to}\ m-1 \tag{279}\]

\[ L_m = L(\infty) = 1,\ i = 1\ \text{to}\ m-1 \tag{280}\]

Then the following is calculated for A$_y$, where $n$ is the base number of user-defined equations already in existence, for state variable:

\[ \mathrm{AD}\_x\_y = \sum_{i=2}^{m}\left[(L_i(\Delta\tau) - L_{i-1}(\Delta\tau))A_i(t) - (1 - L_{m-1})\Delta\tau\right] + (L_m - L_{m-1})A_i(t - m\Delta\tau) + \sum_{i=1}^{m}(L_m - L_{m-1}A_i)A_i(t - \tau_m) \tag{281}\]

which is an approximation to the convolution, because each summand approximates (ALG=1, or ALG=NC(1,T)):

\[ \begin{aligned} L_i(\Delta\tau) - L_{i-1}(\Delta\tau) &\approx l(t - \tau)\Delta\tau \\ &= \int_{(i-1)\Delta\tau}^{i\Delta\tau} l(t - \tau)d\tau \\ &= \int_{\tau_{i-1}}^{\tau_i} A_i(t - \tau)d\tau \end{aligned} \tag{282}\]

The discrete delay solver is able to retrieve $A_i(t - \tau\Delta\tau)$ for any past time $t - i\Delta\tau$. The convolution calculation is done neatly and completely by the built-in function DDEFUNCG, and does not require creating any additional differential equations in the control stream.

The above is based on the simplest form of numerical integration, summing the area of discrete rectangular segments. Even more accurate is the trapezoidal rule equivalent (ALG=2, or ALG=NC(2,T)):

\[ \mathrm{AD}\_x\_y = \sum_{i=2}^{m}\left[\frac{1}{2}A_i(t)\Delta\tau + \frac{1}{2}\sum_{i=1}^{m-1}\left(L_{i+1} - L_{i-1}\right)A_i(t - \tau_i)\right] + \frac{1}{2}(L_m - L_{m-1})A_i(t - \tau_m) \tag{283}\]

Even more accurate than trapezoidal rule, is the Simpson’s rule method for unequal intervals. The basis for evaluating these in Numerical Recipes (Press et al. 1992), p. 127.

Another way of performing the integral is to provide equal intervals of CDF positions, that is, equal intervals of quantiles, in which

\[ \tau_i = L^{-1}(i/m) \tag{284}\]

\[ \begin{aligned} L_i &= L(\tau_i) \\ &= i/m \end{aligned} \tag{285}\]

For the $m$th tau, we have

\[ \tau_m = \tau_e \tag{286}\]

\[ \begin{aligned} L_m &= L(\infty) \\ &= 1 \end{aligned} \tag{287}\]

In so doing, the rectangular integral method (ALG=1, or ALG=NC(1,Q)) is performed as follows (based on equation 281):

\[ \mathrm{AD}\_x\_y = \sum_{i=2}^{m}\left[(L_i - L_{i-1})A_i(t) + (L_m - L_{m-1})\Delta\tau\right] A_i(t - \tau_i) \tag{288}\]

The trapezoidal version (ALG=2, or ALG=NC(2,Q)), based on equation 283:

\[ \mathrm{AD}\_x\_y = \frac{1}{2m}A_i(t)\Delta\tau + \frac{1}{2m}\sum_{i=1}^{m-1}(L_{i+1} - L_{i-1})A_i(t - \tau_i) + \frac{1}{2}(L_m - L_{m-1})A_i(t - \tau_m) \]

Or,

\[ \mathrm{AD}\_x\_y = \frac{1}{2m}A_i(t)\Delta\tau + \frac{1}{m}\sum_{i=1}^{m-1} A_i(t - \tau_i) + \frac{1}{2m}A_i(t - \tau_m) \tag{289}\]

For the Simpson’s rule method (ALG=3, or ALG=NC(3,Q)), we have:

\[ \tau_i = L^{-1}(i/m),\ 0 \leq i \leq 1 \tag{290}\]

\[ \tau_m = \tau_e \tag{291}\]

\[ \begin{aligned} L_i &= L(\infty) \\ &= 1 \end{aligned} \tag{292}\]

\[ l_j = l_i/m,\ j = 0\ \text{to}\ m \tag{293}\]

\[ h_j = L_j - L_{j-1} = l/m = 1/m,\ j = 1\ \text{to}\ m \tag{294}\]

Because we have equal steps $h$, we can use the Simpson’s rule for equal intervals that must be even-valued, consistent with Numerical Recipes (Press et al. 1992), p. 128, equation 4.1.13:

\[ \mathrm{AD}\_x\_y = \frac{1}{m}\left[\frac{1}{3}A_i(t) + \frac{4}{3}\sum_{j=1}^{m/2}A_i(t - \tau_{2j-1}) + \frac{2}{3}\sum_{j=1}^{m/2-1}A_i(t - \tau_{2j}) + \frac{1}{3}A_i(t - \tau_m)\right] \tag{295}\]

The general setup is based on that described in Numerical Recipes, p. 127, third paragraph, which describes that the coefficients to an $M$-point integral provides an exact integration exact for the $M-1$th degree polynomial and lower. With this in mind, let the elemental polynomial functions be

\[ f_i(s) = s^{i-1},\ i = 1\ \text{to}\ M \tag{296}\]

Then the following must be satisfied for the $M$ coefficients $b_i$ of the integral segment, for all $M$ of $l(x)$ functions, at the specific $x$ positions on the grid:

\[ \sum_{j=1}^{M} b_j(x_{kj-1} - x_{ik}) = \int_{x_{ik}}^{x_{ik+1}} f_i(s)ds \tag{297}\]

So,

\[ \begin{aligned} \sum_{j=1}^{M} b_j(x_{kj-1} - x_{ik}) &= \int_{x_{ik}}^{x_{ik+1}} s^{i-1}ds \\ &= (x_{ik+1} - x_{ik})^{i}/i \end{aligned} \tag{298}\]

Let $m$ be the total number of points, then $m$ must be an integral multiple of $(M-1)$ (DDEFUNCG will increase $m$ so this is satisfied, if necessary), so we can replicate this $M$-integral segment process for all $m$ points without residual points left over. So, for each $k = 1, 2, \ldots, m/(M-1)$:

\[ b_j = x_{(M-1)(k-1)+j-1} - x_{(M-1)(k-1)},\ j = 1, M \tag{299}\]

\[ \text{let} \]

\[ h_j = x_{(M-1)(k-1)+j-1} - x_{(M-1)(k-1)},\ j = 1, M \tag{300}\]

Where

\[ \begin{aligned} x_{(M-1)(k-1)+j-1} &= (M-1)(k-1) + j - 1)/m \\ &= x_{(M-1)(k-1)+j-1} \end{aligned} \tag{301}\]

for equally spaced tau values, or

\[ \begin{aligned} x_{(M-1)(k-1)+j-1} &= L_{(M-1)(k-1)+j-1} \\ &= [(M-1)(k-1) + j - 1]/m \end{aligned} \tag{302}\]

For equally spaced CDF positions. Then

\[ C_{ikj} = h_j^{i-1} \tag{303}\]

And vector

\[ y_{ki} = (x_{(M-1)k} - x_{(M-1)(k-1)})^{i}/i \tag{304}\]

Then, the solution vector of coefficients is derived as:

\[ \mathbf{C}_k\mathbf{b}_k = \mathbf{y}_k \tag{305}\]

\[ \mathbf{b}_k = \mathbf{C}_k^{-1}\mathbf{y}_k \tag{306}\]

To produce the composite vector:

\[ \begin{aligned} \beta_0 &= b_{11},\ k=1,\ j=1 \\ \vdots \\ \beta_{(M-1)(k-1)+j-1} &= b_{kj},\ j=2,\ldots M-1,\ k=1,\ldots m/(M-1)-1 \\ \vdots \\ \beta_{(M-1)(k-1)+j-1} &= b_{kM} + b_{k+1,1},\ j=M,\ k=1,\ldots m/(M-1)-1 \\ \vdots \\ \beta_{m-M+j} &= b_{m/(M-1),j},\ j=2,\ldots M-1 \\ \vdots \\ \beta_{m} &= b_{m/(M-1),M} \end{aligned} \tag{307}\]

So the summation is:

\[ S = \sum_{i=0}^{m} \beta_i A(t - \tau_i) \tag{308}\]

Note that when equally spaced CDF positions are used, spacings are equal, and

\[ \begin{aligned} h_j &= x_{(M-1)(k-1)+j-1} - x_{(M-1)(k-1)} \\ &= (j - 1)/m \end{aligned} \tag{309}\]

Then

\[ \begin{aligned} C_{kij} &= [(j-1)/m]^{i-1} \\ &= C_{ij} \end{aligned} \tag{310}\]

And vector

\[ \begin{aligned} y_{ki} &= (x_{(M-1)k} - x_{(M-1)(k-1)})^{i}/i \\ &= [(M-1)/m]^{i}/i \\ &= y_i \end{aligned} \tag{311}\]

So there is one set of vectors $\mathbf{b}$ to be calculated:

\[ \mathbf{b} = \mathbf{C}^{-1}\mathbf{y} \tag{312}\]

which are to be repeatedly used over all $m$ points, $M$ points at a time. And

\[ \begin{aligned} \beta_0 &= b_1 \\ \vdots \\ \beta_{(M-1)(k-1)+j-1} &= b_j,\ j=2,\ldots M-1,\ k=1,\ldots m/(M-1)-1 \\ \vdots \\ \beta_{(M-1)(k-1)+j-1} &= b_M + b_1,\ j=M,\ k=1,\ldots m/(M-1)-1 \\ \vdots \\ \beta_{m-M+j} &= b_j,\ j=2,\ldots M-1 \\ \vdots \\ \beta_m &= b_M \end{aligned} \tag{313}\]

The most efficient method, however, is the Gaussian quadrature method, based on section 4.5 of Numerical Recipes (Press et al. 1992).

Appendix O: Distributed Delay Assessments by the ddexpand Utility Using Transit Compartments

For this method, the ddexpand utility adds distributed delayed differential equations as transit compartment equations in the following manner, and can be solved by a standard ODE solver, not requiring a discrete delay solver. Consider a variable $x$ that the user wishes to have convolved with a distribution. The user identifies its initial value with:

\[ \mathrm{APG\_x\_y} = \{\text{expression equivalent to}\ x(0)\} \tag{314}\]

(not AP_x_y, which is for discrete delay), where $x$ is the variable name, and $y$ is an integer index equal to or greater than 1 referring to the $y$th convolution description (the user may have more than one convolution description). If $x$ is an integer, then $x$ represents $A(x)$, a state variable of the series of differential equations in the user’s model.

In the model, the user may refer to

\[ \mathrm{ADG}\_x\_y \]

To represent the convolution integral

\[ \int_{0}^{t} l(s)x(t-s)ds \]

and it is up to ddexpand to create an approximation to this integral and define $\mathrm{ADG\_x\_y}$ accordingly.

Given the further identification by the user of the following (as comments in the control stream file):

;NNy={number of equations to create for convolution approximation}
;PRCy={imprecision value}
;DISTRIBy={distribname}CDF
...
$DES
APG_x_y=x0
...
{equations that may use ADG_x_y, perhaps such as}
DADT(1)=kin-ADG_1_1

where CDF stands for cumulative distribution function. Then these values are used to represent the following: $m$={NNy}

\[ \alpha = 1 - \{\mathrm{PRC}\ y\} \tag{315}\]

\[ L(s) = \{\mathrm{distribname}\}\mathrm{CDF} \tag{316}\]

\[ \begin{aligned} l(s) &= \{\mathrm{distribname}\} \\ &= \frac{dL(s)}{ds} \end{aligned} \tag{317}\]

\[ \tau_e = L^{-1}(\alpha) \tag{318}\]

\[ \begin{aligned} k &= m/\tau_e \\ \Delta\tau &= \frac{1}{k} \end{aligned} \tag{319}\]

Then the following expansion is made, where $n$ is the base number of user-defined equations already in existence:

\[ A_{i+n}(0) = 0,\ i = 1, m \tag{320}\]

\[ \frac{dA_{1+n}(t)}{dt} = k\left(x(t) - A_{1+n}(t) - \mathrm{APG}\_x\_y\right) \tag{321}\]

\[ \frac{dA_{i+n}(t)}{dt} = k\left(A_{i-1+n} - A_{i+n}\right),\ i = 2, m \tag{322}\]

\[ S = \sum_{i=2}^{m} L((i-1)\Delta\tau)\left(A_{i-1+n} - A_{i+n}\right) + A_{m+n} \tag{323}\]

\[ \mathrm{ADG}\_x\_y = S + \mathrm{APG}\_x\_y \tag{324}\]

Note that an alternative expression is

\[ S = \sum_{i=2}^{m}\left[(L((i-1)\Delta\tau) - L((i-1)\Delta\tau))A_{i-1+n}(t) + (L(m\Delta\tau) - L((m-1)\Delta\tau))A_{m+n}(t)\right] \tag{325}\]

\[ \mathrm{ADG}\_x\_y = S + \mathrm{APG}\_x\_y \tag{326}\]

Notice that the sum of the weights is equal to $L(\infty) = 1$.

An alternative calculation of $S$ and $\mathrm{ADG}\_x\_y$ could be:

\[ S = \sum_{i=1}^{m-1}\left[(L((i+1)\Delta\tau) - L((i-1)\Delta\tau))A_{i+n}(t)\right]/2 + \sum_{i=2}^{m}\left[(L(i\Delta\tau) - L((i-1)\Delta\tau))A_{i+n}(t)\right]/2 \tag{327}\]

\[ L(m\Delta\tau) \approx \sum_{i=2}^{m} h_i/k = W \tag{328}\]

\[ \mathrm{ADG}\_x\_y = S + W \cdot \mathrm{APG}\_x\_y \tag{329}\]

However, the use of the CDF appears to be more accurate.

To show that equations 323 and 324 are a suitable approximation for the convolution of $x$ with the pdf $l$, we start with the derivation in reference (Koch and Schropp 2015), in which it is shown that for an equation type of

\[ \frac{dN(t)}{dt} = k_{in}(t) - \int_{0}^{\infty} l(s)k_{in}(t-s)ds \tag{330}\]

an approximate solution is

\[ \frac{dA_{1+n}(t)}{dt} = k_{in}(t) - kA_{1+n}(t) \tag{331}\]

\[ \frac{dA_{i+n}(t)}{dt} = k\left(A_{i-1+n}(t) - A_{i+n}(t)\right),\ i = 2, m \tag{332}\]

\[ N(t) = \sum_{i=2}^{m} L((i-1)\Delta\tau)\left(A_{i-1+n}(t) - A_{i+n}(t)\right) \tag{333}\]

Comparing the form of equation 331 with that of 321, we note that

\[ k_{in}(t) = kx(t) \tag{334}\]

We prefix the derivative form of $N(t)$ for the sake of generality, so that we may find:

\[ \begin{aligned} \int_{0}^{\infty} l(s)x(t-s)ds &= \frac{1}{k}\int_{0}^{\infty} l(s)k_{in}(t-s)ds \\ &= \frac{k_{in}}{k} - \frac{1}{k}\frac{dN(t)}{dt} \end{aligned} \tag{335}\]

We do so as follows, differentiating equation 333 with respect to $t$:

\[ \begin{aligned} \int_{0}^{\infty} l(s)x(t-s)ds &= \frac{1}{k}\int_{0}^{\infty} l(s)k_{in}(t-s)ds \\ &= \frac{k_{in}}{k} - \frac{1}{k}\frac{dN(t)}{dt} \approx \frac{k_{in}}{k} - \frac{1}{k}\sum_{i=1}^{m}(1 - L((i-1)\Delta\tau))\frac{dA_{i+n}(t)}{dt} = \end{aligned} \tag{336}\]

\[ \frac{k_{in}}{k} - \frac{1}{k}\frac{dA_{1+1+n}(t)}{dt} - \sum_{i=2}^{m}(1 - L((i -1)\Delta\tau))(A_{i-1+n}(t) - A_{i+n}(t)) = \tag{337}\]

\[ \frac{k_{in}}{k} - \left(\frac{k_{in}}{k} - A_{1+n}(t)\right) - (A_{1+n}(t) - A_{m+n}(t)) + \sum_{i=2}^{m} L((i-1)\Delta\tau)(A_{i-1+n}(t) - A_{i+n}(t)) = \tag{338}\]

\[ \sum_{i=2}^{m} L((i-1)\Delta\tau)(A_{i-1+n}(t) - A_{i+n}(t)) + A_{m+n}(t) \tag{339}\]

and we get equation 323. In addition, we incorporate the initial value $x(0)$ of $x(t)$ through $\mathrm{APG}\_x\_y$ rather than set the initial value of all $A_{i+n}(t)$ to $x(0)$.

	LEVWT=0	LEVWT=1
\(w_k\)	\(1/N_g\)	\(N_g/N\)
\(v_k\)	\((N_g-L)/N_g(1+N_g)\)	\((N_g-L)/N_g(N_g/N)\)

Appendix A: Matrix Algebra Tools

Appendix B: Positive Definite Properties

Appendix C: The Fischer Score Matrix for Error Assessment in EM Problems

Appendix D: The Exact Second Derivative Matrix for Error Assessment

Appendix E: Adjustment of Error Matrix for Constraints and Non-Positive Definiteness

Appendix F: Obtaining Analytical Derivatives of Likelihood with Respect to Cholesky of Sigma Parameters

Appendix G: Degrees of Freedom Assessment for OMEGA Priors

Appendix H: Technical Note on NonParametric Analysis

Appendix I: Note on TNPRI

Appendix J: T Distribution Sample Generation

Appendix K: Transformation of Parameters during Classical NONMEM Estimation

Appendix L: Estimating Super ID Level Etas

Appendix M: Simple Distributed Delayed Output Problems and How They May Be Solved with Repetition Parameters

Appendix N: Distributed Delay Assessments by the ddexpand Utility, Using a Discrete Delay Solver

Appendix O: Distributed Delay Assessments by the ddexpand Utility Using Transit Compartments

Appendix G: Degrees of Freedom Assessment for `OMEGA` Priors

Appendix I: Note on `TNPRI`