# Repeated measures correlation as mixed model

I'm writing this question to better understand how to set up a mixed model to estimate a linear relationship between two variables from repeated observations.

Suppose the variables $$X$$ and $$Y$$ are observed on $$n$$ individuals at $$k$$ time points each. The observations $$x_{i,t}$$ change linearly over time according to a random intercept model: $$x_{i,t} = x_{i,0}+\beta\, t+\varepsilon_{i,t}^x,\quad x_{i,0}\sim \mathcal{N}(\mu,\sigma^2),\quad i=1,\ldots,n$$ The observations $$y_{i,t}$$ are linearly related to the $$x_{i,t}$$: $$y_{i,t}=\gamma ( x_{i,0}+\beta\, t)+\varepsilon_{i,t}^y,$$ where, for simplicity, I assumed that $$\gamma$$ is independent of $$t$$.

If I haven't made a mistake a simple simulation of the model can be achieved in R by the following code.

set.seed(1)

n <- 30 #number of subjects
t <- c(0,1,3,12) # time points
b <- .2 # coefficient beta
c <- -2 # coefficient gamma

x0 <- rnorm(n,10,1) # random intercepts

x <- unlist(lapply(t, function(t){x0+b*t}))
y <- c*x

df<-data.frame(x=x+rnorm(4*n,0,1), y=y+rnorm(4*n,0,1), t=rep(t,each=n), ID=rep(seq.int(1,n),4))


Plotting the data manifestly shows the linear relationship at each of the four time points.

require(ggplot2)
ggplot(df,aes(x=x,y=y,col=as.factor(t)))+
geom_point()+
geom_smooth(method='lm', se = F)


Consequently, fitting linear models for each time point separately yields estimates for $$\gamma$$.

lapply(t,function(tp)summary(lm(y~x,d=df, subset = t==tp)))


A blind pooling of the data in a multiple regression

summary(lm(y~x+t,df))


reveals highly significant effects of $$x$$ and time

 Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.18658    1.24738  -6.563 1.51e-09 ***
x           -1.16335    0.12010  -9.686  < 2e-16 ***
t           -0.22201    0.03772  -5.885 3.87e-08 ***


but fails to take into account the dependence between observations from the same subject. A simple mixed-effects model with random effects for subjects

require(nlme)
summary(lme(y~x+t, random=~1|ID, d=df))


yields highly variable estimates for $$\gamma$$

Fixed effects: y ~ x + t
Value Std.Error DF    t-value p-value
(Intercept) -17.333357 1.2615782 88 -13.739424  0.0000
x            -0.270954 0.1183723 88  -2.289001  0.0245
t            -0.364641 0.0278553 88 -13.090550  0.0000


which, despite the larger data basis, are considerably less precise than the single-time estimates.

Question How can one formulate a model to pool the data from the separate time points to obtain an even better estimate for $$\gamma$$?

I am aware of the rmcorr package, which claims to answer the question above. However, having looked the code of the package, it seems to compute a linear model for the pooled data without accounting for correlations. If the package is indeed applicable, I wold be grateful for an explanation.

Edit 16 May 2019:

• Cleared up notational confusion between $$\beta$$ and $$\gamma$$.
• Removed error $$\varepsilon_{i,t}^x$$ from the model equation for $$y_{i,t}$$.

As a alternative you could do a two-stage approach in which you first fit a mixed model for x, and the you plug-in the subject-specific fitted values in a linear regression for y. E.g.,
 fm <- lme(x ~ t, data = df, random = ~ 1 | ID)

• Thanks for your comment pointing out inconsistencies in my question. I would like indeed to recover $\gamma$, the relation between $y$ and $x$. The error $\varepsilon_{i,t}^x$ is thought to be observation error in $x$ but is assumed not to impact on $y_{i,t}$, hence the R code is correct. I edited my question. – Eckhard May 16 at 5:37