2016-05-30 15 views
0

Relativ neue R Benutzer hier, versuchen, einen GLMER für Vogel Nest Erfolg, der Tage der Exposition mit enthält eine binomische Antwortvariable (Erfolg = 1, Fehler = 0). Ich benutze Ben Bolker Code für eine benutzerdefinierte Link-Funktion, erhalten here.Fehler: (Maxstephalfit) PIRLS Schritt-Halvings konnte Deviance in pwrssUpdate nicht reduzieren, wenn versuchen, Binomial GLMER mit benutzerdefinierten Link-Funktion auszuführen

Hier ist der vollständige Code:

NestSuccessExposure<-read.csv("NestSuccessExposure.csv") 
NestSuccessExposure<-na.omit(NestSuccessExposure) 
NestSuccessExposureScaled<-scale(NestSuccessExposure[,c(3:10)]) 
NestSuccessExposureScaled<-cbind(NestSuccessExposureScaled,NestSuccessExposure[,c(1,2,11,12,13)]) 

library(MASS) 
library(lme4) 
logexp <- function(exposure = 1) 
{ 
    linkfun <- function(mu) qlogis(mu^(1/exposure)) 

    linkinv <- function(eta) plogis(eta)^exposure 
    logit_mu_eta <- function(eta) { 
    ifelse(abs(eta)>30,.Machine$double.eps, 
      exp(eta)/(1+exp(eta))^2) 
    } 
    mu.eta <- function(eta) {  
    exposure * plogis(eta)^(exposure-1) * 
     logit_mu_eta(eta) 
    } 
    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") 
} 

#Run GLMER (linear mixed-effects model) 

NSExposure<-glmer(SUCCESS~SLOPE+BEERS+TALL+CANOPY+DISTROAD+DISTSTREAM+NESTHT+ 
QUAL+VINENT+MGMT+(1|YEAR), 
    family=binomial(link=logexp(NestSuccessExposureScaled$EXPOSURE)), 
    data=NestSuccessExposureScaled) 

summary(NSExposure) 

Ich habe kontinuierlich, binomische und kategorischen Prädiktorvariablen und ich neu skaliert zuerst die kontinuierlichen Variablen (einschließlich der Exposition Tage).

Das Problem ist, dass ich die Fehlermeldung erhalte:

"Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate. In addition: Warning message: In qlogis(mu^(1/exposure)) : NaNs produced"

ich gefunden habe, dass NaN bedeutet „keine Zahl“, aber ist nicht sicher, warum dies auftritt. Jede Hilfe wird sehr geschätzt!

Unten ist mein Nest Erfolg Dataset:

dput(NestSuccessExposure) 
structure(list(SUCCESS = c(0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 
0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 
1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 
1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 1L, 0L), YEAR = c(2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2014L, 2014L, 
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
2014L, 2014L, 2014L, 2014L, 2015L, 2015L, 2015L, 2015L, 2015L, 
2015L, 2015L, 2015L, 2015L, 2015L, 2015L), EXPOSURE = c(17, 17, 
9, 9, 10, 2, 9, 9, 7, 1, 9, 7, 15, 3, 17, 8, 8, 1.5, 24, 6, 23, 
2, 22, 5.5, 2.5, 17, 16, 1.5, 0.5, 7.5, 21, 1.5, 17, 3.5, 1, 
0.5, 4, 1, 17, 5, 28, 22, 12, 5, 6, 2, 14, 21, 22, 9, 15, 1, 
11, 15, 0, 14, 9, 6, 4, 14, 23, 3, 2, 4.5), SLOPE = c(12, 35.5, 
48.75, 18, 31, 55, 27.5, 31, 31.25, 24.75, 43.75, 36, 18.75, 
17, 19.75, 45, 25, 42, 17.5, 29.5, 39, 26, 27, 39.5, 47.5, 44.5, 
51.5, 51, 48.5, 6, 45.25, 33.75, 36, 19, 22.75, 57.25, 30.5, 
54, 30.75, 34, 38, 20, 27, 48, 10.25, 5.25, 38.5, 39, 23, 12.75, 
48, 51, 35.5, 30.5, 24, 53.5, 67.5, 4, 68.5, 73, 54.5, 18.5, 
41.5, 7), BEERS = c(1.97, 1.71, 1.42, 0.44, 1.99, 1.8, 1.09, 
1.99, 1.92, 0.05, 1.39, 1.78, 1.71, 1.71, 1.17, 1.66, 1.71, 1.98, 
2, 0.12, 1.85, 1.99, 0.08, 0.67, 1.82, 1.91, 1.97, 1.97, 1.97, 
0.01, 1.57, 0.03, 0.05, 0.74, 2, 1.95, 1.39, 1.71, 1.85, 0.64, 
0.71, 1.97, 0.07, 1.98, 0.03, 0.46, 1.97, 1.07, 1.17, 1.84, 2, 
1.96, 1.97, 1.57, 0.91, 2, 1.98, 0.84, 1.98, 1.91, 1.97, 1.91, 
1.84, 1.6), TALL = c(23.85, 28.9, 22.925, 25, 14.7, 24.6, 17.3, 
24.05, 18.675, 25.6, 20.15, 32.75, 28.5, 20.075, 21.35, 23.425, 
27.2, 25.025, 21.9, 20.85, 22.75, 23.35, 29.05, 21.125, 29.65, 
27.575, 27.175, 27.95, 29.35, 23.225, 26.35, 27.6, 23.2, 22.05, 
22.45, 32, 26.267, 23.9, 21.6, 23.55, 23.9, 26.93, 14.98, 23.55, 
24.8, 17.6, 31.38, 28.05, 22.95, 24.8, 27.2, 31.83, 30.48, 21.63, 
18.3, 26.33, 23, 24.95, 14.75, 21.35, 25.55, 23.2, 26.05, 20.45 
), CANOPY = c(0.9, 0.85, 0.9, 0.95, 0.95, 0.95, 1, 0.85, 0.7, 
0.95, 0.95, 0.8, 0.9, 0.9, 1, 1, 0.95, 0.9, 0.85, 0.8, 0.6, 0.85, 
0.85, 1, 1, 0.95, 0.95, 0.75, 0.9, 0.75, 0.95, 0.9, 0.85, 0.95, 
0.9, 0.85, 0.9, 0.95, 0.85, 0.85, 1, 0.85, 0.9, 0.95, 0.4, 0.75, 
0.8, 0.95, 0.85, 0.9, 0.9, 0.8, 0.75, 0.85, 0.8, 0.95, 0.75, 
0.9, 0.85, 0.65, 0.8, 0.85, 0.85, 0.85), DISTROAD = c(6.81, 19.02, 
158.83, 7.56, 70.87, 263.68, 31.28, 39.32, 181.36, 246.67, 48.58, 
38.63, 75.47, 4.51, 80.56, 362.92, 82.1, 197.36, 361.38, 82.29, 
59.05, 31.32, 81.46, 179.79, 211.74, 238.64, 270.93, 318.72, 
329.96, 14.53, 158.3, 104.38, 94.61, 293.89, 99.84, 178.64, 56.28, 
204.82, 425.32, 4.02, 119.75, 8.31, 1.17, 63.24, 4.62, 119.46, 
65.45, 121.6, 4.38, 17.36, 205.12, 243.77, 349.98, 3.98, 60.29, 
209.79, 247.05, 114.51, 100.86, 331.62, 306.9, 0.95, 27.26, 34.58 
), DISTSTREAM = c(348.37742, 233.394296, 149.503103, 181.305039, 
138.067932, 13.087182, 58.590507, 288.128836, 149.061785, 126.667503, 
220.6535, 194.269426, 115.265352, 275.771326, 78.528016, 118.476224, 
174.095054, 44.340495, 82.174798, 62.745934, 227.662779, 55.671038, 
200.910084, 34.939781, 96.984957, 42.842466, 25.45655, 72.374004, 
32.353038, 134.254137, 184.571521, 58.354152, 78.176574, 22.294154, 
137.078915, 51.516704, 244.946159, 16.681571, 62.556975, 80.092302, 
84.607328, 336.424256, 23.248202, 206.251649, 279.365454, 14.091524, 
198.226118, 534.630654, 102.796308, 217.190835, 15.414713, 43.516136, 
42.080907, 67.19459, 417.021499, 64.101315, 14.218346, 1.932141, 
34.491616, 38.305397, 20.481007, 411.768426, 404.101887, 8.842676 
), NESTHT = c(12.6, 26.6, 18, 10, 14, 18, 15.5, 23, 17, 21.5, 
26, 14.5, 18, 29.5, 12.5, 16, 16.5, 24, 20, 14.2, 16.8, 20.2, 
13.4, 20.4, 25.8, 19.6, 18, 18.4, 15.2, 24, 19.8, 16.8, 20.4, 
20, 15, 18, 24, 17.2, 13.4, 23, 16.8, 9, 20, 26.8, 22.2, 13, 
26.4, 14.6, 11.4, 15, 20.6, 20, 14, 24.5, 21.8, 18.8, 11.2, 21.5, 
12.4, 19.4, 17.2, 15.6, 13, 21), QUAL = c(0L, 0L, 0L, 0L, 1L, 
1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 
0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L), VINENT = c(0L, 0L, 
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 
1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), MGMT = structure(c(2L, 
2L, 3L, 2L, 2L, 2L, 3L, 3L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 2L, 2L, 
2L, 1L, 2L, 2L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 1L, 3L, 3L, 
2L, 3L, 2L, 3L, 3L, 1L, 3L, 1L, 2L, 2L, 2L, 1L, 3L, 2L, 1L, 3L, 
1L, 3L, 2L, 2L, 3L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 3L), .Label = c("C", 
"E", "U"), class = "factor")), .Names = c("SUCCESS", "YEAR", 
"EXPOSURE", "SLOPE", "BEERS", "TALL", "CANOPY", "DISTROAD", "DISTSTREAM", 
"NESTHT", "QUAL", "VINENT", "MGMT"), row.names = c(1L, 2L, 3L, 
4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 15L, 16L, 17L, 19L, 20L, 21L, 
23L, 24L, 26L, 27L, 28L, 30L, 31L, 33L, 34L, 35L, 36L, 37L, 39L, 
40L, 41L, 42L, 43L, 44L, 45L, 47L, 48L, 52L, 53L, 54L, 55L, 57L, 
58L, 59L, 60L, 61L, 62L, 63L, 65L, 67L, 68L, 69L, 70L, 72L, 73L, 
74L, 77L, 78L, 80L, 81L, 82L, 88L, 90L), class = "data.frame", na.action = structure(c(12L, 
13L, 14L, 18L, 22L, 25L, 29L, 32L, 38L, 46L, 49L, 50L, 51L, 56L, 
64L, 66L, 71L, 75L, 76L, 79L, 83L, 84L, 85L, 86L, 87L, 89L), .Names = c("12", 
"13", "14", "18", "22", "25", "29", "32", "38", "46", "49", "50", 
"51", "56", "64", "66", "71", "75", "76", "79", "83", "84", "85", 
"86", "87", "89"), class = "omit")) 
+0

Kurzantwort: 'glmer' ist empfindlich für Rundungsfehler/numerische Fehler, die zu einer Wahrscheinlichkeit von Null von der Funktion des umgekehrten Links führen (zB siehe Kommentare in der Antwort unter http://stackoverflow.com/questions/26152986/ glmer-with-user-defined-link-function-gib-error-maxstephalfit-pirls-step-h? rq = 1). Dies ist wahrscheinlich reparierbar, indem entsprechende "Klemm" -Anweisungen in die Familiendefinition eingefügt werden. Können Sie bitte Daten und/oder Code angeben, der uns ein [reproduzierbares Beispiel] liefert (http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example)? –

+0

Ah, ich hatte diese Frage während der Suche gesehen, war mir aber nicht sicher, ob sie hier angewendet wurde. Bearbeitete meine ursprüngliche Frage, um den Datensatz einzuschließen. Vielen Dank! – Setophaga

Antwort

1

tl; dr Sie können nicht kraft- Belichtungswerte in einem Power-Logit-Link-Funktion haben. Ich habe das zum Teil durch Nachdenken herausgefunden, zum Teil indem ich strategische print() Anweisungen in die Link-Funktionen setzte, um herauszufinden, was vor sich ging. Also:

  • Sie sollten nicht die EXPOSURE Variable
  • auch nach nicht-Skalierung werden Skalierung, gibt es noch eine Beobachtung (Zeile 55) mit einem Null-Belichtungswert - ich weggelassen es (ich glaube, man auch einstellen könnte es auf einen kleinen Wert ungleich Null)
  • danach eingestellt wird, das Modell
  • Sie wahrscheinlich zu laufen OK scheint versuchen zu passen komplexer ein Modell auf diese Daten gesetzt - Sie haben 35 Ausfälle und 28 Erfolge -> eine effektive Stichprobengröße von 28 (sensu Harrel l's Regression Modellierung Strategien Buch) -> Sie können sich wahrscheinlich nicht zuverlässig passen mehr als höchstens 3 Parameter
NestSuccessExposure <- na.omit(NestSuccessExposure) 
noscaleVars <- c("SUCCESS","YEAR","EXPOSURE","QUAL","VINENT","MGMT") 
noscaleCols <- which(names(NestSuccessExposure) %in% noscaleVars) 
NestSuccessExposureScaled<- NestSuccessExposure 
NestSuccessExposureScaled[-noscaleCols] <- 
    scale(NestSuccessExposure[-noscaleCols]) 
NestSuccessExposureScaled <- subset(NestSuccessExposureScaled, 
            EXPOSURE>0) 

... danach das Modell OK zu passen scheint ... (obwohl die unterjährige Variation in der Anpassung auf fast Null reduziert wird. Sie befinden sich auch in der Nähe der unteren praktischen Grenze für die Anzahl der Ebenen eines Gruppierungsfaktors für zufällige Effekte.

+0

Wow, vielen Dank für all die Erkenntnisse und Einsichten!Ich habe eine Menge von NAs für die Variable "BEERS" (Beers 'Aspekt), da flacher Grund keinen Aspektwert hat. Daher sollte die Auslassung dieses Parameters die Stichprobengröße etwas erhöhen. Ich werde auch andere Parameter reduzieren. Ist es schädlich, weiterhin ein Mixed-Effects-Modell zu verwenden (da wir wissen, dass der Erfolg in einigen Jahren höher ist als andere) oder würden Sie stattdessen den Wechsel zu Binomial-GLM empfehlen? Danke noch einmal! – Setophaga