The Birnbaum-Saunders family - Santos-Neto et al. (2012) (P7 based on the mean and a bounded variance)

The function BS9() defines the Birnbaum-Saunders distribution, a two-parameter distribution, for a gamlss.family object to be used in GAMLSS fitting using the function gamlss().

Usage

BS9(mu.link = "log", sigma.link = "log")

Arguments

mu.link: defines the mu.link, with "log" link as the default for the mu parameter (representing the variance).
sigma.link: defines the sigma.link, with "log" link as the default for the sigma parameter (representing the shape).

Value

Returns a gamlss.family object which can be used to fit a BS9 distribution in the gamlss() function.

Details

The Birnbaum-Saunders distribution with parameters mu and sigma (where mu represents the true variance $\sigma^2$ and sigma represents the shape parameter $\alpha$) has density given by

$$ f(x \mid \mu,\sigma) = \frac{ \exp\left(\frac{1}{2(\sigma-1)}\right) (x\sigma+\mu) }{ 4\sqrt{\pi\sigma(\sigma-1)\mu}\,x^{3/2} } \exp\left[ -\frac{1}{4(\sigma-1)} \left( \frac{x\sigma}{\mu}+ \frac{\mu}{x\sigma} \right) \right] $$

for $x>0$, $\mu>0$ and $\sigma>1$. In this parameterization, $E(X) = \mu$ and $Var(X) = \frac{\mu^2(\sigma-1)(5\sigma-3)}{\sigma^2}$.

References

Santos-Neto, M., Cysneiros, F. J. A., Leiva, V., & Ahmed, S. E. (2012). On new parameterizations of the Birnbaum-Saunders distribution. Pakistan Journal of Statistics, 28(1), 1-26.

Author

David Villegas Ceballos, david.villegas1@udea.edu.co

Examples

# Example 1
# Generating some random values with
# known mu and sigma

set.seed(12345)
y <- rBS9(n=300, mu=2, sigma=2.1)

# Fitting the model using default link function
require(gamlss)
mod1 <- gamlss(y~1, sigma.fo=~1, family=BS9,
               control=gamlss.control(n.cyc=1000))
#> GAMLSS-RS iteration 1: Global Deviance = 1057.745 
#> GAMLSS-RS iteration 2: Global Deviance = 1053.973 
#> GAMLSS-RS iteration 3: Global Deviance = 1053.732 
#> GAMLSS-RS iteration 4: Global Deviance = 1053.713 
#> GAMLSS-RS iteration 5: Global Deviance = 1053.711 
#> GAMLSS-RS iteration 6: Global Deviance = 1053.711 

# Extracting the fitted values for mu and sigma
# using the inverse link function
exp(coef(mod1, what="mu"))
#> (Intercept) 
#>    2.279186 
exp(coef(mod1, what="sigma"))
#> (Intercept) 
#>    2.202631 

# Fitting again the model using our own link function.
# As sigma > 1, we can use
# inv link function:    g-1(eta) = exp(eta) + 1
#     link function:    g(sigma) = log(sigma - 1)

# The following function "log_1_to_inf" is defined with the information
# of the link and inv link functions:

if (FALSE) { # \dontrun{

log_1_to_inf <- function() {
  linkfun <- function(mu) { log(mu - 1) }
  linkinv <- function(eta) { 
    thresh <- log(.Machine$double.xmax)
    eta <- pmin(thresh, pmax(eta, -700))
    1 + exp(eta) 
  }
  mu.eta <- function(eta) { 
    thresh <- log(.Machine$double.xmax)
    eta <- pmin(thresh, pmax(eta, -700))
    pmax(exp(eta), .Machine$double.eps)
  }
  valideta <- function(eta) { TRUE }
  link <- "log_1_to_inf"
  structure(list(linkfun = linkfun, linkinv = linkinv, mu.eta = mu.eta,
                 valideta = valideta, name = link), class = "link-gamlss")
}

mod99 <- gamlss(y~1, sigma.fo=~1, 
                family=BS9(sigma.link = log_1_to_inf()),
                control=gamlss.control(n.cyc=1000))

# Extracting the fitted values for mu and sigma
# using the inverse link function
exp(coef(mod99, what="mu"))
exp(coef(mod99, what="sigma")) + 1

# Example 2
# Generating random values for a regression model

# A function to simulate a data set with Y ~ BS9
gendat <- function(n) {
  x1 <- runif(n)
  x2 <- runif(n)
  mu <- exp(1.6 - 2.4 * x1)          # Aprox 1.5 with link function log(mu)
  sigma <- exp(1.6 + 1.5 * x2) + 1   # Aprox 11.5 with link function log(sigma - 1)
  y <- rBS9(n=n, mu=mu, sigma=sigma)
  data.frame(y=y, x1=x1, x2=x2)
}

set.seed(1234)
dat <- gendat(n=500)

mod2 <- gamlss(y~x1, sigma.fo=~x2, 
               family=BS9(sigma.link = log_1_to_inf()), 
               data=dat,
               control=gamlss.control(n.cyc=100))

summary(mod2)

} # }