自定义链接函数适用于GLM，但不适用于mgcv GAM。

Question

自定义链接函数适用于GLM，但不适用于mgcv GAM。

rregressionglmgammgcv

3

如果答案很显然，那么对不起，但我已经花了相当多的时间尝试在mgcv.gam中使用自定义链接函数。

简而言之,

我想要使用来自包psyphy的修改后的Probit链接（我想使用psyphy.probit_2asym，我称之为custom_link）
我可以使用此链接创建{stats}family对象，并在glm的“family”参数中使用它。

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中不行的自定义链接函数缺少什么？

如果您能给我提示我该怎么做，我将不胜感激。

链接函数

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

- user1436340

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Zheyuan Li · Accepted Answer

作为您可能已经知道的，glm 使用迭代加权最小二乘拟合迭代。 gam 的早期版本通过拟合一个迭代惩罚的加权最小二乘来扩展这一点，这是由gam.fit函数完成的。这在某些情况下被称为性能迭代。

自2008年以来（甚至更早），基于所谓的外部迭代的gam.fit3已经取代了gam.fit成为默认值。这种变化确实需要有关家族的一些额外信息，有关此信息可以阅读?fix.family.link。

两个迭代之间的主要区别是系数beta的迭代和平滑参数lambda的迭代是否嵌套。

性能迭代采用嵌套方式，在每次更新beta时，执行单个lambda迭代；
外部迭代完全独立这两个迭代，在每次更新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

我将使用函数probit.2asym中的g = 0.1和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（对于三次样条不能更小），因为我的玩具数据接近线性，不需要大的基函数维度。更重要的是，这似乎可以防止我的玩具数据集在性能迭代中出现收敛失败的情况。