Hi:

On Mon, Nov 1, 2010 at 3:59 PM, Chi Yuan <cyuan@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:cyuan@email.arizona.edu <Email%3Acyuan@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]]

