The function GGEO() defines the Generalized Geometric distribution,
a two parameter distribution,
for a gamlss.family object to be used in GAMLSS fitting
using the function gamlss().
GGEO(mu.link = "logit", sigma.link = "log")defines the mu.link, with "log" link as the default for the mu parameter.
defines the sigma.link, with "logit" link as the default for the sigma. Other links are "probit" and "cloglog"'(complementary log-log)
Returns a gamlss.family object which can be used
to fit a GGEO distribution
in the gamlss() function.
The GGEO distribution with parameters \(\mu\) and \(\sigma\) has a support 0, 1, 2, ... and mass function given by
\(f(x | \mu, \sigma) = \frac{\sigma \mu^x (1-\mu)}{(1-(1-\sigma) \mu^{x+1})(1-(1-\sigma) \mu^{x})}\)
with \(0 < \mu < 1\) and \(\sigma > 0\). If \(\sigma=1\), the GGEO distribution reduces to the geometric distribution with success probability \(1-\mu\).
Gómez-Déniz E (2010). “Another generalization of the geometric distribution.” Test, 19, 399-415.
# Example 1
# Generating some random values with
# known mu and sigma
set.seed(123)
y <- rGGEO(n=200, mu=0.95, sigma=1.5)
# Fitting the model
library(gamlss)
mod1 <- gamlss(y~1, family=GGEO,
control=gamlss.control(n.cyc=500, trace=FALSE))
# Extracting the fitted values for mu and sigma
# using the inverse link function
inv_logit <- function(x) 1/(1 + exp(-x))
inv_logit(coef(mod1, what="mu"))
#> (Intercept)
#> 0.9433538
exp(coef(mod1, what="sigma"))
#> (Intercept)
#> 1.924394
# Example 2
# Generating random values under some model
# A function to simulate a data set with Y ~ GGEO
gendat <- function(n) {
x1 <- runif(n)
x2 <- runif(n)
mu <- inv_logit(1.7 - 2.8*x1)
sigma <- exp(0.73 + 1*x2)
y <- rGGEO(n=n, mu=mu, sigma=sigma)
data.frame(y=y, x1=x1, x2=x2)
}
set.seed(78353)
datos <- gendat(n=100)
mod2 <- gamlss(y~x1, sigma.fo=~x2, family=GGEO, data=datos,
control=gamlss.control(n.cyc=500, trace=FALSE))
summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family: c("GGEO", "Generalized geometric")
#>
#> Call: gamlss(formula = y ~ x1, sigma.formula = ~x2, family = GGEO,
#> data = datos, control = gamlss.control(n.cyc = 500, trace = FALSE))
#>
#> Fitting method: RS()
#>
#> ------------------------------------------------------------------
#> Mu link function: logit
#> Mu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.5645 0.1602 9.766 3.89e-16 ***
#> x1 -3.0030 0.3780 -7.945 3.32e-12 ***
#> ---
#> 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.7825 0.3690 4.831 5.03e-06 ***
#> x2 -0.4900 0.6196 -0.791 0.431
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 100
#> Degrees of Freedom for the fit: 4
#> Residual Deg. of Freedom: 96
#> at cycle: 19
#>
#> Global Deviance: 413.8106
#> AIC: 421.8106
#> SBC: 432.2313
#> ******************************************************************
# Example 3
# Number of accidents to 647 women working on H. E. Shells
# for 5 weeks. Taken from Gomez-Deniz (2010) page 411.
y <- rep(x=0:5, times=c(447, 132, 42, 21, 3, 2))
mod3 <- gamlss(y~1, family=GGEO,
control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 1184.469
# Extracting the fitted values for mu and sigma
# using the inverse link function
inv_logit <- function(x) 1/(1 + exp(-x))
inv_logit(coef(mod3, what="mu"))
#> (Intercept)
#> 0.3399179
exp(coef(mod3, what="sigma"))
#> (Intercept)
#> 0.8738869