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.


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.

  geom_smooth(method='lm', se = F)

enter image description here

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


reveals highly significant effects of $x$ and time

            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

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

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)
  • $\begingroup$ 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. $\endgroup$
    – Eckhard
    May 16 '19 at 5:37
  • $\begingroup$ That's an interesting approach, thank you! How would I set up such a model in STAN, presumable using the rstan interface? $\endgroup$
    – Eckhard
    May 16 '19 at 7:50
  • $\begingroup$ Yes, you could do it with rstan. $\endgroup$ May 16 '19 at 8:21
  • $\begingroup$ Thanks. Do you happen to know of a reference or book where the two-stage approach you suggest is explained and its theoretical properties are analysed? $\endgroup$
    – Eckhard
    May 16 '19 at 22:00
  • $\begingroup$ I don’t have a reference for this particular two-stage combination, but in general this is a solution that is used to fit joint models with standard software procedures. Hence, search for “joint models two-stage approach”. $\endgroup$ May 17 '19 at 5:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.