Estimating and testing correlation of longitudinal random variables

Each patient (indexed by $$i$$) contains multiple measurements of two variables $$X_{i,t}$$ and $$Y_{i,t}$$ over time $$t=1, \dots, T$$. For each time point $$t'$$, assume the correlation $$\mathrm{cor}(X_{i,t'}, Y_{i,t'})$$ does not depend on $$t'$$, making the correlation constant over time. I don't know the exact nature of the dependence within each patient over time, but let's assume it's either (1) a simple AR(1) structure or (2) exchangeable. Note, I don't want to assume any directionality in the relationship between $$X_{i,t}$$ and $$Y_{i,t}$$.

What's an estimator for this between-patient correlation? Is there a package in R to compute this correlation? Further, what's a test for whether this correlation is zero or not?

Of course, one way to estimate the correlation is by throwing out all but the first observation for each patient then calculating the usual Pearson correlation; however, this loses a lot of information.

• Please clarify which type of dependence you are interested in: within patient (over all time periods, i correlations), across patients for a single time period (t correlations) or across patients and time (one correlation). Next, Pearson isn't the only measure of dependence, there are many more possible metrics, e.g., see Clark, A Comparison of Correlation Measures, m-clark.github.io/docs/CorrelationComparison.pdf
– user234562
Commented Aug 12, 2020 at 12:50
• @user332577 I'm happy to assume that Pearson correlation is appropriate. The correlation I mean is between patients and it is constant over time. Commented Aug 12, 2020 at 12:54
• I am still not sure what you want. Your language about constant correlation over time is what is throwing me off. Couldn't you just calculate each patient's mean x and mean y, and then run a correlation on those two variables? Commented Aug 12, 2020 at 15:14
• @ErikRuzek Thanks for the reply! That's possible, but presumably the dependence should be taken into account, especially for later testing (which I inadvertently didn't include in the question--including now) Commented Aug 12, 2020 at 15:15
• Ok, that is helpful. I'll add an answer to see if that gets you what you want. Commented Aug 12, 2020 at 15:17

Since you want the between-patient correlation, you have a few options. The first of which is very simple and as you pointed out in your comment, doesn't take into account the correlation in $$x$$ and $$y$$. The second and third ways do this in a more principled way, but require more work:

1. You could calculate the mean of x and the mean of y and then run the correlation on the means. Sorry for the code below, but the editor is doing weird things.

library(dplyr)

dat <- dat %>% group_by(patient) %>% mutate(mn_x=mean(x), mn_y=mean(y)) %>% ungroup()

with(dat, correlation(mn_x, mn_y)

1. You can get the correlation corrected for the repeated measures sampling by running a multivariate multilevel model with both $$x$$ and $$y$$ as the outcomes. This can be done in nlme and also in Bayesian packages brms and MCMCglmm. This is probably the gold standard approach.

In brms(), such a model would have two outcomes that you identify using the mvbind code, e.g.:

m_mv <- brm(mvbind(y, x) ~ 1 + (1|c|patient), data=dat, chains=2, cores=2)

The c in the random effect term (1|c|patient) tells brm to estimate correlated random effects across all patients. In the output, you will see something like the following:

Group-Level Effects:
~patient (Number of levels: 200)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(y_Intercept)                     0.48      0.05     0.39     0.59 1.00      752
sd(x_Intercept)                       0.25      0.07     0.11     0.39 1.01      409
cor(y_Intercept,x_Intercept)        -0.51   0.22    -0.92    -0.09 1.00      582


The cor parameter is what you are interested. It is the latent correlation between patient values of y and x taking into account repeated measures.

1. You could run separate "empty" lmer models for $$x$$ and $$y$$, get the empirical Bayes estimates for the group means of each, reflate them (see code from Mark Lai here), and then calculate the correlation. See this paper by Leckie et al. from 2018. If done correctly this will reproduce the latent correlation for the intercepts of $$x$$ and $$y$$ reported in brms, MCMCglmm, or nlme. This will be more error prone than option 2 simply because of the coding involved.
• Thanks for the answer! Could you include code for approach 2 so that I can better see what you mean? Commented Aug 12, 2020 at 16:42
• The link to the brms package vignette has all the code necessary to run the model, but I added an example anyway. cran.r-project.org/web/packages/brms/vignettes/… Commented Aug 12, 2020 at 18:27
• Thank you, this is very helpful. This brm approach works for an exchangeable error structure; how can it be adjusted to handle an AR(1) dependence structure within a patient? I'm not sure a multilevel approach can handle that. Commented Aug 13, 2020 at 2:37
• I'm pretty sure you can set up autocorrelation of the residuals in any brms model using ar(p = 1), multilevel or otherwise. Are you sure you need an AR(1) structure? The random intercept helps deal with this problem. Commented Aug 13, 2020 at 17:50
• Thank you, this fully explains how to estimate the correlation in question! Commented Aug 17, 2020 at 14:12