определяемая пользователем "отрицательная экспоненциальная" ссылка glm - PullRequest
Новая
       6

определяемая пользователем "отрицательная экспоненциальная" ссылка glm

0 голосов
/ 14 ноября 2018

Я пытался следовать этому примеру изменить glm ... пользовательскую функцию ссылки в r но я получаю ошибки. У меня есть двоичные данные, и я хотел бы изменить функцию ссылки с «logit» на отрицательную экспоненциальную ссылку. Я хочу предсказать

вероятность успеха (p) = 1-exp (линейный предиктор)

Причина, по которой мне нужна эта ссылка вместо одной из встроенных ссылок, заключается в том, что p увеличивается выпуклым образом между 0 и 0,5, но только для "logit", "cloglog", "probit" и "cauchy" позволяют вогнутую форму. См. Прилагаемое фото для справки: прогнозируемые p против биннинг-наблюдений

Имитация данных 

location<-as.character(LETTERS[rep(seq(from=1,to=23),30)])
success<-rbinom(n=690, size=1, prob=0.15)
df<-data.frame(location,success)
df$random_var<-rnorm(690,5,3)
df$seedling_size<-abs((0.1+df$success)^(1/df$random_var))
df<-df[order(df$location)]

Создание пользовательской функции связи. Примечание: эта = линейный предиктор, мю = вероятность 

negex<-function(){
##link
linkfun<-function(mu) log(-mu+1)
linkinv<-function(eta) 1-exp(eta)
## derivative of inverse link with respect to eta
mu.eta<-function(eta)-exp(eta)
valideta<-function(eta) TRUE
link<-"log(-mu+1)"
structure(list(linkfun=linkfun,linkinv=linkinv,
             mu.eta=mu.eta,valideta=valideta,
             name=link),
        class="link-glm")
}

Успех модели как функция размера рассады 

negexp<-negex()

model1<-glm(success~seedling_size,family=binomial(link=negexp),data=df)

Ошибка: не найден действительный набор коэффициентов: укажите начальные значения

Модель с использованием glmer (Моя конечная цель) 

model2<-glmer(success~seedling_size+ (1|location),family=binomial(link=negexp),data=df)

Ошибка в (функция (fr, X, reTrms, семейство, nAGQ = 1L, подробный = 0L, maxit = 100L,: (maxstephalfit) PIRLS-шаг-половинкам не удалось уменьшить отклонение в pwrssUpdate

Я получаю разные сообщения об ошибках, но я думаю, что проблема одна и та же, независимо от того, используется ли glmer или glm, и это то, что моя функция ссылки как-то не так.

1 Ответ

0 голосов
/ 15 ноября 2018

Я нашел ответ. Наиболее полезной была эта R тема от 2016 . Было 2 вопроса. Во-первых, моя функция связи была неправильной. Я исправил это так: 

negex <- function() 
 { 
 linkfun <- function(mu) -log(1-mu) 
 linkinv <- function(eta) 1-exp(-eta) 
 mu.eta <- function(eta) exp(-eta) 
 valideta <- function(eta) all(is.finite(eta)&eta>0) 
 link <- paste0("negexp") 
 structure(list(linkfun = linkfun, linkinv = linkinv, 
             mu.eta = mu.eta, valideta = valideta, name = link), 
        class = "link-glm") 
}

Во-вторых, модель требует определенных начальных значений. Они будут уникальными для ваших данных. Вот первые несколько строк данных, для которых я действительно нашел решение: 

   site plot sub_plot oak_success oak_o1_gt05ft..1
  0001   10        3           1                0
  0001   12        2           0                0
  0001   12        3           0                0
  0001   12        4           0                0
  0001   13        4           0                0

Я не знаю, как опубликовать полные данные на этом сайте, но если кто-то хочет, чтобы он запустил пример, отправьте мне электронное письмо: lake.graboski@gmail.com

negexp<-negex()

Надеюсь, это поможет кому-то в будущем, потому что я не нашел других примеров решения этой проблемы при переполнении стека или в сети. Используя новые начальные значения, я смог запустить модель: 

starting_values<-c(1,0) #1 for the intercept and 0 for the slope
h_gt05_solo_negex2<-glm(oak_success~ oak_o1_gt05ft..1 , 
                    family=binomial(link=negexp),start=starting_values,data=rocdf)
summary(h_gt05_solo_negex2)
Call:
glm(formula = oak_success ~ oak_o1_gt05ft..1, family = binomial(link = negexp), 
data = lt40, start = starting_values)

Deviance Residuals: 
Min       1Q   Median       3Q      Max  
-1.3808  -0.4174  -0.2637  -0.2637   2.5985  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.034774   0.005484   6.341 2.28e-10 ***
oak_o1_gt05ft..1 0.023253   0.002187  10.635  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1416.9  on 2078  degrees of freedom
Residual deviance: 1213.5  on 2077  degrees of freedom
AIC: 1217.5

Number of Fisher Scoring iterations: 6

Были некоторые проблемы с конвергенцией. Когда высота проростков (oak_o1_gt05ft..1) превысила 40 футов, оценки параметров стали ненадежными проблемы сходимости . У меня было очень мало наблюдений в этом диапазоне, поэтому я ограничил данные наблюдениями, в которых предиктор был <40 футов, и перезапустил модель. Я также включил «сайт» (так же, как «местоположение» в смоделированных данных)). На <a href="https://i.stack.imgur.com/cTTDV.png" rel="nofollow noreferrer"> этой фигуре показаны прогнозы успеха дуба в отношении высоты проростков дуба для каждого участка / места (черные кружки), групповые наблюдения за успехами / образцами (большие зеленые точки) и прогноз вероятность успеха без фактора сайта (синяя линия). Похоже, что наклон переменной размера рассады более точен при учете участка.

К сожалению, мне не удалось заставить эту модель работать с блеском, поэтому сайт должен был быть включен как фиксированный эффект, поэтому стандартные ошибки и оценки наклона для высоты проростков дуба могут быть немного консервативными.

...