Simulation of a hierarchical model: the population parameters
Marc Lavielle April 18th, 2017
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:
\[ y_{ij} | V_i \sim {\cal N} \left( \frac{100}{V_i}e^{-k \, t_{ij}} , a^2 \right) . \]
\[ \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). \]
\[ 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:
- \[ \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.98575print(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.67286print(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