make.link.gamlss: Create a Link for GAMLSS families

Description

The function make.link.gamlss() is used with gamlss.family distributions in package gamlss(). Given a link, it returns a link function, an inverse link function, the derivative dpar/deta where 'par' is the appropriate distribution parameter and a function for checking the domain. It differs from the usual make.link of glm() by having extra links as the logshifto1, and the own. For the use of the own link see the example bellow. show.link provides a way in which the user can identify the link functions available for each gamlss distribution. If your required link function is not available for any of the gamlss distributions you can add it in.

Usage

make.link.gamlss(link)
show.link(family = "NO")

Value

For the make.link.gamlss a list with components

linkfun: Link function function(parameter)

linkinv: Inverse link function function(eta)

mu.eta: Derivative function(eta) dparameter/deta

valideta: function(eta) TRUE if all of eta is in the domain of linkinv.

For the show.link a list with components the available links for the distribution parameters

Arguments

link: character or numeric; one of "logit", "probit", "cloglog", "identity", "log", "sqrt", "1/mu^2", "inverse", "logshifted", "logitshifted", or number, say lambda resulting in power link \(\mu^\lambda\).
family: a gamlss distribution family

Author

Mikis Stasinopoulos and Bob Rigby

Details

The own link function is added to allow the user greater flexibility. In order to used the own link function for any of the parameters of the distribution the own link should appear in the available links for this parameter. You can check this using the function show.link. If the own do not appear in the list you can create a new function for the distribution in which own is added in the list. For example the first line of the code of the binomial distribution, BI, has change from

"mstats <- checklink("mu.link", "Binomial", substitute(mu.link), c("logit", "probit", "cloglog", "log")), in version 1.0-0 of gamlss, to

"mstats <- checklink("mu.link", "Binomial", substitute(mu.link), c("logit", "probit", "cloglog", "log", "own"))

in version 1.0-1. Given that the parameter has own as an option the user needs also to define the following four new functions in order to used an own link.

i) own.linkfun

ii) own.linkinv

iii) own.mu.eta and

iv) own.valideta.

An example is given below.

Only one parameter of the distribution at a time is allowed to have its own link, (unless the same four own functions above are suitable for more that one parameter of the distribution).

Note that from gamlss version 1.9-0 the user can introduce its own link function by define an appropriate function, (see the example below).

References

Rigby, R. A. and Stasinopoulos D. M. (2005). Generalized additive models for location, scale and shape,(with discussion), Appl. Statist., 54, part 3, pp 507-554.

Rigby, R. A., Stasinopoulos, D. M., Heller, G. Z., and De Bastiani, F. (2019) Distributions for modeling location, scale, and shape: Using GAMLSS in R, Chapman and Hall/CRC, tools:::Rd_expr_doi("10.1201/9780429298547"). An older version can be found in https://www.gamlss.com/.

Stasinopoulos D. M. Rigby R.A. (2007) Generalized additive models for location scale and shape (GAMLSS) in R. Journal of Statistical Software, Vol. 23, Issue 7, Dec 2007, tools:::Rd_expr_doi("10.18637/jss.v023.i07").

Stasinopoulos D. M., Rigby R.A., Heller G., Voudouris V., and De Bastiani F., (2017) Flexible Regression and Smoothing: Using GAMLSS in R, Chapman and Hall/CRC. tools:::Rd_expr_doi("10.1201/b21973")

(see also https://www.gamlss.com/).

Examples

Run this code

str(make.link.gamlss("logshiftto1"))
l2<-make.link.gamlss("logshiftto1")
l2$linkfun(2) # should close to zero (Note that 0.00001 is added)
l2$linkfun(1-0.00001) # should be -Inf but it is large negative
#---------------------------------------------------------
# now use the own link function
# first if the distribution allows you
show.link(BI)
# seems OK now define the four own functions
# First try the probit link using the own link function
# 1: the linkfun function
own.linkfun <- function(mu) { qNO(p=mu)}
# 2: the inverse link function 
own.linkinv <- function(eta) { 
              thresh <- -qNO(.Machine$double.eps)
               eta <- pmin(thresh, pmax(eta, -thresh))
              pNO(eta)}
# 3: the dmu/deta function 
own.mu.eta <- function(eta) pmax(dNO(eta), .Machine$double.eps)
# 4: the valideta function 
own.valideta <- function(eta) TRUE

## bring the data
# library(gamlss) 
#data(aep)
# fitting the model using "own"   
# h1<-gamlss(y~ward+loglos+year, family=BI(mu.link="own"), data=aep)  
# model h1 should be identical to the probit 
# h2<-gamlss(y~ward+loglos+year, family=BI(mu.link="probit"), data=aep)
# now using a function instead of "own" 
probittest <- function()
{
linkfun <- function(mu) { qNO(p=mu)}
linkinv <- function(eta) 
            { 
              thresh <- -qNO(.Machine$double.eps)
               eta <- pmin(thresh, pmax(eta, -thresh))
              pNO(eta)
            }
mu.eta <- function(eta) pmax(dNO(eta), .Machine$double.eps) 
valideta <- function(eta) TRUE
link <- "probitTest"
structure(list(linkfun = linkfun, linkinv = linkinv, mu.eta = mu.eta, 
        valideta = valideta, name = link), class = "link-gamlss")
}
# h3<-gamlss(y~ward+loglos+year, family=BI(mu.link=probittest()), data=aep)  
# Second try the complementary log-log 
# using the Gumbel distribution  
own.linkfun <- function(mu) { qGU(p=mu)} 
own.linkinv <- function(eta) { 
              thresh <- -qGU(.Machine$double.eps)
               eta <- pmin(thresh, pmax(eta, -thresh))
              pGU(eta)} 
own.mu.eta <- function(eta) pmax(dGU(eta), .Machine$double.eps)
own.valideta <- function(eta) TRUE
# h1 and h2 should be identical to cloglog
# h1<-gamlss(y~ward+loglos+year, family=BI(mu.link="own"), data=aep)  
# h2<-gamlss(y~ward+loglos+year, family=BI(mu.link="cloglog"), data=aep)
# note that the Gumbel distribution is negatively skew
# for a positively skew link function we can used the Reverse Gumbel 
revloglog  <- function()
{
linkfun <- function(mu) { qRG(p=mu)} 
linkinv <- function(eta) { 
              thresh <- -qRG(.Machine$double.eps)
               eta <- pmin(thresh, pmax(eta, -thresh))
              pRG(eta)}
mu.eta <- function(eta) pmax(dRG(eta), .Machine$double.eps)
valideta <- function(eta) TRUE
link <- "revloglog"
structure(list(linkfun = linkfun, linkinv = linkinv, mu.eta = mu.eta, 
        valideta = valideta, name = link), class = "link-gamlss")
}
# h1<-gamlss(y~ward+loglos+year, family=BI(mu.link=revloglog()), data=aep)  
# a considerable improvement in the deviance
# try a shifted logit link function from -1, 1 
own.linkfun <- function(mu)
             { shift = c(-1,1)           
               log((mu-shift[1])/(shift[2]-mu))
             }
own.linkinv <- function(eta) 
            {
            shift = c(-1,1)  
            thresh <- -log(.Machine$double.eps)
               eta <- pmin(thresh, pmax(eta, -thresh))
                      shift[2]-(shift[2]-shift[1])/(1 + exp(eta))
            } 
own.mu.eta <- function(eta) 
            {
        shift = c(-1,1)  
            thresh <- -log(.Machine$double.eps)
               res <- rep(.Machine$double.eps, length(eta))
            res[abs(eta) < thresh] <- ((shift[2]-shift[1])*exp(eta)/(1 + 
                                 exp(eta))^2)[abs(eta) < thresh]
            res
            }
own.valideta <- function(eta) TRUE       
#----------
str(make.link.gamlss("own"))
l2<-make.link.gamlss("own")
l2$linkfun(0) # should be zero
l2$linkfun(1) # should be Inf
l2$linkinv(-5:5)

Run the code above in your browser using DataLab