[R] nested ANCOVA: still confused
Berton Gunter
gunter.berton at gene.com
Thu Jan 26 01:10:23 CET 2006
Jeff:
> However, I have two remaining questions: (1)how concerned should I be
> with the warning message below and
There was a definitive comment on this just a few days ago on the list
(search the archives), the gist of it was: **very concerned** . "False
convergence" means that you're not truly converged, the details for which
I've forgotten (sigh...). Anyway, this means that your parameter estimates
could be far from the correct minimized values = you could be in trouble. I
can't help you any more than that, but hopefully you'll get responses from
those with suitable expertise who can.
Cheers,
Bert
(2) is there a way to invoke output
> to get an estimate of the effect of purban2 (the proportion of urban
> cover 200 m around a box) on feather color (rtot) and if there is a
> difference between the sexes? I used the summary function and it
> doesn't tell me much (see output below).
>
> I'll read up mixed models when Pinheiro arrives but any
> suggestions for
> diagnostics? I'm going to repeat this study and expand it by doubling
> or tripling the number of birds.
>
> Warning message:
> nlminb returned message false convergence (8)
> in: "LMEoptimize<-"(`*tmp*`, value = list(maxIter = 200, tolerance =
> 1.49011611938477e-08,
>
> > summary(eabl)
> Linear mixed-effects model fit by REML
> Formula: rtot ~ sexv + (purban2 | clutch)
> Data: bb
> AIC BIC logLik MLdeviance REMLdeviance
> 5164.284 6997.864 -2052.142 4128.792 4104.284
> Random effects:
> Groups Name Variance Std.Dev. Corr
>
>
>
>
> clutch (Intercept) 502829 709.10
>
>
>
>
> purban20 1341990 1158.44 -0.477
>
>
>
>
> purban20.006711409 5683957 2384.10 -0.226 0.082
>
>
>
>
> purban20.01342282 1772922 1331.51 -0.386 0.176 0.067
> .
> .
> .
> .
> .
> # of obs: 235, groups: clutch, 74
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 5950.01 241.59 24.628
> sexvm 1509.07 145.73 10.355
>
> Correlation of Fixed Effects:
> (Intr)
> sexvm -0.304
>
> Thanks many time over,
>
> Jeff
>
> ****************************************
> Jeffrey A. Stratford, Ph.D.
> Postdoctoral Associate
> 331 Funchess Hall
> Department of Biological Sciences
> Auburn University
> Auburn, AL 36849
> 334-329-9198
> FAX 334-844-9234
> http://www.auburn.edu/~stratja
> ****************************************
> >>> "Doran, Harold" <HDoran at air.org> 01/25/06 6:37 AM >>>
> OK, we're getting somewhere. First, it looks as though (by the error
> message) that you have a big dataset. My first recommendation
> is to use
> lmer instead of lme, you will see a significant benefit in terms of
> computional speed.
>
> For the model this would be
>
> lmer(rtot ~ sexv +(purban|box:chick) + (purban|box), bb,
> na.action=na.omit)
>
> Now, you have run out of memory. I don't know what operating
> system you
> are using, so go and see the appropriate FAQ for increasing memory for
> your OS.
>
> Second, I made a mistake in my reply. Your random statement should be
> random=~purban|box/chick denoting that chicks are nested in boxes, not
> boxes nested in chicks, sorry about that.
>
> Now, why is it that each chick within box has the same value
> for purban?
> If this is so, why are you fitting that as a random effect?
> It seems not
> to vary across individual chicks, right? It seems there is only an
> effect of box and not an effect for chicks. Why not just fit a random
> effect only for box such as:
>
> rtot.lme <- lme(fixed=rtot~sexv, random=~purban2|box,
> na.action=na.omit,bb)
>
> or in lmer
> lmer(rtot ~ sexv + (purban|box), bb, na.action=na.omit)
>
> Harold
>
>
>
> -----Original Message-----
> From: Jeffrey Stratford [mailto:stratja at auburn.edu]
> Sent: Tue 1/24/2006 8:57 PM
> To: Doran, Harold; r-help at stat.math.ethz.ch
> Cc:
> Subject: RE: [R] nested ANCOVA: still confused
>
> R-users and Harold.
>
> First, thanks for the advice; I'm almost there.
>
> The code I'm using now is
>
> library(nlme)
> bb <- read.csv("E:\\eabl_feather04.csv", header=TRUE)
> bb$sexv <- factor(bb$sexv)
> rtot.lme <- lme(fixed=rtot~sexv, random=~purban2|chick/box,
> na.action=na.omit, data=bb)
>
> A sample of the data looks like this
>
> box chick rtot purban2 sexv
> 1 1 6333.51 0.026846 f
> 1 2 8710.884 0.026846 m
> 2 1 5810.007 0.161074 f
> 2 2 5524.33 0.161074 f
> 2 3 4824.474 0.161074 f
> 2 4 5617.641 0.161074 f
> 2 5 6761.724 0.161074 f
> 4 1 7569.673 0.208054 m
> 4 2 7877.081 0.208054 m
> 4 4 7455.55 0.208054 f
> 7 1 5408.287 0.436242 m
> 10 1 6991.727 0.14094 f
> 12 1 8590.207 0.134228 f
> 12 2 7536.747 0.134228 m
> 12 3 5145.342 0.134228 m
> 12 4 6853.628 0.134228 f
> 15 1 8048.717 0.033557 m
> 15 2 7062.196 0.033557 m
> 15 3 8165.953 0.033557 m
> 15 4 8348.58 0.033557 m
> 16 2 6534.775 0.751678 m
> 16 3 7468.827 0.751678 m
> 16 4 5907.338 0.751678 f
> 21 1 7761.983 0.221477 m
> 21 2 6634.115 0.221477 m
> 21 3 6982.923 0.221477 m
> 21 4 7464.075 0.221477 m
> 22 1 6756.733 0.281879 f
> 23 2 8231.496 0.134228 m
>
> The error I'm getting is
>
> Error in logLik.lmeStructInt(lmeSt, lmePars) :
> Calloc could not allocate (590465568 of 8) memory
> In addition: Warning messages:
> 1: Fewer observations than random effects in all level 2 groups in:
> lme.formula(fixed = rtot ~ sexv, random = ~purban2 | chick/box,
> 2: Reached total allocation of 382Mb: see help(memory.size)
>
> There's nothing "special" about chick 1, 2, etc. These were
> simply the
> order of the birds measured in each box so chick 1 in box 1
> has nothing
> to do with chick 1 in box 2.
>
> Many thanks,
>
> Jeff
>
> ****************************************
> Jeffrey A. Stratford, Ph.D.
> Postdoctoral Associate
> 331 Funchess Hall
> Department of Biological Sciences
> Auburn University
> Auburn, AL 36849
> 334-329-9198
> FAX 334-844-9234
> http://www.auburn.edu/~stratja
> ****************************************
> >>> "Doran, Harold" <HDoran at air.org> 01/24/06 2:04 PM >>>
> Dear Jeff:
>
> I see the issues in your code and have provided what I think
> will solve
> your problem. It is often much easier to get help on this
> list when you
> provide a small bit of data that can be replicated and you state what
> the error messages are that you are receiving. OK, with that
> said, here
> is what I see. First, you do not need to use the syntax bb$sex in your
> model, this can be significantly simplified. Second, you do not have a
> random statement in your model.
>
> Here is your original model:
> lme(bb$rtot~bb$sex, bb$purban|bb$chick/bb$box, na.action=na.omit)
>
> Here is what it should be:
>
> lme(fixed = rtot~sex, random=~purban|chick/box, na.action=na.omit,
> data=bb)
>
> Notice there is a fixed and random call. You can simplify this as
>
> lme(rtot~sex, random=~purban|chick/box, na.action=na.omit, bb)
>
> Note, you can eliminate the "fixed=" portion but not the random
> statement.
>
> Last, if you want to do this in lmer, the newer function for mixed
> models in the Matrix package, you would do
>
> lmer(rtot~sex + (purban|box:chick) + (purban|box), na.action=na.omit,
> data=bb)
>
> Hope this helps.
> Harold
>
>
>
>
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of
> Jeffrey Stratford
> Sent: Tuesday, January 24, 2006 11:34 AM
> To: r-help at stat.math.ethz.ch
> Subject: [R] nested ANCOVA: still confused
>
> Dear R-users,
>
> I did some more research and I'm still not sure how to set up
> an ANCOVA
> with nestedness. Specifically I'm not sure how to express
> chicks nested
> within boxes. I will be getting Pinheiro & Bates (Mixed
> Effects Models
> in S and S-Plus) but it will not arrive for another two weeks from our
> interlibrary loan.
>
> The goal is to determine if there are urbanization (purban) effects on
> chick health (rtot) and if there are differences between
> sexes (sex) and
> the effect of being in the same clutch (box).
>
> The model is rtot = sex + purban + (chick)box.
>
> I've loaded the package lme4. And the code I have so far is
>
> bb <- read.csv("C:\\eabl\\eabl_feather04.csv", header=TRUE) bb$sex <-
> factor(bb$sex) rtot.lme <- lme(bb$rtot~bb$sex,
> bb$purban|bb$chick/bb$box,
> na.action=na.omit)
>
> but this is not working.
>
> Any suggestions would be greatly appreciated.
>
> Thanks,
>
> Jeff
>
>
>
>
>
>
>
>
> ****************************************
> Jeffrey A. Stratford, Ph.D.
> Postdoctoral Associate
> 331 Funchess Hall
> Department of Biological Sciences
> Auburn University
> Auburn, AL 36849
> 334-329-9198
> FAX 334-844-9234
> http://www.auburn.edu/~stratja
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
More information about the R-help
mailing list