# Extension to multivariate distributions : Différence entre versions

## The Gaussian model

We would now like to extend the model defined for a unique individual scalar parameter $\psi_i$ to the case where $\psi_i$ is a vector $(\psi_{i,1},\psi_{i,2}, \ldots,\psi_{i,d})$ of individual parameters.

To begin with, we are going to merely generalize the basic model to each component of $\psi_i$. To this end, we suppose that there exists a vector of covariates $c_i = (c_{i,1}, \ldots, c_{i,L})$ and:

• $d$ monotonic transformations $h_1$, $h_2$, $\ldots$, $h_d$
• $d$ vectors of fixed coefficients $\bbeta_1$, $\bbeta_2, \ldots, \bbeta_d$
• $d$ functions $\hmodel_1$, $\hmodel_2$, $\ldots$, $\hmodel_d$
• a vector of random effects $\beeta_i = (\eta_{i,1},\eta_{i,2},\ldots , \eta_{i,d})$,

such that, for each $\iparam=1,2,\ldots,d$,

$$\begin{eqnarray} \hpsi_{i,\iparam} &=& \hmodel_\iparam(\bbeta_\iparam,c_i) \\ h_\iparam(\psi_{i,\iparam}) & =& h_\iparam(\hpsi_{i,\iparam}) +\eta_{i,\iparam} \\ & =& \mmodel_\iparam(\bbeta_\iparam,c_i) +\eta_{i,\iparam}. \end{eqnarray}$$

For instance, a linear covariate model supposes that for each $\iparam=1,2,\ldots,d$, we have:

$$h_\iparam(\hpsi_{i,\iparam}) = h_\iparam(\psi_{ {\rm pop},\iparam})+ \bbeta_{\iparam,1}(c_{i,1} - c_{\rm pop,1}) + \bbeta_{\iparam,2}(c_{i,2} - c_{\rm pop,2}) + \ldots + \bbeta_{\iparam,L}(c_{i,L} - c_{\rm pop,L}) .$$

Dependency can be introduced between parameters by supposing that the random effects $(\eta_{i,\iparam})$ are not independent. In the special case where the random effects are Gaussian, this means considering them to be correlated, i.e., we suppose there exists a $d\times d$ variance-covariance matrix $\Omega$ such that

$$\begin{eqnarray} \esp{\beeta_i}&=&0 \\ \esp{\beeta_i \beeta_i^\prime} &=& \Omega . \end{eqnarray}$$

Here,

$$\Omega = \left( \begin{array}{cccc} \omega_1^2 & \omega_{1,2} & \ldots & \omega_{1,d} \\ \omega_{1,2} & \omega_2^2 & \ldots & \omega_{2,d} \\ \vdots & \vdots & \ddots & \vdots \\ \omega_{1,d} & \omega_{2,d} & \ldots & \omega_d^2 \end{array} \right),$$

where $\omega_\iparam^2$ is the variance of $\eta_{i,\iparam}$ and $\omega_{\iparam,\iparam^\prime}$ the covariance between $\eta_{i,\iparam}$ and $\eta_{i,\iparam^\prime}$.

It will be useful in the following to have a diagonal decomposition of $\Omega$. To this end, let us define the correlation matrix $R=(R_{\iparam,\iparam^\prime}, 1 \leq \iparam,\iparam^\prime \leq d)$ of the vector $\eta_i$:

$$R_{\iparam,\iparam^\prime} = \left\{ \begin{array}{ll} 1 & {\rm if \quad } \iparam=\iparam^\prime \\ \rho_{\iparam,\iparam^\prime}=\frac{\omega_{\iparam,\iparam^\prime} }{\omega_{\iparam}\omega_{\iparam^\prime} } & \hbox{otherwise,} \end{array} \right.$$

and let $D=(D_{\iparam,\iparam^\prime})$ be a diagonal matrix which contains the standard deviations $(\omega_\iparam)$:

$$D_{\iparam,\iparam^\prime} = \left\{ \begin{array}{ll} \omega_{\iparam} & {\rm if \quad } \iparam=\iparam^\prime \\ 0 & {\rm otherwise.} \end{array} \right.$$

Then we have the diagonal decomposition: $\Omega = D \, R \, D$.

Example:

Consider a PK model with three PK parameters: the absorption rate constant $ka$, the volume $V$ and the clearance $Cl$. Here, $\psi_i=(ka_i,V_i,Cl_i)$.

• If we make the assumption that $\eta_{i,V}$ and $\eta_{i,Cl}$ are correlated, it means that the log-volume and the log-clearance are linearly correlated, with correlation:

$$\begin{eqnarray} \rho_{V,Cl} & = & \corr{\eta_{i,V},\eta_{i,Cl} } \\ & = & \corr{\log(V_i),\log(Cl_i)} . \end{eqnarray}$$

• Assuming that $ka$ is fixed in the population means that $ka_i = ka_{\rm pop}$ for any $i$, which implies that $\eta_{i,ka}=0$, and thus $\omega_{ka}=0$. The correlation matrix $R$ and the variance covariance matrix $\Omega$ of $(\eta_{i,ka}, \eta_{i,V}, \eta_{i,Cl})$ are therefore

$$R = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & \rho_{V,Cl} \\ 0 & \rho_{V,Cl} & 1 \\ \end{array} \right) , \quad \quad \Omega = DRD = \left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & \omega_V^2 & \omega_V\omega_{Cl}\, \rho_{V,Cl} \\ 0 & \omega_V\omega_{Cl}\, \rho_{V,Cl} & \omega_{Cl}^2 \\ \end{array} \right) .$$

• ## The probability distribution function

We have now all the elements needed for computing the pdf of $\psi_i=(\psi_{i,1},\psi_{i,2}, \ldots,\psi_{i,d})$. Here, $\theta = (\psi_{{\rm pop},1}, \ldots, \psi_{{\rm pop},d}, \bbeta_1, \ldots,\bbeta_d,\Omega)$.

• If $\Omega$ is a positive-definite matrix, it can be inverted and a straightforward extension of the pdf proposed in (9) of The covariate model for a scalar variable gives
 $$\ppsii(\psi_i;c_i,\theta )= \left( \prod_{\iparam=1}^d h_\iparam^\prime(\psi_{i,\iparam}) \right) (2 \pi)^{-\frac{d}{2} } |\Omega|^{-\frac{1}{2} } {\rm exp} \left\{-\frac{1}{2} ( h(\psi_i) - \mmodel(\bbeta,c_i) )^\prime \Omega^{-1} ( h(\psi_i) - \mmodel(\bbeta,c_i) ) \right\} ,$$ (1)
where $h(\psi_i)$ is the column vector $(h_1(\psi_{i,1}), h_2(\psi_{i,2}), \ldots, h_d(\psi_{i,d}))^\prime$ and $\mmodel(\bbeta,c_i)$ the column vector $(h_1(\hpsi_{i,1}), h_2(\hpsi_{i,2}), \ldots, h_d(\hpsi_{i,d}))^\prime$.

• If the variance of some of the random effects is null, $\Omega$ is not positive-definite. The pdf in (1) does not apply anymore for the complete $d$-vector $\psi_i$ but only for the $d_1$-vector subset $\psi_i^{(1)}$ of $\psi_i$ whose variance matrix $\Omega_1$ is positive-definite. The distribution of the remaining fixed parameters $\psi_i^{(0)}$ is a Dirac delta distribution. Let $I_0$ be the indices of the parameters $\psi_i^{(0)}$ and $I_1$ those of the parameters $\psi_i^{(1)}$, i.e., $\omega_\iparam =0$ if $\iparam \in I_0$ and $\omega_\iparam >0$ if $\iparam \in I_1$. Then,

$$\ppsii(\psi_i;c_i,\theta ) = \pmacro(\psi_i^{(0)};c_i,\theta )\,\,\pmacro(\psi_i^{(1)};c_i,\theta ) ,$$

where

$$\begin{eqnarray} \pmacro(\psi_i^{(0)};c_i,\theta ) &= & \prod_{\iparam \in I_0} \delta_{ \{ h(\psi_{i,\iparam})=\mmodel_{\iparam}(\bbeta_{\iparam},c_i) \} } \\ \pmacro(\psi_i^{(1)};c_i,\theta )&=& \left( \prod_{\iparam \in I_1} h_\iparam^\prime(\psi_{i,\iparam}) \right) (2 \pi)^{-\frac{d_1}{2} } |\Omega_1|^{-\frac{1}{2} } {\rm exp} \left\{ -\frac{1}{2} ( h(\psi_i) - \mmodel(\bbeta,c_i) )^{(1)^\prime} \Omega_1^{-1} ( h(\psi_i) - \mmodel(\bbeta,c_i) )^{(1)} \right\}, \end{eqnarray}$$

with $( h(\psi_i) - \mmodel(\bbeta,c_i) )^{(1)}$ the same as $( h(\psi_i) - \mmodel(\bbeta,c_i) )$ but with the $I_0$ entries removed.

• There exist other situations where $\Omega$ is not positive-definite. This is the case for instance when two random effects are equal: $\eta_{i,\iparam} \equiv \eta_{i,\iparam^\prime}$. For them, we can calculate a joint distribution

$$\pmacro(\psi_{i,\iparam}, \psi_{i,\iparam^\prime},\eta_{i,\iparam}; \bbeta_{\iparam},\bbeta_{\iparam^\prime},\omega^2_\iparam ,c_i ) = \pmacro(\psi_{i,\iparam} | \eta_{i,\iparam}; \bbeta_{\iparam}, c_i ) \ \pmacro(\psi_{i,\iparam^\prime} | \eta_{i,\iparam};\bbeta_{\iparam^\prime} , c_i ) \ \pmacro(\eta_{i,\iparam} ; \omega^2_\iparam) ,$$

where

$$\begin{eqnarray} \pmacro(\psi_{i,\iparam}| \eta_{i,\iparam}; \bbeta_{\iparam}, c_i ) &=& \delta_{\{h(\psi_{i,\iparam})=\mmodel_{\iparam}(\bbeta_{\iparam},c_i)+\eta_{i,\iparam} \} } \\ \pmacro( \psi_{i,\iparam^\prime} | \eta_{i,\iparam};\bbeta_{\iparam^\prime} , c_i ) &=& \delta_{\{h(\psi_{i,\iparam^\prime})=\mmodel_{\iparam^\prime}(\bbeta_{\iparam^\prime},c_i)+\eta_{i,\iparam} \} } \\ \pmacro(\eta_{i,\iparam} ; \omega^2_\iparam) &=& \displaystyle{ \frac{ 1}{\sqrt{2 \, \pi \omega_\iparam^2 } } }\ \exp\left\{-\displaystyle{ \frac{\eta_{i,\iparam}^2}{2 \omega_\iparam^2} }\right\}. \end{eqnarray}$$

All kinds of combinations are possible, including parameters with and without variability, algebraic relationships between random effects, etc. In all possible cases it is possible to find an adequate decomposition that lets us characterize a pdf. This pdf turns out to play a fundamental role for tasks such as population parameter estimation with maximum likelihood, where we start with the observations $\by = (y_i , 1\leq i \leq N)$ and the individual parameters $(\psi_i)$ are not observed.