[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