自定义链接功能适用于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

我没有得到这个错误的原因，如果我指定标准link=probit，glm和gam都可以工作。

所以我的问题可以概括为：

缺什么，用于GLM的工作原理是自定义链接而不是GAM？

在此先感谢您，如果您能给我提示我该怎么做。

---

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需要迭代重加权最小二乘法拟合迭代。 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和功能probit.2asym的lam = 0.1要使用：

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（对于三次样条不能更小）。更重要的是，这似乎可以防止我的玩具数据集的性能迭代失败。