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:

• Given $$V_i$$, each $$y_{ij}$$ is normally distributed:
1. $y_{ij} | V_i \sim {\cal N} \left( \frac{100}{V_i}e^{-k \, t_{ij}} , a^2 \right) .$
• $$V_i$$ is lognormally distributed and depends on the covariate $$w_i$$:
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