Dear Wilhelm
There is an article, including a part about fitting linear mixed
models. There the problem with the degrees of freedom is
described.
You can have a look to the second link, too, discussing the
problem as well.
http://cran.r-project.org/doc/Rnews/Rnews_2005-1.pdf
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/67414.html
Regards,
Christoph Buser
--------------------------------------------------------------
Christoph Buser <buser at stat.math.ethz.ch>
Seminar fuer Statistik, LEO C13
ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
phone: x-41-44-632-4673 fax: 632-1228
http://stat.ethz.ch/~buser/
--------------------------------------------------------------
Wilhelm B. Kloke writes:
> I have a problem moving from multistratum aov analysis to lmer.
>
> My dataset has observations of ampl at 4 levels of gapf and 2 levels of bl
> on 6 subjects levels VP, with 2 replicates wg each, and is balanced.
>
> Here is the summary of this set with aov:
> >> summary(aov(ampl~gapf*bl+Error(VP/(bl*gapf)),hframe2))
> >
> >Error: VP
> > Df Sum Sq Mean Sq F value Pr(>F)
> >Residuals 5 531 106
> >
> >Error: VP:bl
> > Df Sum Sq Mean Sq F value Pr(>F)
> >bl 1 1700 1700 37.8 0.0017 **
> >Residuals 5 225 45
> >---
> >Signif. codes: 0 '***' 0.001 '**' 0.01 '*'
0.05 '.' 0.1 ' ' 1
> >
> >Error: VP:gapf
> > Df Sum Sq Mean Sq F value Pr(>F)
> >gapf 3 933 311 24.2 5.3e-06 ***
> >Residuals 15 193 13
> >---
> >Signif. codes: 0 '***' 0.001 '**' 0.01 '*'
0.05 '.' 0.1 ' ' 1
> >
> >Error: VP:bl:gapf
> > Df Sum Sq Mean Sq F value Pr(>F)
> >gapf:bl 3 93.9 31.3 3.68 0.036 *
> >Residuals 15 127.6 8.5
> >---
> >Signif. codes: 0 '***' 0.001 '**' 0.01 '*'
0.05 '.' 0.1 ' ' 1
> >
> >Error: Within
> > Df Sum Sq Mean Sq F value Pr(>F)
> >Residuals 48 318 7
> >
> This is mostly identical the analysis by BMDP 4V, except for the
> Greenhouse-Geisser epsilons, which are not estimated this way.
>
> I have to analyse a similar dataset, which is not balanced. So I need to
> change the method. Following Pinheiro/Bates p.90f, I tried
> > hf2.lme <-
lme(ampl~gapf*bl,hframe2,random=list(VP=pdDiag(~gapf*bl),bl=pdDiag(~gapf)))
> and some variations of this to get the same F tests generated. At least,
> I got the F-test on error stratum VP:bl this way, but not the other two:
> >> anova(hf2.lme)
> > numDF denDF F-value p-value
> >(Intercept) 1 78 764.86 <.0001
> >gapf 3 78 17.68 <.0001
> >bl 1 5 37.81 0.0017
> >gapf:bl 3 78 2.99 0.0362
>
> Then I tried to move to lmer.
> I tried to find something equivalent to the above lme call, with no
> success at all.
>
> In case, that the problem is in the data, here is the set:
>
> VP ampl wg bl gapf
> 1 WJ 22 w s 144
> 2 CR 23 w s 144
> 3 MZ 25 w s 144
> 4 MP 34 w s 144
> 5 HJ 36 w s 144
> 6 SJ 26 w s 144
> 7 WJ 34 w s 80
> 8 CR 31 w s 80
> 9 MZ 33 w s 80
> 10 MP 36 w s 80
> 11 HJ 37 w s 80
> 12 SJ 32 w s 80
> 13 WJ 34 w s 48
> 14 CR 37 w s 48
> 15 MZ 38 w s 48
> 16 MP 38 w s 48
> 17 HJ 40 w s 48
> 18 SJ 32 w s 48
> 19 WJ 36 w s 16
> 20 CR 40 w s 16
> 21 MZ 39 w s 16
> 22 MP 40 w s 16
> 23 HJ 40 w s 16
> 24 SJ 38 w s 16
> 25 WJ 16 g s 144
> 26 CR 28 g s 144
> 27 MZ 18 g s 144
> 28 MP 33 g s 144
> 29 HJ 37 g s 144
> 30 SJ 28 g s 144
> 31 WJ 28 g s 80
> 32 CR 33 g s 80
> 33 MZ 24 g s 80
> 34 MP 34 g s 80
> 35 HJ 36 g s 80
> 36 SJ 30 g s 80
> 37 WJ 32 g s 48
> 38 CR 38 g s 48
> 39 MZ 34 g s 48
> 40 MP 37 g s 48
> 41 HJ 39 g s 48
> 42 SJ 30 g s 48
> 43 WJ 36 g s 16
> 44 CR 34 g s 16
> 45 MZ 36 g s 16
> 46 MP 40 g s 16
> 47 HJ 40 g s 16
> 48 SJ 36 g s 16
> 49 WJ 22 w b 144
> 50 CR 24 w b 144
> 51 MZ 20 w b 144
> 52 MP 26 w b 144
> 53 HJ 22 w b 144
> 54 SJ 16 w b 144
> 55 WJ 26 w b 80
> 56 CR 24 w b 80
> 57 MZ 26 w b 80
> 58 MP 27 w b 80
> 59 HJ 26 w b 80
> 60 SJ 18 w b 80
> 61 WJ 28 w b 48
> 62 CR 23 w b 48
> 63 MZ 28 w b 48
> 64 MP 29 w b 48
> 65 HJ 27 w b 48
> 66 SJ 24 w b 48
> 67 WJ 32 w b 16
> 68 CR 26 w b 16
> 69 MZ 30 w b 16
> 70 MP 28 w b 16
> 71 HJ 30 w b 16
> 72 SJ 22 w b 16
> 73 WJ 22 g b 144
> 74 CR 18 g b 144
> 75 MZ 18 g b 144
> 76 MP 26 g b 144
> 77 HJ 22 g b 144
> 78 SJ 18 g b 144
> 79 WJ 24 g b 80
> 80 CR 26 g b 80
> 81 MZ 30 g b 80
> 82 MP 26 g b 80
> 83 HJ 26 g b 80
> 84 SJ 24 g b 80
> 85 WJ 28 g b 48
> 86 CR 28 g b 48
> 87 MZ 27 g b 48
> 88 MP 30 g b 48
> 89 HJ 26 g b 48
> 90 SJ 16 g b 48
> 91 WJ 28 g b 16
> 92 CR 19 g b 16
> 93 MZ 24 g b 16
> 94 MP 32 g b 16
> 95 HJ 30 g b 16
> 96 SJ 22 g b 16
> --
> Dipl.-Math. Wilhelm Bernhard Kloke
> Institut fuer Arbeitsphysiologie an der Universitaet Dortmund
> Ardeystrasse 67, D-44139 Dortmund, Tel. 0231-1084-257
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html