修改glm函数以在R中采用用户指定的链接函数 [英] modify glm function to adopt user-specified link function in R

问题描述

在R的glm中，Gamma系列的默认链接功能是inverse，identity和log.现在，对于我的特定问题，我需要对响应Y使用伽马回归，并使用log(E(Y)-1))形式的修改后的链接函数.因此，我考虑在R中修改一些与glm相关的功能.可能有几个相关的功能，我正在寻求经验的人寻求帮助.

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.

例如，功能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功能?谢谢！

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!

更新和新问题

根据@Ben Bolker的评论，我们需要编写一个名为vlog(真实名称为"log(exp(y)-1)")的新链接函数.我发现make.link函数可能负责这种修改.定义为

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

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会话中，我们都可以使用直接，我们是否应该使用fix(make.link)，然后将vlog的定义添加到其主体中?还是fix()只能在当前的R会话中做到这一点?再次感谢！

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!

还有一件事:我意识到也许需要修改另一个功能.它是Gamma，定义为

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

我认为我们也需要修改

okLinks <- c("inverse", "log", "identity")

到

okLinks <- c("inverse", "log", "identity", "log(exp(y)-1)")

?

推荐答案

我基本上遵循?family中示例的形式，该示例显示了用户指定的qlogis(mu^(1/days))形式的链接.

I'm basically following the form of the example in ?family which shows a user-specified link of the form qlogis(mu^(1/days)).

我们想要一个格式为eta = log(exp(y)-1)的链接(因此，反向链接为y=log(exp(eta)+1)和mu.eta = dy/d(eta) = 1/(1+exp(-eta))

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

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

基本检查:

vv <- vlog()
vv$linkfun(vv$linkinv(27))  ## check invertibility
library("numDeriv")
all.equal(grad(vv$linkinv,2),vv$mu.eta(2))  ## check derivative

示例:

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 
## 

这篇关于修改glm函数以在R中采用用户指定的链接函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！