Logistische Regression gibt Fehler zurück, läuft aber bei reduziertem Datensatz in Ordnung

Question

Ich würde mich sehr über Ihre Eingabe freuen!

ich auf einer logistischen Regression arbeite, aber es ist nicht aus irgendeinem Grund arbeiten:

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

Wenn ich das gleiche Modell mit weniger Daten laufen funktioniert es! Aber mit dem kompletten Datensatz bekomme ich einen Fehler und Warnmeldungen:

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 

Dies sind die Daten: https://www.dropbox.com/s/8ib8m1fh176556h/NSSH1.csv?dl=0

Log Exposition Link-Funktion von 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

tl; dr Sie bekommen Ärger, weil Ihre Prädiktoren yr und yr2 (vermutlich Jahr und Jahr-Quadrat) mit einer ungewöhnlichen Link-Funktion kombiniert werden n, um numerische Probleme zu verursachen; Sie können mit der glm2 package über das hinaus kommen, aber ich würde zumindest ein wenig darüber nachdenken, ob es sinnvoll ist, in diesem Fall das Quadratjahr zu verwenden.

Update: Brute-Force-Ansatz mit mle2 begann unten; habe es noch nicht geschrieben, um das vollständige Modell mit Interaktionen zu machen.

Andrew Gelmans folk theorem wahrscheinlich gilt auch hier:

Wenn Sie rechnerische Probleme haben, oft gibt es ein Problem mit Ihrem Modell.

Ich begann durch eine vereinfachte Version des Modells versucht, ohne die Wechselwirkungen ...

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

... der gut arbeitet. Jetzt versuchen wir zu sehen, wo das Problem ist:

OK, so haben wir Probleme, beide Interaktionen gleichzeitig. Wie hängen diese Prädiktoren tatsächlich zusammen? Mal sehen ...

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

yr und yr2 sind nicht perfekt korreliert sind, aber sie sind perfekt Rang korreliert und es ist sicherlich nicht verwunderlich, dass sie miteinander sind störende numerischUpdate: natürlich "Jahr" und "Jahr-Quadrat" sehen so aus! sogar poly(yr,2) verwenden, die eine orthogonale Polynom konstruiert, helfen in diesem Fall nicht ... noch, es lohnt sich bei den Parametern suchen, falls es einen Anhaltspunkt bietet ...

Wie oben erwähnt, können wir versuchen glm2 (ein Drop-in Ersatz für glm mit einem robusteren Algorithmus) und sehen, was passiert ...

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

Jetzt bekommen wir eine Antwort. Wenn wir cov2cor(vcov(mod5)) überprüfen, sehen wir, dass die yr und yr2 Parameter (und die Parameter für ihre Interaktion mit reLDM2 stark negativ korreliert (etwa -0,97). Lassen Sie uns das vergegenwärtigen ...

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

Was passiert, wenn wir versuchen, dies mit brutalen Gewalt zu tun?

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 *** 

Dies scheint vernünftig (Sie sollen vorhergesagte Werte überprüfen und sehen, dass sie Sinn machen ...), aber ...

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

Dies gibt ein mle2 Objekt, sondern ein mit einer verstümmelten Call-Slot, so ist es hässlich, die Ergebnisse zu drucken.

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? 

Versuchen von diesen zurückgegebenen Werte neu zu starten ...

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 

Jetzt können wir Konfidenzintervall bekommen ...

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 

Dies scheint zu funktionieren, aber ich würde sehr sein verdächtige dieser Passungen. Sie sollten wahrscheinlich andere Optimierer ausprobieren, um sicherzustellen, dass Sie zu den gleichen Ergebnissen zurückkehren können. Ein besseres Optimierungstool wie der AD Model Builder oder der Template Model Builder könnte eine gute Idee sein.

Ich halte nicht mit hirnlos dropping Prädiktoren mit stark korrelierten Parameter Schätzungen, aber dies könnte eine angemessene Zeit, um es zu betrachten.