[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