Hi: On Mon, Nov 1, 2010 at 3:59 PM, Chi Yuan <cy...@email.arizona.edu> wrote:
> Hello: > I need some help about using mixed for model for unbalanced data. I > have an two factorial random block design. It's a ecology > experiment. My two factors are, guild removal and enfa removal. Both > are two levels, 0 (no removal), 1 (removal). I have 5 blocks. But > within each block, it's unbalanced at plot level because I have 5 > plots instead of 4 in each block. Within each block, I have 1 plot > with only guild removal, 1 plot with only enfa removal, 1 plot for > control with no removal, 2 plots for both guild and enfa removal. I am > looking at how these treatment affect the enfa mortality rate. I > decide to use mixed model to treat block as random effect. So I try > both nlme and lme4. But I don't know whether they take the unbalanced > data properly. So my question is, does lme in nlme and lmer in lme4 > take unbalanced data? How do I know it's analysis in a proper way? > Unbalanced data is not a problem in either package. However, five blocks is rather at the boundary of whether or not one can compute reliable variance components and random effects. Given that the variance estimate of blocks in your models was nearly zero, you're probably better off treating them as fixed rather than random and analyzing the data with a fixed effects model instead. Another question is about p values. > I kind of heard the P value does not matter that much in the mixed > model because it's not calculate properly. No. p-values are not calculated in lme4 (as I understand it) because, especially in the case of severely unbalanced data, the true sampling distributions of the test statistics in small to moderate samples are not necessarily close to the asymptotic distributions used to compute the corresponding p-values. It's the (sometimes gross) disparity between the small-sample and asymptotic distributions that makes the reported p-values based on the latter unreliable, not an inability to calculate the p-value properly. I can assure you that Prof. Bates knows how to compute a p-value. Is there any other way I can > tell whether the treatment has a effect not? I know AIC is for model > comparison, > do I report this in formal publication? > As mentioned above, I would suggest analyzing this as a fixed effects problem. Since the imbalance is not too bad, and it is not unusual in field experiments to have more control EUs than treatment EUs within each level of treatment, a fixed effects analysis may be sufficient. It wouldn't hurt to consult with a local statistician to discuss the options. HTH, Dennis > Here is my code and the result for each method. > I first try nlme > library(nlme) > > > m=lme(enfa_mortality~guild_removal*enfa_removal,random=~1|block,data=com_summer) > It gave me the result as following > Linear mixed-effects model fit by REML > Data: com_summer > AIC BIC logLik > 8.552254 14.81939 1.723873 > > Random effects: > Formula: ~1 | block > (Intercept) Residual > StdDev: 9.722548e-07 0.1880945 > > Fixed effects: enfa_mortality ~ guild_removal * enfa_removal > Value Std.Error DF t-value p-value > (Intercept) 0.450 0.0841184 17 5.349603 0.0001 > guild_removal -0.100 0.1189614 17 -0.840609 0.4122 > enfa_removal -0.368 0.1189614 17 -3.093441 0.0066 > guild_removal:enfa_removal 0.197 0.1573711 17 1.251818 0.2276 > Correlation: > (Intr) gld_rm enf_rm > guild_removal -0.707 > enfa_removal -0.707 0.500 > guild_removal:enfa_removal 0.535 -0.756 -0.756 > > Standardized Within-Group Residuals: > Min Q1 Med Q3 Max > -1.7650706 -0.7017751 0.1594943 0.7974717 1.9139320 > > Number of Observations: 25 > Number of Groups: 5 > > > I then try lme4, it give similar result, but won't tell me the p value. > library(lme4) > m<-lmer(enfa_mortality ~ guild_removal*enfa_removal +(1|block), > data=com_summer) > here is the result > Linear mixed model fit by REML > Formula: enfa_mortality ~ guild_removal * enfa_removal + (1 | block) > Data: com_summer > AIC BIC logLik deviance REMLdev > 8.552 15.87 1.724 -16.95 -3.448 > Random effects: > Groups Name Variance Std.Dev. > block (Intercept) 0.000000 0.00000 > Residual 0.035380 0.18809 > Number of obs: 25, groups: block, 5 > > Fixed effects: > Estimate Std. Error t value > (Intercept) 0.45000 0.08412 5.350 > guild_removal -0.10000 0.11896 -0.841 > enfa_removal -0.36800 0.11896 -3.093 > guild_removal:enfa_removal 0.19700 0.15737 1.252 > > Correlation of Fixed Effects: > (Intr) gld_rm enf_rm > guild_remvl -0.707 > enfa_removl -0.707 0.500 > gld_rmvl:n_ 0.535 -0.756 -0.756 > > > I really appreciate any suggestion! > Thank you! > -- > Chi Yuan > Graduate Student > Department of Ecology and Evolutionary Biology > University of Arizona > Room 106 Bioscience West > lab phone: 520-621-1889 > Email:cy...@email.arizona.edu <email%3acy...@email.arizona.edu> > Website: http://www.u.arizona.edu/~cyuan/<http://www.u.arizona.edu/%7Ecyuan/> > > ______________________________________________ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
