# Repeated measures mixed model correlation between measurements

I am interested in looking at the correlation between two types of heart function measurements over time. test1 can be considered the true value. The true value (test1) changes over time, but not much whereas test2 changes a lot.

The hypothesis is that the two tests have a good correlation in the beginning, and then the new type of measurement (test2) loses its ability to follow test1, but I need to quantify it somehow.

I was thinking of some way of measuring if test2 becomes less correlated with test1 at each time point. Maybe using the correlation at baseline as a reference. I have four distinct measuring time points (factor), my testtype (2 types) as factor and my output (value) is continuous on both test 1 and 2 but very different scales. Finally, I would use the subject (ID) for the random effect. So my idea was to do a mixed model like this (using R):

lme(value ~ time * testtype, random = ~ 1 | subject, data = df, na.action = na.omit)


My output would then be estimates for test overall, time1 vs time2, time1 vs time3 etc. And finally test1:time2, test1:time3, etc.

My question is can I use this method to evaluate if the correlation between tests is different at different time points?

And how do I turn the estimates with standard error into something meaningful? Something I can interpret as "the correlation between tests at time2 was xx lower than at time1."

EDIT: Using Dimitris model I get a correlation structure like this:

             1         2         3         4         5         6         7         8
1 1.0000000 0.2868521 0.2868521 0.2868521 0.1973987 0.1973987 0.1973987 0.1973987
2 0.2868521 1.0000000 0.2868521 0.2868521 0.1973987 0.1973987 0.1973987 0.1973987
3 0.2868521 0.2868521 1.0000000 0.2868521 0.1973987 0.1973987 0.1973987 0.1973987
4 0.2868521 0.2868521 0.2868521 1.0000000 0.1973987 0.1973987 0.1973987 0.1973987
5 0.1973987 0.1973987 0.1973987 0.1973987 1.0000000 0.6351789 0.6351789 0.6351789
6 0.1973987 0.1973987 0.1973987 0.1973987 0.6351789 1.0000000 0.6351789 0.6351789
7 0.1973987 0.1973987 0.1973987 0.1973987 0.6351789 0.6351789 1.0000000 0.6351789
8 0.1973987 0.1973987 0.1973987 0.1973987 0.6351789 0.6351789 0.6351789 1.0000000


EDIT 2 For Dimitri - Here is the output from getVarCov(fm)

Marginal variance covariance matrix
1       2       3       4       5       6       7
1 133.240  38.220  38.220  43.692  43.692  43.692  43.692
2  38.220 133.240  38.220  43.692  43.692  43.692  43.692
3  38.220  38.220 133.240  43.692  43.692  43.692  43.692
4  43.692  43.692  43.692 367.700 233.550 233.550 233.550
5  43.692  43.692  43.692 233.550 367.700 233.550 233.550
6  43.692  43.692  43.692 233.550 233.550 367.700 233.550
7  43.692  43.692  43.692 233.550 233.550 233.550 367.700
Standard Deviations: 11.543 11.543 11.543 19.175 19.175 19.175 19.175

• It is doable by fitting mixed model. But it is very complicated. You need to find a person who can write down $Y, X, Z, G, R$ matrices which are used in mixed model before fitting the mixed model. – user158565 Oct 18 '18 at 2:09
• Do the two tests measure the same thing? I.e., is the value from test1 the same outcome as the value from test2? – Dimitris Rizopoulos Oct 18 '18 at 6:23
• The two tests represent the output from the heart, but not directly so one is mililiter blood pumped out and the other is difference in heart chamber size at maximum and minimum - so a percentage of how much the heart contracts. contracts. – User LG Oct 18 '18 at 9:17
• @a_statistician - Would it be better to split my data into 4 so I only analyze 1 timepoint at a time? – User LG Oct 18 '18 at 9:29
• No. Should fit one model to answer all of questions in your case. – user158565 Oct 19 '18 at 4:10

Based on the comments above, it seems that the two tests are measuring two different related outcomes. Hence, it would make more sense to fit a bivariate mixed model, in which you account that the error variance is different for the two tests. You could fit this model using something like this:

fm <- lme(value ~ 0 + testtype + testtype:time, data = df, na.action = na.omit,
random = ~ 0 + testtype + testtype:time | subject,
weights = varIdent(form = ~ 1 | testtype))

summary(fm)


The coefficients you will obtain will be per test-type. Regarding the correlations, you can use

V <- getVarCov(fm, type = "marginal", individual = 1)
R <- cov2cor(V[])
R


This will give you the marginal covariance matrix V of the two outcomes per test (i.e., it will be an $$8 \times 8$$ matrix), which we turn to the corresponding correlation matrix R. The two block-diagonal parts (i.e., R[1:4, 1:4] and R[5:8, 5:8]) will be the correlations over time within the same test. The part R[1:4, 5:8] will be the correlation between the tests over time. Note: you will need to select an individual in the individual argument that has all eight measurements.

However, just my 2c, I don't think correlations are easier to interpret than the coefficients.

• just to be sure, when you write testtype:typeis that testtype:time? – User LG Oct 18 '18 at 14:49
• Yes, my mistake. – Dimitris Rizopoulos Oct 18 '18 at 15:08
• So when I look at the output I first get the fixed effect values for test 1 and test 2 and then for test1:time1, test2:time1, test1:time2etc. So how do you understand these fixed effects? – User LG Oct 18 '18 at 17:33
• I added the correlation structure in the original post. Because when you say that block R[1:4,5:8] is the correlation between tests over time I am puzzled because the correlation is constant over time. So what if the correlation changes between time 3 and 4 but not 1 and 2? – User LG Oct 18 '18 at 17:41
• Is the model fitted with random = ~ 1 | subject or random = ~ 0 + testtype + testtype:time | subject? Could you perhaps also show what the output of getVarCov(fm) is, i.e., without specifying the type argument. – Dimitris Rizopoulos Oct 18 '18 at 18:17