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

Question

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

9

有时候我发现使用lme4包中的glmer函数进行建模后，当调用其summary时，会出现以下警告信息：

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）生成的模型调用。

- kdarras

请问您能否提供数据和/或代码，以便为我们提供一个可重现的示例？这表明您的标准误差估计可能不太准确，但没有示例很难确定。 - Ben Bolker

我现在做了一个可重复的例子。我的最小数据集相当大，所以我更想问一下警告的一般含义和令人担忧程度。 - kdarras

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

这个警告表明您的标准误差估计可能不太准确。但像所有警告一样，很难确定最好的方法是尽可能进行交叉检查。

在这种情况下，我将你的两个拟合结果，从glmer.nb和glmer，保存为g1和g2。您可以看到估计值（点估计、标准误差、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分支上），希望很快就能集成到主分支中：如果你想要安装它，可以使用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 **