Glmer with custom link function giving error: (maxstephalfit) PIRLS step-halvings failed to reduce deviation
While trying to use a random effect user-defined link function, I ran into an error that I don't know how to troubleshoot:
Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate
Does anyone have any tips on how this error could be resolved? It doesn't provide much direction.
I tried to follow the instructions for defining a new link function (in particular with a scaled logit) as described at rpubs.com/bbolker/logregexp , but I wouldn't be surprised if some aspect of my definition is wrong. See what I am missing?
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')
}
There must be something wrong with this implementation, as the estimate works fine with the standard binomial family (assumed boolean reference), but errors occur when I reference this reference with s = 1 (which should be identical). Sample data can be generated as follows:
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)]
Then I tried to evaluate a simple model
logit_link <- glmer(cbind(deg1, indict-deg1) ~ (1|ID) + (1|court), family=binomial, data=test_data[indict > 0])
and the corresponding alternative
scaled_link <- glmer(cbind(deg1, indict-deg1) ~ (1|ID) + (1|court), family=binomial(link=scaled_logit()), data=test_data[indict > 0])
Any ideas would be appreciated! I am using lme4 1.1.6 on R 3.0.3.
source to share
I thought that your problem would not be to "pinch" the backbind function (ie, keep the results strictly from 0 to 1), but it turns out (I think) it is much easier - just confusion max()
and pmax()
. ( max()
has a rather dangerous design!) This works for me:
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')
}
However, it would probably be a good idea for future reliability to do this pmax(epsilon,...)
, rather than pmax(0,...)
restrict the backbind function between epsilon
and 1-epsilon
(where epsilon
something like 1e-6).
PS we (the lme4
maintainers) should probably try to put in even more robust error checking in the PIRLS step - there are a lot of issues with NaN
/ non-final values looking like PIRLS crashes when that's not what they are ( NaN
seems to propagate through C code ++ without triggering any immediate crashes ...)
source to share