[R] Difficulties in trying to do a mixed effects model using the lmer function

Bert Gunter gunter.berton at gene.com
Wed Oct 3 19:08:25 CEST 2012


Please post this on the r-sig-mixed-models list, not here.  You are
waayyy more likely to get useful help there.

Also this is primarily a statistics, not an R issue. You may wish to
consider consulting with your local statistician to help you
understand the statistics.

Cheers,
Bert

On Wed, Oct 3, 2012 at 8:54 AM, Riedel  Judith
<judith.riedel at ipw.agrl.ethz.ch> wrote:
> Dear people of the help list
>
> I am drying to analyze my data using the 'lmer' function and I keep having problems.
>
> This is the model:
>> fm1<-lmer(dbh~spec+scheme+(1|Plot),data=d, REML=FALSE).
>
> I analyse tree size (dbh) of  3 different species (spec) and  3 planting schemes (scheme). I have 5 plots, which I hope to model as a random factor. (However, the subsequent output is based on some simplified dummy data, which is based on only two plots and ha only few observations).
>
> No I do:
>> anova(fm1)
> and I get some output, which I don't understand. Looks like this:
>
> Analysis of Variance Table
>        Df Sum Sq Mean Sq F value
> spec    2 6.098 3.0490  0.6142
> scheme  2 13.161  6.5803  1.3255
>
> The problems I have are:
> (1) How can I get the P-values?
> (2) How can I get the overall model statistic?
>
> Than I do:
>> summary(fm1)
>
> and get:
> Linear mixed model fit by maximum likelihood
> Formula: dbh ~ spec + scheme + (1 | Plot)
>    Data: d
>    AIC BIC logLik deviance REMLdev
>  147.2 157  -66.6    133.2   125.8
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  Plot     (Intercept) 0.0000   0.0000
>  Residual             4.9644   2.2281
> Number of obs: 30, groups: Plot, 2
>
> Fixed effects:
>             Estimate Std. Error t value
> (Intercept)   6.9074     0.9424   7.329
> specCED       0.3859     1.0265   0.376
> specTAB       0.8585     0.9828   0.874
> schemeMON     0.6572     0.9554   0.688
> schemePRO    -1.0344     1.1259  -0.919
>
> Correlation of Fixed Effects:
>           (Intr) spcCED spcTAB schMON
> specCED   -0.537
> specTAB   -0.529  0.500
> schemeMON -0.588  0.002 -0.072
> schemePRO -0.565  0.064  0.063  0.510
>
> What is this? What does it tell me?
>
> The statistics help advised me to do a second model, like this:
>> fm2<-lmer(dbh~scheme+(1|Plot),data=d,REML=FALSE)
>> anova(fm1,fm2)
>
> But why would I compare the two models?
>
>  What I get is:
> Data: d
> Models:
> fm2: dbh ~ scheme + (1 | Plot)
> fm1: dbh ~ spec + scheme + (1 | Plot)
>     Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
> fm2  5 143.96 150.97 -66.982
> fm1  7 147.21 157.01 -66.602 0.7584      2     0.6844
>
> What does this mean? Why Chi?
>
> Finally I would like to do some LSD post hoc tests, but I have no idea how to do it.
>
> In the end I would like to be able to report something like: 'DBH differed significantly between, species, planting schemes, and plots (Fx,xx = X; P = X). DBH of species 1 was significantly larger than DBH of species 2 (LSD post hoc test, P = X)'.
>
> I greatly appreciate any suggestions! Thank You a lot for Your help!
>
> Kind regards,
>
> Judith
>
> PS. the complete output is attached.
>
>
>
> XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
> Judith Riedel
> ETH Zurich
> Institute of Agricultural Sciences
> Applied Entomology
> Schmelzbergstrasse 9/LFO
> 8092 Zurich
> Switzerland
>
> Tel:  ++41 44 632 3923
> Fax: ++41 44 632 1171
> judith.riedel at ipw.agrl.ethz.ch<mailto:judith.riedel at ipw.agrl.ethz.ch>
> http://www.em.ipw.agrl.ethz.ch
>
>
> ______________________________________________
> 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

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm




More information about the R-help mailing list