[R-sig-ME] "varFunc" classes

[Ben Bolker]

  It doesn't generally make sense to include a factor ("study" in your
case) as both a random and a fixed effect in the same model; I
appreciate that you're trying to include only the variability among
treatments within studies, but this is a little bit tricky from both a
practical and a conceptual point of view.  Is there a reason that
you're using

 OUT ~ factor(study)+treat, data = DD, random = ~ treat-1|study

and not

 OUT ~ treat, data=DD, random= ~treat|study

or

  OUT ~ treat, data=DD, random=~1|study/treat

?  The latter estimates the variation of treatments within studies as
compound symmetric, which is a little bit easier computationally.



On Thu, Dec 15, 2016 at 6:42 AM, Dan Jackson
<dan.jackson at mrc-bsu.cam.ac.uk> wrote:
> Thanks to everyone for their help, Ben and Wolfgang in particular. I managed
> to get some results for a toy-warm up example using
>
> fm10 <- lme(OUT ~ factor(study)+treat, data = DD, random = ~
> treat-1|study,weights=varIdent(form = ~ 1 | study),control=list(maxIter =
> 500, msMaxIter = 500))
>
> here DD is a dataframe containing simulated data: 10,000 rows of OUT
> ("outcome", just normally distributed data) treat (0/1 with probability 0.5
> for "success") and study contains the study number from 1 to 10, 1 for the
> first 1000 observations, 2 for the second 1000 and so on. For fun I made the
> first and last groups have larger variances to see if varIndent would pick
> this up -- which it clearly did which is wonderful.
>
> I would not have gotten this far without your help Ben and Wolfgang in
> particular. My only slight residual (no pun intended?!) concern is that I
> seem to be getting some O's for DOF for the which is a bit perplexing (for
> the the factor(study) s in the model). I have done some googling and this
> DOF issue seems to be a special issue in its own right. Since I am happy to
> use normal approximations I can get my inferences but this is something that
> might be worth thinking about, DOF=0 would appear to be just plain wrong to
> me!?!
>
> I am confident I will be able to fit this type of model to my real data in
> the new year. Thanks all, Merry Christmas, Dan
>
>
> On 14/12/2016 19:02, Ben Bolker wrote:
>>
>>    I'm too lazy to see if anyone has answered this one already, so ...
>>
>> On 16-12-13 03:48 AM, Dan Jackson wrote:
>>>
>>> Dear lme4 authors,
>>>
>>> I am sure you are very busy so I will just ask my question very quickly.
>>> I
>>> was reading the book "Mixed-effects models in S and S-plus" by Pinheiro
>>> and
>>> Bates. On the top of page 208 of this book, there is a Table 5.1 that
>>> implements various "varFunc" classes. One of these classes would seem to
>>> be
>>> what I need for my data: varIdent - different variances per stratum. I do
>>> know that different subets in my data have very different variances you
>>> see,
>>> so I would need to include this.
>>>
>>> However this book relates to S-plus and I am not sure if this has been
>>> implemented in R, in the glmer package? My data are continuous so I would
>>> just need this for lmer (and not glmer). If it has not been implemented
>>> is
>>> there any "workaround"?
>>
>>    It has been implemented in R, but in the nlme package rather than the
>> lme4 package (which contains lmer and glmer).
>>
>>     Historical note: nlme is the package associated with Pinheiro and
>> Bates's 2000 book. R's first "stable beta" version (according to
>> Wikipedia) was released in 2000.  Doug Bates has been involved in R
>> since the beginning (or almost?).
>>
>>    If you need to do this in lme4::lmer, you can do it in a sort of
>> clunky way by setting up dummy variables for group differences and
>> adding an individual-level random effect, e.g.
>>
>>
>> data(Orthodont,package="nlme")
>>
>> O2 <- transform(Orthodont,
>>                  obs=seq(nrow(Orthodont))) ## observation-level variance
>>
>>
>> ## since Female var < Male var, we have to use a dummy for Male
>> ## to add extra variance for males (won't work the other way because
>> ## we can't have a negative variance)
>>
>> m1 <-lmer(distance ~ age + (age|Subject) +
>>           (0+dummy(Sex,"Male")|obs),
>>       control=lmerControl(check.nobs.vs.nlev="ignore",
>>                           check.nobs.vs.nRE="ignore"),
>>       O2)
>>
>> library(nlme)
>> m2 <- lme(distance~age,random=~age|Subject,
>>            data=Orthodont,
>>            weights=varIdent(form=~1|Sex))
>> summary(m2)
>>
>> all.equal(c(logLik(m1)),c(logLik(m2)))
>> all.equal(c(fixef(m1)),c(fixef(m2)),tolerance=1e-6)
>>
>>
>>
>>
>>> Thanks in advance for any advice, Dan Jackson
>>>
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>