[R-sig-ME] Specifying correlation structure

[Gavin Simpson]

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
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%