modify glm function to adopt user-specified link function in R

Try gnlm::gnlr(). Using x, y, sh from Ben Bolker's example:

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

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

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

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