Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate when trying to run binomial GLMER with user-defined link function

Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate when trying to run binomial GLMER with user-defined link function

这里是相对较新的 R 用户,正在尝试 运行 将暴露天数与二项式响应变量(成功 = 1,失败 = 0)结合在一起的 GLMER 鸟巢成功。我将 Ben Bolker 的代码用于用户定义的 link 函数,获得 here.

完整代码如下:

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)

我有连续、二项式和分类预测变量,我首先重新调整了连续变量(包括曝光天数)。

问题是我不断收到错误消息:

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

我发现 NaN 表示 "not a number" 但不确定为什么会这样。非常感谢任何帮助!

下面是我的嵌套成功数据集:

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

tl;dr 您不能在 power-logit link 函数中使用非正曝光值。我部分地通过思考解决了这个问题,部分地通过在 link 函数中放置战略性的 print() 语句来弄清楚发生了什么。所以:

  • 你不应该缩放 EXPOSURE 变量
  • 即使在不缩放之后,仍然有一个观察值(第 55 行)具有零曝光值 - 我省略了它(我猜你也可以将它设置为一个小的非零值)
  • 这样调整后模型好像运行OK
  • 您可能正在尝试将过于复杂的模型拟合到此数据集 - 您有 35 次失败和 28 次成功 -> 有效样本量为 28(sensu Harrell 的 回归建模策略 书)-> 你可能无法可靠地适应最多 3 个参数
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)

...在此之后,模型似乎适合 OK...(尽管年内变化在拟合中减少到几乎为零。 您也接近随机效应分组因子水平数量的实际下限)。