[R-sig-ME] [R] updating arguments of formulae

Fri Dec 11 09:27:29 CET 2009

Moreno,

to my understanding, 

(1) depM ~ OS + VR + (1 + OS + VR | Sb2)
(2) depM ~ OS + VR + (1 |Sb2) + (1 + OS | Sb2) + (1 + VR | Sb2)

will not yield the same results. In (1), you model a random intercept
and slopes for OS and VR once for each group of Sb2. In (2), you model
an intercept, then again an intercept together with a slope for OS, and
then the same once more for VR. If anything, you might want to model the
three random effects as independent, then you should use something like 

(3) depM ~ OS + VR + (1 |Sb2) + (0 + OS | Sb2) + (0 + VR | Sb2)

hence not estimating three random intercepts for each group of Sb2 (as
is done by (2)). I would not know how to build the models
semi-automatically, for me, there are a lot of assumptions in it that I
may or may not justify conceptually or based on the data. There are a
lot of choices, and different people might make different ones at times,
so I doubt that this could be automatized. 

Try to look at 

require(lme4)
vignette("Implementation")

In Section 3 and in particular pages 15-16, you find a similar model to
(3), and there is some motivation why this model was chosen.

I hope that I have not completely misinterpreted lme4 now and hope that
some of the experts in the area will correct me in case.

HTH, Michael

> -----Original Message-----
> From: Moreno Ignazio Coco [mailto:M.I.Coco at sms.ed.ac.uk] 
> Sent: Donnerstag, 10. Dezember 2009 18:43
> To: Meyners,Michael,LAUSANNE,AppliedMathematics
> Cc: R-help at r-project.org; r-sig-mixed-models at r-project.org
> Subject: RE: [R] updating arguments of formulae
> 
> Michael,
> 
> Thanks a lot for your reply, I have now understood how to 
> fiddle around with the formulae updates...my question (see my 
> previous e-mail where I was sketching this problem out) about 
> LME models remains open...
> whether:
> 
> depM ~ (1 |Sb2) + OS + (1 + OS | Sb2) + VR + (1 + VR | Sb2)
> 
> is equivalent to:
> 
> depM ~ OS + VR + (1 + OS + VR | Sb2)
> 
> and if probably not what is the best approach to it and where 
> I can find a kind of guideline/rule of thumb list to build 
> "semi-automatically" linear mixed effect models with fixed 
> effects and random intercepts/slopes on it.
> 
> I am putting in copy the group  you suggested me...
> 
> Thanks again,
> 
> Moreno
> 
> Quoting "Meyners,Michael,LAUSANNE,AppliedMathematics"  
> <Michael.Meyners at rdls.nestle.com>:
> 
> > Moreno,
> >
> > I leave the discussion on the mixed models to others (you might 
> > consider the SIG group on mixed models as well for this), but try a 
> > few hints to make your code more accessible:
> >
> > * The "." in updating a formula is substituted with the 
> respective old 
> > formula (depending on the side), but is not mandatory. You 
> could give 
> > the new formula explicitly, i.e. consider something like
> > model1 = update(model, . ~ (1 |Sb2) + OS) if you loose 
> control about 
> > your models. See ?update.formula
> >
> > * I don't see the need for using your construct with 
> > as.formula(paste()), this makes things unnecessarily 
> complicated. See 
> > my above example, which should work as well on your data (and see 
> > ?update)
> >
> > * There is also the "-" operator available in 
> update.formula to remove 
> > terms (because it uses formula, see ?formula). As to your 
> question on 
> > how to move from
> >> depM ~ OS + (1 + OS | Sb2)
> >> to
> >> depM ~ OS + VR + (1 + OS + VR | Sb2)
> > try something like
> > update(model1, .~. - (1 + OS|Sb2) + VR + (1 + OS + VR | 
> Sb2)) while it 
> > goes without saying that in this case, it would be easier 
> to drop the 
> > "." and use something like update(model1, .~ OS + VR + (1 + 
> OS + VR | 
> > Sb2)) directly.
> >
> > * paste accepts more than just two arguments to be pasted: Try 
> > somthing like
> > model2 = update(model1, as.formula(paste(". ~ . + (1 + ", 
> "OS", "|" , 
> > "Sb2", ")")) instead of your construct with several nested calls to 
> > paste, and see ?paste. (Note that I added quotes to "OS" 
> and "Sb2", it 
> > didn't work for me otherwise as I have no object OS, not sure what 
> > happens if you have such an object on our search path, but I would 
> > suspect you encounter problems as well.)
> >
> > If you work yourself through these and thereby simplify 
> your code, you 
> > are more likely to get responses to your questions on which 
> model to 
> > use (which is actually independent from the use of update). 
> As far as 
> > I see it, it doesn't make sense to use a formula like in 
> your model4, 
> > but the mixed model experts might tell me wrong (and I got 
> a bit lost 
> > in your code as well). Please also try to provide 
> commented, minimal, 
> > self-contained, reproducible code for further enquiries 
> (use e.g. one 
> > of the examples on ?lmer to create appropriate examples for your 
> > questions).
> >
> > HTH, Michael
> >
> >
> >> -----Original Message-----
> >> From: r-help-bounces at r-project.org
> >> [mailto:r-help-bounces at r-project.org] On Behalf Of Moreno Ignazio 
> >> Coco
> >> Sent: Donnerstag, 10. Dezember 2009 13:35
> >> To: R-help at r-project.org
> >> Subject: [R] updating arguments of formulae
> >>
> >> Dear R-Community,
> >>
> >> I am relatively new with R, so sorry for things which for 
> you might 
> >> be obvious...
> >> I am trying to automatically update lmer formulae.
> >>
> >> the variables of the model are:
> >>
> >> depM= my dependent measure
> >> Sb2= a random factor
> >> OS = a predictor
> >> VR= another predictor
> >>
> >> So, I am building the first model with random intercept only:
> >>
> >> model = lmer(depM ~ (1 |Sb2))
> >>
> >> then I update the formula adding the first predictor
> >>
> >> model1 = update(model, as.formula(paste(". ~ . + ", OS)))
> >>
> >> the resulting formula will be:
> >>
> >> depM ~ (1 |Sb2) + OS
> >>
> >> let suppose now I want to update the model to have OS both 
> as a fixed 
> >> effect and in the random term, something like:
> >>
> >> depM ~ (1 + OS |Sb2) + OS
> >>
> >> I can do something very ugly (please tell me if there is a more 
> >> elegant way to do it) that looks like:
> >>
> >> model2 = update(model1, as.formula(paste(paste(paste(paste(".
> >> ~ . + (1
> >> + ", OS), "|" ), Sb2), ")")))
> >>
> >> the resulting model2 formula will be:
> >>
> >> depM ~ (1 |Sb2) + OS + (1 + OS | Sb2)
> >>
> >> one first thing I am wondering at this point is whether having
> >> (1 |Sb2) and (1 + OS | Sb2) in the same expression is redundant.
> >> in the output it will obviously tell me that group Sb2 is 
> considered 
> >> twice:
> >>
> >> number of obs: 6514, groups:  Sb2, 23; Sb2, 23
> >>
> >> and i am not sure if am doing it correctly...any advice?
> >>
> >> So let suppose now I want to add the new predictor VR 
> again both in 
> >> the fixed and in the random part of the formula.
> >> If i just repeat the two steps above:
> >>
> >> model3 = update(model2, as.formula(paste(". ~ . + ", VR)))
> >>
> >> and then:
> >>
> >> model4 = update(model3, as.formula(paste(paste(paste(paste(".
> >> ~ . + (1
> >> + ", VR), "|" ), Sb2), ")")))
> >>
> >> the formula I get is:
> >>
> >> depM ~ (1 |Sb2) + OS + (1 + OS | Sb2) + VR + (1 + VR | Sb2)
> >>
> >> so, basically I am adding new stuff on the right side of the 
> >> formula...
> >>
> >> My first question at this point is whether the above formula is 
> >> equivalent to:
> >>
> >> depM ~ OS + VR + (1 + OS + VR | Sb2)
> >>
> >> if is not equivalent, which one of the two is correct?
> >>
> >> obviously in the second case, group Sb2, is considered only once.
> >>
> >> If the second version of the formula is the correct one, I don't 
> >> understand how I can update arguments inside the formula 
> rather than 
> >> adding things on his right side...
> >>
> >> thus, in the ideal case,  how do I go from something like this:
> >>
> >> depM ~ OS + (1 + OS | Sb2)
> >>
> >> to something like this:
> >>
> >> depM ~ OS + VR + (1 + OS + VR | Sb2)
> >>
> >> Thanks a lot for your help,
> >> Best,
> >>
> >> Moreno
> >> --
> >> The University of Edinburgh is a charitable body, registered in 
> >> Scotland, with registration number SC005336.
> >>
> >> ______________________________________________
> >> R-help at r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.