[R-sig-ME] Random effects, autocorrelation and nesting in lme

[Anders Christian Jensen]

Dear R - users

I have a problem with lme when trying to specify both "random" and
"correlation" arguments to my model.

The problem seems to be related to the nesting structure of my data:

For each person (p) I have multiple replicates (r) of a time series of
observations. Observations within a time series is indexed by t. That is,
my outcome (y) is uniquely determined by indices y_{prt}.

I would like to specify a model that includes
1)  a random slope and intercept for each person, and
2)  some sort of correlation on the time series, within each
(person,replicate); possibly an AR(1) process.

In other words, I would like to specify the following model (disregarding
the fixed effects part of the model):

y_{prt} = a_p + b_p * t + e_{prt},

where a_p and b_p are the (Gaussian) random slope and intercept on person
level, and e_{prt} should be such that

cor(e_{prt},e_{p'r't'}) = f(|t-t'|), for some function f

when p=p' and r=r', and 0 otherwise.

When I try to fit the model with lme I get an error message:

"Incompatible formulas for groups in "random" and "correlation""

A very small toy data set (without any autocorrelation) illustrates my
problem:

> Data <- data.frame(Person           = factor(rep(1:3,each=12)),
+                    Person_replicate = factor(rep(1:12,each=3)),
+                    Time             = rep(1:4,9),
+                    response         = round(5+rnorm(36),2))
> str(Data)
'data.frame':   36 obs. of  4 variables:
 $ Person          : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
 $ Person_replicate: Factor w/ 12 levels "1","2","3","4",..: 1 1 1 2 2 2 3
3 3 4 ...
 $ Time            : int  1 2 3 4 1 2 3 4 1 2 ...
 $ response        : num  5.21 4.95 3.32 4.86 6.18 5.68 5.14 3.81 6.17
5.08 ...
> head(Data)
  Person Person_replicate Time response
1      1                1    1     5.21
2      1                1    2     4.95
3      1                1    3     3.32
4      1                2    4     4.86
5      1                2    1     6.18
6      1                2    2     5.68

I then try to fit the model, and get an error:
> library(splines)
> library(nlme)
> lme(response ~ Time,
+     random = list(~Time|Person),
+     correlation=corGaus(form=~1|Person_replicate),
+     data = Data)
Error in lme.formula(response ~ Time, random = list(~Time | Person),
correlation = corGaus(form = ~1 |  :
  Incompatible formulas for groups in "random" and "correlation"

Basically I would like the "random"-argument for lme to work on person
level, but the "correlation"-argument should work on (person,replicate)
level. I found old posts on the mailing list with similar problems but
none of them include an answer.
Any help would be much appreciated! Thanks