2014-10-02 2 views
3

임의 효과 glmer로 사용자 정의 링크 기능을 활용하려고 시도하면서 나는 문제를 해결하는 방법을 모른다 :오류를주는 사용자 정의 링크 기능이있는 glmer : (maxstephalfit) PIRLS step-halvings가 deviance를 줄이지 못함

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

누구든지이 오류를 해결하는 방법에 대한 조언이 있습니까? 그것은 많은 방향을 제공하지 않습니다.

rpubs.com/bbolker/logregexp에 설명 된대로 새 링크 기능 (구체적으로 스케일 된 로그)을 정의하는 지침을 따르려고했지만 내 정의의 일부 측면이 잘못된 경우 놀라지 않을 것입니다. 내가 뭘 놓친 거 보여? 나는 S = 1 (동일해야 함)와이 링크를 참조 할 때 추정 밖으로 표준 이항 가족 (가정 로짓 링크)와 잘 작동하지만 오류 때문에

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') 
} 

이 구현에 문제가 있어야합니다.

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을 사용하고 있습니다.

답변

4

당신의 문제가 역 링크 기능을 "클램핑"하지 못하게 될 것이라고 생각했는데 (즉, 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/PI가 아닌 것처럼 보이지 않는 많은 값의 문제가 PIRLS 실패처럼 보입니다. 그들은 (nan는 즉각적인 실패를 유발하지 않고 C++ 코드를 통해 전파되는 것처럼 보입니다 ...)

+0

Argh. 이렇게 간단한 것을 놓친 것에 대한 사과드립니다. 두 개의 스칼라에 max를 사용하는 것이 좋다고 생각했지만, 분명히이 함수는 glmer 코드의 전체 벡터에 적용될 수 있습니다. 감사합니다. – MGL