[R-sig-ME] Specifying correlation structure
Gang Chen
gangchen at mail.nih.gov
Tue May 15 18:52:26 CEST 2012
The mixing of cond vs. condition was really a silly mistake of mine!
Thanks a lot for spotting that...
Yes, corARMA(c(0.02, 0.03), form=~time|cond, p=1,q=1) seems to work by
imposing the SAME ARMA(1, 1) structure for all the conditions.
May I ask one more question: how can I specify an ARMA model for each
condition separately while also accounting for the fact the three
conditions are correlated? That is why I originally thought about
something like "form=~time:condition". Any thought?
Thanks again,
Gang
On Tue, May 15, 2012 at 11:39 AM, Gavin Simpson <gavin.simpson at ucl.ac.uk> wrote:
> Your variables are looked up in *Dat*, it would be a mistake to try rely
> upon variable being found in various environments as you are doing here.
> That would also force you to consider what you are doing and how you are
> setting up your data in R. Put everything into `Dat`, including the
> `regx` variables.
>
> I didn't realise this at the time but you are doing this:
>
> Dat <- data.frame(time=rep(trials, 3), cond=rep(condition, each=tp),
> res=rnorm(3*tp))
>
> notice you are now calling the condition variable `cond`! So it is no
> wonder that `condition` doesn't have the required length.
>
> If you use the `Dat` as defined, then
>
> (fm <- gls(res ~ 1+reg1+reg2+reg3,
> correlation=corARMA(c(0.02, 0.03), form=~time|cond, p=1,q=1),
> data=Dat))
>
> works fine
>
>> (fm <- gls(res ~ 1+reg1+reg2+reg3, correlation=corARMA(c(0.02, 0.03),
> form=~time|cond, p=1,q=1), data=Dat))
> Generalized least squares fit by REML
> Model: res ~ 1 + reg1 + reg2 + reg3
> Data: Dat
> Log-restricted-likelihood: -425.2665
>
> Coefficients:
> (Intercept) reg1 reg2 reg3
> 0.01305633 -0.03300839 -0.08407831 0.03701266
>
> Correlation Structure: ARMA(1,1)
> Formula: ~time | cond
> Parameter estimate(s):
> Phi1 Theta1
> -0.1197791 0.1345616
> Degrees of freedom: 300 total; 296 residual
> Residual standard error: 0.9848909
>
> Although I would probably clean up my working environment and call `Dat
> $cond` `Dat$condition` but it is up to you. Just make sure you are
> referring to things **in** `Dat` with the correct names.
>
> HTH
>
> G
>
> On Tue, 2012-05-15 at 09:07 -0400, Gang Chen wrote:
>> Many thanks for the suggestion. It seems that gls does not like that either:
>>
>> > (fm <- gls(res ~ 1+reg1+reg2+reg3, correlation=corARMA(c(0.02, 0.03), form=~time|condition, p=1,q=1), data=Dat))
>>
>> Error in model.frame.default(formula = ~time + condition + res + reg1 + :
>> variable lengths differ (found for 'condition')
>>
>>
>> On Tue, May 15, 2012 at 8:11 AM, Gavin Simpson <gavin.simpson at ucl.ac.uk> wrote:
>> > On Tue, 2012-05-15 at 06:27 -0400, Gang Chen wrote:
>> >> I have some time series data (with 100 time points) collected from one
>> >> subject. At each time point three measures were taken under three
>> >> conditions, and there is an explanatory variable for each condition.
>> >
>> > The temporal correlation is, in this sense, *nested* within the three
>> > different conditions. You were close with time:condition, but what I
>> > think you need is
>> >
>> > corARMA(c(0.02, 0.03), form = ~ time | condition)
>> >
>> > HTH
>> >
>> > G
>> >
>> >> Below is just some fake data to demonstrate the data structure:
>> >>
>> >> # 100 time points per condition
>> >> tp <- 100 # time points
>> >> trials <- seq(1, tp, 1) # time counter
>> >> condition <- c('a', 'b', 'c') # three conditions
>> >>
>> >> # fake data
>> >> set.seed(5)
>> >> Dat <- data.frame(time=rep(trials, 3), cond=rep(condition, each=tp),
>> >> res=rnorm(3*tp))
>> >>
>> >> reg1 <- c(rnorm(tp), rep(0, 2*tp)) # explanatory variable
>> >> for condition 1
>> >> reg2 <- c(rep(0, tp), rnorm(tp), rep(0, tp)) # explanatory variable
>> >> for condition 2
>> >> reg3 <- c(rep(0, 2*tp), rnorm(tp)) # explanatory variable
>> >> for condition 3
>> >>
>> >> First I thought that I'd start with gls in nlme package with an ARMA
>> >> model for the correlation structure:
>> >>
>> >> (fm <- gls(res ~ 1+reg1+reg2+reg3, correlation=corARMA(c(0.02, 0.03),
>> >> form=~time, p=1,q=1), data=Dat))
>> >>
>> >> The above model does not work because of the following error:
>> >>
>> >> Error in Initialize.corARMA(X[[1L]], ...) :
>> >> Covariate must have unique values within groups for corARMA objects
>> >>
>> >> Moreover, I would like to account for the fact that the ARMA structure
>> >> should be similar or the same across the three conditions. My
>> >> questions are:
>> >>
>> >> 1) how to impose a same ARMA structure across the three conditions?
>> >> 2) alternatively I'd like to have something like:
>> >>
>> >> (fm <- gls(res ~ 1+reg1+reg2+reg3, correlation=corARMA(c(0.02, 0.03),
>> >> form=~time:condition, p=1,q=1), data=Dat))
>> >>
>> >> but this would not work either, and seems to crash R!
>> >>
>> >> 3) maybe construct a multivariate gls model, but how to do that? what package?
>> >>
>> >> Thanks,
>> >> Gang
>> >>
>> >> _______________________________________________
>> >> R-sig-mixed-models at r-project.org mailing list
>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> >>
>> >
>> > --
>> > %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>> > Dr. Gavin Simpson [t] +44 (0)20 7679 0522
>> > ECRC, UCL Geography, [f] +44 (0)20 7679 0565
>> > Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
>> > Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/
>> > UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
>> > %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>> >
>> >
>>
>
> --
> %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
> Dr. Gavin Simpson [t] +44 (0)20 7679 0522
> ECRC, UCL Geography, [f] +44 (0)20 7679 0565
> Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
> Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/
> UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
> %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>
>
More information about the R-sig-mixed-models
mailing list