Logistic回归返回错误但在减少的数据集上运行正常
MSS*_*MSS 3 statistics regression r glm logistic-regression
非常感谢您对此的投入!
我正在进行逻辑回归,但由于某种原因它不起作用:
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
用户定义链接函数的日志曝光链接功能,用于已知命运生存建模的glmer:
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")
}
文艺青年最爱的你要遇到麻烦,因为你yr和yr2预测(大概一年一年平方)是一个不寻常的链接功能相结合,使数值麻烦; 你可以使用glm2软件包来解决这个问题,但我会至少考虑一下,在这种情况下尝试适合平方年份是否合理.
更新:蛮力方法,mle2从下面开始; 尚未编写它来完成与交互的完整模型.
安德鲁·格尔曼的民间定理可能适用于此:
当您遇到计算问题时,通常会出现模型问题.
我开始尝试模型的简化版本,没有交互......
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没有完全相关,但它们完全是排名相关的,它们在数字上相互干扰当然不足为奇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可能是一个好主意.
我没有盲目地删除具有强相关参数估计的预测变量,但这可能是考虑它的合理时间.