[R-sig-ME] Modeling error covariance in linear mixed models with non-default correlation structures

[Pomponio, Raymond]

Hello,

I'm trying to fit a linear mixed model that accounts for correlated errors at the subject level. This has been described as modeling the "R" matrix in a model where the overall variance of Y can be expressed as Var(Y) = ZGZ' + R.

The nlme package offers the ability to model covariance structures through its corClasses feature. I'm capable of fitting a model with first order auto-regressive covariance structure as in the following example:

  > library(nlme)
  > dat <- data.frame(
     id      = c(rep(1, 18), rep(2, 14), rep(3, 12)),
     rec_num = c(1:18, 1:14, 1:12),
     day     = c(rep(1:9, each=2), rep(1:7, each=2), rep(1:6, each=2)),
     am_pm   = c("am", "pm"),
     y       = c(rnorm(18, mean=2), rnorm(14, mean=1), rnorm(12, mean=0))
     )
  > lme(
     y ~ 1,
     data        = dat,
     random      = ~ 1 | id,
     correlation = corAR1(form = ~ rec_num |id),
     control     = lmeControl(opt="optim")
     )


Each row of the above data.frame corresponds to one recorded observation, with multiple observations per id. Measurements were taken twice daily (one in the AM and one in the PM) over multiple days. The auto-regressive correlation structure makes sense since observations that are closer in time would be more-highly correlated with one another.

However, I want to use a non-default correlation structure to address the specific design of this dataset. I am hoping to implement a "direct product" structure that accounts for doubly-repeated measures, since we have reason to believe that measurements taken in the morning would be correlated with one another, as would measurements taken on nearby days.

Here I have to refer to SAS's "REPEATED" statement since one of the default choices for covariance structures is the direct product called UN using AR(1), the Kronecker product of an unstructured matrix and an auto-regressive matrix.

I am aware that non-default correlation structures can be defined by users in R. In the documentation for nlme, the authors mention:

  "Users may define their own corStruct classes by specifying a constructor function and, at a minimum, methods for the functions corMatrix and coef. For examples of these functions, see the methods for classes corSymm and corAR1."

Unfortunately the examples of code for corSymm and corAR1 are beyond my understanding of object-oriented programming in R. I am wondering if anyone has successfully defined their own correlation structure. More generally, I wonder if anyone has encountered the same issue in linear mixed models where there may be multiple sources of repeated measures, and modeling correlation via AR1 alone is not addressing the specific design of the dataset.

I expect that my user-defined correlation structure would be called as in the following sample code:

  > lme(
     y ~ 1,
     data        = dat,
     random      = ~ 1 | id,
     correlation = corUNxAR1(form = ~ am_pm + day |id),
     control     = lmeControl(opt="optim")
     )

Lastly, for context this is for a graduate school course on longitudinal data analysis. Until this point, we have relied heavily on SAS since we have not found a way to implement Kronecker product covariance structures in R. We are hoping to find a solution here and demonstrate that such complex cases can be handled in R.

Some relevant details below:

* R version 4.2.1
* nlme version 3.1-157

	[[alternative HTML version deleted]]