# Simulating longitudinal lognormal data in R

I want to design a study that will eventually allow me to estimate centile curves for a given longitudinal outcome measured in a sample. I want to simulate data that I can then use to evaluate likely sample sizes needed to estimate the LMS centile curves with a given degree of precision at different centiles and time points.

Unfortunately, there isn't that much info out there on the distribution of this outcome, as the vast majority of studies just report cross-sectional means and SDs (even when the data themselves are longitudinal). However, based on a small number of studies, I can make decent guesses about the mean and variance of the outcome at each time point, and the covariances of measures that are adjacent in time.

I have no information on covariances of non-consecutive measures, and so would like to simulate data under multiple, reasonable autoregressive scenarios. I also have no info on possible skew or kurtosis at any time point, so again would like to evaluate multiple plausible scenarios.

Question: Is there a way to simulate longitudinal data in R when there is possible skewness and/or kutosis? Ignoring the possible skew and kurtosis, I think I can do this with mvrnorm, per this thread: How to simulate multivariate outcomes in R?

Any nudges in the right direction would be appreciated.

• Use mvrnorm and then exponentiate to get log-normal variables?
– Bill
Nov 19, 2013 at 20:53
• Log-normal makes sense (the data are human growth related), so I'm now looking around for tips to simulate multivariate lognormal distributions given the means, variances, and covariances I have at hand. Nov 20, 2013 at 12:05
• You simulate log normal by simulating normal and then exponentiating. You need to figure out how to get the log normal variables to have the variances and covariances you want, though. Wikipedia has the formula to go from the normal means, variances and covariances to the log normal means, variances, and covariances. en.m.wikipedia.org/wiki/Log-normal_distribution. Do you need help inverting the formula there, or is this enough info?
– Bill
Nov 21, 2013 at 0:40
• I see how to take the observed mean and variance of the assumed log-normal variable and calcuate the 2 parameters of the log-normal distribution. I also see how to similarly transform the covariances I have. I still need to get my brain around adding correlations beyond the consecutive measures, but it seems like I can work that out through a multivariate normal first. Thanks Bill. Nov 21, 2013 at 10:07
• @DLDahly Can't you sample your "meanlog" from a trendline or a normal distribution using the weighted average of the previous x points as the mean? Nov 25, 2013 at 14:55

Without more information, it is difficult to guess exactly what kind of data you would like to simulate, but here is an example. In my experience, growth data are often (approximately) linear when transforming (log, square-root, etc.) the predictor (time) variable, so this would be my suggestion.

library(MASS)
library(nlme)

### set number of individuals
n <- 200

### average intercept and slope
beta0 <- 1.0
beta1 <- 6.0

### true autocorrelation
ar.val <- .4

### true error SD, intercept SD, slope SD, and intercept-slope cor
sigma <- 1.5
tau0  <- 2.5
tau1  <- 2.0
tau01 <- 0.3

### maximum number of possible observations
m <- 10

### simulate number of observations for each individual
p <- round(runif(n,4,m))

### simulate observation moments (assume everybody has 1st obs)
obs <- unlist(sapply(p, function(x) c(1, sort(sample(2:m, x-1, replace=FALSE)))))

### set up data frame
dat <- data.frame(id=rep(1:n, times=p), obs=obs)

### simulate (correlated) random effects for intercepts and slopes
mu  <- c(0,0)
S   <- matrix(c(1, tau01, tau01, 1), nrow=2)
tau <- c(tau0, tau1)
S   <- diag(tau) %*% S %*% diag(tau)
U   <- mvrnorm(n, mu=mu, Sigma=S)

### simulate AR(1) errors and then the actual outcomes
dat$eij <- unlist(sapply(p, function(x) arima.sim(model=list(ar=ar.val), n=x) * sqrt(1-ar.val^2) * sigma)) dat$yij <- (beta0 + rep(U[,1], times=p)) + (beta1 + rep(U[,2], times=p)) * log(dat$obs) + dat$eij

### note: use arima.sim(model=list(ar=ar.val), n=x) * sqrt(1-ar.val^2) * sigma
### construction, so that the true error SD is equal to sigma

### create grouped data object
dat <- groupedData(yij ~ obs | id, data=dat)

### profile plots
plot(dat, pch=19, cex=.5)

### fit corresponding growth model
res <- lme(yij ~ log(obs), random = ~ log(obs) | id, correlation = corAR1(form = ~ 1 | id), data=dat)
summary(res)


A single run of this yields the following profile plots:

And the output from the model:

Linear mixed-effects model fit by REML
Data: dat
AIC      BIC    logLik
5726.028 5762.519 -2856.014

Random effects:
Formula: ~log(obs) | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr
(Intercept) 2.611384 (Intr)
log(obs)    2.092532 0.391
Residual    1.509075

Correlation Structure: AR(1)
Formula: ~1 | id
Parameter estimate(s):
Phi
0.3708575
Fixed effects: yij ~ log(obs)
Value Std.Error   DF  t-value p-value
(Intercept) 1.409415 0.2104311 1158  6.69775       0
log(obs)    6.076326 0.1601022 1158 37.95279       0
Correlation:
(Intr)
log(obs) 0.166

Standardized Within-Group Residuals:
Min          Q1         Med          Q3         Max
-2.58849482 -0.53571963  0.04011378  0.52296310  3.11959082

Number of Observations: 1359
Number of Groups: 200


So, the estimates are quite close to the actual parameter values (of course, the larger $n$ is, the better this will work). Maybe this gives you a starting point for tweaking the code to your needs.

• Many thanks. I'll go through this ASAP, and put up the means and a covariance matrix of my observed data as well. Nov 25, 2013 at 16:49
• That would be useful. Also useful would be an idea what kind of model you assume would provide a good fit to your data. That is essentially the basis for simulating the data. I just gave one simple example, but ultimately that needs to be tweaked based on your application. Nov 25, 2013 at 17:01