Simulating correlated data for time points

In this post, we will consider the simulation of one study with correlated data. This process can then fit into the multi-trial simulation of a whole evidence base, as discussed in our previous post.

We will assume there are two arms that participants are randomly allocated to. This is a special case because the randomisation enforces zero correlation between the arms. However, when we add multiple measurements within each participant, there will be non-zero correlation. This mixture means we need to proceed with caution, but that’s how we can learn a great deal from simulation.

Let’s consider the multiple measurements to be time points, like baseline, interim and end-of-intervention. The same logic would apply to situations with “paired” data, e.g. where one kidney, eye, etc is treated at a time with some localised intervention, though some of the numbers will differ. Identical twins are another example (that looms large in the public imagination but is more limited in reality). You might also consider multiple outcomes or subscales of an outcome scale, for example the domains of a quality of life scale will be correlated.

• Population correlation between randomised arms: \( \rho = 0 \)
• Population correlation in the same arm but between time points: \( \rho > 0 \)
• Population correlation between paired arms at the same time point: \( \rho \geq 0 \)

Note that we can comfortably rule out negative correlations in this setting.

Multivariate normal distribution

The multiple measurements present a vector of numbers for each participant. That vector can be modelled as being drawn from a multivariate distribution. (The terms multivariate and multivariable are unfortunately used differently depending on countries, academic fields, software…). The easiest by far to grapple with is the multivariate normal (MVN). In fact, it is so much harder to do anything else that it is quite popular in various aspects of statistics to model a latent MVN vector that is then somehow transformed into the observed data.

The MVN distribution’s parameters are a vector of means and a covariance matrix. Having read Chapter 1, you will now be ready to bore your colleagues with the definition that this is a probability density function.

To get mathematically rigorous about it, $$f_{\mathrm{MVN}} : \mathbb{R}^p \rightarrow \mathbb{R}^{+}$$ which returns the derivative of the cumulative probability with respect to the variable(s).

$$\mathbf{X} \sim \mathrm{MVN}(\mathbf{\mu}, \mathbf{V}) \Leftrightarrow \nabla P(\mathbf{x}) = f_{\mathrm{MVN}}(\mathbf{x} | \mathbf{\mu}, \mathbf{V})$$

Our experience is that it helps to have these very precise definitions in mind, but we exercise a little caution in the level of technical difficulty that we present to the client / colleague / committee. It is unfortunate to lose their interest at the start of a presentation.

The element of the covariance matrix in row \(i\) and column \(j\) is the correlation (between the \(i\)th and the \(j\)th variable), multiplied by the SD of variable \(i\), multiplied by the SD of variable \(j\):

$$\mathbf{V}_{ij} = \rho_{ij} \sigma_i \sigma_j$$

Consider what happens in the diagonal, where we compare a variable with itself. The correlation must be 1, and so:

$$\mathbf{V}_{ii} = \sigma^2_i$$

The diagonal contains the variances of each of the variables. The matrix is sometimes called the variance-covariance matrix for this reason. We don’t like that term because the variances are just covariances of each variable with itself, not something distinct that has been bolted on.

This matrix also has to be square and symmetric. For \(p\) variables, there are \(p(p+1) / 2\) elements to be filled: the diagonal and the upper (or lower) triangle. It’s important to realise that, even within those constraints, only some combinations of numbers will be valid. If variables 1 and 2 are strongly correlated, and 1 and 3 likewise, then 2 and 3 must be quite strongly correlated. We’ll leave aside the details, but it is worth finding out more if you are mathematically inclined.

You can skip this paragraph and the next formula if you are not comfortable with linear algebra. But if you are, it may help you clarify the idea. Suppose you have the SDs in a (column) vector \(\mathbf{\sigma}\), and the correlations in a matrix \(\mathbf{\rho}\). First, create a matrix with the standard deviations on the diagonal, zeros elsewhere: \(\mathbf{D} = \mathrm{diag}(\mathbf{\sigma}) \). Then:

$$\mathbf{V} = \mathbf{D} \mathbf{V} \mathbf{D}$$

Most people seem to find it easier to think about the correlation matrix and then multiply in the standard deviations, so let’s proceed in that way. All the diagonal elements of the correlation matrix are 1 by definition, so there are \(p(p-1) / 2\) elements to be decided (the variance vector will supply the other \(p\) numbers).

If we have only two time points (and two randomised arms), we can’t trip ourselves up with impossible matrices: there is only one correlation to be chosen, so it can’t contradict any other correlations.

Comparing the arms at baseline, the correlation is zero, because t:

Considering only the time points in the control group (where we assume no effect, placebo, nocebo, dropout, or regression to the mean), there will be a positive correlation. Let’s call it 0.7 for our simulation.

In the intervention group, the mean may be different, but the correlation might also be 0.7. When we put it all together, there are four variables: control at baseline, intervention at baseline, control at end, intervention at end:

We have omitted the lower triangular portion of the matrix for clarity. You will spot three correlations that are not defined so far. The randomisation of the arms saves us here because any comparison across arms is arbitrary and has no pairing and will, at a population level, have zero correlation.

We can enter this in R with the matrix() function:

cr <- matrix(c(1, 0, 0.7, 0,
               0, 1, 0, 0.7,
               0.7, 0, 1, 0,
               0, 0.7, 0, 1),
             ncol=4)
eigen(cr)

The eigen() function is a fast way to check that you have a valid correlation or covariance matrix. It outputs eigenvalues and eigenvectors. The eigenvalues should all be positive or zero:

> eigen(cr)
eigen() decomposition
$values
[1] 1.7 1.7 0.3 0.3

$vectors
          [,1]       [,2]       [,3]       [,4]
[1,] 0.0000000 -0.7071068  0.0000000  0.7071068
[2,] 0.7071068  0.0000000 -0.7071068  0.0000000
[3,] 0.0000000 -0.7071068  0.0000000 -0.7071068
[4,] 0.7071068  0.0000000  0.7071068  0.0000000

Now, we can multiply in the standard deviations. We’ll imagine the SD in the intervention arm is 10 (at both time points) and in the control arm is 9, just to help show what’s going on. You can either use a nested loop over the rows and columns, remembering that the order is int_base, ctl_base, int_end, ctl_end:

sds <- c(10, 9, 10, 9)
cv <- matrix(0,nrow=4,ncol=4)
for(i in 1:4) {
  for(j in 1:4) {
    cv[i,j] <- cr[i,j]*sds[i]*sds[j]
  }
}

or by matrix multiplication:

sdiag <- diag(sds)
cv <- sdiag %*% cr %*% sdiag

Either way, you should get this:

> cv
     [,1] [,2] [,3] [,4]
[1,]  100  0.0   70  0.0
[2,]    0 81.0    0 56.7
[3,]   70  0.0  100  0.0
[4,]    0 56.7    0 81.0

We are ready to simulate! In R, the mvtnorm package is useful for MVN (and multivariate t). Once again, the order of columns will be int_base, ctl_base, int_end, ctl_end.

library(mvtnorm)

one_study <- rmvnorm(n=50,
                     mean=c(100, 100, 90, 100),
                     sigma=cv)

We can get means and SDs for this one pseudo-study:

> apply(one_study,2,function(z){c(mean(z),sd(z))})
          [,1]      [,2]     [,3]      [,4]
[1,] 102.18162 101.14376 91.99543 100.00975
[2,]  10.29121  10.06034 11.29704  10.38248

Baseline + change

We can also obtain MVN data by modelling baseline data from normal distributions, and then adding change, also modelled from normal distributions. The change will be (in our model, anyway) unrelated to the baseline value. Bear in mind that this is a special normal (Gaussian) attribute. In general, we can’t add random variates with some distribution together and hope to get the same distribution as output.

The data after intervention end is simply the sum of the baseline and the change for each participant (subject). To see what the mean and variance will be at end, refer back to Formula 1.12 in the book, because change and baseline are uncorrelated.

You might prefer to work with baseline + change rather than simultaneously drawing all the data from a multivariate distribution. In fact, it has some advantages in that it is easier to add in complications, such as specifying a baseline-change correlation. To do so, you need the baseline and the change to come from a bivariate normal distribution. Let’s make them uncorrelated, with equal SDs:

int_data <- rmvnorm(n=50,
                    mean=c(100, -10),
                    sigma=100*diag(2))
ctl_data <- rmvnorm(n=50,
                    mean=c(100, 0),
                    sigma=81*diag(2))
int_base <- int_data[,1]
int_end <- apply(int_data,1,sum)
ctl_base <- ctl_data[,1]
ctl_end <- apply(ctl_data,1,sum)
mydata<- cbind(int_base, ctl_base, int_end, ctl_end)
apply(mydata,2,function(z){c(mean(z),sd(z))})
cor(mydata)

Leading to stats like these:

> apply(mydata,2,function(z){c(mean(z),sd(z))})
      int_base  ctl_base  int_end  ctl_end
[1,] 100.78602 99.026101 89.48192 96.98182
[2,]  10.67414  9.568115 15.01123 13.27851
> cor(mydata)
             int_base   ctl_base     int_end      ctl_end
int_base  1.000000000 -0.1226163  0.70698716 -0.009241648
ctl_base -0.122616345  1.0000000 -0.11991633  0.755156203
int_end   0.706987157 -0.1199163  1.00000000  0.097531493
ctl_end  -0.009241648  0.7551562  0.09753149  1.000000000

Two points are worth thinking about here, but start to go beyond the scope of this post:

• The end-of-intervention SDs are inflated, and they were the same as baseline before, when we used the MVN simulation. If you want them to stay the same, you will need a negative correlation between baseline and change. This is not unrealistic: suppose a high value is a good outcome of treatment; participants who are already doing well may plausibly benefit less.
• Because the SDs of baseline and change are the same, the correlation between baseline and end is around \( \sqrt{2} \).

Paired data

How does this calculation change with paired data? There would be a correlation between intervention baseline and control baseline data. Then, to add a correlated change over the top of that, the MVN approach starts to be more attractive.

Even more time points

You can simulate studies with three or more time points too. You just have to be careful about the correlation structure over all the time points. Here are two classic assumptions that you might want to consider, and which you might recognise if you have studied the modelling method called generalized estimating equations (GEE):

t <- 5

# exchangeable correlation
cr <- 0.3*diag(t)
cr <- cr+0.7
print(cr)
eigen(cr)

# autoregressive AR(1) correlation
rho <- 0.7
cr <- matrix(0, nrow=t, ncol=t)
for(i in 1:t) {
  for(j in 1:t) {
    cr[i,j] <- cr[i,j] + rho^(abs(i-j))
  }
}
print(cr)
eigen(cr)

There are a couple of tricks in the coding here to keep it compact while also ensuring the diagonal elements are ones.

It is probably sensible to base your correlation structure on anything you can glean from other studies with similar population, intervention and time points.

Conclusion

As we saw in the previous post, we cna then extract just the summary stats for the trial in question. This can then be looped over \( m \) trials to provide an “evidence base” of study statistics, and to test any meta-analytic procedure that builds on it.

If you want to read more about matrices in the context of statistical analysis (and who does not?), we recommend the book Matrix Algebra by James Gentle (published by Springer).