41

I have within-subject physiological data from participants (part), who have all looked at stimuli (reading newspapers) on three rounds (round), which each have five papers (paper), and within each there are variable number of visits (visit) in the newspaper. I have two fixed factors (CONDhier and CONDabund) plus interaction to predict the physiological state (e.g., EDA), which is usually autoregressive. I try to take into account individual differences in physiology with random effects (let's settle for intercept only for now), and perhaps fatigue over rounds with another random effect.

Thus, my model that I would like to run in R would be, in SPSS:

MIXED EDA BY CONDhier CONDabund  
 /FIXED=CONDhier CONDabund CONDhier*CONDabund | SSTYPE(3)  
 /RANDOM=INTERCEPT | SUBJECT(part) COVTYPE(VC)  
 /RANDOM=INTERCEPT | SUBJECT(part*round) COVTYPE(VC)  
 /PRINT=SOLUTION  
 /METHOD=REML  
 /REPEATED=visit | SUBJECT(part*round*paper) COVTYPE(AR1).

Now, I have understood that while lme does not do crossed terms well, lmer (that handles crossed terms with no problem) cannot use different covariance structures. I can run simple lme models such as

    lme(EDA ~ factor(CONDhier) * factor(CONDabund), random= ~1
   |part, na.action=na.exclude, data=phys2)

but a more complicated model is beyond me. I have read that crossed terms in lme can be done with random definitions like 

    random=pdBlocked(list(pdCompSymm(~part), pdCompSymm(~round-1), pdCompSymm(~paper-1), 
pdCompSymm(~visit-1)))

but that seems to block the AR1 structure, and the second random intercept for part*round, from me. And I'm not so sure it's the same as my SPSS syntax anyway.

So, any advice? Although there are lots of different writings on lme and lmer, I couldn't find one that would have both crossed terms and AR1.

(Also, the syntax on lme seems quite obscure: from several different sources I have understood that | groups what is on the left under what's on the right, that / makes nested terms, that ~1 is random intercept, ~x is random slope, and ~1+x is both, but there seems to be at least : and -1 definitions that I couldn't find anywhere. Is there a tutorial that would explain all the different definitions?)

Sachu
  • 7,555
  • 7
  • 55
  • 94
RandomMonitor
  • 439
  • 1
  • 8
  • 16
  • 2
    not complete, but for a bit more on R mixed model syntax, see http://glmm.wikidot.com/faq#modelspec – Ben Bolker Nov 21 '13 at 15:27
  • Thanks! (+filler as comments must be at least 15 chars in length...) – RandomMonitor Nov 21 '13 at 16:08
  • 5
    You're right that `lme4` lacks "R-side" (autocorrelation) structures (and probably will for a while, we're swamped). I'm not sure (a reproducible example would be nice), but you *might* want something like `random=pdBlocked(list(pdCompSymm(~part-1), pdCompSymm(~round-1), pdCompSymm(~paper:round), pdCompSymm(~visit:paper:round)))` ... and I don't see quite what you mean by "block the AR structure". You might want `correlation=corAR1()` (although you may be right in saying that doesn't work). AD Model Builder/JAGS/BUGS/Stan (build-your-own) are the only open-source tools I know of for this – Ben Bolker Nov 21 '13 at 19:55
  • Are you sure you mean `pdCompSymm(~paper:round)` and `pdCompSymm(~visit:paper:round)`, and not the other way around - round:paper and round:paper:visit? I mean, all visits occur within papers, and all papers occur within rounds, and looking at the specifications in your link it seems to me that the grouping variable (on the right side of | ) is on left and the smaller thing within that variable on its right. (Or am I reading something wrong? English is not my first language.) – RandomMonitor Nov 22 '13 at 11:56
  • (stupid character limitations) And am I correct that : signifies calculating the intercept or slope only for the rightmost variable while / means calculating both the group variable and the variable within it? Further, by "block the AR structure" I meant that using pdCompSymm would define using compound symmetry, and that I couldn't use AR1 in addition to that. Or does a separate correlation specification override the compound symmetry structure? – RandomMonitor Nov 22 '13 at 11:57

2 Answers2

1

Consider the R package MCMCglmm which allows for complex mixed effects models.

https://cran.r-project.org/web/packages/MCMCglmm/vignettes/CourseNotes.pdf

Although it can be challenging to implement, it may solve the problems you've been having. It allows the fixed and random effects formulas to be given separately, eg.

fixed <- formula(EDA ~ CONDhier * CONDabund)
rand <- formula( ~(us(1+ CONDhier):part + us(1+ CONDhier):round + us(1+ CONDhier):paper + us(1+ CONDhier):visit))

The covariance structure between the random effects are given as coefficients which can be examined using summary() on the MCMCglmm object once the model has been run.

Benjamin Christoffersen
  • 4,703
  • 15
  • 37
FlamingGoose
  • 153
  • 1
  • 9
-2

If you have a cross-covariance matrix use canonical correlation analysis (CCA). There is a documented R package for CCA.

noumenal
  • 1,077
  • 2
  • 16
  • 36