[R] F-test degree of freedoms in lme4 ?

Christoph Buser buser at stat.math.ethz.ch
Thu Jan 12 10:07:28 CET 2006


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




More information about the R-help mailing list