逻辑回归 returns 错误，但在减少的数据集上运行正常
Logistic regression returns error but runs okay on reduced dataset
非常感谢您的意见！
我正在研究逻辑回归，但由于某种原因它不起作用：
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
日志曝光 link 函数来自 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")
}
tl;dr 你遇到麻烦了，因为你的 yr 和 yr2 预测变量（大概是年和年的平方）与不寻常的 link 函数会导致数值问题；您可以使用 glm2 package 解决这个问题，但我至少会考虑一下在这种情况下尝试拟合平方年是否有意义。
更新：下面开始使用 mle2 的暴力方法；还没有编写它来完成带有交互的完整模型。
Andrew Gelman 的 folk theorem 可能适用于此：
When you have computational problems, often there’s a problem with your model.
我开始尝试您模型的简化版本，没有交互...
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 并非 完全 相关，但它们是完全等级相关的，这当然不足为奇他们在数值上相互干扰 更新：当然 "year" 和 "year-squared" 看起来像这样！即使使用构造正交多项式的 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"
这个 returns 是一个 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 可能是个好主意。
我不赞成盲目放弃具有强相关参数估计值的预测变量，但这可能是考虑它的合理时机。
非常感谢您的意见！
我正在研究逻辑回归，但由于某种原因它不起作用：
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
日志曝光 link 函数来自 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")
}
tl;dr 你遇到麻烦了，因为你的 yr 和 yr2 预测变量（大概是年和年的平方）与不寻常的 link 函数会导致数值问题；您可以使用 glm2 package 解决这个问题，但我至少会考虑一下在这种情况下尝试拟合平方年是否有意义。
更新：下面开始使用 mle2 的暴力方法；还没有编写它来完成带有交互的完整模型。
Andrew Gelman 的 folk theorem 可能适用于此：
When you have computational problems, often there’s a problem with your model.
我开始尝试您模型的简化版本，没有交互...
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)
更新：当然 "year" 和 "year-squared" 看起来像这样！即使使用构造正交多项式的 yr 和 yr2 并非 完全 相关，但它们是完全等级相关的，这当然不足为奇他们在数值上相互干扰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"
这个 returns 是一个 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 可能是个好主意。
我不赞成盲目放弃具有强相关参数估计值的预测变量，但这可能是考虑它的合理时机。