GLMER警告：方差 - 协方差矩阵[...]是不是正定或包含NA值

我有时会发现从glmer，包lme4我GLMMs参数，显示下面的警告消息，当他们的总结叫做：GLMER警告：方差 - 协方差矩阵[...]是不是正定或包含NA值

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.nb估算的theta调用glmer确实给了我警告。后者是在相应的glmer.nb完整型号上使用dredge（MuMIn）生成的型号调用。

来源

2016-08-17 kdarras

能否请您包括数据和/或代码，将为我们提供一个[重复的例子（HTTP：/ /stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example）？这表明您的标准误差估计值*可能*不准确，但没有示例就很难知道。 –

我现在做了一个可重现的例子。我的最小数据集非常大，所以我想问一下警告的一般含义和令人担忧的问题。 – kdarras

此警告表明您的标准错误估计值可能是不准确。但正如所有的警告一样，很难确切地知道，最好的办法是尽可能地进行交叉检查。

在这种情况下，我从glmer.nb和glmer保存了两个吻合，分别为g1和g2。你可以看到估计值（点估计，SEs，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 **

我有lme4在Github上（在test_mods分支，有望集成到主分支很快：如果你要安装它，您可以使用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"

的glmmTMB包（在Github上）也给出了几乎相同的结果：

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 **

来源

2016-08-17 21:58:04

GLMER警告：方差 - 协方差矩阵[...]是不是正定或包含NA值

回答

相关问题