Логистическая регрессия возвращает ошибку, но работает нормально при уменьшении набора данных

Question

Буду признателен за ваш вклад в это много!

Я работаю на логистическую регрессию, но он не работает по какой-то причине:

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") 
} 

Answer 1

ТЛ; др вы попадаете в неприятности, потому что ваши предсказатели yr и yr2 (предположительно год и год) объединены с необычной связью functio n вызвать числовые проблемы; вы можете пройти мимо этого, используя glm2 package, но я бы дал хотя бы немного подумать о том, имеет ли смысл попытаться соответствовать квадратному сроку года в этом случае.

Обновление: подход с грубой силой с mle2 запущен ниже; еще не написали его, чтобы сделать полную модель с взаимодействиями.

Эндрю Гельман 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 не совершенно коррелируют, но они совершенно рангом-коррелируют, и это, конечно, не удивительно, что они мешают друг другу численноОбновление: конечно, «год» и «год-квадрат» выглядят так! даже при использовании 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? 

Try перезапуск из этих возвращаемых значений ...

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.

Я не придерживаюсь бессмысленно сбрасывает предикторов с сильно коррелированными оценками параметров, но это может быть разумное время для его рассмотрения.