R scripts: hierarchical2.R
Mlxtran codes: model/hierarchical2.txt
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]
.
We consider the model defined in the previous article:
We assume furthermore that \(w_i\) is 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, 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