The function DMOLBE() defines the Discrete Marshall-Olkin Length Biased
Exponential distribution, a two parameter
distribution, for a gamlss.family object to be used in GAMLSS fitting
using the function gamlss().
DMOLBE(mu.link = "log", sigma.link = "log")Returns a gamlss.family object which can be used
to fit a DMOLBE distribution
in the gamlss() function.
The DMOLBE distribution with parameters \(\mu\) and \(\sigma\) has a support 0, 1, 2, ... and mass function given by
\(f(x | \mu, \sigma) = \frac{\sigma ((1+x/\mu)\exp(-x/\mu)-(1+(x+1)/\mu)\exp(-(x+1)/\mu))}{(1-(1-\sigma)(1+x/\mu)\exp(-x/\mu)) ((1-(1-\sigma)(1+(x+1)/\mu)\exp(-(x+1)/\mu))}\)
with \(\mu > 0\) and \(\sigma > 0\)
Aljohani, H. M., Ahsan-ul-Haq, M., Zafar, J., Almetwally, E. M., Alghamdi, A. S., Hussam, E., & Muse, A. H. (2023). Analysis of Covid-19 data using discrete Marshall–Olkinin length biased exponential: Bayesian and frequentist approach. Scientific Reports, 13(1), 12243.
# Example 1
# Generating some random values with
# known mu and sigma
set.seed(1234)
y <- rDMOLBE(n=100, mu=10, sigma=7)
# Fitting the model
library(gamlss)
mod1 <- gamlss(y~1, sigma.fo=~1, family=DMOLBE,
control=gamlss.control(n.cyc=500, trace=FALSE))
# Extracting the fitted values for mu and sigma
# using the inverse link function
exp(coef(mod1, what="mu"))
#> (Intercept)
#> 9.042679
exp(coef(mod1, what="sigma"))
#> (Intercept)
#> 6.713267
# Example 2
# Generating random values under some model
# A function to simulate a data set with Y ~ DMOLBE
gendat <- function(n) {
x1 <- runif(n, min=0.4, max=0.6)
x2 <- runif(n, min=0.4, max=0.6)
mu <- exp(1.21 - 3 * x1) # 0.75 approximately
sigma <- exp(1.26 - 2 * x2) # 1.30 approximately
y <- rDMOLBE(n=n, mu=mu, sigma=sigma)
data.frame(y=y, x1=x1,x2=x2)
}
set.seed(123)
dat <- gendat(n=350)
# Fitting the model
mod2 <- NULL
mod2 <- gamlss(y~x1, sigma.fo=~x2, family=DMOLBE, data=dat,
control=gamlss.control(n.cyc=500, trace=FALSE))
summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family: c("DMOLBE", "Discrete Marshall-Olkin Length Biased Exponential" )
#>
#> Call: gamlss(formula = y ~ x1, sigma.formula = ~x2, family = DMOLBE,
#> data = dat, control = gamlss.control(n.cyc = 500, trace = FALSE))
#>
#> Fitting method: RS()
#>
#> ------------------------------------------------------------------
#> Mu link function: log
#> Mu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.1253 0.3342 3.367 0.000845 ***
#> x1 -3.0347 0.6771 -4.482 1.01e-05 ***
#> ---
#> 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.2774 0.8489 1.505 0.133
#> x2 -1.6791 1.6745 -1.003 0.317
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 350
#> Degrees of Freedom for the fit: 4
#> Residual Deg. of Freedom: 346
#> at cycle: 9
#>
#> Global Deviance: 952.4395
#> AIC: 960.4395
#> SBC: 975.8712
#> ******************************************************************
# Example 3
# Data Set I (death due to coronavirus in China). The first data set is the number
# of deaths due to coronavirus in China from 23 January to 28 March.
# The data sets used in the paper was collected from 2020 year. The data set
# is reported in https://www.worldometers.info/coronavirus/country/china/.
# The data are:
y <- c(8, 16, 15, 24, 26, 26, 38, 43, 46, 45, 57, 64, 65, 73, 73, 86, 89, 97,
108, 97, 146, 121, 143, 142, 105, 98, 136, 114, 118, 109, 97, 150, 71,
52, 29, 44, 47, 35, 42, 31, 38, 31, 30, 28, 27, 22, 17, 22, 11, 7,
13, 10, 14, 13, 11, 8, 3, 7, 6, 9, 7, 4, 6, 5, 3, 5)
# Fitting the model
mod3 <- gamlss(y~1, sigma.fo=~1, family=DMOLBE,
control=gamlss.control(n.cyc=500, trace=FALSE))
summary(mod3)
#> ******************************************************************
#> Family: c("DMOLBE", "Discrete Marshall-Olkin Length Biased Exponential" )
#>
#> Call: gamlss(formula = y ~ 1, sigma.formula = ~1, family = DMOLBE,
#> control = gamlss.control(n.cyc = 500, trace = FALSE))
#>
#> Fitting method: RS()
#>
#> ------------------------------------------------------------------
#> Mu link function: log
#> Mu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 3.6682 0.2663 13.78 <2e-16 ***
#> ---
#> 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.3031 0.5855 -2.226 0.0296 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 66
#> Degrees of Freedom for the fit: 2
#> Residual Deg. of Freedom: 64
#> at cycle: 2
#>
#> Global Deviance: 654.3575
#> AIC: 658.3575
#> SBC: 662.7368
#> ******************************************************************
# Extracting the fitted values for mu and sigma
# using the inverse link function
mu_hat <- exp(coef(mod3, what="mu"))
mu_hat
#> (Intercept)
#> 39.18313
sigma_hat <- exp(coef(mod3, what="sigma"))
sigma_hat
#> (Intercept)
#> 0.2716752
# Example 4
# Data Set II (daily death due to coronavirus in Pakistan). The second data
# set is the daily deaths due to coronavirus in Pakistan from 18 March
# to 30 June. The data sets used in the paper was collected from 2020 year.
# The data is reported in
# https://www.worldometers.info/coronavirus/country/Pakistan.
# The data are:
y <- c(1, 6, 6, 4, 4, 4, 1, 20, 5, 2, 3, 15, 17, 7, 8, 25, 8, 25, 11,
25, 16, 16, 12, 11, 20, 31, 42, 32, 23, 17, 19, 38, 50, 21, 14,
37, 23, 47, 31, 24, 9, 64, 39, 30, 36, 46, 32, 50, 34, 32, 34,
30, 28, 35, 57, 78, 88, 60, 78, 67, 82, 68, 97, 67, 65, 105,
83, 101, 107, 88, 178, 110, 136, 118, 136, 153, 119, 89, 105,
60, 148, 59, 73, 83, 49, 137, 91)
# Fitting the model
mod4 <- gamlss(y~1, sigma.fo=~1, family=DMOLBE,
control=gamlss.control(n.cyc=500, trace=FALSE))
summary(mod4)
#> ******************************************************************
#> Family: c("DMOLBE", "Discrete Marshall-Olkin Length Biased Exponential" )
#>
#> Call: gamlss(formula = y ~ 1, sigma.formula = ~1, family = DMOLBE,
#> control = gamlss.control(n.cyc = 500, trace = FALSE))
#>
#> Fitting method: RS()
#>
#> ------------------------------------------------------------------
#> Mu link function: log
#> Mu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 3.558 0.222 16.02 <2e-16 ***
#> ---
#> 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) -0.9268 0.4984 -1.86 0.0664 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 87
#> Degrees of Freedom for the fit: 2
#> Residual Deg. of Freedom: 85
#> at cycle: 2
#>
#> Global Deviance: 864.0005
#> AIC: 868.0005
#> SBC: 872.9324
#> ******************************************************************
# Extracting the fitted values for mu and sigma
# using the inverse link function
mu_hat <- exp(coef(mod4, what="mu"))
mu_hat
#> (Intercept)
#> 35.0795
sigma_hat <- exp(coef(mod4, what="sigma"))
sigma_hat
#> (Intercept)
#> 0.3958021