# Simulate longitudinal, curvilinear, convergent data in R

I would like to simulate data that is similar to a few given observations of longitudinal data (28 measurements per unit) in two groups (see below). The distribution of initial values could be normal or lognormal, and the two groups converge on different values (close to 0 for B, about 25 for A with residual variance) with a slightly non-linear "decay".

How can I simulate the group-specific mean time course and the within-group variability in R with a simple parametric model?

Update: I came up with this function here that simulates somewhat similar data. However, the individual data is too variable, the jumps from time-point to time-point should be smaller. In particular, increases should be less frequent and less pronounced. I'd be most grateful for any hints on how to improve the simulation!

Check the following code on how to simulate from a mixed model. The specification below is with a quadratic model, but you could change it to exponential decay in the specification of the design matrices for the fixed effect X and random effects Z, and with a suitable choice for the fixed effects coefficients beta and covariance matrix for the random effects D to resemble your setting:

set.seed(1234)
n <- 100 # number of subjects
K <- 28 # number of measurements per subject
t_max <- 28 # maximum follow-up time

# we constuct a data frame with the design:
# everyone has a baseline measurment, and then measurements at random follow-up times
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
group = rep(gl(2, n/2, labels = c("B", "A")), each = K))

# design matrices for the fixed and random effects
X <- model.matrix(~ group * (time + I(time^2)), data = DF)
Z <- model.matrix(~ (time + I(time^2)), data = DF)

betas <- c(23, 25, -1.74, 0.033, -0.7, 0.01) # fixed effects coefficients
D11 <- 0.0001 # variance of random intercepts
D22 <- 0.0001 # variance of random slopes
D33 <- 0.0001 #  variance quadratic random slopes
# you could instead here define a full covariance matrix D
# to achieve that the variance at the end is smaller than the variance at the start
# you would probably need to consider negative correlations between the random effects
sigma_y <- 1 # measurement error variance

# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)),
rnorm(n, sd = sqrt(D33)))
# you could instead simulate D from a multivariate normal distribution, e.g., via
# b <- MASS::mvrnorm(n, mean = c(0, 0, 0), D)
# linear predictor
eta_y <- drop(X %*% betas + rowSums(Z * b[DF$$id, ])) # we simulate normal longitudinal data DF$$y <- rnorm(n * K, eta_y, sigma_y)

• That looks very promising, thanks! However, as you say, the within-group variance increases with time whereas I'm looking for a decreasing within-group variance. I'll try to follow your suggestion. – caracal Feb 26 '19 at 21:53