The function COMPO2() defines the COM-POISSON 2
distribution, a two parameter
distribution, for a gamlss.family object to be used in GAMLSS fitting
using the function gamlss().
COMPO2(mu.link = "log", sigma.link = "identity")Returns a gamlss.family object which can be used
to fit a COMPO2 distribution
in the gamlss() function.
The COMPO2 distribution with parameters \(\mu\) and \(\sigma\) has a support 0, 1, 2, ... and mass function given by
\(f(x | \mu, \sigma) = \left(\mu + \frac{\exp(\sigma)-1}{2 \exp(\sigma)} \right)^{x \exp(\sigma)} \frac{(x!)^{\exp(\sigma)}}{Z(\mu, \sigma)} \)
with \(\mu > 0\), \(\sigma \in \Re\) and
\(Z(\mu, \sigma)=\sum_{j=0}^{\infty} \frac{\mu^j}{(j!)^\sigma}\).
The proposed functions here are based on the functions from the COMPoissonReg package.
Ribeiro Jr, Eduardo E., et al. "Reparametrization of COM–Poisson regression models with applications in the analysis of experimental data." Statistical Modelling 20.5 (2020): 443-466.
# Example 1
# Generating some random values with
# known mu and sigma
set.seed(1234)
y <- rCOMPO2(n=100, mu=10, sigma=-1)
# Fitting the model
library(gamlss)
mod1 <- gamlss(y~1, sigma.fo=~1, family=COMPO2,
control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 576.6895
#> GAMLSS-RS iteration 2: Global Deviance = 576.6811
#> GAMLSS-RS iteration 3: Global Deviance = 576.6808
# Extracting the fitted values for mu and sigma
# using the inverse link function
exp(coef(mod1, what="mu"))
#> (Intercept)
#> 8.889252
coef(mod1, what="sigma")
#> (Intercept)
#> -0.9100409
# Example 2
# Generating random values under some model
# \donttest{
# A function to simulate a data set with Y ~ COMPO2
gendat <- function(n) {
x1 <- runif(n)
x2 <- runif(n)
mu <- exp(2 + 1 * x1) # 12 approximately
sigma <- 1 - 2 * x2 # 0 approximately
y <- rCOMPO2(n=n, mu=mu, sigma=sigma)
data.frame(y=y, x1=x1, x2=x2)
}
set.seed(123)
dat <- gendat(n=200)
# Fitting the model
mod2 <- NULL
mod2 <- gamlss(y~x1, sigma.fo=~x2, family=COMPO2, data=dat)
#> GAMLSS-RS iteration 1: Global Deviance = 1057.584
#> GAMLSS-RS iteration 2: Global Deviance = 1057.304
#> GAMLSS-RS iteration 3: Global Deviance = 1057.304
summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family: c("COMPO2", "Comway Maxwell Poisson second parameterization")
#>
#> Call: gamlss(formula = y ~ x1, sigma.formula = ~x2, family = COMPO2,
#> data = dat)
#>
#> Fitting method: RS()
#>
#> ------------------------------------------------------------------
#> Mu link function: log
#> Mu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.03088 0.04070 49.9 <2e-16 ***
#> x1 0.95385 0.06534 14.6 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> Sigma link function: identity
#> Sigma Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.8724 0.2034 4.290 2.79e-05 ***
#> x2 -1.7247 0.3734 -4.618 6.95e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 200
#> Degrees of Freedom for the fit: 4
#> Residual Deg. of Freedom: 196
#> at cycle: 3
#>
#> Global Deviance: 1057.304
#> AIC: 1065.304
#> SBC: 1078.498
#> ******************************************************************
# }