Specifying correlation structure
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 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%