[R-sig-ME] R-sig-mixed-models Digest, Vol 62, Issue 29
anthony.sealey at utoronto.ca
anthony.sealey at utoronto.ca
Mon Feb 13 21:40:09 CET 2012
9sbnopoi
-----Original Message-----
From: r-sig-mixed-models-request at r-project.org
Sender: r-sig-mixed-models-bounces at r-project.org
Date: Mon, 13 Feb 2012 20:07:11
To: <r-sig-mixed-models at r-project.org>
Reply-To: r-sig-mixed-models at r-project.org
Subject: R-sig-mixed-models Digest, Vol 62, Issue 29
Send R-sig-mixed-models mailing list submissions to
r-sig-mixed-models at r-project.org
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
or, via email, send a message with subject or body 'help' to
r-sig-mixed-models-request at r-project.org
You can reach the person managing the list at
r-sig-mixed-models-owner at r-project.org
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-mixed-models digest..."
Today's Topics:
1. Comparing against a negative control in an LMM (Masca, Nick)
2. MCMCglmm with cross-classified random effects (Agostino Moro)
3. Re: MCMCglmm with cross-classified random effects
(Jarrod Hadfield)
4. Re: Considerable discrepancies between fixed and random
effect estimates of lme4 (glmer) and glmmADMB (Ben Bolker)
5. Any package for best subset selection for random effects
model (Tao Zhang)
6. Interpretation of nonlinear mixed-effects modeling results
(Gang Chen)
7. Re: Considerable discrepancies between fixed and random
effect estimates of lme4 (glmer) and glmmADMB (Adam Smith)
----------------------------------------------------------------------
Message: 1
Date: Mon, 13 Feb 2012 12:33:42 +0000
From: "Masca, Nick" <Nick.Masca at effem.com>
To: "r-sig-mixed-models at r-project.org"
<r-sig-mixed-models at r-project.org>
Subject: [R-sig-ME] Comparing against a negative control in an LMM
Message-ID:
<8295A4D50D4C644CAC4323DD070D9597106F2D at 034-CH1MPN1-014.034d.mgd.msft.net>
Content-Type: text/plain
Hi all,
I have a problem based on a colleague's experiment that I've been asked to analyse, which is more of a general mixed modelling issue rather than specifically an R issue, and I would be extremely grateful for any help that any readers of this list can provide.
An experiment was conducted in which the aim was to compare 3 concentrations of 2 active treatments (i.e. 6 active treatments in total) to a negative control. Three batches of each of the actives have been tested, and 3 reps tested for each batch. In contrast, 20 replicates have been taken of the negative control - but, by definition, there is no "batch" for this treatment.
Here is some code to reproduce the experimental design:
Treat<- factor(c(rep("NC", 20), rep("A", 27), rep("B", 27)))
Conc<-factor(c(rep(1, 20), rep(1:3, each=9), rep(1:3, each=9)))
Batch<-factor(c(rep(1, 20), rep( rep(1:3, each=3), 6)))
Treatment<-factor(Treat:Conc) #specify new treatment variable (so don't attempt to estimate Conc. 2&3 for NC)
I originally planned to analyses these data in a LMM, with Treat*Conc as a 7 level fixed effect (i.e. 3*2 actives + control), and with Treat:Conc:Batch as random. The following code simulates my response variable assuming this model:
Resp<- rep(9, 74) + #simulate intercept
c( rep(rnorm(1, 0, sd=2.5), 20)^2, rep(rnorm(18, 0, sd=2.5), each=3)^2) + #simulate treat.conc.batch variance
rep(rnorm(74, 0, sd=.2)^2) + #simulate residual variance
c(rep(0,20), rep(c(-4, 0,0,-4, 0,0), each= 9)) #simulate fixed effects
Data<-data.frame(Treatment, Conc, Batch, Resp)
While this code models the data using lme4:
Mod<-lmer(Resp ~ Treatment + (1|Treatment:Batch), data=Data)
I can now obtain and plot treatment means/CIs using glht in the multcomp package:
library(multcomp)
Mean.mat<-diag(rep(1,7))
Mean.mat[,1]<-rep(1,7)
rownames(Mean.mat)<-levels(Data$Treatment)
Est.means<-glht(Mod, Mean.mat)
plot(Est.means)
Hopefully from the above plot you can see what my issue is. The negative control, which I want to compare everything against, has by far the least precision around its estimate, despite the data for the control hardly varying at all. This happens because the greatest source of variability in the model (by far) is the variability between batches, but different batches of the negative control don't exist. As such, I'm not sure that this is a fair way to model the data, because the negative control is unfairly penalised by the variability between the batches of the other treatments.
I imagine that this kind of problem isn't particularly uncommon, but it's the first time I've had to deal with something like this myself. The only potential solution I've come up with so far is to scrap the negative control from the model, and simply subtract the negative control's mean "count" from all other values (either by specifying this mean as an offset or by subtracting it from all data-points). But this will probably give "anti-conservative" results, as it would assume the mean for the negative control doesn't vary.
I would be extremely grateful if anyone would care to share their thoughts on possible solutions to this problem - and whether anyone has dealt with this kind of issue before. I feel that I may well be missing something obvious - but can't see at the moment how else to get around it!
Many thanks for any help you can provide.
Cheers,
Nick
[[alternative HTML version deleted]]
------------------------------
Message: 2
Date: Mon, 13 Feb 2012 15:02:20 +0000
From: Agostino Moro <agostino.moro99 at gmail.com>
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] MCMCglmm with cross-classified random effects
Message-ID:
<CAMS_pxvdSZVhe_qSFgqsnkMofyTGxL6eLAY4g114VzS=k9HFpA at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Dear R-users,
I would like to fit ?a glmm with cross-classified random effects with
the function MCMCglmm. Something along the lines:
model1<-MCMCglmm(response~pred1, random=~re1+re2, data=data)
where re1 and re2 should be crossed random effects. I was wondering
whether you could tell me specifying cross-classified random effects
in MCMCglmm requires a particular syntax? Are there any examples
somewhere? I have had a look at the manual and the package vignette,
but I have not been able to find any examples relevant to what I want
to do.
Thanks,
Agostino
------------------------------
Message: 3
Date: Mon, 13 Feb 2012 15:19:07 +0000
From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
To: Agostino Moro <agostino.moro99 at gmail.com>
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] MCMCglmm with cross-classified random effects
Message-ID: <20120213151907.57957fqy2tmlxcg8 at www.staffmail.ed.ac.uk>
Content-Type: text/plain; charset=ISO-8859-1; DelSp="Yes";
format="flowed"
Hi,
As long as the levels of re1 and re2 are uniquely labelled any cross
classification will be dealt with appropriately.
Cheers,
Jarrod
Quoting Agostino Moro <agostino.moro99 at gmail.com> on Mon, 13 Feb 2012
15:02:20 +0000:
> Dear R-users,
>
> I would like to fit ?a glmm with cross-classified random effects with
> the function MCMCglmm. Something along the lines:
>
> model1<-MCMCglmm(response~pred1, random=~re1+re2, data=data)
>
> where re1 and re2 should be crossed random effects. I was wondering
> whether you could tell me specifying cross-classified random effects
> in MCMCglmm requires a particular syntax? Are there any examples
> somewhere? I have had a look at the manual and the package vignette,
> but I have not been able to find any examples relevant to what I want
> to do.
>
> Thanks,
>
> Agostino
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
------------------------------
Message: 4
Date: Mon, 13 Feb 2012 16:43:51 +0000 (UTC)
From: Ben Bolker <bbolker at gmail.com>
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Considerable discrepancies between fixed and
random effect estimates of lme4 (glmer) and glmmADMB
Message-ID: <loom.20120213T045329-772 at post.gmane.org>
Content-Type: text/plain; charset=utf-8
Adam Smith <raptorbio at ...> writes:
> I hope I'm not overlooking something elementary here, but estimated
> fixed and random effects are considerably different from the
> following Poisson model in lme4 and glmmADMB.? The fixed effects
> seem to differ most considerably...? Thanks for any thoughts...
> Adam Smith
> > sessionInfo()
> R version 2.14.1 (2011-12-22)
> Platform: x86_64-pc-mingw32/x64 (64-bit)
>
> (SNIP)
>
> other attached packages:
> [1] glmmADMB_0.7.2.5?? lme4_0.999375-42?? Matrix_1.0-3??????
bbmle_1.0.4.1????? numDeriv_2010.11-1
> [6] lattice_0.20-0???? R2admb_0.7.5?????? MASS_7.3-16??????
> > str(cons09) # The dataset
> 'data.frame':?? 394 obs. of? 15 variables:
> ?$ plot?????? : Factor w/ 16 levels "n_10","n_2","n_3",..:
> 2 2 2 2 2 2 2 2 2 2 ...
> ?$ plot_trt?? : Factor w/ 32 levels
> "cont_n_10","cont_n_2",..: 2 2 2 2 2 2 2 2 2 2 ...
> ?$ geog?????? : Factor w/ 2 levels "n","s": 1 1 1 1 1 1 1 1 1 1 ...
> ?$ trt??????? : Factor w/ 2 levels "cont","trt": 1 1 1 1 1 1 1 1 1 1 ...
> ?$ count????? : Factor w/ 14 levels "1","2","3","4",..:
> 1 2 3 4 5 6 7 8 9 10 ...
> ?$ total????? : int? 341 326 257 244 185 141 128 121 115 84 ...
> ?$ cons?????? : int? 12 52 8 57 36 8 0 1 20 27 ...
> ?$ dt???????? : int? 4 3 3 3 3 3 3 3 3 3 ...
> ?$ obs??????? : Factor w/ 394 levels "1","2","3","4",..:
> 1 2 3 4 5 6 7 8 9 10 ...
> ?$ logtotal?? : num? 5.83 5.79 5.55 5.5 5.22 ...
> ?$ logdt????? : num? 1.39 1.1 1.1 1.1 1.1 ...
> # The models
> > Poiss <- glmmadmb(cons ~ count + geog + trt + count:trt +
> count:geog + trt:geog + offset(logtotal) +
> ??? ??? offset(logdt) + (1|plot) + (1|plot_trt),
> zeroInflation=FALSE, family="poisson", data=cons09)
>
> > P_glmer <- glmer(cons ~ count + geog + trt +
> count:trt + count:geog + trt:geog + offset(log(total)) +
> ??? ? offset(log(dt)) + (1|plot) +
> (1|plot_trt), family="poisson", data=cons09)
By the way, you can specify these fixed effects
more compactly as (count+geog+trt)^2 ...
> > fixef(Poiss)
> 0.77333 0.58885 -0.21547 0.79101 0.24924 -0.02565 0.36004
> -0.94852 0.39442 0.29032 0.13658 -0.39564
> -1.37470 -0.88496 -0.53259 -0.25434
> -0.41033 1.18840 0.25678 0.69734
> 0.59823 1.11430 -0.17330 -0.15364
> 0.98442 0.51533 -0.99431 -1.63920 -0.40455 -1.55670 0.15891
> 0.35378 0.45748 0.92293 0.78373 1.05440
> 0.39598 0.32019 1.41600
> 0.86656 2.45790 1.07350 -0.31149
> # After detaching glmmADMB and lme4, then
> re-requiring lme4 to avoid masking of lme4's fixef function
> > fixef(P_glmer)
> -5.00326871 0.69554781 0.17105622 1.29467717
> 0.97268949 0.87707266 1.39399967 0.22923370
> 1.61829807 1.84546906 1.99506279 1.37859852
> 0.53987085 1.13282524 -0.12901214 0.07021103
> -0.52669841 0.87378323 0.05397134 0.46901053
> 0.39863837 0.94486828 -0.18146287 -0.20340655
> 0.77866847 0.49935778 -0.65101999 -1.05381279
> -0.05290984 -1.54103716 -0.03887633 0.04570463 0.02076874
> 0.47883964 0.23057821 0.49408747
> -0.15116087 -0.32683907 0.50098187 0.14670125
> 1.63961096 0.40944384 -0.22555478
>
> I juxtapose ranef() estimates here for comparison's sake...
>
[snip]
There's nothing obviously wrong here. It's not a full solution,
but I wonder how wide the confidence intervals are ... if they
are very wide, then the practical answer is that these are poorly
determined estimates. You're probably overfitting the model --
2 random effects plus 43 fixed-effect coefficients is
quite a lot for 394 observations (the general rule of thumb
is N/(# params)>10), especially if the Poisson data are sparse
(although they don't look that way from the first few
values listed in str() -- is there anyway you can allow
'count' to be continuous, or ordinal, rather than insisting
on it being categorical?) But it would be best to try to answer
more definitively ... can you send data?
Ben Bolker
------------------------------
Message: 5
Date: Mon, 13 Feb 2012 09:22:24 -0800
From: Tao Zhang <zt020200 at gmail.com>
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] Any package for best subset selection for random
effects model
Message-ID:
<CADMWRionUd97x-oUuUqj0nUStGiSZohudJLVh4kctEY0tC8=cw at mail.gmail.com>
Content-Type: text/plain
Hi Pros,
I know leaps() computes the best subset selection for linear model,
and
the bestglm() computes the best subset selection for generalized linear
model. Is there any package for best subset selection on random effects
model, or mixed effects model?
Thank you!
Tao
[[alternative HTML version deleted]]
------------------------------
Message: 6
Date: Mon, 13 Feb 2012 13:41:43 -0500
From: Gang Chen <gangchen6 at gmail.com>
To: r-sig-mixed-models <r-sig-mixed-models at r-project.org>
Subject: [R-sig-ME] Interpretation of nonlinear mixed-effects modeling
results
Message-ID:
<CAHmzXO6hydefDRtcZAgbuHBhiJXGfM6zDxDfPkP=k8LouY5suQ at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
I'm fitting a nonlinear mixed-effects model to some data with two
groups (controls and patients) with something like
fm <- nlme(response ~ myFunc(time, a, b), data=myData, fixed = a + b
~ group, start=...)
myFunc is a nonlinear function defined with two parameters a and b.
I'm very confused with the results between summary(fm) and anova(fm):
> summary(fm)
...
Fixed effects: a + b ~ group
Value Std.Error DF t-value p-value
a.(Intercept) 29.905889 10.532769 2196 2.839319 0.0046
a.groupPat 6.437218 16.045223 2196 0.401192 0.6883
b.(Intercept) 0.290943 0.072544 2196 4.010559 0.0001
b.groupPat -0.138361 0.077339 2196 -1.789010 0.0738
...
> anova(fm)
numDF denDF F-value p-value
a.(Intercept) 1 2196 497.8594 <.0001
a.group 1 2196 12.6109 0.0004
b.(Intercept) 1 2196 45.2787 <.0001
b.group 1 2196 3.2006 0.0738
If I understand it correctly, the last row in the fixed effects table
of summary(fm) is the difference in parameter b between the two
groups, and the t-statistic (and p-value) matches the F-statistic (and
p-value) from the last row of anova(fm): (-1.789010)^2 = 3.2006.
However, I'm totally at a loss for the other three rows in the two
tables? For example, I thought a.groupPat (2nd row) in the summary(fm)
table is the amount in parameter a in Patient group that is more than
parameter a in the Control group (1st row); but this interpretation is
not consistent with what is shown in the 2nd row of anova(fm) table.
What am I missing here?
Thanks,
Gang
------------------------------
Message: 7
Date: Mon, 13 Feb 2012 14:06:57 -0500
From: Adam Smith <raptorbio at hotmail.com>
To: <bbolker at gmail.com>, <r-sig-mixed-models at r-project.org>
Subject: Re: [R-sig-ME] Considerable discrepancies between fixed and
random effect estimates of lme4 (glmer) and glmmADMB
Message-ID: <BAY170-W29F87F0CC98BFF5E9B7F8EA17F0 at phx.gbl>
Content-Type: text/plain; charset="iso-8859-1"
> > I hope I'm not overlooking something elementary here, but estimated
> > fixed and random effects are considerably different from the
> > following Poisson model in lme4 and glmmADMB. The fixed effects
> > seem to differ most considerably... Thanks for any thoughts...
> > Adam Smith
>
> > > sessionInfo()
> > R version 2.14.1 (2011-12-22)
> > Platform: x86_64-pc-mingw32/x64 (64-bit)
> >
> > (SNIP)
> >
> > other attached packages:
> > [1] glmmADMB_0.7.2.5 lme4_0.999375-42 Matrix_1.0-3
> bbmle_1.0.4.1 numDeriv_2010.11-1
> > [6] lattice_0.20-0 R2admb_0.7.5 MASS_7.3-16
>
>
> > > str(cons09) # The dataset
> > 'data.frame': 394 obs. of 15 variables:
> > $ plot : Factor w/ 16 levels "n_10","n_2","n_3",..:
> > 2 2 2 2 2 2 2 2 2 2 ...
> > $ plot_trt : Factor w/ 32 levels
> > "cont_n_10","cont_n_2",..: 2 2 2 2 2 2 2 2 2 2 ...
> > $ geog : Factor w/ 2 levels "n","s": 1 1 1 1 1 1 1 1 1 1 ...
> > $ trt : Factor w/ 2 levels "cont","trt": 1 1 1 1 1 1 1 1 1 1 ...
> > $ count : Factor w/ 14 levels "1","2","3","4",..:
> > 1 2 3 4 5 6 7 8 9 10 ...
> > $ total : int 341 326 257 244 185 141 128 121 115 84 ...
> > $ cons : int 12 52 8 57 36 8 0 1 20 27 ...
> > $ dt : int 4 3 3 3 3 3 3 3 3 3 ...
> > $ obs : Factor w/ 394 levels "1","2","3","4",..:
> > 1 2 3 4 5 6 7 8 9 10 ...
> > $ logtotal : num 5.83 5.79 5.55 5.5 5.22 ...
> > $ logdt : num 1.39 1.1 1.1 1.1 1.1 ...
>
>
> > # The models
> > > Poiss <- glmmadmb(cons ~ count + geog + trt + count:trt +
> > count:geog + trt:geog + offset(logtotal) +
> > offset(logdt) + (1|plot) + (1|plot_trt),
> > zeroInflation=FALSE, family="poisson", data=cons09)
> >
> > > P_glmer <- glmer(cons ~ count + geog + trt +
> > count:trt + count:geog + trt:geog + offset(log(total)) +
> > offset(log(dt)) + (1|plot) +
> > (1|plot_trt), family="poisson", data=cons09)
>
> By the way, you can specify these fixed effects
> more compactly as (count+geog+trt)^2 ...
Indeed, I was being explicit for explicitness' sake...
>
> > > fixef(Poiss)
>
> > 0.77333 0.58885 -0.21547 0.79101 0.24924 -0.02565 0.36004
> > -0.94852 0.39442 0.29032 0.13658 -0.39564
> > -1.37470 -0.88496 -0.53259 -0.25434
> > -0.41033 1.18840 0.25678 0.69734
> > 0.59823 1.11430 -0.17330 -0.15364
> > 0.98442 0.51533 -0.99431 -1.63920 -0.40455 -1.55670 0.15891
> > 0.35378 0.45748 0.92293 0.78373 1.05440
> > 0.39598 0.32019 1.41600
> > 0.86656 2.45790 1.07350 -0.31149
>
> > # After detaching glmmADMB and lme4, then
> > re-requiring lme4 to avoid masking of lme4's fixef function
> > > fixef(P_glmer)
> > -5.00326871 0.69554781 0.17105622 1.29467717
> > 0.97268949 0.87707266 1.39399967 0.22923370
> > 1.61829807 1.84546906 1.99506279 1.37859852
> > 0.53987085 1.13282524 -0.12901214 0.07021103
> > -0.52669841 0.87378323 0.05397134 0.46901053
> > 0.39863837 0.94486828 -0.18146287 -0.20340655
> > 0.77866847 0.49935778 -0.65101999 -1.05381279
> > -0.05290984 -1.54103716 -0.03887633 0.04570463 0.02076874
> > 0.47883964 0.23057821 0.49408747
> > -0.15116087 -0.32683907 0.50098187 0.14670125
> > 1.63961096 0.40944384 -0.22555478
> >
> > I juxtapose ranef() estimates here for comparison's sake...
> >
> [snip]
>
> There's nothing obviously wrong here. It's not a full solution,
> but I wonder how wide the confidence intervals are ... if they
> are very wide, then the practical answer is that these are poorly
> determined estimates. You're probably overfitting the model --
> 2 random effects plus 43 fixed-effect coefficients is
> quite a lot for 394 observations (the general rule of thumb
> is N/(# params)>10), especially if the Poisson data are sparse
> (although they don't look that way from the first few
> values listed in str() -- is there anyway you can allow
> 'count' to be continuous, or ordinal, rather than insisting
> on it being categorical?) But it would be best to try to answer
> more definitively ... can you send data?
I started with a general specification of count expecting to model
it with an additive term when I start comparing fixed effects.? I
suppose I could model it as a lower order polynomial initially...
Doing so should drastically reduce the number of fixed effects in
the model.? I suppose I should center it before creating the
polynomial terms...
However, I'll still send the data off-list.?
Thanks for looking this over.
>
> Ben Bolker
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
------------------------------
_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
End of R-sig-mixed-models Digest, Vol 62, Issue 29
**************************************************
More information about the R-sig-mixed-models
mailing list