自定义 Link 函数适用于 GLM 但不适用于 mgcv GAM
Custom Link function works for GLM but not mgcv GAM
抱歉,如果答案很明显,但我花了很多时间尝试在 mgcv.gam
中使用自定义 link 函数
简而言之,
- 我想使用包 psyphy ( I want to use psyphy.probit_2asym 中的修改概率 link,我称之为
custom_link
)
我可以用这个 link 创建一个 {stats} 系列对象并在 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 都可以工作 link=probit
。
所以我的问题可以概括为:
这个适用于 glm 但不适用于 gam 的自定义 link 缺少什么?
如果你能给我一个提示我应该做什么,在此先感谢你。
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")
}
如您所知,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"
,它会进行性能迭代。
所以,经过上面的概述,我们有两个选择:
- 仍然使用外部迭代,但扩展您自定义的 link 函数;
- 使用性能迭代与
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)
s(x)
我用的是自然三次样条基cr
,至于单变量平滑不需要薄板样条的默认设置。我还设置了一个小的基础维度 k = 3
(对于三次样条曲线不能更小)因为我的玩具数据接近线性并且不需要大的基础维度。更重要的是,这似乎可以防止我的玩具数据集的性能迭代收敛失败。
抱歉,如果答案很明显,但我花了很多时间尝试在 mgcv.gam
中使用自定义 link 函数简而言之,
- 我想使用包 psyphy ( I want to use psyphy.probit_2asym 中的修改概率 link,我称之为
custom_link
) 我可以用这个 link 创建一个 {stats} 系列对象并在 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 都可以工作 link=probit
。
所以我的问题可以概括为:
这个适用于 glm 但不适用于 gam 的自定义 link 缺少什么?
如果你能给我一个提示我应该做什么,在此先感谢你。
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")
}
如您所知,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"
,它会进行性能迭代。
所以,经过上面的概述,我们有两个选择:
- 仍然使用外部迭代,但扩展您自定义的 link 函数;
- 使用性能迭代与
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)
s(x)
我用的是自然三次样条基cr
,至于单变量平滑不需要薄板样条的默认设置。我还设置了一个小的基础维度 k = 3
(对于三次样条曲线不能更小)因为我的玩具数据接近线性并且不需要大的基础维度。更重要的是,这似乎可以防止我的玩具数据集的性能迭代收敛失败。