Anders,

I think that the correct sintax or model specification is different. I think you should use form=~1|Person instead of form=~1|Person_replicate at correlation argument, this because repeated measures on the same person along time induces correlation. More about it in Pinheiro & Bates (2000). The code is the following:

> 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))
> 
> lme(response~Time, random=list(~Time|Person),
+     correlation=corGaus(form=~1|Person), data=Data)
Linear mixed-effects model fit by REML
  Data: Data 
  Log-restricted-likelihood: -41.12866
  Fixed: response ~ Time 
(Intercept)        Time 
  5.3595834  -0.2110061 

Random effects:
 Formula: ~Time | Person
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev       Corr  
(Intercept) 2.270560e-05 (Intr)
Time        1.236894e-08 0     
Residual    7.507016e-01       

Correlation Structure: Gaussian spatial correlation
 Formula: ~1 | Person 
 Parameter estimate(s):
    range 
0.8230017 
Number of Observations: 36
Number of Groups: 3 
> 

Walmes Zeviani. Brasil.




________________________________
De: Anders Christian Jensen <m03acj@math.ku.dk>
Para: r-sig-mixed-models@r-project.org
Enviadas: Terça-feira, 16 de Fevereiro de 2010 13:57:15
Assunto: [R-sig-ME] Random effects, autocorrelation and nesting in lme

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

_______________________________________________
R-sig-mixed-models@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



      ____________________________________________________________________________________
[[elided Yahoo spam]]

	[[alternative HTML version deleted]]


