110 Three Stage EM Analysis
NONMEM 7 Technical Guide
There are times when one desires to use information from a previous analysis and incorporate it into the present analysis. This would be in the form of prior information for thetas and/or omegas. The principle on which this is based is as follows. Let \(\boldsymbol{\theta}_0\) be the priors to the thetas (theta priors, which could be estimates of theta from a previous analysis). Let the matrix \(\boldsymbol{\Omega}_\theta^{-1}\) be the information matrix (which could be the theta portion of the inverse of the standard error variance matrix of a previous analysis) of the theta priors. Then \(\boldsymbol{\Omega}_\theta\) may be called variance to theta priors, or theta variance priors. Let \(\boldsymbol{\Omega}_\Omega\) be the prior to the omegas of the population inter-subject variance-covariance matrix, of dimension \(p\) (Omega priors, could be Omega estimates of a previous analysis), let \(\rho\) be the degrees of freedom of \(\boldsymbol{\Omega}_\Omega\) (degrees of freedom priors, could be the number of subjects in that prior analysis). The contribution to the objective function that incorporates this prior information is
\[ L_p = -\log\bigl(N(\boldsymbol{\theta} \mid \boldsymbol{\theta}_0, \boldsymbol{\Omega}_\theta)\bigr) - \log\bigl(W^{-1}(\rho \boldsymbol{\Omega}_\Omega \mid \boldsymbol{\Omega}_\Omega, \rho + p + 1)\bigr) \tag{1}\]
And is then added to equation 11:
\[ L = -\sum_{i=1}^{m} \log\left(\int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\mu}, \boldsymbol{\Omega}) d\boldsymbol{\phi}\right) + L_p \tag{2}\]
where
\[ -\log(N(\boldsymbol{\theta} \mid \boldsymbol{\theta}_0, \boldsymbol{\Omega}_\theta)) = \frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\theta}_0)^{\prime} \boldsymbol{\Omega}_\theta^{-1}(\boldsymbol{\theta} - \boldsymbol{\theta}_0) + \frac{1}{2}\log(\det(\boldsymbol{\Omega}_\theta)) \tag{3}\]
\[ \begin{aligned} & -\log(W^{-1}(\rho \boldsymbol{\Omega}_\Omega \mid \boldsymbol{\Omega}_\Omega, d_w)) = \\ &\quad \frac{1}{2}\left(\rho\operatorname{tr}(\boldsymbol{\Omega}_\Omega \boldsymbol{\Omega}^{-1}) + (d_w - n - 1)\ln(|\boldsymbol{\Omega}|) - d_w\bigl[\ln(|\boldsymbol{\Omega}_\Omega|) + n\ln(\rho)\bigr]\right) \end{aligned} \tag{4}\]
(not including constants) where \(n\) is the dimension of \(\boldsymbol{\Omega}\). The degrees of freedom \(d_w\) will be described later. It follows that the partial derivatives to \(L\) contributed by this prior information are:
\[ \frac{\partial L_p}{\partial \boldsymbol{\theta}} = -\boldsymbol{\Omega}_\theta^{-1}(\boldsymbol{\theta} - \boldsymbol{\theta}_0) \tag{5}\]
\[ \frac{\partial L_p}{\partial \omega_{j_1, j_2}} = \rho c(j_1, j_2) \mathbf{I}_{j_1}^{\prime} \boldsymbol{\Omega}^{-1}(z_w \boldsymbol{\Omega} - \boldsymbol{\Omega}_\Omega) \boldsymbol{\Omega}^{-1} \mathbf{I}_{j_2} \tag{6}\]
where
\[ \begin{aligned} c(j_1, j_2) &= 1 \text{ for } j_1 \ne j_2 \\ &= 1/2 \text{ for } j_1 = j_2 \end{aligned} \tag{7}\]
and
\[ z_w = (d_w - n - 1)/\rho \tag{8}\]
Which suggest the following updates. To determine the \(\boldsymbol{\mu}\) modeled theta that minimize the objective function, we must solve adding the contribution from the prior:
\[ \begin{aligned} \frac{\partial L}{\partial \boldsymbol{\theta}_\mu} & = -\sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \boldsymbol{\Omega}^{-1}(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) - \boldsymbol{\Omega}_\theta^{-1}(\boldsymbol{\theta} - \boldsymbol{\theta}_0) \\ &= \mathbf{g}_{\boldsymbol{\theta}_\mu} \end{aligned} \tag{9}\]
and
\[ \mathbf{H}_{\boldsymbol{\theta}_\mu} = \sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta} } \boldsymbol{\Omega}^{-1} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta} } + \boldsymbol{\Omega}_\theta^{-1} \tag{10}\]
And performing the following update As before:
\[ \hat{\boldsymbol{\theta}}_\mu = \boldsymbol{\theta}_\mu + \alpha \mathbf{H}_{\boldsymbol{\theta}_\mu}^{-1} \mathbf{g}_{\boldsymbol{\theta}_\mu} \tag{11}\]
For non-MU-modeled thetas,
\[ \begin{aligned} \frac{ \partial L_i }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } &= E_{\boldsymbol{\phi}, i}\left(\frac{ -\partial \log(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})) }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} }\right) \\ &= \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i} \end{aligned} \tag{12}\]
\[ \begin{aligned} \frac{ \partial L }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } & = \sum_{i=1}^{m} \frac{ \partial L_i }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } - \boldsymbol{\Omega}_\theta^{-1}(\boldsymbol{\theta} - \boldsymbol{\theta}_0) \\ &= \sum_{i=1}^{m} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i} - \boldsymbol{\Omega}_\theta^{-1}(\boldsymbol{\theta} - \boldsymbol{\theta}_0) \\ &= \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} \\ &= \mathbf{0} \end{aligned} \tag{13}\]
\[ \mathbf{H}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} = \sum_{i=1}^{m} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i}^{\prime} + \boldsymbol{\Omega}_\theta^{-1} \tag{14}\]
\[ \hat{\boldsymbol{\theta}}_{\enclose{horizontalstrike}{\mu}} = \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \mathbf{H}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}}^{-1} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} \tag{15}\]
The inter-subject variances are updated as
\[ \hat{\boldsymbol{\Omega}} = \frac{1}{m + d_w - n - 1}\left[\sum_{i=1}^{m} (\bar{\boldsymbol{\phi}}_i - \hat{\boldsymbol{\mu}})(\bar{\boldsymbol{\phi}}_i - \hat{\boldsymbol{\mu}})^{\prime} + \sum_{i=1}^{m} \bar{\mathbf{B}}_i + \rho \boldsymbol{\Omega}_\Omega\right] \tag{16}\]
These were derived from setting partial derivatives of the objective function to 0 and using the appropriate “inverse” densities and particular degrees of freedom in the objective function.
For maximization methods (all methods except BAYES), the degrees of freedom of the inverse Wishart (actually, omega inverse as Wishart distributed), is selected as
\[ d_w = \rho + n + 1 \tag{17}\]
so that the maximization of these densities leads to a centering about the prior inter-subject variances, weighted according to the number of subjects from that previous analysis, and with a denominator term of \(m + \rho\), yielding an intuitive update. That is, the density of omega inverse whose mode is at \(\boldsymbol{\Omega}_\Omega^{-1}\) is \(W(\boldsymbol{\Omega}^{-1} \mid [\rho \boldsymbol{\Omega}_\Omega]^{-1}, \rho + n + 1)\). We shall call this the “modal” or maximization version of adding the prior information. The density of the omega inverse whose mean is at \(\boldsymbol{\Omega}_\Omega^{-1}\) is \(W(\boldsymbol{\Omega}^{-1} \mid [\rho \boldsymbol{\Omega}_\Omega]^{-1}, \rho)\). A BAYES analysis is concerned with obtaining average population parameters rather than best fit or modal population parameters, so it utilizes the degrees of freedom
\[ d_w = \rho \tag{18}\]
which we shall call the “mean” version of adding the prior information.
The above equations are also suggested by the sample distribution equations listed on page 341 of (Bennett et al. 1996).
The priors to \(\boldsymbol{\Sigma}\) are also inverse Wishart distributed with prior parameters \((\boldsymbol{\Sigma}_L, \rho_L)\) so similar relationships hold, as for \(\boldsymbol{\Omega}\) priors. However, \(\boldsymbol{\Sigma}\) parameters are embedded in the data likelihood portion of the total likelihood in a non-linear manner, so updates need to be performed by extending the first and second derivatives of the total likelihood with respect to \(\boldsymbol{\Sigma}\), and then using them in the Gauss-Newton update process. With this in mind, we need to find the partial derivative of the prior portion of the objective function with respect to each of the cholesky elements to \(\boldsymbol{\Sigma}\), since this is how we vary the parameters in \(\boldsymbol{\Sigma}\). Let \(\boldsymbol{\Lambda}\) be the cholesky matrix to \(\boldsymbol{\Sigma}\):
\[ \boldsymbol{\Sigma} = \boldsymbol{\Lambda}\boldsymbol{\Lambda}^{\prime} \tag{19}\]
So,
\[ L_{\boldsymbol{\Sigma}_L} = \frac{1}{2}\left(\rho_L \operatorname{tr}(\boldsymbol{\Sigma}_L \boldsymbol{\Sigma}^{-1}) + (d_\Sigma - n_\Sigma - 1)\ln(|\boldsymbol{\Sigma}|) - d_\Sigma\bigl[\ln(|\boldsymbol{\Sigma}_L|) + n_\Sigma\ln(\rho_L)\bigr]\right) \tag{20}\]
\[ \frac{1}{2}\left(\rho_L \operatorname{tr}(\boldsymbol{\Sigma}_L \boldsymbol{\Lambda}^{\prime-1} \boldsymbol{\Lambda}^{-1}) + (d_\Sigma - n_\Sigma - 1)\ln(|\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{\prime}|) - d_\Sigma\bigl[\ln(|\boldsymbol{\Sigma}_L|) + n_\Sigma\ln(\rho_L)\bigr]\right) \tag{21}\]
It is the cholesky elements in \(\boldsymbol{\Lambda}\) that are varied to optimize the likelihood, so,
\[ \frac{ \partial L_{\boldsymbol{\Sigma}_L} }{ \partial \boldsymbol{\Lambda} } = -\rho_L \boldsymbol{\Lambda}^{\prime-1} \boldsymbol{\Lambda}^{-1} \boldsymbol{\Sigma}_L \boldsymbol{\Lambda}^{-1} + (d_\Sigma - n_\Sigma - 1) \boldsymbol{\Lambda}^{-1} \tag{22}\]
\[ \frac{ \partial L_{\boldsymbol{\Sigma}_L} }{ \partial \lambda_{j_1, j_2} \partial \lambda_{j_3, j_4} } = \rho_L \mathbf{I}_{j_3}^{\prime} \boldsymbol{\Lambda}^{\prime-1} \mathbf{I}_{k_2, k_1} \boldsymbol{\Sigma}_L \boldsymbol{\Lambda}^{-1} \mathbf{I}_{j_1} + \rho_L \mathbf{I}_{j_3}^{\prime} \boldsymbol{\Lambda}^{\prime-1} \boldsymbol{\Sigma}_L \mathbf{I}_{k_2, k_1} \boldsymbol{\Lambda}^{-1} \mathbf{I}_{j_1} - \mathbf{I}_{j_3} \frac{ \partial L_{\boldsymbol{\Sigma}_L} }{ \partial \boldsymbol{\Lambda} } \mathbf{I}_{k_1, j_2} \boldsymbol{\Lambda}^{-1} \mathbf{I}_{j_1} \tag{23}\]
Where
\(\mathbf{I}_j = 1\) for vector element \(j\)
\(= 0\) otherwise
\(\mathbf{I}_{j,k} = 1\) for matrix element \(j, k\)
\(= 0\) otherwise