R scripts: hierarchical3.R

Mlxtran codes: model/hierarchical3.txt


$$ \newcommand{\esp}[1]{\mathbb{E}\left(#1\right)} \newcommand{\var}[1]{\mbox{Var}\left(#1\right)} \newcommand{\deriv}[1]{\dot{#1}(t)} \newcommand{\prob}[1]{ \mathbb{P}\!(#1)} \newcommand{\eqdef}{\mathop{=}\limits^{\mathrm{def}}} \newcommand{\by}{\boldsymbol{y}} \newcommand{\bc}{\boldsymbol{c}} \newcommand{\bpsi}{\boldsymbol{\psi}} \def\pmacro{\texttt{p}} \def\like{{\cal L}} \def\llike{{\cal LL}} \def\logit{{\rm logit}} \def\probit{{\rm probit}} \def\one{{\rm 1\!I}} \def\iid{\mathop{\sim}_{\rm i.i.d.}} \def\simh0{\mathop{\sim}_{H_0}} \def\df{\texttt{df}} \def\res{e} \def\xomega{x} \newcommand{\argmin}[1]{{\rm arg}\min_{#1}} \newcommand{\argmax}[1]{{\rm arg}\max_{#1}} \newcommand{\Rset}{\mbox{$\mathbb{R}$}} \def\param{\theta} \def\setparam{\Theta} \def\xnew{x_{\rm new}} \def\fnew{f_{\rm new}} \def\ynew{y_{\rm new}} \def\nnew{n_{\rm new}} \def\enew{e_{\rm new}} \def\Xnew{X_{\rm new}} \def\hfnew{\widehat{\fnew}} \def\degree{m} \def\nbeta{d} \newcommand{\limite}[1]{\mathop{\longrightarrow}\limits_{#1}} \def\ka{k{\scriptstyle a}} \def\ska{k{\scriptscriptstyle a}} \def\kel{k{\scriptstyle e}} \def\skel{k{\scriptscriptstyle e}} \def\cl{C{\small l}} \def\Tlag{T\hspace{-0.1em}{\scriptstyle lag}} \def\sTlag{T\hspace{-0.07em}{\scriptscriptstyle lag}} \def\Tk{T\hspace{-0.1em}{\scriptstyle k0}} \def\sTk{T\hspace{-0.07em}{\scriptscriptstyle k0}} \def\thalf{t{\scriptstyle 1/2}} \newcommand{\Dphi}[1]{\partial_\pphi #1} \def\asigma{a} \def\pphi{\psi} \newcommand{\stheta}{{\theta^\star}} \newcommand{\htheta}{{\widehat{\theta}}}$$

1 Introduction

We have seen in the previous article that the vector of individual parameters \(\psi_i\) and the vector of individual covariates \(c_i\) can be treated as random vectors in a population context.

The probability distributions of \(\psi_i\) and \(c_i\) depend usually on a vector of population parameters \(\theta\), that can also be treated as a random vector with probability distribution \(\pmacro(\theta)\). This distribution can be either used as a prior distribution in a Bayesian context, as a description of the uncertainty about the population parameters, or a a description of some inter-population variability.

In this context, the model for individual \(i\) is the joint distribution of the longitudinal data \(y_i\), the individual parameters \(\psi_i\), the individual covariates \(c_i\) and the population parameters \(\theta\): \[\pmacro(y_i , \psi_i, c_i; \theta)=\pmacro(y_i | \psi_i ) \, \pmacro(\psi_i | c_i , \theta) \, \pmacro(c_i | \theta) \, \pmacro(\theta).\]

The model \(\pmacro(\theta)\) for the population parameters is implemented in the section [POPULATION].


2 Example

We consider the model defined in the previous article:

  1. \[ y_{ij} | V_i \sim {\cal N} \left( \frac{100}{V_i}e^{-k \, t_{ij}} , a^2 \right) . \]

  2. \[ \log(V_i) \iid {\cal N}\left(\log \left(V_{\rm pop} \left({w_i}/{w_{\rm pop}}\right)^{\beta}\right) \ ,\ \omega_V^2\right). \]

  3. \[ w_i \iid {\cal N}\left(w_{\rm pop} \ ,\ \omega_w^2\right). \]

We furthermore assume that \(w_{\rm pop}\) and \(V_{\rm pop}\) are, respectively, normally and log-normally distributed:

  1. \[ \begin{align} w_{\rm pop} & \iid {\cal N}\left(w_{s} \ ,\ \gamma_w^2\right) \\ \log(V_{\rm pop}) & \iid {\cal N}\left(\log(V_{s}) \ ,\ \gamma_V^2\right) \end{align} \]

This model is implemented in hierarchical3.txt.

[LONGITUDINAL]
input = {V, k, a}
EQUATION:
f = 100/V*exp(-k*t)
DEFINITION:
y = {distribution = normal, prediction = f, sd = a}

;----------------------------------------------
[INDIVIDUAL]
input = {V_pop, omega_V, beta, w, w_pop}
EQUATION:
V_pred = V_pop*(w/w_pop)^beta
DEFINITION:
V = {distribution = lognormal, prediction = V_pred, sd = omega_V}

;----------------------------------------------
[COVARIATE]
input = {w_pop, omega_w}
DEFINITION:
w = {distribution = normal, mean = w_pop, sd = omega_w}

;----------------------------------------------
[POPULATION]
input = {ws, gw, Vs, gV}
EQUATION:
lV = log(Vs)
DEFINITION:
w_pop = {distribution = normal,    mean = ws, sd = gw}
V_pop = {distribution = lognormal, mean = lV, sd = gV}

We will use simulx for sampling one set of population parameters (\(w_{\rm pop}\), \(V_{\rm pop}\)), one individual covariate \(w_i\), one individual parameter \(V_i\) and one sequence \(y_i=(y_{ij})\) for one individual \(i\).

p <- c(Vs=10, gV=0.1, ws=70, gw=10, omega_w=12, omega_V=0.15, beta=1, k=0.15, a=0.5)
f <- list(name='f', time=seq(0, 30, by=0.1))
y <- list(name='y', time=seq(1, 30, by=3))
ind <- list(name=c('w','V'))
pop <- list(name=c('w_pop','V_pop'))

res <- simulx(model     = 'model/hierarchical3.txt', 
              parameter = p, 
              output    = list(pop, ind, f, y),
              settings  = list(seed = 123456))

print(res$parameter)
##      w_pop    V_pop        w        V
## 1 63.70506 12.80637 70.11372 12.98575
print(ggplot() + geom_line(data=res$f, aes(x=time, y=f), size=1) + 
geom_point(data=res$y, aes(x=time, y=y), colour='red'))

Instead of only one individual, we can easily simulate a group of several individuals combining several levels of randomization.

We can draw, for instance, 2 populations and 3 individual covariates per population,

g <- list(size=c(2,3),
          level=c('population','covariate'))

res <- simulx(model     = 'model/hierarchical3.txt', 
              parameter = p, 
              output    = list(pop, ind, f, y),
              group     = g,
              settings  = list(seed = 123456))

print(res$parameter)
##   id    w_pop    V_pop        w        V
## 1  1 63.70506 12.80637 70.11372 12.98575
## 2  2 63.70506 12.80637 57.14946 13.26632
## 3  3 63.70506 12.80637 75.21543 13.62218
## 4  4 67.24433 10.54857 58.89745 10.90190
## 5  5 67.24433 10.54857 80.48285 14.91550
## 6  6 67.24433 10.54857 80.58027 12.67286
print(ggplot() + geom_line(data=res$f, aes(x=time, y=f, colour=id), size=1) + 
geom_point(data=res$y, aes(x=time, y=y, colour=id)))

or 2 populations and 3 individual parameters per population (in such case, the three individuals of a same population share the same covariates).

g <- list(size=c(2,3),
          level=c('population','individual'))

res <- simulx(model     = 'model/hierarchical3.txt', 
              parameter = p, 
              output    = list(pop, ind, f, y),
              group     = g,
              settings  = list(seed = 123456))

print(res$parameter)
##   id    w_pop    V_pop        w         V
## 1  1 63.70506 12.80637 70.11372 12.985748
## 2  2 63.70506 12.80637 70.11372 16.275764
## 3  3 63.70506 12.80637 70.11372 12.698216
## 4  4 67.24433 10.54857 60.68873 11.233467
## 5  5 67.24433 10.54857 60.68873 11.247154
## 6  6 67.24433 10.54857 60.68873  9.544515