首页 文章

glmer用户定义的链接函数给出错误：（maxstephalfit）PIRLS步长减半未能减少偏差

提问于
浏览 次
3

在尝试利用具有随机效应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 .

r lme4

1 回答

  • 4

    我认为你的问题会变成没有"clamp"反向链接函数（即严格保持结果在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,...) 而是约束 epsilon 和 1-epsilon 之间的反向链接函数（其中 epsilon 类似于1e-6） .

    PS我们（ lme4 维护者）应该尝试在PIRLS步骤中插入一些更强大的错误检查 - 很多问题与 NaN /非有限值弹出看起来像PIRLS失败，当它不是真正的它们时（ nan s似乎通过C代码传播而不会引发任何直接的故障...）

    回复于

相关问题