R scripts: hierarchical2.R

Mlxtran codes: model/hierarchical2.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\) is treated as a random vector in a population context.

When the probability distribution of \(\psi_i\) depends on a vector of individual covariates \(c_i\), we can also consider that \(c_i\) is a random vector sampled from a population distribution \(\pmacro(\, c_i \,;\theta)\).

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

The model \(\pmacro(c_i ; \theta)\) for the individual covariates is implemented in the section [COVARIATE].


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) . \]
  1. \[ \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). \]

We assume furthermore that \(w_i\) is normally distributed:

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

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, w, w_pop}

EQUATION:
V_pred = V_pop*(w/w_pop)

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}

We will use simulx for sampling one individual covariate \(w_i\), one individual parameter \(V_i\) and one sequence \(y_i=(y_{ij})\) for one individual \(i=1\), with some given values of the population parameters \(V_{pop}\), \(\omega_V\), \(w_{pop}\), \(\omega_w\), \(k\) and \(a\).

p <- c(V_pop=10, omega_V=0.1, beta=1, w_pop=70, omega_w=12, 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'))
out <- list(ind, f, y)

res1 <- simulx(model     = 'model/hierarchical2.txt', 
               parameter = p, 
               output    = out,
               settings  = list(seed = 12345))

When several individual parameters and/or covariates are defined in the list of outputs, simulx also returns them as an additional element parameter.

print(res1$parameter)
##          w        V
## 1 70.98498 11.04154
plot(ggplot() + 
  geom_line( data=res1$f, aes(x=time, y=f), colour="black") + 
  geom_point(data=res1$y, aes(x=time, y=y), colour="red"))

Instead of only one individual, let us now simulate a group of 5 individuals with individual covariates \(w_i\) drawn for each individual \(i=1,2,\ldots,5\). The level of randomization is therefore ‘covariate’.

g <- list( size  = 5,
           level = 'covariate')

res2 <- simulx(model     = 'model/hierarchical2.txt', 
               parameter = p, 
               output    = out,
               group     = g,
               settings  = list(seed = 12345))

print(res2$parameter)
##   id        w         V
## 1  1 70.98498 11.041537
## 2  2 80.21273 11.424322
## 3  3 69.63671  9.498118
## 4  4 64.44544  9.627048
## 5  5 75.36014  8.524051
plot(ggplot() + geom_line(data=res2$f, aes(x=time, y=f, colour=id), size=0.75) + 
       geom_point(data=res2$y, aes(x=time, y=y, colour=id), size=2))

We can combine several levels of randomization, and draw, for instance, 2 individual covariates \(w_1\), \(w_2\) and 3 individual parameters \(V_i\) per covariate

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

res3 <- simulx(model     = 'model/hierarchical2.txt', 
               parameter = p, 
               output    = out,
               group     = g,
               settings  = list(seed = 12345))

print(res3$parameter)
##   id        w         V
## 1  1 70.98498 11.041537
## 2  2 70.98498 10.110058
## 3  3 70.98498  9.682016
## 4  4 80.21273 11.982412
## 5  5 80.21273  9.072930
## 6  6 80.21273 10.050901

In the next example, we draw 2 individual covariates, 2 individual parameters per covariate and 2 sequences of longitudinal data per individual parameters.

g <- list( size  = c(2,2,2),
           level = c('covariate','individual','longitudinal'))

res4 <- simulx(model     = 'model/hierarchical2.txt', 
               parameter = p, 
               output    = out,
               group     = g,
               settings  = list(seed = 123123))

print(res4$parameter)
##   id         w        V
## 1  1  69.97862 12.85809
## 2  2  69.97862 12.85809
## 3  3  69.97862 10.12943
## 4  4  69.97862 10.12943
## 5  5 100.20319 13.28821
## 6  6 100.20319 13.28821
## 7  7 100.20319 16.36998
## 8  8 100.20319 16.36998
plot(ggplot() + geom_line(data=res4$f, aes(x=time, y=f, colour=id), size=0.75) + 
  geom_point(data=res4$y, aes(x=time, y=y, colour=id), size=2))

pa <- c(V_pop=10, omega_V=0.1, w_pop=70, beta=1, k=0.15, a=0.5)
pw <- data.frame(id=1:4, w=c(60,70,80,90))
res5a <- simulx(model     = 'model/hierarchical1b.txt', 
               parameter = list(pa,pw), 
               output    = out,
               settings  = list(seed = 123123))

print(res5a$parameter)
##   id  w         V
## 1  1 60  8.569901
## 2  2 70 12.862014
## 3  3 80 11.580032
## 4  4 90 11.935142
res5b <- simulx(model     = 'model/hierarchical1b.txt', 
               parameter = list(pa,pw), 
               output    = out,
               group     = list(size = 7),
               settings  = list(seed = 1231))

print(res5b$parameter)
##   id  w         V
## 1  1 60  7.774017
## 2  2 70  9.277483
## 3  3 80 10.296795
## 4  4 90 13.118423
## 5  5 60  8.206328
## 6  6 70  9.853925
## 7  7 90 14.364631
res5c <- simulx(model     = 'model/hierarchical1b.txt', 
               parameter = list(pa,pw), 
               output    = out,
               group     = list(size = 2),
               settings  = list(seed = 123123))

print(res5c$parameter)
##   id  w         V
## 1  1 70  9.998218
## 2  2 80 14.699444