Функция пользовательской Link работает для GLM, но не mgcv GAM

Question

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

Короче говоря,

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

• Я могу создать {статистики} семейного объект по этой ссылке и использовать его в «семейном» аргументе GLM.

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

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

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

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

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

Так что мой вопрос можно резюмировать следующим образом:

чего не хватает в этой пользовательской ссылке, которая работает на ГЖС, но не гам?

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

---

функция Link

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

Answer 1

Как вы знаете, glm занимает итеративно reweighted наименьших квадратов фитингов итераций. Ранняя версия 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"), то есть метод внешней итерации 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 (не может быть меньше для кубического сплайна), так как мои данные о игрушке близки к линейной и большой размерности базиса не требуется. Что еще более важно, это, по-видимому, предотвращает схождение производительности итерации производительности для моего набора данных игрушек.