Hi Tobias,

Thank your reply.

On 06/10/2013 10:18 AM, Tobias Heed wrote:
> Gabriel,
>
> did you see the "add-on" to the paper you are citing:
> Barr, D. J. (2013). Random effects structure for testing interactions 
> in linear mixed-effects models. /Frontiers in Psychology/, /4/. 
> doi:10.3389/fpsyg.2013.00328
> It deals with exactly your case.

Yes, it does :-). Thank you.

>
> The random effect structure for your experiment should look like this:
>> (1 + task * vf | su0)

In my dataset, I had no replication within cell (the reason is that I 
wanted to compare lmer with old-fashion repeated measure ANOVA which is 
done on the mean value for each subject and condition). In this case, 
including the interaction term led to convergence problems. As noted in 
Barr's add-on, interaction term should be included for design with 
replications.

When I use the dataset with the replications (5 per cell), I am able to 
get the desired results (with bobyqa optimizer):

fit<-lmer(log(e.area)~group*task*vf + (1+task*vf|su0), CPvf,
     control = lmerControl(optimizer ="bobyqa" ))
VarCorr(fit)
  Groups   Name        Std.Dev. Corr
  su0      (Intercept) 0.461895
           task1       0.120315 0.170
           task2       0.156975  0.421 -0.626
           vf1         0.109569 -0.104 -0.301 -0.060
           task1:vf1   0.024895 -0.715 -0.474  0.098 -0.394
           task2:vf1   0.067276  0.197  0.339  0.101  0.452 -0.594
  Residual 0.373385

I still have some questions about lmer specifications if one wanted to 
remove some correlations.

Note 2 things:
> 1, all within-subjects effects relating to your cluster variable (su0) 
> are listed just like they are as fixed effects - i.e., note the * 
> between task and vf, which includes the main effects and the 
> interaction as random slopes in your random effects.
> 2, all effects relating to your cluster variable are included in 1 
> bracket. Like this, everything is allowed to be correlated. This would 
> be the maximum RE structure for your experiment.
>
> The way you coded it, you did not include the interaction term.
> I think you got several intercepts because you had several random 
> terms for su0. If you write (vf | su0), then this creates an intercept 
> and a slope for vf. Then, by writing another (task | su0), you create 
> another intercept. If you code (task * vf | su0), then you create 1 
> intercept and slopes for all effects.

Yes.

> If you wanted the different random effects to not be correlated, you 
> would write: + (1 | su0) + (0 + vf | su0) + (0 + task | su0) + (0 + 
> vf:task | su0)
> (I am not sure whether the way you coded, with -1, does the same 
> thing). This way, you should also get only 1 intercept.
>

Writing 0+task gives the same thing as task-1.

fit<-lmer(log(e.area)~group*task*vf + (1 | su0) + (0 + vf | su0) + (0 + 
task | su0) + (0 + vf:task | su0), CPvf,
     control = lmerControl(optimizer ="bobyqa",optCtrl=list(maxfun=100000)))
  VarCorr(fit)
#  Groups   Name              Std.Dev. Corr
# su0      (Intercept) 1.5432e-07
# su0.1    vfVF 4.1476e-08
#          vfnoVF            3.3620e-07 -1.000
# su0.2    taskstand 1.4901e-01
#          taskltstat        3.5674e-01 1.000
#          taskltdyn         1.5498e-01 1.000 1.000
# su0.3    vfVF:taskstand 4.4405e-01
#          vfnoVF:taskstand  5.2225e-01 0.923
#          vfVF:taskltstat   4.4352e-01 0.950 0.815
#          vfnoVF:taskltstat 4.3955e-01 0.826  0.941 0.762
#          vfVF:taskltdyn    3.8040e-01 0.889  0.740 0.840 0.709
#          vfnoVF:taskltdyn  3.5549e-01 0.944  0.981 0.849 0.954 0.845
# Residual                   3.7338e-01

As you can see, factor + 0 (or factor - 1) yields a design matrix for 
random effects which appears  quite
redundant. In fact,  Variance for (1|Su0) is now very small and 
correlation coefficients are
very high (this is why convergence is not easy and required increase 
maximum number of iteration). Also,
  this design matrix does not take into account the constrasts 
associated with each factor.

This is something I find surprising in lmer coding of random effects 
(see my original
post below). I would have expected that the contrasts would be used in 
the random
design matrix (like for the fixed effect design matrix).

In other words, what I hoping to get with the above syntax was 
VarianceCovariance matrix like that

  Groups   Name        Std.Dev. Corr

  su0      (Intercept) 0.461895
su0.1    task1       0.120315
           task2       0.156975  0.421
su0.2    vf1 0.109569
su0.3    task1:vf1 0.024895
           task2:vf1   0.067276  0.197
  Residual             0.373385

where task1 and task2 are the two contrasts for the task factor, vf1 the 
contrast for the vf factor, etc.

A possible variant would be to also include correlation with the intercept.

  Groups   Name        Std.Dev. Corr

  su0      (Intercept) 0.461895
su0.1    task1       0.120315  cor(task1,intercept)
           task2       0.156975 cor(task2,intercept) cor(task1,task2)
su0.2    vf1         0.109569 cor(vf1,intercept)
su0.3    task1:vf1 0.024895 cor(task1:vf1,intercept)
           task2:vf1   0.067276 cor(task2:vf1,intercept) 
cor(task1:vf1,task2:vf1)
  Residual             0.373385

but this would definitively be incompatible with the way random effects 
are specified with lmer.

> Finally, about your last question concerning task: In your output, I 
> see two terms for the random effect of task, taskltstat and taskltdyn. 
> If you have 3 levels, you code them with 2 contrasts, so this is 
> exactly what you should expect. The intercept codes for the baseline, 
> and each random effect relates to one of your two contrasts. 
> Interpretation depends on how you code your "task" contrasts. You are 
> probably familiar with this, but if not, then you can google for "sum 
> coding", "deviation coding" and "treatment coding". There is also a 
> bit of this in the Barr article you cited, if I remember correctly 
> (maybe in the supplementary info).
>

For the random design matrix, we also get the same number of column as 
there are levels in the factor.

 > fit<-lmer(log(e.area)~group*task*vf + (1 | su0) + (vf | su0) + (task 
| su0) + (vf:task | su0), CPvf,
+     control = lmerControl(optimizer 
="bobyqa",optCtrl=list(maxfun=100000)))
 > VarCorr(fit)
#  Groups   Name              Std.Dev. Corr
#  su0      (Intercept) 0.0000e+00
#  su0.1    (Intercept) 1.6116e-05
#          vf1               5.4197e-06 1.000
#  su0.2    (Intercept) 2.1474e-01
#           task1             6.9445e-02 -1.000
#           task2             1.3308e-01  1.000 -1.000
#  su0.3    (Intercept) 3.8718e-01
#           vfVF:taskstand    1.4249e-01 0.254
#           vfnoVF:taskstand  1.9936e-01  0.545 0.320
#           vfVF:taskltstat   1.9431e-01  0.102  0.653 -0.212
#           vfnoVF:taskltstat 1.4865e-01  0.239 -0.650 0.439 -0.570
#           vfVF:taskltdyn    1.9090e-01 -0.274  0.280 -0.620  0.191 -0.779
#           vfnoVF:taskltdyn  6.9589e-02 -0.505 -0.019 0.347 -0.542  
0.163 -0.049
#  Residual                   3.7338e-01

As you can see, contrasts are being used for the main effects. For the 
interaction, I get 7 variance components!
As in the previous example, this appears very redundant.

I am really not sure how to interpret that and I would feel much better 
if the design matrix was build in a way
that is similar to the fixed effects.


I would appreciate any feedback on why things have been setup in this 
manner.

> Should you have problems getting the full model to converge, you can 
> try the bobyqa optimizer (passed to lmer/gmler using g(l)merControl) 
> -- in my experience, this lets models converge better than the 
> standard optimizer.
>

Yes, I agree.  Convergence appears easier with bobyqa.

Thanks a lot.

Gabriel

> Hope this helps.
> Best,
> Tobias
>
> -- 
> --------------------------------------------------------------------------------------------------------------
> *Tobias Heed, PhD*
> Biological Psychology and Neuropsychology  |  University of Hamburg
> Von-Melle-Park 11, Room 208  |  D-20146 Hamburg, Germany
> Phone: (49) 40 - 42838 5831  |  Fax:   (49) 40 - 42838 6591
> tobias.heed@uni-hamburg.de <mailto:tobias.heed@uni-hamburg.de>  | 
> Website <http://www.epb.uni-hamburg.de/de/personen/heed>  | Google 
> Scholar <http://tinyurl.com/9gh83ea>  | ResearcherID 
> <http://www.researcherid.com/rid/C-5330-2008>
> --------------------------------------------------------------------------------------------------------------
>
> On 6 Oct 2013, at 02:53, Gabriel Baud-Bovy <baud-bovy.gabriel@hsr.it 
> <mailto:baud-bovy.gabriel@hsr.it>> wrote:
>
>> Hi r-sig-me,
>>
>> I am analyzing the results of an experiment with one between-subject 
>> factor (group, 3 levels: yo, yf, ef)
>> and two within-subject factors (vf, 2 levels: vf and novf, and task, 
>> 3 levels: stand, ltstat, ltmov).
>> Each group includes uniquely identified 15 subjects. All factors are 
>> coded with sum contrasts (contr.sum).
>>
>> Barr et al. 2013 advise to keep the random effects structure maximal. 
>> Moreover, they write
>> that "In general, within-unit treatments require both the by-unit 
>> intercepts and slopes in the
>> random effects specification, whereas between-unit treatments require 
>> only the by-unit
>> random intercepts."
>>
>> Barr et al. (2013) Random effects structure for confirmatory 
>> hypothesis testing:
>> Keep it maximal. Journal of Memory and Language 68, 255–278
>>
>> I have two questions:
>>
>> First, how would you specify the model if you were trying to follow 
>> Barr et al.
>> (2013) advice to have a maximal random effect covariance structure 
>> (granted that it still has
>> to be identifiable) for such a design.
>>
>> Second, I could not find a discussion discussing design matrix for 
>> random effects corresponding
>> to two or more within-subject factors. It seems that it would make 
>> sense in cases like that
>> to define one intercept for each subject plus one random effect for 
>> each treatment relative
>> to the baseline but the syntax to achieve it is ackward. It seems 
>> that the default leads
>> to design matrices with two many random effects. See the examples below.
>>
>> Following Barr, my initial attempt was:
>>
>> > fit<-lmer(log(e.area)~group*task*vf + (1|su0) + (task|su0) + 
>> (vf|su0), CP0vf)
>> In checkZrank(reTrms$Zt, n = n, control, nonSmall = 1e+06) :
>> number of observations <= rank(Z); variance-covariance matrix will be 
>> unidentifiable.
>> > VarCorr(fit)
>> Groups Name Std.Dev. Corr
>> su0 (Intercept) 0.21461
>> su0.1 (Intercept) 0.36124
>> taskltstat-stand 0.16478 0.501
>> taskltdyn-stand 0.12364 -0.949 -0.676
>> su0.2 (Intercept) 0.19043
>> vfnovf-vf 0.10934 -0.224
>> Residual 0.17567
>>
>> This model yields a warning message and includes more random effects 
>> that I wanted (three "intercepts" !)
>>
>> My first attempt to remove intercept to simplify the covariance 
>> structure was not successful
>>
>> > fit<-lmer(log(e.area)~group*task*vf + (1|su0) + (task-1|su0) + 
>> (vf-1|su0),CP0vf)
>> which yielded the warning message In checkZrank(reTrms$Zt, n = n, 
>> control, nonSmall = 1e+06) :
>> number of observations <= rank(Z); variance-covariance matrix will be 
>> unidentifiable.
>> > VarCorr(fit)
>> Groups Name Std.Dev. Corr
>> su0 (Intercept) 0.30722
>> su0.1 taskstand 0.38166
>> taskltstat 0.43823 0.804
>> taskltdyn 0.22527 0.949 0.888
>> su0.2 vfVF 0.11602
>> vfnoVF 0.10147 -1.000
>> Residual 0.17508
>>
>> This yields the same error message and looking at the structure of 
>> the covariance matrix
>> I see that the number of random effects has not changed (only that 
>> now the coding is
>> made with a dummy matrix instead of the contrasts).
>>
>> The following model however does not give any warning message
>>
>> > fit<-lmer(log(e.area)~group*task*vf + (task|su0) + (vf|su0),CP0vf)
>> > VarCorr(fit)
>> Groups Name Std.Dev. Corr
>> su0 (Intercept) 0.39050
>> taskltstat-stand 0.16468 0.492
>> taskltdyn-stand 0.12581 -0.891 -0.659
>> su0.1 (Intercept) 0.27859
>> vfnovf-vf 0.11002 -0.129
>> Residual 0.17198
>>
>> However, I find difficult to interpret the fact that there are "two 
>> intercept".
>> Any suggestion would be appreciated.
>>
>> It seems to me that a model with one random intercept for
>> the subject + one random effects for each treatment relative to the
>> baseline is easire to interpret. The syntax to get it is however not
>> straightforward and requires coding dummy variables by hand.
>>
>> > tmp<-CP0vf
>> > tmp$vf.novf<-ifelse(tmp$vf=="noVF",1,0)
>> > tmp$task.ltstat<-ifelse(tmp$task=="ltstat",1,0)
>> > tmp$task.ltmov<-ifelse(tmp$task=="ltsmov",1,0)
>> > fit<-lmer(log(e.area)~group*task*vf + (1|su0) + (task.ltstat-1|su0) 
>> + (task.ltmov-1|su0) + (vf.novf-1|su0),tmp)
>> > VarCorr(fit)
>> Groups Name Std.Dev.
>> su0 (Intercept) 0.43475
>> su0.1 task.ltstat 0.23228
>> su0.2 task.ltmov 0.37742
>> su0.3 vf.novf 0.19578
>> Residual 0.20493
>>
>> or, to allow for some correlations,
>>
>> > fit<-lmer(log(e.area)~group*task*vf + (task.ltstat-1|su0) + 
>> (task.ltmov-1|su0) + (vf|su0),tmp)
>> > VarCorr(fit)
>> Groups Name Std.Dev. Corr
>> su0 task.ltstat 0.23284
>> su0.1 task.ltmov 0.20882
>> su0.2 (Intercept) 0.43813
>> vfnovf-vf 0.10022 0.090
>> Residual 0.20410
>>
>> > fit<-lmer(log(e.area)~group*task*vf + (task|su0) + (vf.novf-1|su0), 
>> tmp)
>> > VarCorr(fit)
>> Groups Name Std.Dev. Corr
>> su0 (Intercept) 0.45837
>> taskltstat-stand 0.16109 0.413
>> taskltdyn-stand 0.12072 -0.698 -0.674
>> su0.1 vf.novf 0.20669
>> Residual 0.17932
>>
>> I am a bit confused about how to reason about these different 
>> possibilities and
>> I wonder what to do in the spirit of Barr et al. (2013) who suggest 
>> not to try
>> to model the random structure but to keep it maximal. To quote them:
>>
>> "The generalization performance of LMEMs including data-driven random 
>> effects structures
>> strongly depends upon modeling criteria and sample size, yielding 
>> reasonable results on
>> moderately-sized samples when conservative criteria are used, but 
>> with little or no power
>> advantage over maximal models."
>>
>> Gabriel
>>
>>
>>
>>
>> -- 
>> ---------------------------------------------------------------------
>> Gabriel Baud-Bovy               tel.: (+39) 02 2643 4839 (office)
>> UHSR University                       (+39) 02 2643 3429 (laboratory)
>> via Olgettina, 58                     (+39) 02 2643 4891 (secretary)
>> 20132 Milan, Italy               fax: (+39) 02 2643 4892
>>
>> _______________________________________________
>> R-sig-mixed-models@r-project.org 
>> <mailto:R-sig-mixed-models@r-project.org> mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>


-- 
---------------------------------------------------------------------
Gabriel Baud-Bovy               tel.: (+39) 02 2643 4839 (office)
UHSR University                       (+39) 02 2643 3429 (laboratory)
via Olgettina, 58                     (+39) 02 2643 4891 (secretary)
20132 Milan, Italy               fax: (+39) 02 2643 4892
---------------------------------------------------------------------


	[[alternative HTML version deleted]]

