在R中的
例如,
此外,为了使用命令
更新和新问题 根据@Ben Bolker的评论,我们需要编写一个名为
我的问题是:如果我们想要将这个链接函数vlog永久地添加到glm中,以便在每个R会话中,我们可以直接使用glm(y~x,family=Gamma(link="log(exp(y)-1)")),那么我们应该使用fix(make.link),然后将vlog的定义添加到其主体中吗?或者fix()只能在当前R会话中执行此操作?再次感谢!
还有一件事:我意识到可能需要修改另一个函数。它是Gamma,定义为
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的评论,我们需要编写一个名为
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会话中,我们可以直接使用glm(y~x,family=Gamma(link="log(exp(y)-1)")),那么我们应该使用fix(make.link),然后将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()
不需要修改(请参见我的答案)。 - Ben Bolkerfamily=Gamma(link=vlog())
... - Ben Bolkervlog
函数。再次感谢你的帮助;-) - alittleboy