[R-sig-ME] "varFunc" classes

[Ben Bolker]

  That Simmonds and Higgins's example uses a numeric treatment (0/1)
variable is what makes this work.  If you have a categorical treatment
variable, if it's binary you should turn it into a numeric dummy
variable (you can use dummy() as shown in one of my earlier examples);
if it has >2 categories, this is going to be a little tricky. It makes
perfect sense that you want to allow for variation in the intercept
(baseline) among studies, but I don't quite know why S&H made the main
effect of study a fixed effect, which makes all of this harder ... ?

  Ben Bolker

On 16-12-15 12:39 PM, dan jackson wrote:
> Hi thanks for this Ben I am still at the playing stage (and I am new
> to all this and so am at the bewildered and learning stage) but I am
> trying to fit an individual patient data meta-analysis with
> continuous data for the first time.
> 
> In the paper "A general framework for the use of logistic regression
> models in meta-analysis", Simmonds and Higgins suggest fitting the
> model for binary outcome data
> 
> glmer(cbind(event,n-event) ∼ factor(study) + factor(treat) +
> (treat-1|study), data=HRT2, family=binomial(link="logit”))
> 
> (see their appendix) so I thought
> 
> fm10 <- lme(OUT ~ factor(study)+treat, data = DD, random = ~
> treat-1|study,weights=varIdent(form = ~ 1 |
> study),control=list(maxIter = 500, msMaxIter = 500))
> 
> Would just be a natural continuous version of this (allowing
> different residual variances in each study also, this is usually
> considered important, and also upping the maxits to get convergence).
> In Higgins and Simmonds treat is binary so +treat and +factor(treat)
> are the same thing in the fixed effect part of the model, the
> intuition is (as I understand it) you want a different baseline
> average in each study (so factor(study)) because this is bound to
> vary and also the treatment effect to differ randomly from one study
> to the next, to relax the strong assumption that the treatment effect
> is the same for each study.
> 
> I hope that helps explain my motivation, essentially I am "just"
> trying to "copy" what Simmonds and Higgins did for binary data in the
> continuous case.
> 
> Dan
> 
> 
> -----Original Message----- From: Ben Bolker
> [mailto:bbolker at gmail.com] Sent: 15 December 2016 16:59 To: Dan
> Jackson Cc: Dan Jackson; r-sig-mixed-models at r-project.org;
> daniel.jackson at mrc-bsu.cam.ac.uk Subject: Re: [R-sig-ME] "varFunc"
> classes
> 
> 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 | 
>> treat-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
>>>> 
>> 
>