I'm trying to use lmerTest to have p-values for my fixed effects. I have 4 different random intercepts, 3 crossed and one nested :
test.reml <- lmerTest::lmer(y ~ s1 + min + cot + min:cot + ge
+ vis + dur + mo + nps + dist + st1 + st2 + di1 + s1:cot
+ s1:min + s1:cot:min + s1:ge + s1:vis + s1:dur + s1:mo
+ s1:nps + s1:dist + s1:st1 + s1:st2 + s1:di1 + (1|Unique_key)
+ (s1-1|object) + (ns1-1|object)
+ (1|region), bdr, REML=1)
The objects are observed two times and the correlation between the two measures is introduced by the random effect on Unique_key, an unique identifier of the object i in a region j. Each object could be observed in any region. S1 is a binary variable that takes the value 1 if the observation is observed the first time period and 0 either. There is one random intercept for the first periode and one random intercept for the second period for each object. ns1 is in fact a binary variable that is the complement of s1 and s1 + ns1 = 1 for each observation.
I can fit the model and get the estimates and p-values with summary() :
summary( test.reml)
Linear mixed model fit by REML ['merModLmerTest']
Formula: y ~ s1 + min + cot + min:cot + ge
+ vis + dur + mo + nps + dist + st1 + st2 + di1 + s1:cot
+ s1:min + s1:cot:min + s1:ge + s1:vis + s1:dur + s1:mo
+ s1:nps + s1:dist + s1:st1 + s1:st2 + s1:di1 + (1|Unique_key)
+ (s1-1|object) + (ns1-1|object)
+ (1|region), bdr, REML=1)
Data: bdr
REML criterion at convergence: 204569.1
Random effects:
Groups Name Variance Std.Dev.
Unique_key (Intercept) 0.2023 0.4497
object s1 0.3528 0.5940
object.1 ns1 0.5954 0.7716
Region (Intercept) 0.7563 0.8697
Residual 0.1795 0.4237
Number of obs: 113396, groups: Unique_key , 58541; object, 1065; Region, 87
Fixed effects:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.7341569 0.2382673 28.263 < 2e-16 ***
s1 0.7391924 0.2004413 3.688 0.000233 ***
min -0.0067606 0.0171385 -0.394 0.694205
cot 0.1235093 0.0353693 3.492 0.000499 ***
ge2 -0.1535452 0.0800998 -1.917 0.055525 .
ge3 -0.2131246 0.0986559 -2.160 0.030982 *
ge4 -0.1032694 0.1115603 -0.926 0.354830
ge5 -0.1769347 0.1296558 -1.365 0.172663
ge6 0.0117401 0.1115897 0.105 0.916231
ge7 -0.2692483 0.1022565 -2.633 0.008589 **
vis2 -0.0928661 0.0607950 -1.528 0.126938
vis3 -0.3026112 0.1246595 -2.428 0.015375 *
dur2 0.1479195 0.0786369 1.881 0.060249 .
dur3 0.1406340 0.0809379 1.738 0.082590 .
dur4 0.2742243 0.0884301 3.101 0.001981 **
dur5 0.1946761 0.1065815 1.827 0.068059 .
mo2 -0.1168591 0.1256017 -0.930 0.352386
mo3 -0.0611162 0.1267657 -0.482 0.629824
mo4 -0.2725720 0.1263740 -2.157 0.031248 *
mo5 -0.6107000 0.1379264 -4.428 1.05e-05 ***
mo6 -0.3635142 0.1299799 -2.797 0.005260 **
mo7 -0.0899233 0.1275164 -0.705 0.480846
mo8 -0.2349548 0.1253422 -1.875 0.061140 .
mo9 -0.2624888 0.1263051 -2.078 0.037934 *
mo10 -0.2882749 0.1244404 -2.317 0.020724 *
mo11 -0.1702823 0.1356031 -1.256 0.209497
mo12 0.1989155 0.1322339 1.504 0.132819
nps 0.0278418 0.0010393 26.790 < 2e-16 ***
dist2 0.4065093 0.1118916 3.633 0.000294 ***
dist3 0.0155691 0.0906664 0.172 0.863693
dist4 -0.2910960 0.1595805 -1.824 0.068424 .
dist5 -0.1316553 0.0913394 -1.441 0.149782
dist6 0.0477956 0.0995679 0.480 0.631308
dist7 0.1383000 0.0981247 1.409 0.159011
dist8 -0.3985620 0.0886316 -4.497 7.69e-06 ***
dist9 -0.2036683 0.0799584 -2.547 0.011005 *
st11 -0.0258775 0.0591631 -0.437 0.661919
st21 0.0089230 0.0573352 0.156 0.876356
di11 -0.0910207 0.0838321 -1.086 0.277846
min:cot 0.0066210 0.0006195 10.688 < 2e-16 ***
s1:cot -0.1505670 0.0443186 -3.397 0.000694 ***
s1:min 0.0079478 0.0015051 5.280 1.29e-07 ***
s1:ge2 0.0329272 0.1007943 0.327 0.743948
s1:ge3 0.2150927 0.1241590 1.732 0.083367 .
s1:ge4 0.1786057 0.1404119 1.272 0.203526
s1:ge5 -0.0422380 0.1631757 -0.259 0.795780
s1:ge6 0.1372051 0.1404415 0.977 0.328717
s1:ge7 0.1343314 0.1287059 1.044 0.296755
s1:vis2 0.1354091 0.0765084 1.770 0.076913 .
s1:vis3 0.2449180 0.1568745 1.561 0.118637
s1:dur2 -0.0888179 0.0989573 -0.898 0.369547
s1:dur3 -0.0532473 0.1018481 -0.523 0.601167
s1:dur4 -0.1239068 0.1112907 -1.113 0.265696
s1:dur5 -0.1191069 0.1341435 -0.888 0.374705
s1:mo2 -0.1357615 0.1574365 -0.862 0.388618
s1:mo3 0.0130976 0.1588743 0.082 0.934306
s1:mo4 0.0343900 0.1579532 0.218 0.827669
s1:mo5 0.2257241 0.1732449 1.303 0.192761
s1:mo6 0.0500347 0.1628755 0.307 0.758728
s1:mo7 -0.0451271 0.1596277 -0.283 0.777435
s1:mo8 -0.0200467 0.1572383 -0.127 0.898564
s1:mo9 0.0394005 0.1584268 0.249 0.803620
s1:mo10 0.0641038 0.1562518 0.410 0.681662
s1:mo11 -0.3136235 0.1703456 -1.841 0.065764 .
s1:mo12 -0.7003775 0.1660455 -4.218 2.58e-05 ***
s1:nps -0.0095428 0.0013077 -7.297 4.31e-13 ***
s1:dist2 -0.3867962 0.1407463 -2.748 0.006050 **
s1:dist3 -0.0516400 0.1140519 -0.453 0.650762
s1:dist4 -0.0567491 0.2008542 -0.283 0.777562
s1:dist5 0.0025780 0.1147143 0.022 0.982073
s1:dist6 -0.1456445 0.1252219 -1.163 0.244940
s1:dist7 -0.0452712 0.1234110 -0.367 0.713785
s1:dist8 0.0546400 0.1114865 0.490 0.624117
s1:dist9 0.0540697 0.1000415 0.540 0.588934
s1:st11 0.0784027 0.0744677 1.053 0.292549
s1:st21 -0.0394419 0.0721720 -0.546 0.584788
s1:di11 0.0463040 0.1055326 0.439 0.660882
s1:min:cot -0.0012850 0.0006004 -2.140 0.032344 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
but with anova() I get :
type3.bonmodele <- lmerTest::anova(test.reml, ddf="Satterthwaite")
Analysis of Variance Table
Df Sum Sq Mean Sq F value
s1 1 7.385 7.385 41.1448
min 1 0.081 0.081 0.4536
cot 1 29.384 29.384 163.7026
ge 6 25.198 4.200 23.3968
vis 2 0.464 0.232 1.2929
dur 4 22.763 5.691 31.7042
mo 11 15.581 1.416 7.8914
nps 1 234.535 234.535 1306.6487
dist 8 18.547 2.318 12.9162
st1 1 0.034 0.034 0.1879
st2 1 0.058 0.058 0.3220
di1 1 0.261 0.261 1.4549
min:cot 1 22.537 22.537 125.5611
s1:cot 1 9.146 9.146 50.9555
s1:min 1 18.383 18.383 102.4171
s1:ge 6 5.152 0.859 4.7843
s1:vis 2 1.698 0.849 4.7311
s1:dur 4 2.829 0.707 3.9404
s1:mo 11 8.157 0.742 4.1312
s1:nps 1 10.102 10.102 56.2803
s1:dist 8 2.233 0.279 1.5550
s1:st1 1 0.188 0.188 1.0481
s1:st2 1 0.046 0.046 0.2560
s1:di1 1 0.035 0.035 0.1927
s1:min:cot 1 0.822 0.822 4.5804
When I try to remove the triple interaction the anova() function returns the p-values...I have also tried to split my data frame and to fit the model on half the data and anova() works well to.
There is no warning when I use the functions and I have also tried to change the ddf option and the method but nothing seems to work.
Here is my session info :
R version 3.0.0 (2013-04-03)
Platform: x86_64-w64-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=French_Canada.1252 LC_CTYPE=French_Canada.1252 LC_MONETARY=French_Canada.1252 LC_NUMERIC=C
[5] LC_TIME=French_Canada.1252
attached base packages:
[1] parallel splines stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggplot2_0.9.3.1 snow_0.3-13 Snowball_0.0-10 xtable_1.7-1 lmerTest_2.0-0 pbkrtest_0.3-7 MASS_7.3-29
[8] papeR_0.3 gmodels_2.15.4 survival_2.37-4 nlme_3.1-111 car_2.0-19 lme4_1.1-1 Matrix_1.1-0
[15] lattice_0.20-15
loaded via a namespace (and not attached):
[1] bitops_1.0-6 caTools_1.16 cluster_1.14.4 colorspace_1.2-4 dichromat_2.0-0 digest_0.6.3
[7] gdata_2.13.2 gplots_2.12.1 grid_3.0.0 gtable_0.1.2 gtools_3.0.0 Hmisc_3.12-2
[13] KernSmooth_2.23-10 labeling_0.2 minqa_1.2.1 munsell_0.4.2 nnet_7.3-7 numDeriv_2012.9-1
[19] plyr_1.8 proto_0.3-10 RColorBrewer_1.0-5 RCurl_1.95-4.1 reshape2_1.2.2 rJava_0.9-4
[25] ROAuth_0.9.3 rpart_4.1-3 scales_0.2.3 stringr_0.6.2 tools_3.0.0 twitteR_1.1.7
I can't share the data but I can add more infos if needed! I wanted to use the Satterthwaite approximation for the degrees of freedom but if you have any other suggestions to get p-values please share! Thank you very much!