有时我发现glmer包lme4中的GLMM在调用其摘要时显示以下警告消息:
Warning messages:
1: In vcov.merMod(object, use.hessian = use.hessian) :
variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RX
2: In vcov.merMod(object, correlation = correlation, sigm = sig) :
variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RX
我在Stackoverflow上发现的类似问题是指其他功能，而不是glmer，而且LME4 Wiki也没有对此进行详细说明。在this问题中，问题是在解决此类错误消息之前解决的，而here的讨论重点是特定模型而不是警告消息的含义。
所以问题是:我应该担心该消息，还是可以，因为它只是警告而不是错误，所以就可以了，正如它所说的，它“回落到RX估计的var-cov”(无论RX是什么)反正。
有趣的是，尽管摘要指出模型无法收敛，但是我没有得到通常的红色收敛警告。
这是一个(最小的)数据集:
testdata=structure(list(Site = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("EO1", "EO2",
"EO3", "EO4", "EO5", "EO6"), class = "factor"), Treatment = structure(c(1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L), .Label = c("control",
"no ants", "no birds", "no birds no ants"), class = "factor"),
Tree = structure(c(2L, 3L, 4L, 16L, 12L, 13L, 14L, 15L, 5L,
6L, 7L, 8L, 1L, 9L, 10L, 11L, 28L, 29L, 30L, 31L, 17L, 25L,
26L, 27L, 18L, 19L, 20L, 32L, 21L, 22L, 23L, 24L, 33L, 41L,
42L, 43L, 37L, 38L, 39L, 40L, 44L, 45L, 46L, 47L, 34L, 35L,
36L, 48L, 49L, 57L, 58L, 59L, 50L, 51L, 52L, 64L, 53L, 54L,
55L, 56L, 60L, 61L, 62L, 63L, 66L, 67L, 68L, 80L, 69L, 70L,
71L, 72L, 76L, 77L, 78L, 79L, 65L, 73L, 74L, 75L, 82L, 83L,
84L, 96L, 92L, 93L, 94L, 95L, 85L, 86L, 87L, 88L, 81L, 89L,
90L, 91L), .Label = c("EO1 1", "EO1 10", "EO1 11", "EO1 12",
"EO1 13", "EO1 14", "EO1 15", "EO1 16", "EO1 2", "EO1 3",
"EO1 4", "EO1 5", "EO1 6", "EO1 7", "EO1 8", "EO1 9", "EO2 1",
"EO2 10", "EO2 11", "EO2 12", "EO2 13", "EO2 14", "EO2 15",
"EO2 16", "EO2 2", "EO2 3", "EO2 4", "EO2 5", "EO2 6", "EO2 7",
"EO2 8", "EO2 9", "EO3 1", "EO3 10", "EO3 11", "EO3 12",
"EO3 13", "EO3 14", "EO3 15", "EO3 16", "EO3 2", "EO3 3",
"EO3 4", "EO3 5", "EO3 6", "EO3 7", "EO3 8", "EO3 9", "EO4 1",
"EO4 10", "EO4 11", "EO4 12", "EO4 13", "EO4 14", "EO4 15",
"EO4 16", "EO4 2", "EO4 3", "EO4 4", "EO4 5", "EO4 6", "EO4 7",
"EO4 8", "EO4 9", "EO5 1", "EO5 10", "EO5 11", "EO5 12",
"EO5 13", "EO5 14", "EO5 15", "EO5 16", "EO5 2", "EO5 3",
"EO5 4", "EO5 5", "EO5 6", "EO5 7", "EO5 8", "EO5 9", "EO6 1",
"EO6 10", "EO6 11", "EO6 12", "EO6 13", "EO6 14", "EO6 15",
"EO6 16", "EO6 2", "EO6 3", "EO6 4", "EO6 5", "EO6 6", "EO6 7",
"EO6 8", "EO6 9"), class = "factor"), predators_trunk = c(7L,
10L, 9L, 15L, 18L, 11L, 5L, 7L, 15L, 12L, 6L, 12L, 7L, 13L,
24L, 17L, 3L, 0L, 0L, 2L, 4L, 3L, 0L, 6L, 2L, 3L, 5L, 1L,
5L, 12L, 18L, 15L, 7L, 0L, 5L, 1L, 17L, 7L, 13L, 19L, 7L,
3L, 5L, 10L, 11L, 7L, 13L, 7L, 7L, 0L, 4L, 2L, 5L, 7L, 4L,
7L, 8L, 7L, 9L, 20L, 13L, 2L, 12L, 7L, 0L, 7L, 2L, 2L, 2L,
4L, 17L, 2L, 3L, 1L, 1L, 1L, 11L, 1L, 1L, 8L, 8L, 18L, 5L,
6L, 6L, 5L, 6L, 5L, 9L, 2L, 8L, 13L, 13L, 5L, 3L, 5L), pH_H2O = c(4.145,
4.145, 4.145, 4.145, 4.1825, 4.1825, 4.1825, 4.1825, 4.1325,
4.1325, 4.1325, 4.1325, 4.14125, 4.14125, 4.14125, 4.14125,
4.265, 4.265, 4.265, 4.265, 4.21, 4.21, 4.21, 4.21, 4.18375,
4.18375, 4.18375, 4.18375, 4.09625, 4.09625, 4.09625, 4.09625,
4.1575, 4.1575, 4.1575, 4.1575, 4.1125, 4.1125, 4.1125, 4.1125,
4.20875, 4.20875, 4.20875, 4.20875, 3.97125, 3.97125, 3.97125,
3.97125, 4.025, 4.025, 4.025, 4.025, 4.005, 4.005, 4.005,
4.005, 4.04, 4.04, 4.04, 4.04, 4.03125, 4.03125, 4.03125,
4.03125, 4.4575, 4.4575, 4.4575, 4.4575, 4.52, 4.52, 4.52,
4.52, 4.505, 4.505, 4.505, 4.505, 4.34875, 4.34875, 4.34875,
4.34875, 4.305, 4.305, 4.305, 4.305, 4.32, 4.32, 4.32, 4.32,
4.35, 4.35, 4.35, 4.35, 4.445, 4.445, 4.445, 4.445), ant_mean_abundance = c(53.85714,
53.85714, 53.85714, 53.85714, 24.28571, 24.28571, 24.28571,
24.28571, 45.5, 45.5, 45.5, 45.5, 51.14286, 51.14286, 51.14286,
51.14286, 66.28571, 66.28571, 66.28571, 66.28571, 76.5, 76.5,
76.5, 76.5, 65.71429, 65.71429, 65.71429, 65.71429, 8.642857,
8.642857, 8.642857, 8.642857, 109.3571, 109.3571, 109.3571,
109.3571, 25.14286, 25.14286, 25.14286, 25.14286, 101.3571,
101.3571, 101.3571, 101.3571, 31.78571, 31.78571, 31.78571,
31.78571, 78.64286, 78.64286, 78.64286, 78.64286, 93.28571,
93.28571, 93.28571, 93.28571, 63.14286, 63.14286, 63.14286,
63.14286, 67.14286, 67.14286, 67.14286, 67.14286, 44.0625,
44.0625, 44.0625, 44.0625, 23.875, 23.875, 23.875, 23.875,
95.8125, 95.8125, 95.8125, 95.8125, 49.125, 49.125, 49.125,
49.125, 57, 57, 57, 57, 38.125, 38.125, 38.125, 38.125, 40.6875,
40.6875, 40.6875, 40.6875, 22, 22, 22, 22), bird_activity = c(153.24,
153.24, 153.24, 153.24, 153.24, 153.24, 153.24, 153.24, 0,
0, 0, 0, 0, 0, 0, 0, 240.96, 240.96, 240.96, 240.96, 240.96,
240.96, 240.96, 240.96, 0, 0, 0, 0, 0, 0, 0, 0, 154.54, 154.54,
154.54, 154.54, 154.54, 154.54, 154.54, 154.54, 0, 0, 0,
0, 0, 0, 0, 0, 107.68, 107.68, 107.68, 107.68, 107.68, 107.68,
107.68, 107.68, 0, 0, 0, 0, 0, 0, 0, 0, 172.42, 172.42, 172.42,
172.42, 172.42, 172.42, 172.42, 172.42, 0, 0, 0, 0, 0, 0,
0, 0, 113.8, 113.8, 113.8, 113.8, 113.8, 113.8, 113.8, 113.8,
0, 0, 0, 0, 0, 0, 0, 0)), .Names = c("Site", "Treatment",
"Tree", "predators_trunk", "pH_H2O", "ant_mean_abundance", "bird_activity"
), class = "data.frame", row.names = c(NA, -96L))
这是导致警告的代码:
library(lme4)
summary(glmer.nb(predators_trunk ~ scale(ant_mean_abundance) + scale(bird_activity) + scale(pH_H2O) + (1 | Site/Treatment), testdata, na.action = na.fail))
summary(glmer(predators_trunk ~ scale(ant_mean_abundance) + scale(bird_activity) + scale(pH_H2O) + (1 | Site/Treatment), testdata, family = negative.binomial(theta = 4.06643400243645), na.action = na.fail))
对我来说有趣的是，
glmer.nb的摘要未产生任何警告，但是使用glmer估计的theta调用glmer.nb的确给了我警告。后者是通过在相应的dredge完整模型上使用glmer.nb(MuMIn)生成的模型调用。 最佳答案
此警告表明您的标准错误估算值可能不太准确。但是与所有警告一样，很难确定，最好的办法是尝试进行交叉检查。
在这种情况下，我将glmer.nb和glmer中的两个拟合保存为g1和g2。您可以看到估算值(点估算值，SE，Z值...)已经有所变化，但变化不大，因此至少可以使您放心。
printCoefmat(coef(summary(g1)),digits=2)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.844 0.111 16.7 <2e-16 ***
scale(ant_mean_abundance) -0.347 0.077 -4.5 7e-06 ***
scale(bird_activity) -0.122 0.076 -1.6 0.107
scale(pH_H2O) -0.275 0.104 -2.6 0.008 **
> printCoefmat(coef(summary(g2)),digits=2)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.846 0.108 17.1 <2e-16 ***
scale(ant_mean_abundance) -0.347 0.077 -4.5 6e-06 ***
scale(bird_activity) -0.122 0.075 -1.6 0.102
scale(pH_H2O) -0.275 0.102 -2.7 0.007 **
我在Github上有一个
lme4的开发版本(test_mods分支，希望很快会集成到master分支中:如果您要安装它，则可以使用devtools::install_github("lme4/lme4",ref="test_mods"))，它使您可以为标准选择更准确(但更慢)的计算错误:这使我们回到(几乎)与glmer.nb相同的标准错误。g3 <- update(g2, control=glmerControl(deriv.method="Richardson"))
printCoefmat(coef(summary(g3)),digits=2)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.846 0.111 16.7 <2e-16 ***
scale(ant_mean_abundance) -0.347 0.077 -4.5 6e-06 ***
scale(bird_activity) -0.122 0.076 -1.6 0.106
scale(pH_H2O) -0.275 0.104 -2.6 0.008 **
all.equal(coef(summary(g1))[,"Std. Error"],
coef(summary(g3))[,"Std. Error"])
[1] "Mean relative difference: 0.001597978"
(在Github上)
glmmTMB包也给出了几乎相同的结果:library(glmmTMB)
g5 <- glmmTMB(predators_trunk ~ scale(ant_mean_abundance) +
scale(bird_activity) + scale(pH_H2O) +
(1 | Site/Treatment), testdata,
family=nbinom2)
printCoefmat(coef(summary(g5))[["cond"]],digits=2)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.852 0.110 16.8 <2e-16 ***
scale(ant_mean_abundance) -0.348 0.077 -4.5 7e-06 ***
scale(bird_activity) -0.123 0.076 -1.6 0.106
scale(pH_H2O) -0.276 0.105 -2.6 0.008 **
关于r - GLMER警告: variance-covariance matrix [. ..]不是肯定的或包含NA值，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/38997371/