R help - Oct 2010 - question in using nlme and lme4 for unbalanced data

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
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/

Chi:

I would say that these are not really R questions, but about
statistical procedures in general. However, there are some important
issues that you touch on, so I cannot resist commenting. Hopefully my
comments might pique others to add theirs -- and perhaps to disagree
with me.

1. Unbalanced data?
Yes that's what they're for. However, you need to  recognize that with
unbalanced data, interpretability of coefficients may be lost (the
estimates depend on what factors you include in your model).

2.Proper analysis?
TRUST in DOUG. (Doug Bates, the author of the packages and an
acknowledged expert in the field).

3. Without P values, how do I I tell what's significant?
Well, first of all, P values/statistical significance are only
meaningful when hypotheses are pre-specified; and even then,
statistical significance is a function of effect magnitude, sample
size, and experimental noise, so usually doesn't accord with
scientific importance anyway unless the experiments have been powered
appropriately, which is rarely the case outside clinical trials. So
statistical significance is probably irrelevant to the scientific
questions iin any case.

I would say that you need to focus on effect (coefficient) sizes and
perhaps do sensitivity analyses to see how much predicted results
change as you change them. Good graphs of your data are also always a
good idea.

4, With only 5 blocks, you have practically no information to estimate
variance anyway. You're probably better off treating them as fixed
effects and just estimating the model via lm. You might find a sqrt or
log transformation of the enfa numbers to be useful, though...

Cheers,
Bert

On Mon, Oct 25, 2010 at 11:44 AM, Chi Yuan <cyuan at email.arizona.edu>
wrote:
> 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
> 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/
>
> ______________________________________________
> R-help at 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.
>


-- 
Bert Gunter
Genentech Nonclinical Biostatistics

On 10-10-25 04:59 PM, Bert Gunter wrote:
>  ...ignore the block variation entirely, "
> 
> If the between block variability is large, this will lose precision;
> with imbalance, it could also result in bias (prhaps not in this
> study...). The mixed or fixed effects choice is arbitrary; this is not
> -- blocks should be modeled, not ignored.
> 
> -- Bert
   I agree that in general blocking factors should not be ignored!  I
don't even think they should be discarded if they are
small/non-significant (in the ecological literature this is called
'sacrificial pseudoreplication' sensu Hurlbert (1984)).

 But the point is that for this particular case the *estimated* block
variation is *exactly* zero for one function and trivially small for the
other (I think it's exactly zero for lmer, but not sure), so in fact the
original poster will (I think?) get exactly (or almost exactly) the same
answer if they simply use lm() and ignore blocks entirely ...

  Or am I missing something?

  (Taking the liberty of cc'ing back to the list so I can look like an
idiot in public if necessary ...)

  cheers
    Ben

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/