Simulating correlated variables in R

Simulate correlated variables in R by hand, with chol() and with mvrnorm(), and why a correlation set on the log scale is not the correlation of the counts.

A pilot survey counted two ground beetle species in pitfall traps at 40 grassland sites. The first species averaged a few beetles per trap, the second several more, and sites that were good for one tended to be good for the other: the Pearson correlation between the two columns of counts was about 0.6. The next step is a power analysis for a larger survey, and that needs simulated surveys that look like the pilot, with the two species correlated as they were in the field. How do you simulate two correlated variables, and how do you make the counts come out correlated at 0.6?

Power analysis by simulation in R builds the simulation loop, the approach Bolker (2008) teaches for ecological models, and Probability distributions in R: d, p, q and r explains rnorm() and rpois(). Neither makes two variables move together. Dozens of posts on this site do, usually with one line built around chol() or MASS::mvrnorm() and no comment. This post explains that line, then measures the trap that is easy to walk into the first time you simulate correlated counts.

The short answer. Draw independent standard normals and mix them: for two variables, x2 <- rho * z1 + sqrt(1 - rho^2) * z2; for any number, Z %*% chol(Sigma), with the Cholesky factor on the right of the matrix of draws, or let MASS::mvrnorm() do it. For counts, the usual route puts the correlation on a hidden log scale and draws Poisson counts from it. The counts then correlate less than the number you typed, so decide first which scale your target belongs to. If it is the correlation of the counts, solve for the log-scale value that produces it, and check that it can be produced at all.

All data here are simulated, and every seed is in the code.

library(ggplot2)

te_paper  <- "#f5f4ee"
te_ink    <- "#16241d"
te_body   <- "#2c3a31"
te_forest <- "#275139"
te_rust   <- "#b5534e"
te_gold   <- "#c9b458"
te_line   <- "#dad9ca"

theme_datasheet <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.background  = element_rect(fill = te_paper, colour = NA),
          panel.background = element_rect(fill = te_paper, colour = NA),
          panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
          panel.grid.minor = element_blank(),
          text             = element_text(colour = te_body),
          axis.text        = element_text(colour = te_body),
          strip.text       = element_text(colour = te_ink, face = "bold"))
}

Two correlated variables by hand

Start with two independent standard normals, z1 and z2. Keep z1 as the first variable and build the second from both: a share rho of z1 plus enough of z2 to bring the variance back to one.

set.seed(1015)
n_big <- 1e5
rho   <- 0.6

z1 <- rnorm(n_big)
z2 <- rnorm(n_big)
x1 <- z1
x2 <- rho * z1 + sqrt(1 - rho^2) * z2

round(c(sd_x1 = sd(x1), sd_x2 = sd(x2), correlation = cor(x1, x2)), 3)

      sd_x1       sd_x2 correlation 
      1.000       1.001       0.598 

The algebra is two lines. The variance of x2 is rho^2 + (1 - rho^2) = 1, because z1 and z2 are independent with variance one, and the covariance of x1 and x2 is rho, because only the z1 part is shared. Both standard deviations are one, so the covariance is also the correlation, and a hundred thousand draws give 0.598 against the 0.6 we asked for. To give the variables other means and spreads, multiply and add afterwards: 10 + 2 * x2 has correlation rho with x1 too, because shifting and rescaling does not change a correlation.

Any number of variables: chol()

With three or more variables the mixing weights are no longer obvious, and the Cholesky decomposition supplies them. Say three species have log-abundances that correlate at 0.6, 0.3 and 0.5 in pairs. Write the targets into a correlation matrix and hand it to chol():

sigma3 <- matrix(c(1.0, 0.6, 0.3,
                   0.6, 1.0, 0.5,
                   0.3, 0.5, 1.0), nrow = 3)
chol_upper <- chol(sigma3)
round(chol_upper, 3)

     [,1] [,2]  [,3]
[1,]    1  0.6 0.300
[2,]    0  0.8 0.400
[3,]    0  0.0 0.866

chol() returns an upper triangular matrix, zeros below the diagonal, and it is the factor R for which t(R) %*% R gives back sigma3. That orientation decides which side of the draws it goes on. Put independent standard normals in a matrix with one row per site and one column per variable, and multiply by R on the right:

set.seed(2015)
z3 <- matrix(rnorm(3 * n_big), ncol = 3)
x3 <- z3 %*% chol_upper
round(cor(x3), 3)

      [,1]  [,2]  [,3]
[1,] 1.000 0.602 0.302
[2,] 0.602 1.000 0.504
[3,] 0.302 0.504 1.000

The covariance of z3 %*% R is t(R) %*% R, which is sigma3, and the simulated correlations match the targets. The first column of R is (1, 0, 0), so the first variable is simply the first column of draws, and the second column of R holds the same two weights as the hand-made x2 above: the by-hand recipe is the two-variable Cholesky factor written out.

Put the factor on the other side, z3 %*% t(R), and the covariance becomes R %*% t(R), which is a different matrix:

x3_wrong <- z3 %*% t(chol_upper)
round(cor(x3_wrong), 3)

      [,1]  [,2]  [,3]
[1,] 1.000 0.559 0.253
[2,] 0.559 1.000 0.452
[3,] 0.253 0.452 1.000

round(diag(var(x3_wrong)), 3)

[1] 1.458 0.803 0.750

The targets 0.6, 0.3 and 0.5 come out as 0.559, 0.253 and 0.452, and the variances, which should all be one, run from 0.75 to 1.46. Nothing warns you. The same mistake with only two variables is no safer: R %*% t(R) then has correlation rho / sqrt(1 + rho^2), which is 0.514 for a target of 0.6. Some posts on this site draw one realisation at a time as a column, t(chol(Sigma)) %*% z; that is the correct product for a column, because it is the transpose of z %*% chol(Sigma) for a row.

The shortcut is mvrnorm() from the MASS package that accompanies Venables and Ripley (2002). It takes the means and the covariance matrix and returns one row per draw:

set.seed(3015)
x3_mass <- MASS::mvrnorm(n_big, mu = c(0, 0, 0), Sigma = sigma3)
round(cor(x3_mass), 3)

      [,1]  [,2]  [,3]
[1,] 1.000 0.599 0.303
[2,] 0.599 1.000 0.502
[3,] 0.303 0.502 1.000

Internally mvrnorm() uses an eigen decomposition of Sigma rather than chol(), so the individual draws differ from the Cholesky ones, but the covariance they are drawn from is the same. It is the line to use in practice; writing the Cholesky version once is what makes it readable.

Correlated counts: the correlation you set is not the one you get

Now the beetles. Counts are not normal, so the usual route, and the one behind most joint species models for counts, puts the correlation on a hidden scale. Each site gets a pair of correlated normal log-abundances with standard deviation s; each species’ expected count is its site mean times the exponential of its log-abundance; and the trap catches a Poisson count around that expectation. Subtracting s^2 / 2 inside the exponential keeps the average count at the site mean, because the mean of exp(x) for a normal x with mean zero and standard deviation s is exp(s^2 / 2).

simulate_counts <- function(n_sites, rho_latent, s, mu) {
  sigma  <- s^2 * matrix(c(1, rho_latent, rho_latent, 1), nrow = 2)
  latent <- matrix(rnorm(2 * n_sites), ncol = 2) %*% chol(sigma)
  expected <- cbind(mu[1] * exp(latent[, 1] - s^2 / 2),
                    mu[2] * exp(latent[, 2] - s^2 / 2))
  counts <- matrix(rpois(2 * n_sites, expected), ncol = 2)
  list(latent = latent, counts = counts)
}

mu_beetle <- c(3, 8)   # mean beetles per trap, species A and B
s_log     <- 1         # standard deviation of log-abundance between sites
n_sites   <- 40

set.seed(4015)
survey <- simulate_counts(n_sites, rho_latent = 0.6, s = s_log, mu = mu_beetle)
head(survey$counts)

     [,1] [,2]
[1,]    4    6
[2,]    2    3
[3,]    0    1
[4,]    0    4
[5,]   10   26
[6,]    2    3

round(c(latent = cor(survey$latent)[1, 2], counts = cor(survey$counts)[1, 2]), 3)

latent counts 
 0.593  0.416 

In this one simulated survey the log-abundances correlate at 0.593 and the counts at 0.416. One survey of 40 sites is noisy, so the chunk below asks what the counts correlate at in the long run: forty batches of fifty thousand sites each, with the correlation computed in every batch and averaged, so that the spread between batches gives a Monte Carlo standard error.

long_run <- function(rho_latent, s, mu, n_batch = 40, batch_size = 5e4) {
  per_batch <- replicate(n_batch, {
    sim <- simulate_counts(batch_size, rho_latent, s, mu)
    c(latent    = cor(sim$latent)[1, 2],
      lognormal = cor(exp(sim$latent))[1, 2],
      counts    = cor(sim$counts)[1, 2],
      var_a     = var(sim$counts[, 1]),
      var_b     = var(sim$counts[, 2]),
      cov_ab    = cov(sim$counts)[1, 2])
  })
  rbind(mean = rowMeans(per_batch),
        mc_se = apply(per_batch, 1, sd) / sqrt(n_batch))
}

set.seed(5015)
run_06 <- long_run(0.6, s_log, mu_beetle)
round(run_06, 3)

      latent lognormal counts  var_a   var_b cov_ab
mean   0.600     0.481  0.424 18.642 117.763 19.883
mc_se  0.001     0.002  0.002  0.102   0.847  0.137

The latent correlation is the 0.60 we set. The exponentials of the log-abundances already correlate at only 0.481, and the counts at 0.424. A power analysis run with rho_latent = 0.6 would simulate beetles that are much less tied to each other than the pilot’s, and every power figure from it would answer a different question.

None of this needs a simulation to predict. Two short results give the long-run values. The count formulas are in Aitchison and Ho (1989); the first result is the same calculation without the Poisson term. For the exponentials: the sum of the two log-abundances is normal with variance 2 s^2 (1 + rho), so the mean of the product of the exponentials is exp(s^2 (1 + rho)), and dividing the covariance by the variance gives

cor(exp(x1), exp(x2)) = (exp(rho * s^2) - 1) / (exp(s^2) - 1)

For the counts: given the expected counts, the two Poisson draws are independent, so the Poisson noise adds its mean to each variance and nothing to the covariance,

var(y)       = mu + mu^2 * (exp(s^2) - 1)
cov(y1, y2)  = mu1 * mu2 * (exp(rho * s^2) - 1)

and the count correlation is the covariance over the product of the two standard deviations. The Poisson term mu in each variance is the second loss: it is largest, relative to the rest, for the rarer species.

count_cor <- function(rho_latent, s, mu) {
  cov_ab <- mu[1] * mu[2] * (exp(rho_latent * s^2) - 1)
  var_ab <- mu + mu^2 * (exp(s^2) - 1)
  cov_ab / sqrt(var_ab[1] * var_ab[2])
}
closed <- c(lognormal = (exp(0.6 * s_log^2) - 1) / (exp(s_log^2) - 1),
            counts    = count_cor(0.6, s_log, mu_beetle),
            var_a     = mu_beetle[1] + mu_beetle[1]^2 * (exp(s_log^2) - 1),
            var_b     = mu_beetle[2] + mu_beetle[2]^2 * (exp(s_log^2) - 1),
            cov_ab    = prod(mu_beetle) * (exp(0.6 * s_log^2) - 1))
z_scores <- (run_06["mean", names(closed)] - closed) / run_06["mc_se", names(closed)]
round(rbind(closed = closed, simulated = run_06["mean", names(closed)], z = z_scores), 3)

          lognormal counts  var_a   var_b cov_ab
closed        0.478  0.423 18.465 117.970 19.731
simulated     0.481  0.424 18.642 117.763 19.883
z             1.332  0.952  1.737  -0.244  1.108

Each simulated value is within 1.8 Monte Carlo standard errors of its formula: the loss from 0.60 to 0.423 is arithmetic, not a finding of the simulation. It depends only on rho, s and the two means.

The pilot’s own 0.6 has a second weakness: it comes from 40 sites. Here is the spread of the count correlation across four thousand simulated 40-site surveys, all with the same log-scale correlation of 0.6:

set.seed(6015)
cor_40 <- replicate(4000, {
  cor(simulate_counts(n_sites, 0.6, s_log, mu_beetle)$counts)[1, 2]
})
round(quantile(cor_40, c(0.10, 0.50, 0.90)), 3)

  10%   50%   90% 
0.174 0.420 0.687 

The long-run count correlation of this design is 0.423, yet one survey in ten reports less than 0.174 and one in ten more than 0.687, with the middle survey at 0.420; 19.9 per cent of these surveys report 0.6 or more. A pilot value of 0.6 from 40 sites can therefore come from a process whose long-run count correlation is well below 0.6, which is worth remembering whichever scale you settle on.

The fix: decide which scale the target lives on

If your 0.6 is a correlation on the log scale, for example a residual correlation from a fitted Poisson or negative binomial mixed model or a joint species model, then rho_latent = 0.6 is already right and the counts should correlate less. If it is the Pearson correlation of the raw counts, as in the pilot, solve for the log-scale value that produces it. count_cor() rises steadily with rho_latent, so uniroot() can search between 0 and 1:

target <- 0.6
rho_needed <- uniroot(function(r) count_cor(r, s_log, mu_beetle) - target,
                      interval = c(0, 1))$root
rho_needed

[1] 0.7732494

set.seed(7015)
run_fixed <- long_run(rho_needed, s_log, mu_beetle)
round(run_fixed[, c("latent", "counts")], 3)

      latent counts
mean   0.773  0.597
mc_se  0.000  0.002

A log-scale correlation of 0.773 gives counts that correlate at 0.597 in the long run (Monte Carlo standard error 0.002), 1.5 standard errors from the target. That is the number to put into the power analysis.

There is a limit, and it is better to see it before calling uniroot(). Even with rho_latent = 1, when both species follow exactly the same log-abundance, the Poisson noise keeps the counts apart, so the count correlation has a ceiling. The ceiling depends on how much the sites differ, s, and on how large the counts are: what matters is the between-site variance mu^2 * (exp(s^2) - 1) against the Poisson variance mu, so rarer species and more similar sites both lower it:

ceiling_s1  <- count_cor(1, s = 1,   mu = mu_beetle)
ceiling_s05 <- count_cor(1, s = 0.5, mu = mu_beetle)
ceiling_s05_x10 <- count_cor(1, s = 0.5, mu = 10 * mu_beetle)
round(c(ceiling_s1 = ceiling_s1, ceiling_s05 = ceiling_s05, ceiling_s05_x10 = ceiling_s05_x10), 3)

     ceiling_s1     ceiling_s05 ceiling_s05_x10 
          0.884           0.565           0.926 

no_root <- tryCatch(
  uniroot(function(r) count_cor(r, 0.5, mu_beetle) - target, interval = c(0, 1)),
  error = function(e) conditionMessage(e))
no_root

[1] "f() values at end points not of opposite sign"

With s = 1 the counts can reach 0.884, so 0.6 is attainable. With s = 0.5, sites that differ less, and these small means, the ceiling is 0.565: no log-scale correlation gives counts correlated at 0.6, and uniroot() stops with the message above, because the function it searches is below zero at both ends of the interval. That error is the model telling you the target is impossible, not a numerical problem to work around. If you meet it, either the pilot correlation is higher than this kind of process can produce (and with 40 sites, the previous section shows how easily that happens), or the counts vary between sites more than your s says. The same s = 0.5 with means ten times larger, 30 and 80 beetles per trap, lifts the ceiling to 0.926, because the Poisson noise then matters much less.

Figure 1: Long-run correlation of two simulated Poisson counts (site means 3 and 8) against the correlation set on the log scale, for a between-site log-scale standard deviation of 1 and of 0.5, from the closed form. The dashed line is where the two would be equal. Points are the simulated long-run values at a log-scale correlation of 0.6 and at the value solved for a count correlation of 0.6 (dotted line). Each curve ends at its ceiling, the count correlation when the log-scale correlation is 1.

What to check in your own data

Check the orientation once, on a large draw: round(cor(Z %*% chol(Sigma)), 2) should reproduce your target matrix, and the variances should be the diagonal of Sigma. If they are not, the factor is on the wrong side.

Write down which scale your target correlation was measured on. A residual correlation from a mixed model or a joint species model lives on the link scale; a correlation of raw counts lives on the count scale. Simulating one as if it were the other is the mistake measured above. Checking a joint species distribution model goes the other way for presence-absence data, reading a latent correlation back from a two by two table with the tetrachoric correlation.

Estimate s from the pilot before you solve for the log-scale correlation. For each species, the variance formula above rearranges to s^2 = log(1 + (var(y) - mean(y)) / mean(y)^2); if the variance is not above the mean, the counts show no extra spread between sites, s is near zero and so is any count correlation this model can make. If the two species give different values, the covariance becomes mu1 * mu2 * (exp(rho * s1 * s2) - 1) with s1^2 and s2^2 in the two variances.

Compare the ceiling, count_cor(1, s, mu), with your target before calling uniroot(). Then simulate a large batch at the solved value and check the counts hit the target. Testing your analysis code shows how to keep a check like that as a test that runs every time the simulator changes.

Honest limits

The count model here is Poisson on a lognormal, with one s for both species; a negative binomial simulator, or a zero-inflated one, has other formulas and another ceiling, although the same uniroot() approach works once the formula is written down, or with a simulated long run in its place. Only positive correlations were shown, and a negative target meets a floor as well: at s = 1 and these two means, count_cor(-1, 1, mu) is -0.325, so the counts cannot correlate more negatively than that. The whole post is about the Pearson correlation of two counts, which is only one summary of how two species go together; Copulas for dependent ecological data measures the limits that any pair of marginal distributions puts on it, and Co-occurrence is not interaction explains why a correlation between two species, on either scale, says nothing on its own about why they go together. Finally, the pilot correlation is itself uncertain, as the 40-site spread showed, so a power analysis that is sensitive to it should be run over a range of plausible correlations, not at one solved value, and over a range of s, which a 40-site pilot also pins down only loosely.

References

Aitchison J, Ho CH 1989 Biometrika 76(4):643-653 (10.1093/biomet/76.4.643)

Bolker BM 2008 Ecological Models and Data in R (ISBN 9780691125220)

Venables WN, Ripley BD 2002 Modern Applied Statistics with S, 4th edition (ISBN 9780387954578)

Related tutorials

• Power analysis by simulation in R
• Probability distributions in R: d, p, q and r
• Copulas for dependent ecological data
• Checking a joint species distribution model