逻辑回归返回错误，但在减少的数据集上运行正常答案

【问题标题】：Logistic regression returns error but runs okay on reduced dataset逻辑回归返回错误，但在减少的数据集上运行正常
【发布时间】：2016-09-17 01:35:57
【问题描述】：

非常感谢您对此的意见！

我正在研究逻辑回归，但由于某种原因它不起作用：

mod1<-glm(survive~reLDM2+yr+yr2+reLDM2:yr +reLDM2:yr2+NestAge0,
         family=binomial(link=logexp(NSSH1$exposure)),
                       data=NSSH1, control = list(maxit = 50))

当我用更少的数据运行相同的模型时，它可以工作！但是对于完整的数据集，我会收到错误和警告消息：

Error: inner loop 1; cannot correct step size
In addition: Warning messages:
1: step size truncated due to divergence 
2: step size truncated due to divergence

这是数据：https://www.dropbox.com/s/8ib8m1fh176556h/NSSH1.csv?dl=0

来自User-defined link function for glmer for known-fate survival modeling的日志曝光链接功能：

library(MASS)
logexp <- function(exposure = 1) {
    linkfun <- function(mu) qlogis(mu^(1/exposure))
    ## FIXME: is there some trick we can play here to allow
    ##   evaluation in the context of the 'data' argument?
    linkinv <- function(eta)  plogis(eta)^exposure
    mu.eta <- function(eta) exposure * plogis(eta)^(exposure-1) *
      .Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats")
    valideta <- function(eta) TRUE
    link <- paste("logexp(", deparse(substitute(exposure)), ")",
               sep="")
    structure(list(linkfun = linkfun, linkinv = linkinv,
               mu.eta = mu.eta, valideta = valideta, 
               name = link),
          class = "link-glm")
}

【问题讨论】：

尽可能内嵌代码会更好。另请注意，链接问题具有此类更新的链接功能...

标签： r statistics regression glm logistic-regression

【解决方案1】：

tl;dr 你遇到麻烦了，因为你的 yr 和 yr2 预测器（大概是年和年平方）与不寻常的链接函数结合会导致数值问题；你可以使用glm2 package 来解决这个问题，但我至少会考虑一下在这种情况下尝试拟合平方年份是否有意义。

更新：下面开始使用mle2 的蛮力方法；还没有写它来做完整的交互模型。

Andrew Gelman 的folk theorem 可能适用于这里：

当您遇到计算问题时，您的模型通常会出现问题。

我首先尝试了您的模型的简化版本，没有交互...

NSSH1 <- read.csv("NSSH1.csv")
source("logexpfun.R")  ## for logexp link
mod1 <- glm(survive~reLDM2+yr+yr2+NestAge0,
          family=binomial(link=logexp(NSSH1$exposure)),
          data=NSSH1, control = list(maxit = 50))

... 效果很好。现在让我们试着看看问题出在哪里：

mod2 <- update(mod1,.~.+reLDM2:yr)  ## OK
mod3 <- update(mod1,.~.+reLDM2:yr2) ## OK
mod4 <- update(mod1,.~.+reLDM2:yr2+reLDM2:yr)  ## bad

好的，所以我们无法同时包含两个交互。这些预测变量实际上是如何相互关联的？让我们看看...

pairs(NSSH1[,c("reLDM2","yr","yr2")],gap=0)

yr 和 yr2 并非完全相关，但它们完全是排名相关的，它们在数字上相互干扰当然不足为奇 s> 更新：当然“年”和“年平方”看起来像这样！即使使用构造正交多项式的poly(yr,2)，在这种情况下也无济于事......仍然，值得查看参数以防提供线索......

如上所述，我们可以尝试glm2（使用更强大的算法替代glm）看看会发生什么...

library(glm2)
mod5 <- glm2(survive~reLDM2+yr+yr2+reLDM2:yr +reLDM2:yr2+NestAge0,
          family=binomial(link=logexp(NSSH1$exposure)),
          data=NSSH1, control = list(maxit = 50))

现在我们确实得到了答案。如果我们检查cov2cor(vcov(mod5))，我们会看到yr 和yr2 参数（以及它们与reLDM2 交互的参数是强烈负相关的（大约-0.97）。让我们想象一下......

library(corrplot)
corrplot(cov2cor(vcov(mod5)),method="ellipse")

如果我们尝试通过蛮力来做到这一点呢？

library(bbmle)
link <- logexp(NSSH1$exposure)
fit <- mle2(survive~dbinom(prob=link$linkinv(eta),size=1),
     parameters=list(eta~reLDM2+yr+yr2+NestAge0),
     start=list(eta=0),
     data=NSSH1,
     method="Nelder-Mead")  ## more robust than default BFGS
summary(fit)
##                   Estimate Std. Error  z value   Pr(z)    
## eta.(Intercept)  4.3627816  0.0402640 108.3545 < 2e-16 ***
## eta.reLDM2      -0.0019682  0.0011738  -1.6767 0.09359 .  
## eta.yr          -6.0852108  0.2068159 -29.4233 < 2e-16 ***
## eta.yr2          5.7332780  0.1950289  29.3971 < 2e-16 ***
## eta.NestAge0     0.0612248  0.0051272  11.9411 < 2e-16 ***

这似乎是合理的（您应该检查预测值并查看它们是否有意义......），但是......

cc <- confint(fit)  ## "profiling has found a better solution"

这会返回一个mle2 对象，但它的调用槽被破坏了，所以打印结果很难看。

coef(cc)
## eta.(Intercept)                      eta.reLDM2 
##     4.329824508                    -0.011996582 
##       eta.yr                         eta.yr2 
##     0.101221970                     0.093377127 
##     eta.NestAge0 
##      0.003460453 
##
vcov(cc) ## all NA values! problem?

尝试从这些返回值重新开始...

fit2 <- update(fit,start=list(eta=unname(coef(cc))))
coef(summary(fit2))
##                     Estimate  Std. Error    z value        Pr(z)
## eta.(Intercept)  4.452345889 0.033864818 131.474082 0.000000e+00
## eta.reLDM2      -0.013246977 0.001076194 -12.309102 8.091828e-35
## eta.yr           0.103013607 0.094643420   1.088439 2.764013e-01
## eta.yr2          0.109709373 0.098109924   1.118229 2.634692e-01
## eta.NestAge0    -0.006428657 0.004519983  -1.422274 1.549466e-01

现在我们可以得到置信区间...

ci2 <- confint(fit2)
##                       2.5 %       97.5 %
## eta.(Intercept)  4.38644052  4.519116156
## eta.reLDM2      -0.01531437 -0.011092655
## eta.yr          -0.08477933  0.286279919
## eta.yr2         -0.08041548  0.304251382
## eta.NestAge0    -0.01522353  0.002496006

这似乎可行，但我会非常怀疑这些配合。您可能应该尝试其他优化器以确保您可以得到相同的结果。一些更好的优化工具，例如 AD Model Builder 或 Template Model Builder 可能是个好主意。

我不赞成无意识地放弃具有强相关参数估计的预测变量，但这可能是考虑它的合理时机。

【讨论】：

我在调用 glm 时遇到错误：Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, : NAs in d(mu)/d(eta) 在 Mac El Cap 上使用 R 3.3.0。
第一个glm 电话？我可以相信这一切都足够敏感，以至于它会在不同平台上的不同点出现问题。如果您从公式中删除yr2，glm 是否有效？ glm2 适合你吗？
glm2 在有或没有 yr2 的情况下也会失败。同样的错误信息。
嗯。在那种情况下，如果我需要这样做，我可能会使用蛮力 MLE。也许我会等着听 OP 的回复是否对他们有用。 cloglog 链接（与 power-logistic 非常相似）也可能会更好地工作。
感谢您的回答 @Ben Bolker 和 42，是的，我收到一条关于新链接功能和 glm2 包的错误消息。使用yr 和yr2 作为连续变量（即1,2,3 和1^2, 2^2, 3^2 之一）会得到相同的错误（在我与您分享的数据中，它是用函数尺度标准化的，它们都不起作用）。分析的目的是寻找季节随时间的变化（由于捕食，这是非线性的，因此将年和相对产卵日期（reLDM）交互作为平方项包括在内）。你能帮我实现这个MLE 方法吗？