DiscreteDists (version 1.1.3)

NPGL2: Mean-parametrized New Poisson-generalised Lindley distribution (NPGL2)

Description

The function NPGL2() defines the mean-parametrized version of the Poisson-generalised Lindley distribution, a two-parameter distribution, for a gamlss.family object to be used in GAMLSS fitting using the function gamlss().

Usage

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

Value

Returns a gamlss.family object which can be used to fit a mean-parametrized NPGL distribution in the gamlss() function.

Arguments

mu.link

defines the mu.link, with "log" link as the default for the mu parameter (the mean of the distribution).

sigma.link

defines the sigma.link, with "log" link as the default for the sigma parameter (parameter alpha of the original NPGL).

Author

Tomas Mesa,

Details

This family uses the mean-parametrized version of the NPGL distribution proposed in Section 6 of Altun (2021). The new parameters are:

• \(\mu > 0\): the mean of the distribution, i.e. \(E[X]\).

• \(\sigma > 0\)

The reparametrization links the mean and \(\sigma\) to the original parameter \(\mu\) via the mean equation \(E[X] = (\sigma + \mu) / (\mu (1 + \mu))\)

References

Altun, E. A new two-parameter discrete poisson-generalized Lindley distribution with properties and applications to healthcare data sets. Comput Stat 36, 2841-2861 (2021). https://doi.org/10.1007/s00180-021-01097-0

See Also

dNPGL2, NPGL.

Examples

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

set.seed(123)
y <- rNPGL2(n=5000, mu=4, sigma=2)

# Fitting the model
library(gamlss)
mod1 <- gamlss(y~1, sigma.fo=~1, family=NPGL2,
               control=gamlss.control(n.cyc=500, trace=TRUE))

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

# Example 2
# Generating random values under some model

# A function to simulate a data set with Y ~ NPGL2
gendat <- function(n) {
  x1 <- runif(n)
  x2 <- runif(n)
  mu    <- exp(1.5 - 1.8 * x1)  # mean of the distribution
  sigma <- exp(0.5 + 0.8 * x2)  # shape parameter alpha
  y <- rNPGL2(n=n, mu=mu, sigma=sigma)
  data.frame(y=y, x1=x1, x2=x2)
}

set.seed(123)
dat <- gendat(n=1000)

# Fitting the model
mod2 <- gamlss(y~x1, sigma.fo=~x2, family=NPGL2, data=dat,
               control=gamlss.control(n.cyc=800, trace=TRUE))

summary(mod2)

# Example 3
# Using the data from Altun (2021), Section 8.1.
# The dataset is about the 1991 Arizona cardiovascular patient.
# We model the length of hospital stay (los) with cardiovascular
# procedure, sex, type of admission and age as covariates.
# The response variable has a dispersion index of 5.43,
# confirming over-dispersion.
# Data available in the azpro dataset of the COUNT package.

mod3 <- gamlss(
  los ~ procedure + sex + admit + age75,
  sigma.fo = ~1,
  family = NPGL2(mu.link = "log"),
  data = azpro,
  control = gamlss.control(n.cyc = 500, trace = TRUE)
  )

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

# Note: coef(mod3, what="mu") is the estimated mean length
# of hospital stay for the reference group (PTCA procedure,
# female, elective admission, age <= 75), directly interpretable
# because mu is the mean of the NPGL2 distribution.

# Example 4
# Using the data from Altun (2021), Section 8.2.
# The dataset originates from the US National Medical Expenditure
# Survey (NMES) conducted in 1987 and 1988. We model the number
# of physician office visits with hospital stays, chronic conditions,
# activity limitations, age, gender, marital status, income,
# employment status, Medicaid, private insurance and self-perceived
# health status as covariates.
# The response variable has a dispersion index of 7.91,
# confirming over-dispersion.
# Data available in the NMES1988 dataset of the AER package.
# Set "average" as the reference category for health status,
# so that healthexcellent and healthpoor are the two dummies,
# matching x11 and x12 in Altun (2021).

NMES1988$health <- relevel(factor(NMES1988$health), ref = "average")

mod4 <- gamlss(
  visits ~ hospital + chronic + adl + age + gender + married +
    income + employed + medicaid + insurance + health,
  sigma.fo = ~1,
  family = NPGL2(mu.link = "log"),
  data = NMES1988,
  control = gamlss.control(n.cyc = 500, trace = TRUE)
)

coef(mod4, what="mu")
coef(mod4, what="sigma")