Repeated measures correlation as mixed model

Question

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}$.

Answer 1

AFAIK a full (restricted) maximum likelihood estimation approach for this model is not possible with neither nlme nor lme4. But you could, for example, do it under a Bayesian approach using either STAN or JAGS.

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)
 lm(y ~ 0 + fitted(fm), data = df)