The Discrete Perks family

The function DPERKS() defines the Discrete Perks distribution, a two parameter distribution, for a gamlss.family object to be used in GAMLSS fitting using the function gamlss().

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

Arguments

mu.link: defines the mu.link, with "log" link as the default for the mu parameter.
sigma.link: defines the sigma.link, with "log" link as the default for the sigma.

Details

The discrete Perks distribution with parameters \(\mu > 0\) and \(\sigma > 0\) has a support 0, 1, 2, ... and density given by

\(f(x | \mu, \sigma) = \frac{\mu(1+\mu)(e^\sigma-1)e^{\sigma x}}{(1+\mu e^{\sigma x})(1+\mu e^{\sigma(x+1)})} \)

Note: in this implementation we changed the original parameters \(\lambda\) for \(\mu\) and \(\beta\) for \(\sigma\), we did it to implement this distribution within gamlss framework.

References

Tyagi, A., Choudhary, N., & Singh, B. (2020). A new discrete distribution: Theory and applications to discrete failure lifetime and count data. J. Appl. Probab. Statist, 15, 117-143.

Author

Veronica Seguro Varela, vseguro@unal.edu.co

Examples

# Example 1
# Generating some random values with
# known mu and sigma
set.seed(123)
y <- rDPERKS(n=1000, mu=2.5, sigma=0.4)

# Fitting the model
library(gamlss)
mod1 <- gamlss(y~1, family=DPERKS,
               control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 4443.326 
#> GAMLSS-RS iteration 2: Global Deviance = 4269.683 
#> GAMLSS-RS iteration 3: Global Deviance = 4197.172 
#> GAMLSS-RS iteration 4: Global Deviance = 4164.998 
#> GAMLSS-RS iteration 5: Global Deviance = 4149.949 
#> GAMLSS-RS iteration 6: Global Deviance = 4142.595 
#> GAMLSS-RS iteration 7: Global Deviance = 4138.868 
#> GAMLSS-RS iteration 8: Global Deviance = 4136.94 
#> GAMLSS-RS iteration 9: Global Deviance = 4135.921 
#> GAMLSS-RS iteration 10: Global Deviance = 4135.375 
#> GAMLSS-RS iteration 11: Global Deviance = 4135.081 
#> GAMLSS-RS iteration 12: Global Deviance = 4134.92 
#> GAMLSS-RS iteration 13: Global Deviance = 4134.832 
#> GAMLSS-RS iteration 14: Global Deviance = 4134.784 
#> GAMLSS-RS iteration 15: Global Deviance = 4134.757 
#> GAMLSS-RS iteration 16: Global Deviance = 4134.742 
#> GAMLSS-RS iteration 17: Global Deviance = 4134.734 
#> GAMLSS-RS iteration 18: Global Deviance = 4134.729 
#> GAMLSS-RS iteration 19: Global Deviance = 4134.727 
#> GAMLSS-RS iteration 20: Global Deviance = 4134.725 
#> GAMLSS-RS iteration 21: Global Deviance = 4134.724 

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

# Example 2
# Generating random values under some model

# A function to simulate a data set with Y ~ DPERKS
gendat <- function(n) {
  x1 <- runif(n)
  x2 <- runif(n)
  mu    <- exp(-1.6 + 5 * x1)
  sigma <- exp(1.7 - 5 * x2)
  y <- rDPERKS(n=n, mu=mu, sigma=sigma)
  data.frame(y=y, x1=x1, x2=x2)
}

set.seed(12345)
datos <- gendat(n=1000)

mod2 <- NULL
mod2 <- gamlss(y~x1, sigma.fo=~x2, family=DPERKS, data=datos,
                 control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 4510.373 
#> GAMLSS-RS iteration 2: Global Deviance = 4371.655 
#> GAMLSS-RS iteration 3: Global Deviance = 4318.232 
#> GAMLSS-RS iteration 4: Global Deviance = 4292.156 
#> GAMLSS-RS iteration 5: Global Deviance = 4278.969 
#> GAMLSS-RS iteration 6: Global Deviance = 4272.119 
#> GAMLSS-RS iteration 7: Global Deviance = 4268.541 
#> GAMLSS-RS iteration 8: Global Deviance = 4266.698 
#> GAMLSS-RS iteration 9: Global Deviance = 4265.771 
#> GAMLSS-RS iteration 10: Global Deviance = 4265.32 
#> GAMLSS-RS iteration 11: Global Deviance = 4265.108 
#> GAMLSS-RS iteration 12: Global Deviance = 4265.01 
#> GAMLSS-RS iteration 13: Global Deviance = 4264.966 
#> GAMLSS-RS iteration 14: Global Deviance = 4264.946 
#> GAMLSS-RS iteration 15: Global Deviance = 4264.937 
#> GAMLSS-RS iteration 16: Global Deviance = 4264.933 
#> GAMLSS-RS iteration 17: Global Deviance = 4264.931 
#> GAMLSS-RS iteration 18: Global Deviance = 4264.931 

summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family:  c("DPERKS", "discrete-Perks") 
#> 
#> Call:  gamlss(formula = y ~ x1, sigma.formula = ~x2, family = DPERKS,  
#>     data = datos, control = gamlss.control(n.cyc = 500, trace = TRUE)) 
#> 
#> Fitting method: RS() 
#> 
#> ------------------------------------------------------------------
#> Mu link function:  log
#> Mu Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  -1.3044     0.2423  -5.382 9.18e-08 ***
#> x1            4.7674     0.9593   4.969 7.90e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> Sigma link function:  log
#> Sigma Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  1.77139    0.07243   24.45   <2e-16 ***
#> x2          -5.15249    0.11661  -44.19   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> No. of observations in the fit:  1000 
#> Degrees of Freedom for the fit:  4
#>       Residual Deg. of Freedom:  996 
#>                       at cycle:  18 
#>  
#> Global Deviance:     4264.931 
#>             AIC:     4272.931 
#>             SBC:     4292.562 
#> ******************************************************************

# Example 3
# The dataset comes from Tyagi et al. (2020) page 21
# The dataset contains the number of outbreaks of strikes in the
# UK coal mining industries in four successive week periods
# in the year 1948-59.
y <- rep(0:4, times=c(46, 76, 24, 9, 1))

# Fitting the model
library(gamlss)
mod3 <- gamlss(y~1, family=DPERKS,
               control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 396.3619 
#> GAMLSS-RS iteration 2: Global Deviance = 390.0067 
#> GAMLSS-RS iteration 3: Global Deviance = 385.4729 
#> GAMLSS-RS iteration 4: Global Deviance = 382.3254 
#> GAMLSS-RS iteration 5: Global Deviance = 380.219 
#> GAMLSS-RS iteration 6: Global Deviance = 378.8279 
#> GAMLSS-RS iteration 7: Global Deviance = 377.8972 
#> GAMLSS-RS iteration 8: Global Deviance = 377.2971 
#> GAMLSS-RS iteration 9: Global Deviance = 376.9094 
#> GAMLSS-RS iteration 10: Global Deviance = 376.659 
#> GAMLSS-RS iteration 11: Global Deviance = 376.4991 
#> GAMLSS-RS iteration 12: Global Deviance = 376.395 
#> GAMLSS-RS iteration 13: Global Deviance = 376.3274 
#> GAMLSS-RS iteration 14: Global Deviance = 376.2838 
#> GAMLSS-RS iteration 15: Global Deviance = 376.2556 
#> GAMLSS-RS iteration 16: Global Deviance = 376.2374 
#> GAMLSS-RS iteration 17: Global Deviance = 376.2257 
#> GAMLSS-RS iteration 18: Global Deviance = 376.2181 
#> GAMLSS-RS iteration 19: Global Deviance = 376.2133 
#> GAMLSS-RS iteration 20: Global Deviance = 376.2102 
#> GAMLSS-RS iteration 21: Global Deviance = 376.2082 
#> GAMLSS-RS iteration 22: Global Deviance = 376.2069 
#> GAMLSS-RS iteration 23: Global Deviance = 376.206 

# Extracting the fitted values for mu and sigma
# using the inverse link function
exp(coef(mod3, what="mu"))
#> (Intercept) 
#>  0.09176146 
exp(coef(mod3, what="sigma"))
#> (Intercept) 
#>    1.835224 

# Example 4
# The dataset comes from Tyagi et al. (2020) page 21
# The dataset contains the number fishes caught
# in traps (David, 1971, pg. 168).
y <- rep(0:9, times=c(1, 2, 11, 20, 29, 23, 10, 3, 1, 0))

# Fitting the model
library(gamlss)
mod4 <- gamlss(y~1, family=DPERKS,
               control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 359.9905 
#> GAMLSS-RS iteration 2: Global Deviance = 358.9466 
#> GAMLSS-RS iteration 3: Global Deviance = 358.1495 
#> GAMLSS-RS iteration 4: Global Deviance = 357.5214 
#> GAMLSS-RS iteration 5: Global Deviance = 357.0343 
#> GAMLSS-RS iteration 6: Global Deviance = 356.6513 
#> GAMLSS-RS iteration 7: Global Deviance = 356.3544 
#> GAMLSS-RS iteration 8: Global Deviance = 356.1232 
#> GAMLSS-RS iteration 9: Global Deviance = 355.9441 
#> GAMLSS-RS iteration 10: Global Deviance = 355.803 
#> GAMLSS-RS iteration 11: Global Deviance = 355.6919 
#> GAMLSS-RS iteration 12: Global Deviance = 355.6054 
#> GAMLSS-RS iteration 13: Global Deviance = 355.5368 
#> GAMLSS-RS iteration 14: Global Deviance = 355.4823 
#> GAMLSS-RS iteration 15: Global Deviance = 355.4392 
#> GAMLSS-RS iteration 16: Global Deviance = 355.4049 
#> GAMLSS-RS iteration 17: Global Deviance = 355.3778 
#> GAMLSS-RS iteration 18: Global Deviance = 355.3563 
#> GAMLSS-RS iteration 19: Global Deviance = 355.3392 
#> GAMLSS-RS iteration 20: Global Deviance = 355.3257 
#> GAMLSS-RS iteration 21: Global Deviance = 355.315 
#> GAMLSS-RS iteration 22: Global Deviance = 355.3064 
#> GAMLSS-RS iteration 23: Global Deviance = 355.2997 
#> GAMLSS-RS iteration 24: Global Deviance = 355.2943 
#> GAMLSS-RS iteration 25: Global Deviance = 355.29 
#> GAMLSS-RS iteration 26: Global Deviance = 355.2867 
#> GAMLSS-RS iteration 27: Global Deviance = 355.284 
#> GAMLSS-RS iteration 28: Global Deviance = 355.2818 
#> GAMLSS-RS iteration 29: Global Deviance = 355.2802 
#> GAMLSS-RS iteration 30: Global Deviance = 355.2788 
#> GAMLSS-RS iteration 31: Global Deviance = 355.2777 
#> GAMLSS-RS iteration 32: Global Deviance = 355.2768 

# Extracting the fitted values for mu and sigma
# using the inverse link function
exp(coef(mod4, what="mu"))
#> (Intercept) 
#> 0.003652202 
exp(coef(mod4, what="sigma"))
#> (Intercept) 
#>    1.238806 

# Example 5
# The dataset comes from Tyagi et al. (2020) page 24
# This dataset consists of remission times in weeks
# for 20 leukemia patients randomly assigned to a certain treatment
y <- c(1, 3, 3, 6, 7, 7, 10, 12, 14, 15, 18,
       19, 22, 26, 28, 29, 34, 40, 48, 49)

# Fitting the model
library(gamlss)
mod4 <- gamlss(y~1, family=DPERKS)
#> GAMLSS-RS iteration 1: Global Deviance = 358.4939 
#> GAMLSS-RS iteration 2: Global Deviance = 253.1631 
#> GAMLSS-RS iteration 3: Global Deviance = 200.0494 
#> GAMLSS-RS iteration 4: Global Deviance = 175.8855 
#> GAMLSS-RS iteration 5: Global Deviance = 165.5142 
#> GAMLSS-RS iteration 6: Global Deviance = 161.1295 
#> GAMLSS-RS iteration 7: Global Deviance = 159.2337 
#> GAMLSS-RS iteration 8: Global Deviance = 158.389 
#> GAMLSS-RS iteration 9: Global Deviance = 157.9998 
#> GAMLSS-RS iteration 10: Global Deviance = 157.8111 
#> GAMLSS-RS iteration 11: Global Deviance = 157.7192 
#> GAMLSS-RS iteration 12: Global Deviance = 157.6744 
#> GAMLSS-RS iteration 13: Global Deviance = 157.6517 
#> GAMLSS-RS iteration 14: Global Deviance = 157.6401 
#> GAMLSS-RS iteration 15: Global Deviance = 157.6341 
#> GAMLSS-RS iteration 16: Global Deviance = 157.631 
#> GAMLSS-RS iteration 17: Global Deviance = 157.6296 
#> GAMLSS-RS iteration 18: Global Deviance = 157.6286 

# Extracting the fitted values for mu and sigma
# using the inverse link function
exp(coef(mod4, what="mu"))
#> (Intercept) 
#>   0.4293858 
exp(coef(mod4, what="sigma"))
#> (Intercept) 
#>  0.08486852

Arguments

Details

References

See also

Author

Examples