R scripts: hierarchical3.R
Mlxtran codes: model/hierarchical3.txt
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]
.
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:
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