# R/plinks.R In glmx: Generalized Linear Models Extended

#### Documented in angularao1ao2foldexpgjgossetloglognblogitpregibonrocketalpha

```## angular transformation (aka arcsine transformation)

angular <- function(verbose = FALSE)
{
linkfun <- function(mu) asin(sqrt(mu))
linkinv <- function(eta) {
etastar <- pmin(asin(1) - .Machine\$double.eps, pmax(.Machine\$double.eps, eta))
if(verbose && !isTRUE(all.equal(as.vector(eta), as.vector(etastar))))
warning("truncation in inverse link function")
return((sin(etastar))^2)
}
mu.eta <- function(eta) 2 * cos(eta) * sin(eta)
valideta <- function(eta) {
if(verbose && (!all(eta <= asin(1) & 0 <= eta)))
warning("some of the current etas are out of range")
TRUE
}
name <- "angular"

## return link-glm object
valideta = valideta, name = name), class = "link-glm")
}

## just for completeness: log-log link in addition to complementary log-log

loglog <- function() {
structure(list(
linkfun = function(mu) -log(-log(mu)),
linkinv = function(eta) pmax(pmin(exp(-exp(-eta)), 1 - .Machine\$double.eps), .Machine\$double.eps),
mu.eta = function(eta) {
eta <- pmin(eta, 700)
pmax(exp(-eta - exp(-eta)), .Machine\$double.eps)
},
valideta = function(eta) TRUE,
name = "loglog"),
}

## transformations from Aranda-Ordaz (1981)
## ao1 is the symmetric link from Section 2.1

ao1 <- function(phi, verbose = FALSE)
{
## parameter processing
if(verbose && phi < 0) warning("sign of phi ignored as ao1(phi) = ao1(-phi)")
phi <- abs(phi)
if(phi == 0) {
rval\$name <- "ao1"
return(rval)
}

linkfun <- function(mu) (2/phi) * (mu^phi - (1 - mu)^phi)/(mu^phi + (1 - mu)^phi)

linkinv <- function(eta) {
etastar <- pmin(2/phi - .Machine\$double.eps, pmax(-2/phi + .Machine\$double.eps, eta))
if(verbose && !isTRUE(all.equal(as.vector(eta), as.vector(etastar))))
warning("truncation in inverse link function")
((1 + 0.5 * phi * etastar)^(1/phi))/((1 + 0.5 * phi * etastar)^(1/phi) + (1 - 0.5 * phi * etastar)^(1/phi))
}

mu.eta <- function(eta) {
phieta1 <- 1 + 0.5 * phi * eta
phieta2 <- 1 - 0.5 * phi * eta
phieta1^((1/phi) - 1) * 0.5/(phieta1^(1/phi) + phieta2^(1/phi)) - (phieta1^(1/phi)) *
(phieta1^((1/phi) - 1) * 0.5 - phieta2^((1/phi) - 1) * (0.5))/(phieta1^(1/phi) + phieta2^(1/phi))^2
}

valideta <- function(eta) {
if(verbose && !all(abs(0.5 * phi * eta) < 1)) warning("some of the current etas are out of range")
TRUE
}

name <- "ao1" ## "Aranda-Ordaz symmetric"

## return link-glm object
valideta = valideta, name = name), class = "link-glm")
}

## transformations from Aranda-Ordaz (1981)
## ao2 is the assymmetric link from Section 2.2

ao2 <- function(phi, verbose = FALSE)
{
## parameter processing
if(phi == 1) {
rval\$name <- "ao2"
return(rval)
}
if(phi == 0) {
rval\$name <- "ao2"
return(rval)
}

linkfun <- function(mu) log(((1 - mu)^(-phi) - 1)/phi)

if(phi < 0) {
etastar <- pmin(eta, log(-1/phi) - .Machine\$double.eps)
if(verbose && !isTRUE(all.equal(as.vector(eta), as.vector(etastar))))
warning("truncation in inverse link function")
eta <- etastar
}
1- (1 + phi * exp(eta))^(-1/phi)
}

mu.eta <- function(eta) exp(eta) * (1 + phi * exp(eta))^(-(1 + phi)/phi)

valideta <- function(eta) {
if(verbose && !all(phi * exp(eta) > -1)) warning("some of the current etas are out of range")
TRUE
}

name <- "ao2" ## "Aranda-Ordaz asymmetric"

## return link-glm object
valideta = valideta, name = name), class = "link-glm")
}

## folded exponential transformation (Piepho 2003)

foldexp <- function(phi, verbose = FALSE)
{
## parameter processing
if(phi == 0) {
structure(list(
linkfun = function(mu) mu - 0.5,
linkinv = function(eta) {
etastar <- pmin(0.5 - .Machine\$double.eps, pmax(-0.5 + .Machine\$double.eps, eta))
if(verbose && !isTRUE(all.equal(as.vector(eta), as.vector(etastar))))
warning("truncation in inverse link function")
etastar + 0.5
},
mu.eta = function(eta) rep(1, length(eta)),
valideta = function(eta) TRUE,
name = "foldexp"),
}

upper <- (exp(phi) - 1)/(2 * phi)
lower <- (1 - exp(phi))/(2 * phi)

linkfun <- function(mu) (exp(phi * mu) - exp(phi * (1 - mu)))/(2 * phi)

dfoldexp <- function(mu) (exp(phi * mu) + exp(phi * (1 - mu)))/2

linkinv <- function(eta) {
etastar <- pmin(upper - .Machine\$double.eps, pmax(lower + .Machine\$double.eps, eta))
if(verbose && !isTRUE(all.equal(as.vector(eta), as.vector(etastar))))
warning("truncation in inverse link function")
0.5 + asinh((phi * etastar)/exp(phi/2))/phi
}

mu.eta <- function(eta) 1/dfoldexp(linkinv(eta))

valideta <- function(eta) TRUE

name <- "foldexp" ## "folded exponential"

## return link-glm object
valideta = valideta, name = name), class = "link-glm")
}

## Guerrero-Johnson (1982) transformation

gj <- function(phi, verbose = FALSE)
{
## parameter processing
if(phi == 0) {
rval\$name <- "gj"
return(rval)
}

linkfun <- function(mu) ((mu/(1 - mu))^phi - 1)/phi

etastar <- pmax((-sign(phi)/phi) + .Machine\$double.eps, eta)
if(verbose && !isTRUE(all.equal(as.vector(eta), as.vector(etastar))))
warning("truncation in inverse link function")
((phi * etastar + 1)^(1/phi))/(1 + ((phi * etastar + 1)^(1/phi)))
}

mu.eta <- function(eta) {
pe1 <- (phi * eta + 1)^((1/phi) - 1)
pe2 <- (phi * eta + 1)^(1/phi)
pe1/(1 + pe2) - pe2 * (pe1/(1 + pe2)^2)
}

valideta <- function(eta) TRUE
## old version:
##  all(-sign(phi)/phi < eta)

name <- "gj" ## "Guerrero-Johnson"

## return link-glm object
valideta = valideta, name = name), class = "link-glm")
}

## t_alpha link function

talpha <- function(alpha, verbose = FALSE,
splineinv = TRUE, eps = 2 * .Machine\$double.eps, maxit = 100)
{
## parameter processing
if(alpha > 2 | alpha < 0) stop("alpha must be between 0 and 2")
if(alpha == 1) {
rval\$name <- "talpha"
return(rval)
}

linkfun <- function(mu) alpha * log(mu) - (2 - alpha) * log(1 - mu)

dtrafo <- function(alpha, x) alpha/x + (2 - alpha)/(1 - x)

linkinv <- function(eta) {
## inverse link is easy to compute for alpha 0, 1, 2
## alpha == 1 already special-cased above
if(alpha == 0){
etastar <- pmax(eta, .Machine\$double.eps)
if(verbose && !isTRUE(all.equal(as.vector(eta), as.vector(etastar))))
warning("truncation in inverse link function")
return(1 - exp(-etastar/2))
}
if(alpha == 2) {
etastar <- pmin(eta, -.Machine\$double.eps)
if(verbose && !isTRUE(all.equal(as.vector(eta), as.vector(etastar))))
warning("truncation in inverse link function")
return(exp(etastar/2))
}

## otherwise some numeric technique is required:
## either via splines (quick & dirty) or via Newton (cleaner but somewhat slower)
if(splineinv) {

## this is somewhat less precise than the Newton solution
## but faster for moderately to large data sets
mu0 <- plogis(seq(-30, 30, by = 0.01))
mu <- splinefun(eta0, mu0)(eta)

} else {

## simple newton algorithm (maxit is rarely reached, only when alpha close to 0 or 2)
## 8 times faster than uniroot
mu <- plogis(eta)
ok <- FALSE
n <- 0
while(!ok) {
mu.prev <- mu
mu <- pmax(pmin(mu - (linkfun(mu) - eta)/(dtrafo(alpha, mu) - 1), 1 - eps), eps)
ok <- max(abs(mu.prev - mu)) < eps
if(n > maxit) {
warning("maxit reached")
ok <- TRUE
}
n <- n + 1
}

}
mu
}

mu.eta <- function(eta) 1/dtrafo(alpha, linkinv(eta))

valideta <- function(eta) TRUE
## old version:
##    if(alpha == 0){  if(any(eta <= 0)){FALSE}
##            if(any(eta > 0)){TRUE}}
##    if(alpha == 2){  if(any(eta >= 0)){FALSE}
##            if(any(eta < 0)){TRUE}}

name <- "talpha"

## return link-glm object
valideta = valideta, name = name), class = "link-glm")
}

## Rocke (1993) beta transformation as link function

rocke <- function(shape1, shape2, verbose = FALSE)
{
## parameter processing
if(missing(shape2)) {
shape2 <- rep(shape1, length.out = 2)[2]
shape1 <- shape1[1]
}
if(shape1 == 0 & shape2 == 0) {
rval\$name <- "rocke"
return(rval)
}
if(shape1 < 0 | shape2 < 0) stop("parameters of Rocke transformation need to be both strictly positive or both equal to 0")
# what to do if shape1 == 0, and shape2 > 0? nothing: then the integral does not exist!

shift <- pbeta(0.5, shape1, shape2)
betafactor <- beta(shape1, shape2)

linkfun <- function(mu) betafactor * (pbeta(mu, shape1, shape2) - shift)

etastar <- pmin((1-shift)*betafactor, pmax(eta, -shift*betafactor))
if(verbose && !isTRUE(all.equal(as.vector(eta), as.vector(etastar))))
warning("truncation in inverse link function")
pmin(1 - .Machine\$double.eps, pmax(.Machine\$double.eps,
qbeta(etastar/betafactor + shift, shape1, shape2)))
}

mu.eta <- function(eta)
1/(betafactor * (dbeta(linkinv(eta), shape1, shape2)))

valideta <- function(eta) {
if(verbose && any( eta < -shift * betafactor | eta > (1 - shift) * betafactor ))
warning("some of the current etas are out of range")
TRUE
}

name <- "rocke" ## "Rocke's beta"

## return link-glm object
valideta = valideta, name = name), class = "link-glm")
}

## some form of generalized logistic link
## constructed as right-censored negative
## binomial distribution with parameter theta

nblogit <- function(theta) {
## parameter processing
if(theta <= 0) stop("theta must be positive")
if(theta == 1) {
rval\$name <- "nblogit"
return(rval)
}
if(!is.finite(theta)) {
rval\$name <- "nblogit"
return(rval)
}

linkfun <- function(mu) log(1 - (1 - mu)^(1/theta)) - log(1 - mu)/theta + log(theta)
linkinv <- function(eta) 1 - (1 / (1 + exp(eta)/theta))^theta
mu.eta <- function(eta) exp( -(theta + 1) * log(1 + exp(eta)/theta) + eta )
valideta <- function(eta) TRUE
name <- "nblogit"

## return link-glm object
valideta = valideta, name = name), class = "link-glm")
}

gosset <- function(nu, verbose = FALSE)
{
if(nu <= 0) stop("nu must be positive")
if(verbose && nu < 0.2) warning("gosset link might be unreliable for nu < 0.2")
qqt <- function(p, nu) sign(p - 0.5) * sqrt(qf(1 - 2 * pmin(p, 1 - p), 1, nu))

linkfun <- function(mu) qqt(mu, nu)
linkinv <- function(eta) {
thresh <- -qqt(.Machine\$double.eps, nu)
eta <- pmin(thresh, pmax(eta, -thresh))
pt(eta, nu)
}
mu.eta <- function(eta) pmax(dt(eta, nu), .Machine\$double.eps)
valideta <- function(eta) TRUE
name <- "gosset"

## return link-glm object
valideta = valideta, name = name), class = "link-glm")
}

pregibon <- function(a, b) {

## allow specification via vector as first argument
if(missing(b)) {
b <- rep(a, length.out = 2)[2]
a <- a[1]
}

linkfun <- function(mu) -qpregibon(1 - mu, a = a, b = b)

linkinv <- function(eta) {
eps <- .Machine\$double.eps^0.5
tlo <- qpregibon(eps, a = a, b = b)
thi <- qpregibon(1 - eps, a = a, b = b)
eta <- -pmin(thi, pmax(-eta, tlo))
1 - ppregibon(-eta, a = a, b = b)
}

mu.eta <- function(eta) pmax(dpregibon(-eta, a = a, b = b), .Machine\$double.eps^0.5)
valideta <- function(eta) TRUE
name <- "pregibon"

## return link-glm object