Функция настраиваемой ссылки работает для GLM, но не для MGCV GAM

Приносим извинения, если ответ очевиден, но я потратил некоторое время, пытаясь использовать настраиваемую функцию ссылки в mgcv.gam

Суммируя,

  • Я хочу использовать измененную пробит-ссылку из пакета psyphy (я хотите использовать psyphy.probit_2asym, я называю это custom_link )

  • С помощью этой ссылки я могу создать семейный объект {stats} и использовать его в аргументе "семья" glm.

    m <- glm(y~x, family=binomial(link=custom_link), ... )

  • Не работает при использовании в качестве аргумента для gam {mgcv}.

    m <- gam(y~s(x), family=binomial(link=custom_link), ... )

    Я получаю сообщение об ошибке Error in fix.family.link.family(family) : link not recognised

Я не понимаю причину этой ошибки, и glm, и gam работают, если я укажу стандартный link=probit.

Итак, мой вопрос можно резюмировать так:

чего не хватает в этой настраиваемой ссылке, которая работает для glm, но не работает?

Заранее спасибо, если вы подскажете, что мне делать.

Функция ссылки

probit.2asym <- function(g, lam) {
    if ((g < 0 ) || (g > 1))
        stop("g must in (0, 1)")
    if ((lam < 0) || (lam > 1))
        stop("lam outside (0, 1)")
    linkfun <- function(mu) {
        mu <- pmin(mu, 1 - (lam + .Machine$double.eps))
        mu <- pmax(mu, g + .Machine$double.eps)
        qnorm((mu - g)/(1 - g - lam))
        }
    linkinv <- function(eta) {
        g + (1 - g - lam) * 
         pnorm(eta)
        }
    mu.eta <- function(eta) {
        (1 - g - lam) * dnorm(eta)      }
    valideta <- function(eta) TRUE
    link <- paste("probit.2asym(", g, ", ", lam, ")", sep = "")
    structure(list(linkfun = linkfun, linkinv = linkinv, 
    mu.eta = mu.eta, valideta = valideta, name = link), 
    class = "link-glm")
}

r regression glm gam mgcv
person user1436340    schedule 30.09.2016    source источник
comment
Большое спасибо за редактирование вопроса.   -  person user1436340    schedule 03.10.2016

Ответы (1)


arrow_upward
4
arrow_downward

Как вы, возможно, знаете, glm выполняет подгонку итеративно пересмотренным методом наименьших квадратов. Ранняя версия gam расширяет это, устанавливая итеративно штрафованные пересмотренные наименьшие квадраты, что выполняется функцией gam.fit. В некотором контексте это называется итерацией производительности.

С 2008 года (или, может быть, немного раньше) gam.fit3 на основе так называемой внешней итерации заменил gam.fit как gam по умолчанию. Такое изменение требует дополнительной информации о семействе, о которой вы можете прочитать о ?fix.family.link.

Основное различие между двумя итерациями заключается в том, вложены ли итерация коэффициентов beta и итерация параметров сглаживания lambda.

  • Итерация производительности осуществляется вложенным способом, где для каждого обновления beta выполняется одна итерация lambda;
  • Внешняя итерация полностью разделяет эти 2 итерации, при этом для каждого обновления beta итерация lambda доводится до конца до сходимости.

Очевидно, что внешняя итерация более стабильна и с меньшей вероятностью пострадает от сбоя сходимости.

gam имеет аргумент optimizer. По умолчанию он принимает optimizer = c("outer", "newton"), то есть метод ньютона внешней итерации; но если вы установите optimizer = "perf", потребуется итерация производительности.

Итак, после приведенного выше обзора у нас есть два варианта:

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

Я ленился, поэтому продемонстрирую второй вариант (на самом деле я не слишком уверен, чтобы использовать первый подход).

Воспроизводимый пример

Вы не предоставили воспроизводимый пример, поэтому я готовлю его, как показано ниже.

set.seed(0)
x <- sort(runif(500, 0, 1))    ## covariates (sorted to make plotting easier)
eta <- -4 + 3 * x * exp(x) - 2 * log(x) * sqrt(x)   ## true linear predictor
p <- binomial(link = "logit")$linkinv(eta)    ## true probability (response)
y <- rbinom(500, 1, p)    ## binary observations

table(y)    ## a quick check that data are not skewed
#  0   1 
#271 229

Я возьму g = 0.1 и lam = 0.1 функции probit.2asym, которую вы собираетесь использовать:

probit2 <- probit.2asym(0.1, 0.1)

par(mfrow = c(1,3))

## fit a glm with logit link
glm_logit <- glm(y ~ x, family = binomial(link = "logit"))
plot(x, eta, type = "l", main = "glm with logit link")
lines(x, glm_logit$linear.predictors, col = 2)

## glm with probit.2asym
glm_probit2 <- glm(y ~ x, family = binomial(link = probit2))
plot(x, eta, type = "l", main = "glm with probit2")
lines(x, glm_probit2$linear.predictors, col = 2)

## gam with probit.2aysm
library(mgcv)
gam_probit2 <- gam(y ~ s(x, bs = 'cr', k = 3), family = binomial(link = probit2),
                   optimizer = "perf")
plot(x, eta, type = "l", main = "gam with probit2")
lines(x, gam_probit2$linear.predictors, col = 2)

введите описание изображения здесь

Я использовал естественный кубический сплайн cr для s(x), так как для одномерного сглаживания настройка по умолчанию с тонкопластинчатым сплайном не нужна. Я также установил небольшой базовый размер k = 3 (не может быть меньше для кубического сплайна), так как мои данные о игрушке близки к линейным, и большой базовый размер не нужен. Что еще более важно, это, похоже, предотвращает сбой сходимости итерации производительности для моего набора данных игрушек.

person Zheyuan Li    schedule 01.10.2016
comment
Большое спасибо за ответ. Использование оптимизатора производительности - простой ответ на вопрос. Я изучу эту внешнюю итерацию, о которой я не знал, что это новый стандарт (очевидно, пропустил этот раздел при чтении книги Саймона Вуда). На данный момент я эмпирически исследую, не ухудшает ли этот оптимизатор по умолчанию больше производительность согласования и сходимость для моей проблемы, и сообщу о своих выводах здесь. Я мог бы расширить и поделиться ссылками для работы с внешней итерацией. - person user1436340; 03.10.2016