[R-sig-ME] lme or lmer?

[Antoine Tremblay]

Giorgio, another thing regarding the Interval variable that came to my 
attention after reading Ben Bolker's reply to your question.

Were you planning on adding Interval to your model specification? 
Because it *shouldn't* be in the model. From what I can see, you'd be 
performing an analysis on the 200-300 time window, then another one on 
the 400-600 time window, and so forth, and in each analysis you'd 
include all electrodes (which, as I mentioned in my reply, is probably 
best re-coded as ROI with 9 levels).

I'm assuming you chose these time windows because you were expecting to 
find amplitude differences between your two conditions with respect to 
ERPs found in these time windows (e.g., maybe a late N400 in the 400-600 
window?).

So, for *each* time window the model would look something like this 
(after recoding Electrode to ROI and your datya frame is named eeg):

m1 <- lmer(Ampl ~ Condition * ROI + (1|Subject) + (1|Item) + 
(ROI|Subject), data = eeg)

Would check model assumptions using e.g., function mcp.fnc from package 
LMERConvenienceFunctions.

mcp.fnc(m1, trim = 2.5)

Typically the qq plot top right isn't good and there are a lot of data 
points with undue influence (bottom right panel). But after removing 
outliers +/- 2.5 std. dev. below and above the m1 residuals mean (could 
use function romr.fnc from package LMERConvenienceFunctions) the plots 
look much, much better.

eeg<- romr.fnc(model = m1, data = eeg, trim = 2.5)$data

Typically 1.5% to 2.5% of the data is removed this way (the function 
tells you the percentage). Then refit the model:

m1 <- update(m1)

Then look at the model criticism plot:

mcp.fnc(m1, trim = 2.5)


In addition, you could start with a simple model with only (1|Subject) 
then fit a more complex model with (1|Subject) + (1|Item) then a more 
complex one with (1|Subject) + (1|Item) + (ROI|Subject) ... and test 
whether the inclusion of the random effects is warranted using a log 
likelihood ratio test (anova(m1,m2) ...). I'd bet that (1|Subject) + 
(1|Item) + (ROI|Subject) will be warranted, maybe also (Condition|Subject).

Sincerely,

Antoine Tremblay
NeuroCognitive Imaging Laboratory
Dalhousie university
Halifax, Canada

> ----------------------------------------------------------------------
>
> Message: 1
> Date: Sat, 21 Jan 2012 15:54:00 +0000 (UTC)
> From: Ben Bolker<bbolker at gmail.com>
> To: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] lme or lmer?
> Message-ID:<loom.20120121T041825-883 at post.gmane.org>
> Content-Type: text/plain; charset=us-ascii
>
> Giorgio Arcara<giorgio.arcara at ...>  writes:
>
>> I would like to use mixed models in R to analyze EEG data, but I don't
>> know if it is more correct to use lme or lmer.
>> My data have the following structure
>>
>> Subject	Electrode	 Interval	Trial 	Condition	Ampl
>> 1		Fp1		200-300		1	A		3.5
>> 1		Fp1		200-300		2	B		4.2
>> 1		Fp2		400-600		1	A		6.5
>> 1		Fp2		400-600		2	B		3.3
>> 2		Fp1		200-300		1	A		2.1
>> 2		Fp1		200-300		2	B		-5.4
>> 2		Fp2		400-600		1	A		-5.6
>> 2		Fp2		400-600		2	B		-3.2
>
>   [snip]
>
>    In general if you *can* use either lme or lmer (i.e. the random effects are
> not nested, there is no heteroscedasticity or correlation structure at the
> level of residuals ...) then either is OK.  I generally recommend lme in
> these cases (with apologies to Doug Bates) because it is better documented,
> although it's slower.
>
>> For sake of simplicity, here I include only 2 Trials but in the real
>> dataset they are many more.
>> In this hypothetical dataset Ampl is the depentent variable. Electrode
>> and Interval are two predictors. I expect that levels of Electrode
>> will be highly correlated as well the levels of Interval.
>
>    ?? Meaning that most instances of FP1 have Interval 200-300, most
> of Fp2 have 400-500, and there are few 'crossover' instances where you
> have FP1 with 400-500 and FP2 with 200-300?  This will in general make
> it hard to estimate to distinguish the two effects ... but presumably
> you know that.
>
>> My goal is to study if Condition influence Ampl and if interact with
>> Electrode variable and Interval Variable.
>>
>> I would fit a model on these data with lmer with the following structure
>>
>> mod=lmer(Ampl~Electrode*Interval*Condition+(1+Electrode|Subject)
>> +(1+Interval|Subject))
>> If I'm correct the corresponding lme model would be
>>
>> mod=lme(Ampl~Electrode*Interval*Condition, random=list(~1+condition|
>> Subject, ~1+Interval|Subject))
>
>    Why not (Electrode+Interval|Subject)?  That (a) avoids a duplicated intercept
> term and (b) estimates correlations between the intercept, electrode, and
> Interval RE.
>
>> Any suggestion for covariance matrix specification in lme?
>
>    Do you need one other than the default -- i.e. do you expect
> any particular correlation structure?
>
>
>
> ------------------------------
>
> Message: 2
> Date: Sat, 21 Jan 2012 12:55:50 -0400
> From: Antoine Tremblay<trea26 at gmail.com>
> To: r-sig-mixed-models at r-project.org
> Cc: r-sig-mixed-models-request at r-project.org
> Subject: Re: [R-sig-ME] lme or lmer?
> Message-ID:<4F1AEE16.3090801 at gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> Hey Giorgio!
>
> I have one published paper in JoCN where we analyzed ERP data using
> LMER, another one is on the way, and a methods paper is almost finished
> and submitted to NeuroImage regarding the addition of by-item random
> intercepts and slopes to the model. I'm sending you a copy of the JoCN
> paper (proofs for now, it'll be officially out any time soon) and a copy
> of the almost finished NeuroImage paper. If others are interested,
> please let me know and I can send them to you (trea26 at gmail dot com).
>
> So LMER-wise, adding (1+Electtrode|Subject), which you could simply
> write (Electrode|Subject) is correct, but fitting a model with such a
> random effect will take forever.
>
> In NewmanTremblay2011 paper we collapsed Electrode into 9 ROIs: left
> anterior, midline anterior, right anterior, left central, midline
> central, right central, left posterior, midline posterior, and right
> posterior and use that new variable (ROI) instead of Electrode. Note
> that we are not averaging at all! The main reason for this is
> computation time: With 9 levels of ROI in the model as (ROI|Subject) it
> takes A LOT of time to fit, tried once with (Electrode|Subject) and
> would run for sooooo long, actually killed it after like two or three
> days. Even the model with ROI takes what seems to be forever.
>
> Then, once fitted, you'll see with print(model) that in the random
> effects portion of the summary there's a table of correlations between
> levels of ROI, something like this from one of our models (incomplete):
>
> Linear mixed model fit by REML
> Formula: Ampl ~ rProficiency * Condition * ROI * Group + (Condition |
>    Subject) + (1 | Subject) + (ROI | Subject)
>      Data: dat
>        AIC     BIC  logLik deviance REMLdev
>    1011929 1013147 -505843  1011340 1011685
> Random effects:
>    Groups   Name          Variance   Std.Dev. Corr
>    Subject  (Intercept)   2.5852e-01 0.508447
>             ConditionGood 1.4045e+00 1.185133 -0.681
>    Subject  (Intercept)   1.8638e-04 0.013652
>    Subject  (Intercept)   2.8105e+00 1.676449
>             ROILcent      9.4985e-01 0.974604 -0.423
>             ROILpost      3.2123e+00 1.792281 -0.644  0.833
>             ROIMant       3.3728e-01 0.580755 -0.240  0.150  0.258
>             ROIMcent      1.8586e+00 1.363301 -0.249  0.687  ...
>             ROIMpost      2.8172e+00 1.678445 -0.577  0.746  ...
>             ROIRant       7.5682e-01 0.869957 -0.421 -0.171  ...
>             ROIRcent      1.4031e+00 1.184519 -0.840  0.531  ...
>             ROIRpost      3.0983e+00 1.760201 -0.797  0.578  ...
>    Residual               3.3207e+01 5.762530
> Number of obs: 159374, groups: Subject, 44
>
> I'm attaching a complete report (using Sweave) for your reference.
>
> There's one other paper where they used lme with a corSpher function in
> Davidson2007 (attached here, see page 90). The problem i see with using
> lme is that you can't really add crossed by-item random effect, which
> you should as demonstrated in Tremblay2012. Basically, adding by-item
> random effects substantially decreases the amount of (partial)
> autocorrelation in the model residuals (i.e., better approximation of
> the assumption of independence of errors. Note that the data used in
> that paper is available on CRAN (data package LCFdata) and I can send
> you the .Rnw file that contains the R code used for data manipulation,
> analysis and plotting in Tremblay2012, so it's fully replicable. Note
> that for simplicity, analyses in Tremblay2012 are on a single electrode,
> and that I'm not comparing (ROI|Subject) in LMER to corSpher in LME.
> Would like to eventually look a variograms to see how well the two take
> care of spatial correlation and also compare them with a model that
> doesn't account for this correlation.
>
> Note that, as I demonstrate in Tremblay2012, by-item random effect
> should be added if warranted (by Log-likelihood Ratio Test or other; I
> suspect it will always be).
>
>
> Cheers,
>
> Antoine
>
>
> End of R-sig-mixed-models Digest, Vol 61, Issue 31
> **************************************************
>