# [R] question in using nlme and lme4 for unbalanced data

Chi Yuan cyuan at email.arizona.edu
Mon Oct 25 20:44:49 CEST 2010

```Hello:
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?
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 kind of heard the P value does not matter that much in the mixed
model. Is there any other way I can tell whether the treatment has a
significant effect or not?
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