make.link.gamlss | R Documentation |
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.
make.link.gamlss(link)
show.link(family = "NO")
link |
character or numeric; one of |
family |
a gamlss distribution |
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).
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
For the links involving parameters as in logshifted
and logitshifted
the parameters
can be passed in the definition of the distribution by calling the checklink
function, for example in the definition of
the tau
parameter in BCPE distribution the following call is made:
tstats <- checklink( "tau.link", "Box Cox Power Exponential", substitute(tau.link),
c("logshifted", "log", "identity"), par.link = c(1))
Mikis Stasinopoulos and Bob Rigby
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, \Sexpr[results=rd]{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, \Sexpr[results=rd]{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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1201/b21973")}
(see also https://www.gamlss.com/).
gamlss.family
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.