机器算法验证 - “错误：未找到有效的系数集：请提供起始值”在尝试获取 R 中的置信区间时 - 吾爱随笔录

“错误：未找到有效的系数集：请提供起始值”在尝试获取 R 中的置信区间时

机器算法验证 r 置信区间泊松回归

2022-04-19 09:16:14

我正在尝试分析一些物种的一些计数数据（“Tetab”表示以下代码中的物种）。咨询了一位比我更懂统计的朋友，他建议用泊松回归分析数据，然后利用置信区间来确定哪些处理导致计数反应显着不同。这对于我分析的其他两个物种来说效果很好，但我得到了标题中列出的错误。比较代码，不同物种的分析中的一切都是相同的，所以我假设它与数据有关——也因为这个物种是唯一不能运行零膨胀泊松回归的物种。其他两个物种的总计数数据是 33 和 47，但 Tetab 只有 22。这可能与错误有关吗？有什么解决方法吗？数据的方差是异质的，所以我不能使用 Kruskal-Wallis 或多重比较。

> Tetab.pglm <- glm(Count ~ Treatment, data = spond.spp.list[['Tetab']], family = poisson)
> Tetab.zpglm <- zeroinfl(Count ~ Treatment, data = spond.spp.list[['Tetab']], dist = "poisson")
Error in solve.default(as.matrix(fit$hessian)) : 
  system is computationally singular: reciprocal condition number = 1.63511e-19

> summary(Tetab.pglm)

Call:
glm(formula = Count ~ Treatment, family = poisson, data = spond.spp.list[["Tetab"]])

Deviance Residuals: 
 Min        1Q    Median        3Q       Max  
-1.35873  -0.00006  -0.00006   0.07899   2.36154  

Coefficients:
          Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.030e+01  4.311e+03  -0.005    0.996
Treatment2   2.004e+01  4.311e+03   0.005    0.996
Treatment3   9.922e-09  6.096e+03   0.000    1.000
Treatment4   2.022e+01  4.311e+03   0.005    0.996

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 54.484  on 51  degrees of freedom
Residual deviance: 23.804  on 48  degrees of freedom
AIC: 68.297

Number of Fisher Scoring iterations: 18

> exp(coef(Tetab.zpglm))
Error in coef(Tetab.zpglm) : object 'Tetab.zpglm' not found
> exp(coef(Tetab.pglm))
 (Intercept)   Treatment2   Treatment3   Treatment4 
1.522998e-09 5.050767e+08 1.000000e+00 6.060920e+08 
> exp(confint(Tetab.zpglm))
Error in confint(Tetab.zpglm) : object 'Tetab.zpglm' not found
> exp(confint(Tetab.pglm))
Waiting for profiling to be done...
Error: no valid set of coefficients has been found: please supply starting values
In addition: Warning messages:
1: glm.fit: fitted rates numerically 0 occurred 
2: glm.fit: fitted rates numerically 0 occurred 
3: glm.fit: fitted rates numerically 0 occurred 
4: glm.fit: fitted rates numerically 0 occurred 
5: glm.fit: fitted rates numerically 0 occurred 
6: glm.fit: fitted rates numerically 0 occurred 
7: glm.fit: fitted rates numerically 0 occurred 
8: glm.fit: fitted rates numerically 0 occurred 
9: glm.fit: fitted rates numerically 0 occurred

感谢您的任何帮助，您可以提供！最大限度

这是数据集：

Species Date    Site    Treatment   Count
Tetab   20160602    2   1   0
Tetab   20160602    2   2   1
Tetab   20160602    2   3   0
Tetab   20160602    2   4   1
Tetab   20160606    1   1   0
Tetab   20160606    1   2   1
Tetab   20160606    1   3   0
Tetab   20160606    1   4   0
Tetab   20160606    2   1   0
Tetab   20160606    2   2   1
Tetab   20160606    2   3   0
Tetab   20160606    2   4   0
Tetab   20160607    2   1   0
Tetab   20160607    2   2   0
Tetab   20160607    2   3   0
Tetab   20160607    2   4   1
Tetab   20160609    1   1   0
Tetab   20160609    1   2   0
Tetab   20160609    1   3   0
Tetab   20160609    1   4   2
Tetab   20160609    2   1   0
Tetab   20160609    2   2   0
Tetab   20160609    2   3   0
Tetab   20160609    2   4   1
Tetab   20160610    1   1   0
Tetab   20160610    1   2   1
Tetab   20160610    1   3   0
Tetab   20160610    1   4   0
Tetab   20160610    2   1   0
Tetab   20160610    2   2   1
Tetab   20160610    2   3   0
Tetab   20160610    2   4   0
Tetab   20160620    1   1   0
Tetab   20160620    1   2   1
Tetab   20160620    1   3   0
Tetab   20160620    1   4   1
Tetab   20160620    2   1   0
Tetab   20160620    2   2   1
Tetab   20160620    2   3   0
Tetab   20160620    2   4   4
Tetab   20160622    1   1   0
Tetab   20160622    1   2   0
Tetab   20160622    1   3   0
Tetab   20160622    1   4   1
Tetab   20160622    2   1   0
Tetab   20160622    2   2   2
Tetab   20160622    2   3   0
Tetab   20160622    2   4   1
Tetab   20160624    2   1   0
Tetab   20160624    2   2   1
Tetab   20160624    2   3   0
Tetab   20160624    2   4   0

1个回答

我（还）不能解决这个问题，但我可以诊断出一些问题。

这是您的数据的摘要：

 with(dd,table(Treatment,Count))
         Count
Treatment  0  1  2  4
        1 13  0  0  0
        2  4  8  1  0
        3 13  0  0  0
        4  5  6  1  1

您可以看到处理 1 和 3 的所有值都为零。当我们对此进行泊松 GLM 拟合时，我们将参数拟合到对数尺度上- 也就是说，截距是第一个处理的对数密度，其他参数是其他处理的对数密度与第一的。如果我们看一下系数表：

printCoefmat(coef(summary(Tetab.pglm)))
               Estimate  Std. Error z value Pr(>|z|)
(Intercept) -2.0303e+01  4.3105e+03 -0.0047   0.9962
Treatment2   2.0040e+01  4.3105e+03  0.0046   0.9963
Treatment3   9.9225e-09  6.0960e+03  0.0000   1.0000
Treatment4   2.0223e+01  4.3105e+03  0.0047   0.9963

我们看到所有参数都很大（+20/-20），除了处理 3 基本上为零；标准误差很大；并且p值基本为1。荒谬的标准误现象称为Hauck-Donner效应，它发生在这种极端情况下。

零通胀的东西在这里似乎完全没有必要，这将使已经很困难的情况变得更加困难。

在这种情况下，处理 1 和 3 的对数密度的最大似然估计实际上是 $-\infty$ ，这将使生活变得更加艰难。原则上可以计算治疗 1 的有限置信区间上限和治疗 1 和 (2,4) 之间差异的下置信区间，但它在数值上会很难看。

可能（？）最好的解决方案是某种减少偏差或惩罚的估计，它将解决方案推离 $-\infty$ . 我以为arm::bayesglm()会这样做，但我仍然遇到了麻烦：

b1 <- bayesglm(Count ~ Treatment, data= dd, family=poisson)

给出了合理的答案，但confint(b1)仍然失败......

其它你可能感兴趣的问题

上一篇如何在知道样本均值和某个分位数的情况下拟合对数正态分布的参数？下一篇numer.ai 如何预测未来？