[R] question in using nlme and lme4 for unbalanced data
Chi Yuan
cyuan at email.arizona.edu
Mon Nov 1 23:59:42 CET 2010
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?
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. 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?
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:cyuan at email.arizona.edu
Website: http://www.u.arizona.edu/~cyuan/
More information about the R-help
mailing list