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