У меня высоко параметризованная модель, подходящая для одного случайного эффекта. Я хотел бы сравнить его соответствие (с lmer) с фиксированной моделью (с lm). Если я не ошибаюсь, мне нужно заменить модель Lmer на REML = FALSE. Это дает мне подгонку ML, вероятность которой (и, следовательно, AI C) сравнима с вероятностью lm.
Однако, когда я пытаюсь подобрать эту модель с ML, результат является единственным. Я понимаю, что это не странно, учитывая высокую параметризацию моей модели, но я не знаю, как действовать в этом случае.
Ниже я приведу воспроизводимый пример того, что я сделал. Любая подсказка или направление для поиска будут полезны.
Большое спасибо.
Asier
y<- c(0.26, 0.0972, 0.0844444444444444, 7.915, 2.01857142857143, 4.36166666666667, 2.111,
0.5725, 0.1615, 5.0015, 1.01833333333333, 0.114, 0.054, 0.59, 0.958, 2.01, 0.102,
3.26844444444444, 1.295, 3.84666666666667, 0.126428571428571, 0.525, 3.627875,
2.63285714285714, 0.0516666666666667, 4.35444444444444, 3.562, 4.396, 1.07,
2.21333333333333, 0.13, 0.14575, 3.4425, 6.588, 2.1545, 1.49833333333333, 3.717,
3.59833333333333, 0.047, 1.46366666666667, 0.981428571428571, 4.188, 3.33428571428571,
0.0892727272727273, 1.997, 0, 0.977777777777778, 2.779, 0.816636363636364, 1.5575,
0.10125, 0.854, 1.34025, 0.056, 3.20714285714286, 8.078, 0.375555555555556,
0.0813333333333333, 0, 2.506, 4.08, 0.053, 0.073, 0, 1.11233333333333, 2.39433333333333,
0.30625, 5.23025, 0.21, 0.806166666666667, 5.342, 5.673375, 0.329, 1.265,
1.36288888888889, 5.2264, 0.127, 1.01, 1.8325, 1.924, 3.3206, 3.2375, 3.6725,
0.193111111111111, 2.893, 1.71, 3.47277777777778, 4.08575, 7.04818181818182, 0.686,
0.124, 5.5065, 0.2178, 0, 0, 9.8975, 0.16375, 6.03166666666667, 0, 3.775,
2.30333333333333, 0.717142857142857, 0.23425, 0)
id <- c(1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7, 7, 8, 8, 8, 9, 9, 10, 12, 12, 13, 13, 14,
14, 14, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 19, 20, 21, 22,
22, 23, 23, 23, 24, 24, 24, 24, 25, 25, 26, 27, 27, 27, 28, 28, 28, 28, 29, 29, 30,
31, 31, 31, 32, 32, 32, 32, 33, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 37, 37, 38,
38, 38, 39, 39, 40, 41, 41, 41, 42, 43, 43, 44, 44, 45)
G <- c("Gen2062","Gen2062","Gen1043","Gen1043","Gen2062","Gen2062","Gen2062","Gen2062",
"Gen2062","Gen1002", "Gen1002","Gen1002","Gen1043","Gen1043","Gen2062","Gen2062",
"Gen1002","Gen1002","Gen1002","Gen1043", "Gen1043","Gen1043","Gen1002","Gen1002",
"Gen1043","Gen1043","Gen1002","Gen1002","Gen1002","Gen2062", "Gen2062","Gen2062",
"Gen2062","Gen2062","Gen2062","Gen1002","Gen1002","Gen1002","Gen2062","Gen2062",
"Gen2062","Gen1043","Gen1043","Gen1043","Gen1043","Gen1002","Gen1043","Gen2062",
"Gen2062","Gen1043", "Gen1043","Gen1043","Gen2062","Gen2062","Gen2062","Gen2062",
"Gen2062","Gen2062","Gen1002","Gen1043", "Gen1043","Gen1043","Gen1002","Gen1002",
"Gen1002","Gen1002","Gen1002","Gen1002","Gen1043","Gen1043", "Gen1043","Gen1043",
"Gen2062","Gen2062","Gen2062","Gen2062","Gen1002","Gen1002","Gen1002","Gen1002",
"Gen1002","Gen1002","Gen1002","Gen2062","Gen2062","Gen2062","Gen1002","Gen1043",
"Gen1043","Gen2062", "Gen2062","Gen2062","Gen2062","Gen2062","Gen1002","Gen2062",
"Gen2062","Gen2062","Gen1002","Gen1043", "Gen1043","Gen1043","Gen1043","Gen1043")
treat <- factor(c("CLIP","CLIP","MWh","MWh","CTRL","CTRL","MJ","MJ","MJ","MWl","MWl","MWl",
"CLIP","CLIP", "MWh","MWh","CTRL","CTRL","CTRL","MJ","MJ","MWl","MWh","MWh","CTRL",
"CTRL","MJ","MJ", "MJ","MWl","MWl","MJ","MJ","MJ","MJ","CTRL","CTRL","CTRL","MWl",
"MWl","MWl","MWh", "MWh","MWh","MWh","CLIP","MJ","CTRL","CTRL","MWl","MWl","MWl",
"MWh","MWh","MWh","MWh", "CLIP","CLIP","MJ","CTRL","CTRL","CTRL","MWl","MWl","MWl",
"MWl","MWh","MWh","CLIP","MWh", "MWh","MWh","MWl","MWl","MWl","MWl","MJ","MJ","MJ",
"MJ","CLIP","CLIP","CLIP","CTRL", "CTRL","CTRL","MWh","MWl","MWl","MJ","MJ","MJ",
"CLIP","CLIP","CTRL","MWh","MWh","MWh", "MWl","MJ","MJ","CLIP","CLIP","CTRL"),
levels= c("CTRL","MJ","CLIP","MWl","MWh"))
year <- factor(c(3, 6, 3, 4, 3, 4, 4, 3, 6, 5, 3, 6, 6, 3, 5, 3, 6, 5, 3, 3, 6, 3, 5, 3, 6, 3,
5, 6, 3, 3, 6, 6, 3, 4, 5, 3, 4, 5, 6, 5, 3, 4, 3, 6, 5, 3, 3, 5, 6, 5, 6, 3,
5, 6, 3, 4, 3, 6, 3, 3, 5, 6, 6, 3, 5, 4, 3, 5, 3, 6, 5, 4, 6, 3, 5, 4, 6, 3,
5, 4, 5, 3, 4, 6, 4, 3, 5, 5, 4, 3, 6, 4, 6, 3, 3, 3, 6, 4, 3, 5, 3, 3, 6, 3))
yearsc <- c(-1.07978316353608, 1.3728671650673, -1.07978316353608, -0.262233054001619,
-1.07978316353608, -0.262233054001619, -0.262233054001619, -1.07978316353608,
1.3728671650673, 0.555317055532841, -1.07978316353608, 1.3728671650673,
1.3728671650673, -1.07978316353608, 0.555317055532841, -1.07978316353608,
1.3728671650673, 0.555317055532841, -1.07978316353608, -1.07978316353608,
1.3728671650673, -1.07978316353608, 0.555317055532841, -1.07978316353608,
1.3728671650673, -1.07978316353608, 0.555317055532841, 1.3728671650673,
-1.07978316353608, -1.07978316353608, 1.3728671650673, 1.3728671650673,
-1.07978316353608, -0.262233054001619, 0.555317055532841, -1.07978316353608,
-0.262233054001619, 0.555317055532841, 1.3728671650673, 0.555317055532841,
-1.07978316353608, -0.262233054001619, -1.07978316353608, 1.3728671650673,
0.555317055532841, -1.07978316353608, -1.07978316353608, 0.555317055532841,
1.3728671650673, 0.555317055532841, 1.3728671650673, -1.07978316353608,
0.555317055532841, 1.3728671650673, -1.07978316353608, -0.262233054001619,
-1.07978316353608, 1.3728671650673, -1.07978316353608, -1.07978316353608,
0.555317055532841, 1.3728671650673, 1.3728671650673, -1.07978316353608,
0.555317055532841, -0.262233054001619, -1.07978316353608, 0.555317055532841,
-1.07978316353608, 1.3728671650673, 0.555317055532841, -0.262233054001619,
1.3728671650673, -1.07978316353608, 0.555317055532841, -0.262233054001619,
1.3728671650673, -1.07978316353608, 0.555317055532841, -0.262233054001619,
0.555317055532841, -1.07978316353608, -0.262233054001619, 1.3728671650673,
-0.262233054001619, -1.07978316353608, 0.555317055532841, 0.555317055532841,
-0.262233054001619, -1.07978316353608, 1.3728671650673, -0.262233054001619,
1.3728671650673, -1.07978316353608, -1.07978316353608, -1.07978316353608,
1.3728671650673, -0.262233054001619, -1.07978316353608, 0.555317055532841,
-1.07978316353608, -1.07978316353608, 1.3728671650673, -1.07978316353608)
size <- c(-0.787432564310817, -0.152782818327027, -0.699224001486484, -0.550015536772591,
-0.781582651675141, -0.610089505479677, -0.577270197592965, -0.741751305963175,
0.356078702136875, -0.13601511806665, -0.638026310946081, 0.471603198642645,
1.15814515868844, -0.789767139940315, -0.0493058764568589, -0.830867368236643,
0.282001152184494, 0.0005151484178005, -0.706132205926861, -0.758959877669081,
0.635067465847522, -0.639996726846006, 0.338851364938649, -0.687759226677031,
1.50611203577326, -0.617838405181706, -0.431216398523762, -0.095943463138333,
-0.828615165337114, -0.894932740870706, -0.190534573435163, 2.09068336019719,
-0.728105080373099, -0.428488978079942, 0.371238758871314, -0.795793311954492,
-0.534287991283384, -0.121219875531701, 1.58721804976744, 0.336882526401706,
-0.704783835032158, -0.288243542229793, -0.599717273194699, 2.3037811783131,
0.628369071795738, -0.931816759347182, -0.68268847216755, -0.395481979078294,
-0.0979933714973827, 1.32533791387265, 3.1571568700452, -0.404143493349221,
0.590877174420973, 1.81784576923909, -0.703980791394687, -0.416364696648416,
-0.768567091923648, 0.23505593448162, -0.939399200506726, -0.617601830894264,
0.63681119086057, 2.55610810427755, 3.8807918928945, -0.469139992701266,
1.20366369228973, -0.000593770934626139, -0.634399579467078, 0.689872113268812,
-0.808520055565538, 2.23682619130318, 0.390822911612295, -0.394403171374082,
2.39599810096244, -0.500738681014898, 0.806525796582534, -0.203580589790564,
1.86137157765943, -0.735217373868222, 0.398014962973706, -0.437491202755023,
-0.517043281980417, -0.819560480524912, -0.793550730064359, 0.0959488057944091,
-0.852469670957697, -0.927580574622887, -0.327342793595416, 0.520510016505482,
-0.182641774798893, -0.653124769092475, 1.61024043157109, -0.420258943460452,
-0.430386880522537, -0.871015659481176, -0.883994250355051, -0.82909061812865,
0.731026629392447, -0.721724460118443, -0.815520300277671, 0.115266586262742,
-0.632005356595319, -0.771733007212575, 0.453258256039623, -0.7880133716792)
my.df2<- data.frame(y, id, G, treat, year, yearsc, size)
contrasts(my.df2$treat) <- contr.treatment
contrasts(my.df2$G) <- contr.sum
contrasts(my.df2$year)<- contr.treatment
library("lme4")
#> Loading required package: Matrix
mod1 <- lmer(y~treat*G + year + year:treat +size +(1|id), my.df2)
summary(mod1)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: y ~ treat * G + year + year:treat + size + (1 | id)
#> Data: my.df2
#>
#> REML criterion at convergence: 330.3
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -1.7475 -0.4820 -0.0221 0.2716 4.1059
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> id (Intercept) 0.06866 0.262
#> Residual 2.43441 1.560
#> Number of obs: 104, groups: id, 44
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 1.32960 0.62538 2.126
#> treat2 -0.09680 0.77170 -0.125
#> treat3 -0.89724 0.80078 -1.120
#> treat4 -0.66361 0.79923 -0.830
#> treat5 1.41161 0.82483 1.711
#> G1 -0.42156 0.51853 -0.813
#> G2 0.58878 0.57904 1.017
#> year2 2.28864 1.13245 2.021
#> year3 2.20964 1.01149 2.185
#> year4 -0.81387 1.06034 -0.768
#> size -0.41995 0.33247 -1.263
#> treat2:G1 0.05079 0.71043 0.071
#> treat3:G1 1.24021 1.00465 1.234
#> treat4:G1 0.03188 0.70573 0.045
#> treat5:G1 0.34787 0.79340 0.438
#> treat2:G2 -0.10445 0.81277 -0.129
#> treat3:G2 -0.82656 0.85744 -0.964
#> treat4:G2 -0.06169 0.75898 -0.081
#> treat5:G2 -0.78486 0.75113 -1.045
#> treat2:year2 0.51063 1.50463 0.339
#> treat3:year2 -0.24544 2.22703 -0.110
#> treat4:year2 1.87065 1.54810 1.208
#> treat5:year2 1.13673 1.47973 0.768
#> treat2:year3 -0.50833 1.35564 -0.375
#> treat3:year3 -0.40223 2.15699 -0.186
#> treat4:year3 -0.16616 1.29332 -0.128
#> treat5:year3 -1.65140 1.28763 -1.283
#> treat2:year4 0.97226 1.23395 0.788
#> treat3:year4 1.06639 1.32601 0.804
#> treat4:year4 1.18690 1.27731 0.929
#> treat5:year4 -0.97460 1.38438 -0.704
#>
#> Correlation matrix not shown by default, as p = 31 > 12.
#> Use print(x, correlation=TRUE) or
#> vcov(x) if you need it
mod2 <- lm(y~treat*G + year + year:treat +size, my.df2)
summary(mod2)
#>
#> Call:
#> lm(formula = y ~ treat * G + year + year:treat + size, data = my.df2)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.7330 -0.7436 -0.0521 0.4474 6.5494
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.30986 0.62317 2.102 0.0390 *
#> treat2 -0.09464 0.76991 -0.123 0.9025
#> treat3 -0.89698 0.79892 -1.123 0.2652
#> treat4 -0.66726 0.79755 -0.837 0.4055
#> treat5 1.40364 0.82312 1.705 0.0924 .
#> G1 -0.42105 0.50724 -0.830 0.4092
#> G2 0.61061 0.56940 1.072 0.2871
#> year2 2.31885 1.14117 2.032 0.0458 *
#> year3 2.24022 1.01549 2.206 0.0305 *
#> year4 -0.76560 1.06268 -0.720 0.4736
#> size -0.44355 0.32974 -1.345 0.1827
#> treat2:G1 0.03985 0.69321 0.057 0.9543
#> treat3:G1 1.23850 0.99694 1.242 0.2181
#> treat4:G1 0.03143 0.68797 0.046 0.9637
#> treat5:G1 0.36039 0.78027 0.462 0.6455
#> treat2:G2 -0.13073 0.80002 -0.163 0.8707
#> treat3:G2 -0.84184 0.84639 -0.995 0.3232
#> treat4:G2 -0.07929 0.74368 -0.107 0.9154
#> treat5:G2 -0.81681 0.73590 -1.110 0.2707
#> treat2:year2 0.46098 1.51702 0.304 0.7621
#> treat3:year2 -0.22867 2.24265 -0.102 0.9191
#> treat4:year2 1.87100 1.55630 1.202 0.2332
#> treat5:year2 1.15153 1.49203 0.772 0.4427
#> treat2:year3 -0.47888 1.36376 -0.351 0.7265
#> treat3:year3 -0.37929 2.17060 -0.175 0.8618
#> treat4:year3 -0.15231 1.30132 -0.117 0.9071
#> treat5:year3 -1.66341 1.29562 -1.284 0.2032
#> treat2:year4 0.97095 1.24366 0.781 0.4375
#> treat3:year4 1.04592 1.33667 0.782 0.4365
#> treat4:year4 1.20306 1.28681 0.935 0.3529
#> treat5:year4 -0.91330 1.39371 -0.655 0.5143
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.578 on 73 degrees of freedom
#> Multiple R-squared: 0.6184, Adjusted R-squared: 0.4616
#> F-statistic: 3.943 on 30 and 73 DF, p-value: 9.092e-07
mod2_ML <- update(mod1, REML=FALSE)
#> boundary (singular) fit: see ?isSingular
summary(mod2_ML)
#> Linear mixed model fit by maximum likelihood ['lmerMod']
#> Formula: y ~ treat * G + year + year:treat + size + (1 | id)
#> Data: my.df2
#>
#> AIC BIC logLik deviance df.resid
#> 419.2 506.5 -176.6 353.2 71
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.0672 -0.5625 -0.0394 0.3384 4.9539
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> id (Intercept) 0.000 0.000
#> Residual 1.748 1.322
#> Number of obs: 104, groups: id, 44
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 1.30986 0.52209 2.509
#> treat2 -0.09464 0.64503 -0.147
#> treat3 -0.89698 0.66934 -1.340
#> treat4 -0.66726 0.66819 -0.999
#> treat5 1.40364 0.68962 2.035
#> G1 -0.42105 0.42497 -0.991
#> G2 0.61061 0.47705 1.280
#> year2 2.31885 0.95608 2.425
#> year3 2.24022 0.85078 2.633
#> year4 -0.76560 0.89033 -0.860
#> size -0.44355 0.27626 -1.606
#> treat2:G1 0.03985 0.58078 0.069
#> treat3:G1 1.23850 0.83524 1.483
#> treat4:G1 0.03143 0.57639 0.055
#> treat5:G1 0.36039 0.65372 0.551
#> treat2:G2 -0.13073 0.67027 -0.195
#> treat3:G2 -0.84184 0.70911 -1.187
#> treat4:G2 -0.07929 0.62306 -0.127
#> treat5:G2 -0.81681 0.61655 -1.325
#> treat2:year2 0.46098 1.27097 0.363
#> treat3:year2 -0.22867 1.87891 -0.122
#> treat4:year2 1.87100 1.30388 1.435
#> treat5:year2 1.15153 1.25003 0.921
#> treat2:year3 -0.47888 1.14257 -0.419
#> treat3:year3 -0.37929 1.81855 -0.209
#> treat4:year3 -0.15231 1.09026 -0.140
#> treat5:year3 -1.66341 1.08548 -1.532
#> treat2:year4 0.97095 1.04195 0.932
#> treat3:year4 1.04592 1.11988 0.934
#> treat4:year4 1.20306 1.07810 1.116
#> treat5:year4 -0.91330 1.16767 -0.782
#>
#> Correlation matrix not shown by default, as p = 31 > 12.
#> Use print(x, correlation=TRUE) or
#> vcov(x) if you need it
#> convergence code: 0
#> boundary (singular) fit: see ?isSingular
detach("package:lme4", unload=TRUE)
Создано в 2020-03-27 пакетом представлены (v0.3.0)
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 3.6.2 (2019-12-12)
#> os Ubuntu 16.04.5 LTS
#> system x86_64, linux-gnu
#> ui X11
#> language eu
#> collate eu_ES.UTF-8
#> ctype eu_ES.UTF-8
#> tz Europe/Madrid
#> date 2020-03-27
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date lib source
#> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.2)
#> backports 1.1.5 2019-10-02 [1] CRAN (R 3.6.2)
#> boot 1.3-23 2019-07-05 [4] CRAN (R 3.6.1)
#> callr 3.4.1 2020-01-24 [1] CRAN (R 3.6.2)
#> cli 2.0.1 2020-01-08 [1] CRAN (R 3.6.2)
#> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.2)
#> desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.2)
#> devtools 2.2.2 2020-02-17 [1] CRAN (R 3.6.2)
#> digest 0.6.23 2019-11-23 [1] CRAN (R 3.6.2)
#> ellipsis 0.3.0 2019-09-20 [1] CRAN (R 3.6.2)
#> evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.2)
#> fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.2)
#> fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.2)
#> glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.2)
#> highr 0.8 2019-03-20 [1] CRAN (R 3.6.2)
#> htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.6.2)
#> knitr 1.28 2020-02-06 [1] CRAN (R 3.6.2)
#> lattice 0.20-40 2020-02-19 [1] CRAN (R 3.6.2)
#> magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.2)
#> MASS 7.3-51.5 2019-12-20 [1] CRAN (R 3.6.2)
#> Matrix * 1.2-18 2019-11-27 [4] CRAN (R 3.6.1)
#> memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.2)
#> minqa 1.2.4 2014-10-09 [1] CRAN (R 3.6.2)
#> nlme 3.1-144 2020-02-06 [4] CRAN (R 3.6.2)
#> nloptr 1.2.2.1 2020-03-11 [1] CRAN (R 3.6.2)
#> pkgbuild 1.0.6 2019-10-09 [1] CRAN (R 3.6.2)
#> pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.6.2)
#> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 3.6.2)
#> processx 3.4.2 2020-02-09 [1] CRAN (R 3.6.2)
#> ps 1.3.0 2018-12-21 [1] CRAN (R 3.6.2)
#> R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.2)
#> Rcpp 1.0.3 2019-11-08 [1] CRAN (R 3.6.2)
#> remotes 2.1.1 2020-02-15 [1] CRAN (R 3.6.2)
#> rlang 0.4.4 2020-01-28 [1] CRAN (R 3.6.2)
#> rmarkdown 2.1 2020-01-20 [1] CRAN (R 3.6.2)
#> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.2)
#> rstudioapi 0.11 2020-02-07 [1] CRAN (R 3.6.2)
#> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.2)
#> stringi 1.4.5 2020-01-11 [1] CRAN (R 3.6.2)
#> stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.2)
#> testthat 2.3.2 2020-03-02 [1] CRAN (R 3.6.2)
#> usethis 1.5.1 2019-07-04 [1] CRAN (R 3.6.2)
#> withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.2)
#> xfun 0.12 2020-01-13 [1] CRAN (R 3.6.2)
#> yaml 2.2.1 2020-02-01 [1] CRAN (R 3.6.2)
#>
#> [1] /home/asier/R/x86_64-pc-linux-gnu-library/3.6
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/library