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

[Anders Christian Jensen]

Walmes,

Thank you for your fast reply!

If I had only one time series per person (no r index) then I would use the
model you suggest, but I dont think this is the solution to my problem:

If I understand correctly, your model would imply that for all r and r',

e_{prt} and e_{pr't'}

would have the same correlation as

e_{prt} and e_{prt'}.

That is, for a given person you would have a correlation between
measurements from different time series. But this is not what Im trying to
model. I want measurements from different replications to be correlated
only through the random person effect.

Thanks!

Anders

> 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 at math.ku.dk>
> Para: r-sig-mixed-models at 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 at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
>
>       ____________________________________________________________________________________
> Veja quais são os assuntos do momento no Yahoo! +Buscados
> http://br.maisbuscados.yahoo.com