【发布时间】：2013-04-02 14:56:03
【问题描述】：
在 R 中的 glm 中，Gamma 系列的默认链接函数是 inverse、identity 和 log。现在对于我的特定问题，我需要使用带有响应Y 的伽玛回归和log(E(Y)-1)) 形式的修改链接函数。因此，我考虑在 R 中修改一些与glm 相关的函数。有几个函数可能是相关的，我正在为以前有此经验的人寻求帮助。
例如，函数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")
}
另外，为了使用命令glm(y ~ log(mu), family = Gamma(link = MyLink))，是否还需要修改glm.fit函数？谢谢！
更新和新问题
根据@Ben Bolker 的cmets，我们需要编写一个名为vlog（实名"log(exp(y)-1)"）的新链接函数。我发现make.link 函数可能是造成这种修改的原因。它被定义为
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")
}
我的问题是：如果我们想永久地添加这个链接函数vlog到glm，这样在每个R会话中，我们可以使用@987654339直接@，我们是不是用fix(make.link)，然后在它的body里面加上vlog的定义？或者fix() 只能在当前的 R 会话中执行此操作？再次感谢！
还有一件事：我意识到可能需要修改另一个函数。是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")
}
我觉得我们也需要修改
okLinks <- c("inverse", "log", "identity")
到
okLinks <- c("inverse", "log", "identity", "log(exp(y)-1)")
?
【问题讨论】：
-
我不明白所有这些额外的复杂性。我展示了下面的示例，只要定义了
vlog，就可以通过glm(...,family=Gamma(link=vlog())拟合备用链接模型。您可以将vlog放在.R文件中，并将source()放在每个会话中，或者创建一个定义函数的小包。如果您愿意，您也可以将它放在您的 R 配置文件中，但在您将要使用它的每个 R 脚本中，它可能对source("vlog.R")更透明。我认为Gamma()不需要修改（再次，请参阅我的回答）。 -
我想如果你坚持按名称调用链接函数，你将不得不做你上面描述的所有额外的黑客攻击，但我看不出@有什么问题987654354@ ...
-
@BenBolker：是的，我尝试了您的代码，它们运行良好！也许我的额外问题是关于
fixing 一个 R 函数以永久包含用户定义的选项。我将在我的包中包含vlog函数。再次感谢您的帮助 ;-) -
我会说您应该从 R 源代码中复制该函数（以便您获得包含任何相关的 cmets）并将其合并到您加载的包中，这将掩盖基本版本。这是一项完全不同的任务，您可能应该将其作为一个单独的问题提出。
-
@BenBolker：是的——我将把它作为一个单独的问题发布；-)