In glm
in R, the default link functions for the Gamma
family are inverse
,identity
and log
. Now for my particular question, I need to use gamma regression with response Y
and a modified link function in the form of log(E(Y)-1))
. Thus, I consider modifying some glm
-related functions in R. There are several functions that may be relevant, and I am seeking help for anyone who had previous experience in doing this.
在glm In R中,Gamma家族的默认链接函数为逆、identity和log。现在,对于我的特殊问题,我需要使用响应Y和一个修改过的链接函数,以log(E(Y)-1)的形式使用gamma回归。因此,我考虑在r中修改一些与glm相关的功能,有一些可能是相关的功能,我正在寻求帮助那些曾经有过这样做的人。
For example, the functions Gamma
is defined as
例如,函数Gamma被定义为。
function (link = "inverse")
{
linktemp <- substitute(link)
if (!is.character(linktemp))
linktemp <- deparse(linktemp)
okLinks <- c("inverse", "log", "identity")
if (linktemp %in% okLinks)
stats <- make.link(linktemp)
else if (is.character(link))
stats <- make.link(link)
else {
if (inherits(link, "link-glm")) {
stats <- link
if (!is.null(stats$name))
linktemp <- stats$name
}
else {
stop(gettextf("link \"%s\" not available for gamma family; available links are %s",
linktemp, paste(sQuote(okLinks), collapse = ", ")),
domain = NA)
}
}
variance <- function(mu) mu^2
validmu <- function(mu) all(mu > 0)
dev.resids <- function(y, mu, wt) -2 * wt * (log(ifelse(y ==
0, 1, y/mu)) - (y - mu)/mu)
aic <- function(y, n, mu, wt, dev) {
n <- sum(wt)
disp <- dev/n
-2 * sum(dgamma(y, 1/disp, scale = mu * disp, log = TRUE) *
wt) + 2
}
initialize <- expression({
if (any(y <= 0)) stop("non-positive values not allowed for the 'gamma' family")
n <- rep.int(1, nobs)
mustart <- y
})
simfun <- function(object, nsim) {
wts <- object$prior.weights
if (any(wts != 1))
message("using weights as shape parameters")
ftd <- fitted(object)
shape <- MASS::gamma.shape(object)$alpha * wts
rgamma(nsim * length(ftd), shape = shape, rate = shape/ftd)
}
structure(list(family = "Gamma", link = linktemp, linkfun = stats$linkfun,
linkinv = stats$linkinv, variance = variance, dev.resids = dev.resids,
aic = aic, mu.eta = stats$mu.eta, initialize = initialize,
validmu = validmu, valideta = stats$valideta, simulate = simfun),
class = "family")
}
Also, in order to use the command glm(y ~ log(mu), family = Gamma(link = MyLink))
, do I also need to modify the glm.fit
function? Thank you!
另外,为了使用命令glm(y ~ log(mu), family = Gamma(link = MyLink)),我也需要修改glm。适应函数?谢谢你们!
Updates and New Question
更新和新问题
According to @Ben Bolker's comments, we need to write a new link function called vlog
(with real name "log(exp(y)-1)"
). I find that the make.link
function might be responsible for such a modification. It is defined as
根据@Ben Bolker的评论,我们需要编写一个名为vlog的新链接函数(使用实名“log(exp(y)-1)”)。我发现制作。链接函数可能负责这样的修改。它被定义为
function (link)
{
switch(link, logit = {
linkfun <- function(mu) .Call(C_logit_link, mu)
linkinv <- function(eta) .Call(C_logit_linkinv, eta)
mu.eta <- function(eta) .Call(C_logit_mu_eta, eta)
valideta <- function(eta) TRUE
},
...
}, log = {
linkfun <- function(mu) log(mu)
linkinv <- function(eta) pmax(exp(eta), .Machine$double.eps)
mu.eta <- function(eta) pmax(exp(eta), .Machine$double.eps)
valideta <- function(eta) TRUE
},
...
structure(list(linkfun = linkfun, linkinv = linkinv, mu.eta = mu.eta,
valideta = valideta, name = link), class = "link-glm")
}
My question is: if we want to permanently add this link function vlog
to glm
, so that in each R session, we can use glm(y~x,family=Gamma(link="log(exp(y)-1)"))
directly, shall we use the fix(make.link)
and then add the definition of vlog
to its body? Or fix()
can only do that in current R session? Thanks again!
我的问题是:如果我们想永久地将这个链接函数vlog添加到glm,那么在每个R会话中,我们可以直接使用glm(y~x,family=Gamma(link="log(exp(y)-1))),我们是否应该使用fix(make.link),然后将vlog的定义添加到它的body中?或者fix()只能在当前的R会话中实现吗?再次感谢!
One more thing: I realize that maybe another function needs to be modified. It is Gamma
, defined as
还有一件事:我意识到可能需要修改另一个函数。它是,定义为。
function (link = "inverse")
{
linktemp <- substitute(link)
if (!is.character(linktemp))
linktemp <- deparse(linktemp)
okLinks <- c("inverse", "log", "identity")
if (linktemp %in% okLinks)
stats <- make.link(linktemp)
else if (is.character(link))
stats <- make.link(link)
else {
if (inherits(link, "link-glm")) {
stats <- link
if (!is.null(stats$name))
linktemp <- stats$name
}
else {
stop(gettextf("link \"%s\" not available for gamma family; available links are %s",
linktemp, paste(sQuote(okLinks), collapse = ", ")),
domain = NA)
}
}
variance <- function(mu) mu^2
validmu <- function(mu) all(mu > 0)
dev.resids <- function(y, mu, wt) -2 * wt * (log(ifelse(y ==
0, 1, y/mu)) - (y - mu)/mu)
aic <- function(y, n, mu, wt, dev) {
n <- sum(wt)
disp <- dev/n
-2 * sum(dgamma(y, 1/disp, scale = mu * disp, log = TRUE) *
wt) + 2
}
initialize <- expression({
if (any(y <= 0)) stop("non-positive values not allowed for the 'gamma' family")
n <- rep.int(1, nobs)
mustart <- y
})
simfun <- function(object, nsim) {
wts <- object$prior.weights
if (any(wts != 1))
message("using weights as shape parameters")
ftd <- fitted(object)
shape <- MASS::gamma.shape(object)$alpha * wts
rgamma(nsim * length(ftd), shape = shape, rate = shape/ftd)
}
structure(list(family = "Gamma", link = linktemp, linkfun = stats$linkfun,
linkinv = stats$linkinv, variance = variance, dev.resids = dev.resids,
aic = aic, mu.eta = stats$mu.eta, initialize = initialize,
validmu = validmu, valideta = stats$valideta, simulate = simfun),
class = "family")
}
I think we also need to revise
我想我们还需要修改。
okLinks <- c("inverse", "log", "identity")
to
来
okLinks <- c("inverse", "log", "identity", "log(exp(y)-1)")
?
吗?
2 个解决方案
#1
12
I'm basically following the form of the example in ?family
which shows a user-specified link of the form qlogis(mu^(1/days))
.
我基本上遵循中的示例的形式?家庭显示用户指定链接的形式qlogis(μ^(1 /天))。
We want a link of the form eta = log(exp(y)-1)
(so the inverse link is y=log(exp(eta)+1)
, and mu.eta = dy/d(eta) = 1/(1+exp(-eta))
我们想要的是form =log(exp(y)-1)的一个链接,所以反向链接是y=log(exp(eta)+1)和mu。= dy/d(eta) = 1/(1+exp(-eta))
vlog <- function() {
## link
linkfun <- function(y) log(exp(y)-1)
## inverse link
linkinv <- function(eta) log(exp(eta)+1)
## derivative of invlink wrt eta
mu.eta <- function(eta) { 1/(exp(-eta) + 1) }
valideta <- function(eta) TRUE
link <- "log(exp(y)-1)"
structure(list(linkfun = linkfun, linkinv = linkinv,
mu.eta = mu.eta, valideta = valideta,
name = link),
class = "link-glm")
}
Basic checks:
基本的检查:
vv <- vlog()
vv$linkfun(vv$linkinv(27)) ## check invertibility
library("numDeriv")
all.equal(grad(vv$linkinv,2),vv$mu.eta(2)) ## check derivative
Example:
例子:
set.seed(101)
n <- 1000
x <- runif(n)
sh <- 2
y <- rgamma(n,scale=vv$linkinv(2+3*x)/sh,shape=sh)
glm(y~x,family=Gamma(link=vv))
##
## Call: glm(formula = y ~ x, family = Gamma(link = vv))
##
## Coefficients:
## (Intercept) x
## 1.956 3.083
##
## Degrees of Freedom: 999 Total (i.e. Null); 998 Residual
## Null Deviance: 642.2
## Residual Deviance: 581.8 AIC: 4268
##
#2
2
Try gnlm::gnlr()
. Using x
, y
, sh
from Ben Bolker's example:
尝试gnlm:gnlr()。用Ben Bolker的例子,用x, y, sh:
library(gnlm)
# custom link / inverse
custom_inv <- function(eta) log(exp(eta)+1)
library(gnlm)
gnlr(y=y,
distribution = "gamma",
mu = ~ custom_inv(beta0 + beta1*x),
pmu = list(beta0=0, beta1=0),
pshape=sh
)
# Location parameters:
# estimate se
# beta0 1.956 0.1334
# beta1 3.083 0.2919
#
# Shape parameters:
# estimate se
# p[1] 0.625 0.04133
#1
12
I'm basically following the form of the example in ?family
which shows a user-specified link of the form qlogis(mu^(1/days))
.
我基本上遵循中的示例的形式?家庭显示用户指定链接的形式qlogis(μ^(1 /天))。
We want a link of the form eta = log(exp(y)-1)
(so the inverse link is y=log(exp(eta)+1)
, and mu.eta = dy/d(eta) = 1/(1+exp(-eta))
我们想要的是form =log(exp(y)-1)的一个链接,所以反向链接是y=log(exp(eta)+1)和mu。= dy/d(eta) = 1/(1+exp(-eta))
vlog <- function() {
## link
linkfun <- function(y) log(exp(y)-1)
## inverse link
linkinv <- function(eta) log(exp(eta)+1)
## derivative of invlink wrt eta
mu.eta <- function(eta) { 1/(exp(-eta) + 1) }
valideta <- function(eta) TRUE
link <- "log(exp(y)-1)"
structure(list(linkfun = linkfun, linkinv = linkinv,
mu.eta = mu.eta, valideta = valideta,
name = link),
class = "link-glm")
}
Basic checks:
基本的检查:
vv <- vlog()
vv$linkfun(vv$linkinv(27)) ## check invertibility
library("numDeriv")
all.equal(grad(vv$linkinv,2),vv$mu.eta(2)) ## check derivative
Example:
例子:
set.seed(101)
n <- 1000
x <- runif(n)
sh <- 2
y <- rgamma(n,scale=vv$linkinv(2+3*x)/sh,shape=sh)
glm(y~x,family=Gamma(link=vv))
##
## Call: glm(formula = y ~ x, family = Gamma(link = vv))
##
## Coefficients:
## (Intercept) x
## 1.956 3.083
##
## Degrees of Freedom: 999 Total (i.e. Null); 998 Residual
## Null Deviance: 642.2
## Residual Deviance: 581.8 AIC: 4268
##
#2
2
Try gnlm::gnlr()
. Using x
, y
, sh
from Ben Bolker's example:
尝试gnlm:gnlr()。用Ben Bolker的例子,用x, y, sh:
library(gnlm)
# custom link / inverse
custom_inv <- function(eta) log(exp(eta)+1)
library(gnlm)
gnlr(y=y,
distribution = "gamma",
mu = ~ custom_inv(beta0 + beta1*x),
pmu = list(beta0=0, beta1=0),
pshape=sh
)
# Location parameters:
# estimate se
# beta0 1.956 0.1334
# beta1 3.083 0.2919
#
# Shape parameters:
# estimate se
# p[1] 0.625 0.04133