Estimating and testing correlation of longitudinal random variables

Question

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.

Answer 1

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)

2. 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.

3. 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.