r - 带有用户定义链接函数的 glmer 给出错误 : (maxstephalfit) PIRLS step-halvings failed to reduce deviance

Question

在尝试将用户定义的链接函数与随机效果 glmer 结合使用时，我遇到了一个我不知道如何解决的错误:

Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

有人对如何解决此错误有任何建议吗？它没有提供太多方向。

我已尝试按照 rpubs.com/bbolker/logregexp 中概述的说明定义新的链接函数(特别是缩放的 logit) ，但如果我的定义的某些方面不正确，我不会感到惊讶。看到我遗漏了什么吗？

scaled_logit <- function(s = 1) {
  linkfun <- function(mu) log( max(0, mu / (s-mu)) )
  linkinv <- function(eta) s / (1 + exp(-eta))
  mu.eta <- function(eta) s * exp(-eta) / (1 + exp(-eta))^2
  valideta <- function(eta) TRUE
  link <- paste0('scaled_logit(',s,')')
  structure(list(linkfun = linkfun, linkinv = linkinv, mu.eta = mu.eta, valideta = valideta, name = link), class = 'link-glm')
}

这个实现一定有问题，因为估计在标准二项式族(假定的 logit 链接)上工作正常，但是当我用 s=1(应该是相同的)引用这个链接时出错。示例数据可以生成如下:

library(data.table)

courts <- 50
test_courts <- data.table(court = 1:courts,
                     court_factor = pmax(0, rnorm(courts, mean=1, sd=0.25)))
setkey(test_courts, court)
pros <- 100
test_pros <- data.table(ID = 1:pros,
                       deg1_rate = pmax(0, rnorm(pros, mean=0.02, sd=0.0075)))
setkey(test_pros, ID)

test_data <- data.table(expand.grid(ID = 1:pros, court = 1:courts))
setkeyv(test_data, c('ID','court'))
test_data <- merge(test_data, test_courts, by='court', all.x=TRUE)
test_data <- merge(test_data, test_pros  , by='ID'   , all.x=TRUE)

test_data[ , indict := sample(0:20, nrow(test_data), replace=TRUE)]
test_data[ , deg1 := rbinom(pros*courts, size=indict, prob=court_factor*deg1_rate)]

然后我一直在尝试估计简单模型

logit_link <- glmer(cbind(deg1, indict-deg1) ~ (1|ID) + (1|court), family=binomial, data=test_data[indict > 0])

和相应的选择

scaled_link <- glmer(cbind(deg1, indict-deg1) ~ (1|ID) + (1|court), family=binomial(link=scaled_logit()), data=test_data[indict > 0])

如有任何见解，我们将不胜感激!我在 R 3.0.3 上使用 lme4 1.1.6。

Answer 1

我原以为您的问题是未能“钳制”反向链接函数(即将结果严格保持在 0 和 1 之间)，但事实证明(我认为)它比这简单得多 - - 只是混淆了 max() 和 pmax()。 (max() 有一个非常危险的设计!)这对我有用:

scaled_logit <- function(s = 1) {
  linkfun <- function(mu) log( pmax(0, mu / (s-mu)) )
  linkinv <- function(eta) s / (1 + exp(-eta))
  mu.eta <- function(eta) s * exp(-eta) / (1 + exp(-eta))^2
  valideta <- function(eta) TRUE
  link <- paste0('scaled_logit(',s,')')
  structure(list(linkfun = linkfun, linkinv = linkinv, mu.eta = mu.eta, valideta = valideta, name = link), class = 'link-glm')
}

也就是说，为了将来的稳健性，使用 pmax(epsilon,...) 而不是 pmax(0,...) 可能是个好主意> 并限制 epsilon 和 1-epsilon 之间的反向链接函数(其中 epsilon 类似于 1e-6)。

PS 我们(lme4 维护者)可能应该尝试在 PIRLS 步骤中插入一些更强大的错误检查——NaN/non-finite 有很多问题弹出的值看起来像 PIRLS 故障，但实际上并非如此(nan 似乎通过 C++ 代码传播而不会立即触发任何故障...)