Package 'flexCountReg'

Conway-Maxwell-Poisson (COM) Distribution

Description

These functions provide the density function, distribution function, quantile function, and random number generation for the Conway-Maxwell-Poisson (COM) Distribution

Usage

dcom(x, mu = NULL, lambda = 1, nu = 1, log = FALSE)

pcom(q, mu = NULL, lambda = 1, nu = 1, lower.tail = TRUE, log.p = FALSE)

qcom(p, mu = NULL, lambda = 1, nu = 1)

rcom(n, mu = NULL, lambda = 1, nu = 1)

Arguments

x	numeric value or a vector of values.
mu	optional. Numeric value or vector of mean values for the distribution (the values have to be greater than 0).
lambda	optional. Numeric value or vector of values for the rate parameter of the distribution (the values have to be greater than 0). If mu is provided, lambda is ignored.
nu	optional. Numeric value or vector of values for the decay parameter of the distribution ((the values have to be greater than 0).
log	logical; if TRUE, probabilities p are given as log(p).
q	quantile or a vector of quantiles.
lower.tail	logical; if TRUE, probabilities p are $P[X\leq x]$ otherwise, $P[X>x]$.
log.p	logical; if TRUE, probabilities p are given as log(p).
p	probability or a vector of probabilities.
n	the number of random numbers to generate.

Details

dcom computes the density (PDF) of the COM Distribution.

pcom computes the CDF of the COM Distribution.

qcom computes the quantile function of the COM Distribution.

rcom generates random numbers from the COM Distribution.

The Probability Mass Function (PMF) for the Conway-Maxwell-Poisson distribution is:

$f(x|\lambda, \nu) = \frac{\lambda^x}{(x!)^\nu Z(\lambda,\nu)}$

Where $\lambda$ and $\nu$ are distribution parameters with $\lambda>0$ and $\nu>0$, and $Z(\lambda,\nu)$ is the normalizing constant.

The normalizing constant is given by:

$Z(\lambda,\nu)=\sum_{n=0}^{\infty}\frac{\lambda^n}{(n!)^\nu}$

The mean and variance of the distribution are given by:

$E[x]=\mu=\lambda \frac{\delta}{\delta \lambda} \log(Z(\lambda,\nu))$

$Var(x)=\lambda \frac{\delta}{\delta \lambda} \mu$

When the mean value is given, the rate parameter ($\lambda$) is computed using the mean and the decay parameter ($\nu$). This is useful to allow the calculation of the rate parameter when the mean is known (e.g., in regression))

Value

dcom gives the density, pcom gives the distribution function, qcom gives the quantile function, and rcom generates random deviates.

The length of the result is determined by n for rcom, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

Examples

dcom(1, mu=0.75, nu=3)
pcom(c(0,1,2,3,5,7,9,10), lambda=0.75, nu=0.75)
qcom(c(0.1,0.3,0.5,0.9,0.95), mu=0.75, nu=0.75)
rcom(30, mu=0.75, nu=0.5)

---

Generate a covariance matrix using a correlation matrix and vector of standard deviations

Description

Generate a covariance matrix using a correlation matrix and vector of standard deviations

Usage

cor2cov(C, S)

Arguments

C	A correlation matrix.
S	A vector of standard deviations.

Value

A covariance matrix

Examples

C <- matrix(c(1,-0.3,0.7,-0.3,1,-0.2,0.7,-0.2,1), 3, 3)
S <- c(0.5, 2, 1.25)
cor2cov(C,S)

---

Generate Correlated Random Variables Using Halton or Scrambled Halton Draws

Description

This function generates N correlated random variables using Halton or scrambled Halton draws. The function supports normal and truncated normal distributions.

Usage

corr_haltons(
  means,
  cholesky = NULL,
  stdev = NULL,
  correlations = NULL,
  hdraws = NULL,
  ndraws = 500,
  scrambled = FALSE,
  dist = "normal",
  lower = -Inf,
  upper = Inf
)

Arguments

means	A numeric vector of means for each variable.
cholesky	A Cholesky decomposition matrix to introduce correlation.
stdev	A numeric vector of standard deviations for each variable. If provided, the function will use these values instead of the Cholesky decomposition matrix (must also provide a correlation matrix if providing standard deviations). Default is NULL.
correlations	A correlation matrix to introduce correlation. If provided, the function will use these values instead of the Cholesky decomposition matrix (must also provide standard deviations). Default is NULL.
hdraws	A matrix of Halton or scrambled Halton draws. If provided, the function will use these draws instead of generating new ones. Default is NULL.
ndraws	An integer specifying the number of values to simulate for each variable. Default is 500.
scrambled	A logical value indicating whether to use scrambled Halton draws. Default is FALSE.
dist	A character string specifying the distribution type. Options are "normal" and "truncated_normal". Default is "normal".
lower	A numeric value specifying the lower bound for truncated normal distribution. Default is -Inf.
upper	A numeric value specifying the upper bound for truncated normal distribution. Default is Inf.

Value

A matrix with N columns and ndraws rows containing the simulated values for the correlated random variables.

Examples

# Define mean, correlation, and standard deviations
means <- c(3, 2, 0.9)
sdevs <- c(0.25,1.5,0.8)
CORR <- matrix(c(1, -0.3, 0.5, -0.3, 1, -0.2, 0.5, -0.2, 1), 3, 3)

# Create the Cholesky decomposition matrix and set values for ndraws, etc.
ndraws <- 5000
scrambled <- TRUE
dist <- "normal"

# simulated the data
simulated_data <- corr_haltons(means, stdev=sdevs, correlations=CORR,
                                ndraws=ndraws, scrambled=scrambled,
                                dist=dist)

# look at the mean, standard deviation, and correlation of the simulated data
apply(simulated_data, 2, mean)
apply(simulated_data, 2, sd)
cor(simulated_data)

# providing a cholesky decomposition matrix
dist <- "normal"
cholesky <- chol(cor2cov(CORR, sdevs))
simulated_data <- corr_haltons(means, cholesky=cholesky, ndraws=ndraws,
                                scrambled=scrambled, dist=dist)
apply(simulated_data, 2, mean)
apply(simulated_data, 2, sd)
cor(simulated_data)

# Truncated normal
dist <- "truncated_normal"
lower <- 0
upper <- 30
simulated_data <- corr_haltons(means, cholesky=cholesky, ndraws=ndraws,
                                scrambled=scrambled, dist=dist,
                                lower=lower, upper=upper)
apply(simulated_data, 2, mean)
apply(simulated_data, 2, sd)
cor(simulated_data)

---

Count regression models

Description

The purpose of this function is to estimate count regression models using maximum likelihood estimation (MLE) or Maximum Simulated Likelihood Estimation (MSLE). The function can estimate the following models:

• Poisson (Poisson)

• Generalized Poisson Version 1 (GP1)

• Generalized Poisson Version 2 (GP2)

• Negative Binomial 1 (NB1)

• Negative Binomial 2 (NB2)

• Negative Binomial P (NBP)

• Poisson-Lognormal (PLN)

• Poisson-Generalized-Exponential (PGE)

• Poisson-Inverse-Gaussian Type 1 (PIG1)

• Poisson-Inverse-Gaussian Type 2 (PIG2)

• Poisson-Inverse-Gamma (PIG)

• Poisson-Lindley (PL)

• Poisson-Lindley-Gamma (PLG), also known as the Negative Binomial-Lindley (NBL)

• Poisson-Lindley-Lognormal (PLL)

• Poisson-Weibull (PW)

• Sichel (SI)

• Generalized Waring (GW)

• Conway-Maxwell-Poisson (COM)

Usage

countreg(
  formula,
  data,
  family = "NB2",
  offset = NULL,
  weights = NULL,
  verbose = FALSE,
  dis_param_formula_1 = NULL,
  dis_param_formula_2 = NULL,
  underreport_formula = NULL,
  underreport_family = "logit",
  ndraws = 1500,
  method = "NM",
  max.iters = 1000,
  start.vals = NULL,
  stderr = "normal",
  bootstraps = NULL,
  engine = "halton"
)

Arguments

formula	a symbolic description of the model to be fitted.
data	a data frame containing the variables in the model.
family	the name of the distribution/model type to estimate. The default "NB2" is the standard negative binomial distribution with a log link. other options are listed below.
offset	the name of a variable, or vector of variable names, in the data frame that should be used as an offset (i.e., included but forced to have a coefficient of 1). The normal method of setting an offset in the equation can also be used (overrides the offset option).
weights	the name of a variable in the data frame that should be used as a frequency weight.
verbose	an optional parameter. If TRUE, the function will print out the progress of the model fitting. Default is FALSE.
dis_param_formula_1	a symbolic description of the model for the natural log of the dispersion parameter or first parameter of the count distribution used. Further details are provided below.
dis_param_formula_2	a symbolic description of the model for the second parameter of the count distribution used. Further details are provided below.
underreport_formula	an optional formula to estimate the underreporting for any of the count model options. The underreporting is estimated as a function (logit or probit) of the predictors in the model. For the model to be tractable, the independent variables cannot be the exact same as the count model. The default is NULL.
underreport_family	the name of the distribution/model type to estimate the underreporting portion of the model when underreport_formula is specified. The default is "logit" for a binary logistic regression model. The other option is "probit" for a probit model.
ndraws	The number of Halton draws for integrating the distribution being compounded with the Poisson distribution when there is not a closed-form solution. Default is 1500. It is recommended to test different numbers of draws to determine if the model is stable (i.e., doesn't change or has minimal change as the number of draws changes within a reasonable range).
method	Optimization method to be used for maximum likelihood estimation. See maxLik documentation for options. The default is "NM" for the Nelder-Mead method.
max.iters	Maximum number of iterations for the optimization method.
start.vals	Optional vector of starting values for the optimization.
stderr	Type of standard errors to use. The default is "normal". Other options include "boot" for bootstrapped standard errors, or "robust" for robust standard errors.
bootstraps	Optional integer specifying the number of bootstrap samples to be used for estimating standard errors when stderr= "boot". Note that this currently does not work when an offset variable is used.
engine	the engine to use for the integration if the Poisson Lognormal (PLN) distribution is used. Options include "halton" which uses Halton draws or "poilog" which uses the poilog package. The default is "halton".

Details

For the family argument, the following options are available:

• "POISSON" for Poisson distribution with a log link.

• "GP1" for Generalized Poisson Version 1 Distribution with a log link

• "GP2" for Generalized Poisson Version 2 Distribution with a log link

• "NB1" for Negative Binomial 1 distribution with a log link.

• "NB2" for Negative Binomial 2 distribution with a log link (i.e., the standard negative binomial model).

• "NBP" for Negative Binomial P distribution with a log link.

• "PLN" for Poisson-Lognormal distribution with a log link.

• "PGE" for Poisson-Generalized-Exponential distribution with a log link.

• "PIG1" for Poisson-Inverse-Gaussian Type-1 distribution with a log link.

• "PIG2" for Poisson-Inverse-Gaussian Type-2 distribution with a log link.

• "PIG" for Poisson-Inverse-Gamma distribution with a log link.

• "PL" for Poisson-Lindley distribution with a log link.

• "PLG" for Poisson-Lindley-Gamma distribution with a log link.

• "PLL" for Poisson-Lindley-Lognormal distribution with a log link.

• "PW" for Poisson-Weibull distribution with a log link.

• "SI" for Sichel distribution with a log link.

• "GW" for Generalized Waring distribution with a log link.

• "COM" for Conway-Maxwell-Poisson (COM) distribution with a log link.

The dis_param_formula_1 and dis_param_formula_2 parameters are used to estimate the dispersion parameter or other parameters of the count distribution used. This leads to the distributions parameters being functions rather than constants in the model. For example, if the user wants to estimate the overdispersion parameter of the Negative Binomial 2 distribution as a function of the variable x1 and x2, the user would specify dis_param_formula_1 = ~ x1 + x2. In the case of the Negative Binomial distributions, the model is known as a Generalized Negative Binomial model when the overdispersion parameter is specified as a function.

The function linking the distribution parameters to the predictors is:

$Param = \exp((Intercept) + \sum \beta X)$

The parameters for the different models are as follows:

For dis_param_formula_1, the models are for the parameters:

• $ln(\phi+1)$ for the Generlized Poisson Version 1.

• $ln(\alpha+1)$ for the Generlized Poisson Version 2.

• $ln(\alpha)$ for the Negative Binomial 1 model.

• $ln(\alpha)$ for the Negative Binomial 2 model.

• $ln(\alpha)$ for the Negative Binomial P model.

• $ln(\sigma)$ for the Poisson-Lognormal model.

• shape parameter for the Poisson-Generalized-Exponential model.

• $ln(\eta)$ for the Poisson-Inverse-Gaussian model.

• $ln(\eta)$ for the Poisson-Inverse-Gamma model.

• $ln(\theta)$ for the Poisson-Lindley model.

• $ln(\theta)$ for the Poisson-Lindley-Gamma model.

• $ln(\theta)$ for the Poisson-Lindley-Lognormal model.

• $ln(\alpha)$ for the Poisson-Weibull model.

• $\gamma$ for the Sichel model.

• $k$ for the Generalized Waring model.

• $ln(\nu)$ for the Conway-Maxwell-Poisson model.

For dis_param_formula_2, the models are for the parameters:

• Not Applicable for the Generalized Poisson Version 1 model.

• Not Applicable for the Negative Binomial 1 model.

• Not Applicable for the Negative Binomial 2 model.

• p for the Negative Binomial P model.

• Not Applicable for the Poisson-Lognormal model.

• scale parameter for the Poisson-Generalized-Exponential model.

• Not Applicable for the Poisson-Inverse-Gaussian model.

• Not Applicable for the Poisson-Inverse-Gamma model.

• Not Applicable for the Poisson-Lindley model.

• $ln(\alpha)$ for the Poisson-Lindley-Gamma model.

• $ln(\sigma)$ for the Poisson-Lindley-Lognormal model.

• $ln(\sigma)$ for the Poisson-Weibull model.

• $ln(\sigma)$ for the Sichel model.

• $ln(\rho)$ for the Generalized Waring model.

• Not Applicable for the Conway-Maxwell-Poisson model.

The ndraws parameter is used to estimate the distribution when there is not a closed-form solution. This uses Halton draws to integrate the distribution being compounded with the Poisson distribution. The default is 1500. The models this is applicable for include:

• Poisson-Lognormal

• Poisson-Generalized-Exponential

• Poisson-Lindley-Gamma (more efficient than using hypergeometric functions)

• Poisson-Lindley-Lognormal

• Poisson-Weibull

Value

An object of class countreg which is a list with the following components:

• model: the fitted model object.

• data: the data frame used to fit the model.

• call: the matched call.

• formula: the formula used to fit the model.

Model Details

Poisson Model

This implements the Poisson regression model using Maximum Likelihood Estimation (MLE), as opposed to the Iteratively Reweighted Least Squares (IRLS) method used in the glm function.

The PMF and log-likelihood functions are:

$P(Y = y) = \frac{e^{-\mu} \mu^y}{y!}$

$LL_{\text{Poisson}}(\beta) = \sum_{i=1}^n \left[ -\mu_i + y_i \ln(\mu_i) - \ln(y_i!) \right]$

The mean is:

$\mu = exp(X\beta)$

The variance is:

$\text{Var}(Y) = \mu$

Generalized Poisson Version 1 Model

This implements the Generalized Poisson Version 1 model using MLE.

The PMF is:

$f(y|\phi,\mu)=\frac{\mu(\mu+\phi y)^{y-1} exp\left(-\frac{\mu+\phi y} {1+\phi}\right)}{(1+\phi)^y y!}$

Where $\phi$ is a scale parameter with the restriction that $\eta>0$, $\mu>0$ is the mean value, and $y$ is a non-negative integer. The mean is modeled as:

$\mu = exp(X\beta)$

The variance of the GP-1 distribution is:

$\sigma^2=(1+\phi)^2 \mu$

If $\phi>0$, the distribution is overdispersed. If $\phi=0$, the distribution is equidispersed. If $\phi<0$, the distribution is underdispersed.

Furthermore, $\phi>-1$ is required for this distribution. When $\phi<0$, there is also a maximum value of support for the integar value $y$. This is $y_max=\left\lfloor -\frac{\mu}{\phi} \right\rfloor$.

Generalized Poisson Version 2 Model

This implements the Generalized Poisson Version 2 model using MLE.

The PMF is:

$f(y|\mu,\alpha)= \frac{\mu(\mu+\alpha\mu y)^{y-1} \exp\left(-\frac{\mu+\alpha\mu y}{1+\alpha\mu}\right)} {(1+\alpha\mu)^y y!}$

Where $\phi$ is a scale parameter with the restriction that $\eta>0$, $\mu>0$ is the mean value, and $y$ is a non-negative integer. The mean is modeled as:

$\mu = exp(X\beta)$

The variance of the GP-2 distribution is:

$\sigma^2=\mu(1+\alpha\mu)^2$

If $\alpha > 0$, the distribution is overdispersed. If $\alpha = 0$, the distribution is equidispersed and reduces to the ordinary Poisson distribution. If $\alpha < 0$ and $1 + \alpha\mu > 0$, the distribution is underdispersed and has finite support.

Under the underdispersed case, the largest possible value of $y$ is bounded above by $\left\lceil -1/\alpha \right\rceil - 1$, i.e., the largest integer satisfying $1 + \alpha y > 0$.

Negative Binomial Models**

The NB-1, NB-2, and NB-P versions of the negative binomial distribution are based on Greene (2008). The details of each of these are provided below.

NB-1 Model

The PMF and log-likelihood functions are:

$P(Y = y) = \frac{\Gamma(y + \frac{\mu}{\alpha})}{y! \, \Gamma(\frac{\mu}{\alpha})} \left( \frac{\frac{\mu}{\alpha}}{\frac{\mu}{\alpha} + \mu} \right) ^{\frac{\mu}{\alpha}} \left( \frac{\mu}{\frac{\mu}{\alpha} + \mu} \right)^y$

$LL_{\text{NB1}}(\beta, \alpha) = \sum_{i=1}^n \left[ \ln \Gamma\left( y_i + \frac{\mu_i}{\alpha} \right) - \ln \Gamma\left( \frac{\mu_i}{\alpha} \right) - \ln y_i! + \frac{\mu_i}{\alpha} \ln \left( \frac{\frac{\mu_i}{\alpha}}{\frac{\mu_i}{\alpha} + \mu_i} \right) + y_i \ln \left( \frac{\mu_i}{\frac{\mu_i}{\alpha} + \mu_i} \right) \right]$

The mean is:

$\mu = exp(X\beta)$

The variance is:

$\text{Var}(Y) = \mu + \alpha\mu$

NB-2 Model

The PMF and log-likelihood functions are:

$P(Y = y) = \frac{\Gamma(y + \alpha)}{y! \, \Gamma(\alpha)} \left( \frac{\alpha}{\alpha + \mu} \right)^\alpha \left( \frac{\mu}{\alpha + \mu} \right)^y$

$LL_{\text{NB2}} = \sum_{i=1}^n \left[ \ln \Gamma(y_i + \alpha) - \ln \Gamma(\alpha) - \ln y_i! + \alpha \ln \left( \frac{\alpha}{\alpha + \mu_i} \right) + y_i \ln \left( \frac{\mu_i}{\alpha + \mu_i} \right) \right]$

The mean is:

$\mu = exp(X\beta)$

The variance is:

$\text{Var}(Y) = \mu + \alpha\mu^2$

NB-P Model

The PMF and log-likelihood functions are:

$P(Y = y) = \frac{\Gamma(y + \frac{\mu^{2-p}}{\alpha})}{y! \, \Gamma(\frac{\mu^{2-p}}{\alpha})} \left( \frac{\frac{\mu^{2-p}}{\alpha}} {\frac{\mu^{2-p}}{\alpha} + \mu} \right)^{\frac{\mu^{2-p}}{\alpha}} \left( \frac{\mu}{\frac{\mu^{2-p}}{\alpha} + \mu} \right)^y$

$LL_{\text{NBP}}(\beta, \alpha, p) = \sum_{i=1}^n \left[ \ln \Gamma\left( y_i + \frac{\mu_i^{2-p}}{\alpha} \right) - \ln \Gamma\left( \frac{\mu_i^{2-p}}{\alpha} \right) - \ln y_i! + \frac{\mu_i^{2-p}}{\alpha} \ln \left( \frac{\frac{\mu_i^{2-p}}{\alpha}}{\frac{\mu_i^{2-p}}{\alpha} + \mu_i} \right) + y_i \ln \left( \frac{\mu_i}{\frac{\mu_i^{2-p}}{\alpha} + \mu_i} \right) \right]$

The mean is:

$\mu = exp(X\beta)$

The variance is:

$\text{Var}(Y) = \mu + \alpha\mu^P$

Poisson-Lognormal (PLN) Model

The compound Probability Mass Function(PMF) for the Poisson-Lognormal distribution is:

$f(y|\lambda,\sigma)=\int_0^\infty \frac{\lambda^y x^y e^{-\lambda x}} {y!}\frac{exp\left(-\frac{ln^2(x)}{2\sigma^2} \right)}{x\sigma\sqrt{2\pi}}dx$

Where $\sigma$ is a parameter for the lognormal distribution with the restriction $\sigma>0$, and $y$ is a non-negative integer.

The expected value of the distribution is:

$E[y]=e^{X\beta+\sigma^2/2} = \mu e^{\sigma^2/2}$

When ln.sigma.formula is used, the parameter $\sigma$ is modeled as:

$ln(\sigma)=\beta_0+\beta_1 x_1 + \cdots + \beta_n x_n$

Thus, the resulting value for the parameter $\sigma$ is:

$\sigma=e^{\beta_0+\beta_1 x_1 + \cdots + \beta_n x_n}$

The t-statistics and p-values for the coefficients related to ln(sigma) are, by default, testing if the coefficients are different from a value of 0. This has little practical meaning given that they are coefficients for ln(sigma). They are not testing if the coefficients have statistical significance in terms of improvement over a Poisson model. The Likelihood-Ratio test results provided in the output provide a test comparing if the Poisson-Lognormal model provides a statistically significant improvement in model fit over the Poisson model.

Poisson Generalized-Exponential (PGE) Model

The Generalized Exponential distribution can be written as a function with a shape parameter $\alpha>0$ and scale parameter $\gamma>0$. The distribution has strictly positive continuous values. The PDF of the distribution is:

$f(x|\alpha,\gamma)=\frac{\alpha}{\gamma}\left(1-e^{-\frac{x}{\gamma}} \right)^{\alpha-1}e^{-\frac{x}{\gamma}}$

Thus, the compound Probability Mass Function(PMF) for the PGE distribution is:

$f(y|\lambda,\alpha,\beta)=\int_0^\infty \frac{\lambda^y x^y e^{-\lambda x}}{y!}\frac{\alpha}{\gamma}\left(1-e^{-\frac{x}{\gamma}} \right)^{\alpha-1}e^{-\frac{x}{\gamma}} dx$

The expected value of the distribution is:

$E[y]=\mu=\lambda \left(\frac{\psi(\alpha+1)-\psi(1)}{\gamma}\right)$

Where $\psi(\cdot)$ is the digamma function.

The variance is:

$\sigma^2=\lambda \left(\frac{\psi(\alpha+1)-\psi(1)}{\gamma}\right) + \left(\frac{-\psi'(\alpha+1)+\psi'(1)}{\gamma^2}\right)\lambda^2$

Where $\psi'(\cdot)$ is the trigamma function.

To ensure that $\mu=e^{X\beta}$, $\lambda$ is replaced with:

$\lambda=\frac{\gamma e^{X\beta}}{\psi(\alpha+1)-\psi(1)}$

This results in:

$f(y|\mu,\alpha,\beta)=\int_0^\infty \frac{\left(\frac{\gamma e^{X\beta}}{\psi(\alpha+1)-\psi(1)}\right)^y x^y e^{-\left(\frac{\gamma e^{X\beta}}{\psi(\alpha+1)-\psi(1)}\right) x}}{y!}\frac{\alpha}{\gamma} \left(1-e^{-\frac{x}{\gamma}}\right)^{\alpha-1}e^{-\frac{x}{\gamma}} dx$

Halton draws are used to perform simulation over the lognormal distribution to solve the integral.

Poisson-Inverse-Gaussian Type 1 (PIG1) and Type 2 (PIG2) Models

The Poisson-Inverse-Gaussian regression model is based on the Poisson-Inverse-Gaussian Distribution.

The expected value of the distribution in the regression utilizes a log-link function. Thus, the mean is:

$\mu=e^{X\beta}$

The variance function for the Type 1 distribution (which is the default) is:

$\sigma^2=\mu+\eta\mu$

While the variance for the Type 2 distribution is:

$\sigma^2=\mu+\eta\mu^2$

The parameter $\eta$ is estimated as the natural logarithm transformed value, $\ln(\eta)$, to ensure that $\eta>0$.

Poisson-Inverse-Gamma (PIG) Model

The PDF of the distribution is:

$f(x|\eta,\mu)=\frac{2\left(\mu\left(\frac{1}{\eta}+1\right)\right) ^{\frac{x+\frac{1}{\eta}+2}{2}}}{x!\Gamma\left(\frac{1}{\eta}+2\right)} K_{x-\frac{1}{\eta}-2}\left(2\sqrt{\mu\left(\frac{1}{\eta}+1\right)}\right)$

Where $\eta$ is a shape parameter with the restriction that $\eta>0$, $\mu>0$ is the mean value, $y$ is a non-negative integer, and $K_i(z)$ is the modified Bessel function of the second kind. This formulation uses the mean directly.

The variance of the distribution is:

$\sigma^2=\mu+\eta\mu^2$

Poisson-Lindley (PL) Model

The Poisson-Lindley regression is based on a compound Poisson-Lindley distribution. It handles count outcomes with high levels of zero observations (or other high densities at low outcome values) that standard count regression methods, including the negative binomial, may struggle to adequately capture or model.

The compound Probability Mass Function(PMF) for the Poisson-Lindley (PL) distribution is:

$f(y|\theta,\lambda)=\frac{\theta^2\lambda^y(\theta+\lambda+ y+1)}{(\theta+1)(\theta+\lambda)^{y+2}}$

Where $\theta$ and $\lambda$ are distribution parameters with the restrictions that $\theta>0$ and $\lambda>0$, and $y$ is a non-negative integer.

The expected value of the distribution is:

$\mu=\frac{\lambda(\theta+2)}{\theta(\theta+1)}$

If a log-link function is used, the mean is:

$\mu=e^{X\beta}=\frac{\lambda(\theta+2)}{\theta(\theta+1)}$

Thus, the parameter $\lambda$ in the PL distribution when applied to regression analysis is:

$\lambda=\frac{\mu\theta(\theta+1)}{\theta+2}$

Using the replacement and simplifying results in:

$f(y \mid \theta, \mu) = \\ \frac{\theta^2 (\mu \theta (\theta+1))^y (\theta^2 (1+\mu) + \theta (2+\mu) + (\theta+2) (y+1))}{(\theta+1) (\theta+2)^{y+1} (\theta^2 (1+\mu) + \theta (2+\mu))^{y+2}}$

And

$LL=2 \log(\theta) + y (\log(\mu) + \log(\theta) + \log(\theta+1)) + \log(\theta^2 (1+\mu) + \theta (2+\mu) + (\theta+2) (y+1)) - \log(\theta+1) - (y+1) \log(\theta+2) - (y+2) \log(\theta^2 (1+\mu) + \theta (2+\mu))$

The variance function is defined as:

$\sigma^2=\mu+\left(1-\frac{2}{(\theta+2)^2}\right)\mu^2$

It should be noted that the p-value for the parameter ln(theta) in the model summary is testing if the parameter theta is equal to a value of 1. This has no practical meaning. The Likelihood-Ratio (LR) test compares the Poisson-Lindley regression with a Poisson regression with the same independent variables. Thus, the PR test result indicates the statistical significance for the improvement in how well the model fits the data over a Poisson regression. This indicates the statistical significance of the theta parameter.

Poisson-Lindley-Gamma (PLG) Model

The Poisson-Lindley-Gamma regression is based on a compound Poisson-Lindley- Gamma distribution. Details of the distribution can be seen at dplindGamma.

The mean for the regression model is:

$\mu=e^{X\beta}$

The variance function is defined as:

$\sigma^2=\mu+\left(2\alpha+1-\frac{2(1+\alpha)} {(\theta+2)^2}\right)\mu^2$

It should be noted that the p-value for the parameters ln(theta) and ln(alpha) in the model summary are testing if the parameter theta and alpha are equal to a value of 1.

Poisson-Lindley-Lognormal (PLL) Model

The Poisson-Lindley-Lognormal regression is based on a compound Poisson- Lindley-Lognormal distribution. Details of the distribution can be seen at dplindLnorm.

The mean for the regression model is:

$\mu=e^{X\beta}$

The variance function is defined as:

$\sigma^2=\mu+\left(\frac{1-\frac{2}{(\theta+2)^2}}{e^{\frac{\sigma^2} {2}}}+e^{\sigma^2}-1\right)\mu^2$

It should be noted that the p-value for the parameters ln(theta) and ln(sigma) in the model summary are testing if the parameter theta and sigma are equal to a value of 1.

Poisson-Weibull (PW) Model

The Poisson-Weibull distribution uses the Weibull distribution as a mixing distribution for a Poisson process. It is useful for modeling overdispersed count data. The density function (probability mass function) for the Poisson-Weibull distribution is given by:

$P(y|\lambda,\alpha,\sigma) = \int_0^\infty \frac{e^{-\lambda x} \lambda^y x^y }{y!} \left(\frac{\alpha}{\sigma}\right)\left(\frac{x}{\sigma} \right)^{\alpha-1}e^{-\left(\frac{x}{\sigma}\right)^\alpha} dx$

where $f(x| \alpha, \sigma)$ is the PDF of the Weibull distribution and $\lambda$ is the mean of the Poisson distribution.

For the Poisson-Weibull Regression model, the expected values is:

$E[Y]=\lambda\sigma\Gamma\left(1+\frac{1}{\alpha}\right)$

Where $\lambda$ is the mean of the Poisson distribution, $\alpha$ is the shape parameter, and $\sigma$ is the scale parameter.

To ensure that the regression model predicts the mean value, the regression utilizes:

$\mu=\exp{X\gamma}=\lambda\sigma\Gamma\left(1+\frac{1}{\alpha}\right)$

Where $X$ is a matrix of independent variables and $\gamma$ is a vector of coefficients.

This leads to:

$\lambda=\frac{\mu}{\sigma\Gamma\left(1+\frac{1}{\alpha}\right)}$

The variance for the Poisson-Weibull regression is:

$V[Y]=\mu+\left(\frac{\Gamma\left(1+\frac{2}{\alpha}\right)} {\Gamma\left(1+\frac{1}{\alpha}\right)^2}-1\right)\mu^2$

Sichel (SI) Model

The compound Probability Mass Function (PMF) for the Sichel distribution uses the formulation from Zhou et al. (2011) and Rigby et al. (2008):

$f(y|\mu, \sigma, \gamma)=\frac{\left(\frac{\mu}{c}\right)^y K_{y+\gamma}(\alpha)}{K_\gamma(1/\sigma)y!(\alpha\sigma)^{y+\gamma}}$

Where $\sigma$ and $\gamma$ are distribution parameters with $-\infty < \gamma < \infty$ and $\sigma>0$, $c=\frac{K_{\gamma+1}(1/\sigma)}{K_\gamma(1/\sigma)}$, $\alpha^2=\sigma^{-2}+2\mu(c\sigma)^{-1}$, a mean value of $\mu$, $y$ is a non-negative integer, and $K_j(x)$ is a modified Bessel function of the third kind with order $j$ and argument $x$.

The variance of the distribution is:

$\sigma^2=\mu+\left(\frac{2\sigma(\gamma+1)}{c}+\frac{1}{c^2}-1\right) \mu^2$

Generalized Waring (GW) Model

The following are the versions of the PMF, mean, and variance used for the Generalized Waring model. This is adjusted from the typical formulation by replacing parameter k with $\mu$

$PMF=\frac{\Gamma(\alpha+y)\Gamma(k+y)\Gamma(\rho+k)\Gamma(\alpha+\rho)} {y!\Gamma(\alpha)\Gamma(k)\Gamma(\rho)\Gamma(\alpha+k+\rho+y)}$

$\mu=e^{X\beta}=\frac{\alpha k}{\rho-1}$

$\sigma^2=\frac{\alpha k(\alpha+k+\rho-1)}{(\rho-1)^2(\rho-2)}$

The distribution parameters are often considered to capture the randomness (parameter

$\alpha$

), proneness (parameter k), and liability (parameter

$\rho$

) of the data.

If we use:

$\alpha=\frac{\mu k}{\rho-1}$

The PMF becomes:

$PMF=\frac{\Gamma\left(\frac{\mu k}{\rho-1}+y\right)\Gamma(k+y) \Gamma(\rho+k)\Gamma\left(\frac{\mu k}{\rho-1}+\rho\right)}{y! \Gamma\left(\frac{\mu k}{\rho-1}\right)\Gamma(k)\Gamma(\rho) \Gamma\left(\frac{\mu k}{\rho-1}+k+\rho+y\right)}$

This results in a regression model where:

$\mu=e^{X\beta}$

$\sigma^2=\frac{\mu k^2\left(\frac{\mu k}{\rho-1}+k+\rho-1\right)} {(\rho-1)^3(\rho-2)}=\left(\frac{k^3+\rho k^2- k^2}{(\rho-1)^3(\rho-2)} \right)\mu+\left(\frac{k^3}{(\rho-1)^4(\rho-2)}\right)\mu^2$

Note that when

$p=1$

or

$p=2$

, the distribution is undefined.

Conway-Maxwell-Poisson (COM) Model

The following is the the PMF for the COM model.

$f(x|\lambda, \nu)=\frac{\lambda^x}{(x!)^\nu Z(\lambda,\nu)}$

Where $\lambda$ and $\nu$ are distribution parameters with $\lambda>0$ and $\nu>0$, and $Z(\lambda,\nu)$ is the normalizing constant.

The normalizing constant is given by:

$Z(\lambda,\nu)=\sum_{n=0}^{\infty}\frac{\lambda^n}{(n!)^\nu}$

The mean and variance are:

$\mu=e^{X\beta}=\lambda \frac{\delta}{\delta \lambda} \log(Z(\lambda, \nu))$

$\sigma^2=\lambda \frac{\delta}{\delta \lambda} \mu$

Note that the COM distribution parameter $\lambda$ is solved for using $\mu$ and $\nu$, so the regression model provides direct predictions for the mean.

Underreporting

Models for underreporting combine a binary probability model (logit or probit) with a count model. This is accomplished using a model for the probability of crashes being reported multiplied by the estimated mean for the count model, based on the observed data. This is discussed in Wood et. al. (2016), Pararai et. al., (2006), and Pararai et. al., (2010). The underreporting model is based on:

$\mu_{true}=\mu_{observed}\cdot P(\text{event is reported})$

This allows the inference of both the true event count and the probability of the event being reported as a function of independent variables.

References

Greene, W. (2008). Functional forms for the negative binomial model for count data. Economics Letters, 99(3), 585-590.

Pararai, M., Famoye, F., & Lee, C. (2006). Generalized Poisson regression model for underreported counts. Advances and applications in Statistics, 6(3), 305-322.

Pararai, M., Famoye, F., & Lee, C. (2010). Generalized Poisson-Poisson mixture model for misreported counts with an application to smoking data. Journal of Data Science, 8(4), 607-617.

Rigby, R. A., Stasinopoulos, D. M., & Akantziliotou, C. (2008). A framework for modelling overdispersed count data, including the Poisson-shifted generalized inverse Gaussian distribution. Computational Statistics & Data Analysis, 53(2), 381-393.

Wood, J.S., Eric T. Donnell, & Christopher J. Fariss. "A method to account for and estimate underreporting in crash frequency research." Accident Analysis & Prevention 95 (2016): 57-66.

Zou, Y., Lord, D., & Zhang, Y. (2012). Analyzing highly dispersed crash data using the Sichel generalized additive models for location, scale and shape.

Examples

# Load the Washington data
data("washington_roads")
washington_roads$AADT10kplus <- ifelse(washington_roads$AADT > 10000, 1, 0)
# Estimate an NB2 model with a dispersion parameter as a function of the 
# variable `speed50` (i.e., generalized NB2), verbose output, and use the 
# BFGS optimization method
nb2 <- countreg(Total_crashes ~ lnaadt + lnlength + speed50 + AADT10kplus,
               data = washington_roads, family = "NB2",
               dis_param_formula_1 = ~ speed50, verbose = TRUE, 
               method='BFGS')
summary(nb2)


# Estimate a Poisson-Lognormal model (a low number of draws is used to speed 
# up the estimation for examples - not recommended in practice)
pln <- countreg(Total_crashes ~ lnaadt + lnlength + speed50 + AADT10kplus,
              data = washington_roads, family = "PLN", engine="poilog")
summary(pln)  

# Estimate an Poisson-Lognormal with underreporting (probit)
plogn_underreport <- countreg(Total_crashes ~ lnaadt + lnlength + speed50 + 
               AADT10kplus,
              data = washington_roads, family = "NB2",
              underreport_formula = ~ speed50 + AADT10kplus, 
              underreport_family = "probit")
summary(plogn_underreport)

# Estimate a Conway-Maxwell-Poisson model
com_model <- countreg(Total_crashes ~ lnaadt + lnlength + speed50 + 
                      AADT10kplus,
                      data = washington_roads, family = "COM", method="BHHH")
summary(com_model)

# Estimate a CGenerlized-Poisson Version 1 model
gp1_model <- countreg(Total_crashes ~ lnaadt + lnlength + speed50 + 
                      AADT10kplus,
                      data = washington_roads, family = "GP1", method="BHHH")
summary(gp1_model)

#' # Estimate a CGenerlized-Poisson Version 2 model
gp2_model <- countreg(Total_crashes ~ lnaadt + lnlength + speed50 + 
                      AADT10kplus,
                      data = washington_roads, family = "GP2", method="BHHH")
summary(gp2_model)

---

Random Parameters Count Regression Models

Description

Random Parameters Count Regression Models

Usage

countreg.rp(
  formula,
  rpar_formula,
  data,
  family = "NB2",
  rpardists = NULL,
  dis_param_formula_1 = NULL,
  dis_param_formula_2 = NULL,
  het_mean_formula = NULL,
  het_var_formula = NULL,
  ndraws = 500,
  scrambled = FALSE,
  correlated = FALSE,
  panel_id = NULL,
  weights = NULL,
  offset = NULL,
  method = "BHHH",
  max.iters = 1000,
  start.vals = NULL,
  verbose = FALSE,
  engine = "halton"
)

Arguments

formula	an R formula. This formula should specify the outcome and the independent variables that have fixed parameters.
rpar_formula	a symbolic description of the model related specifically to the random parameters. This should not include an outcome variable. If the intercept is random, include it in this formula. If the intercept is fixed, include it in formula but not in rpar_formula.
data	a dataframe that has all of the variables in the formula and rpar_formula.
family	the name of the distribution/model type to estimate. Default is "NB2". Options include "Poisson", "NB1", "NB2", "NBP", "PIG", "Sichel", etc. (See countreg for full list).
rpardists	an optional named vector whose names are the random parameters and values the distribution. The distribution options include normal ("n"), lognormal ("ln"), triangular ("t"), uniform ("u"), and gamma ("g"). If not provided, normal is used.
dis_param_formula_1	a symbolic description of the model for the first parameter of the count distribution (e.g., ln(alpha) for NB2).
dis_param_formula_2	a symbolic description of the model for the second parameter of the count distribution (if applicable).
het_mean_formula	an optional symbolic description of the model for heterogeneity in the means of the random parameters.
het_var_formula	an optional symbolic description of the model for heterogeneity in the variances of the random parameters.
ndraws	the number of Halton draws to use for estimating the random parameters.
scrambled	if the Halton draws should be scrambled.
correlated	if the random parameters should be correlated. If TRUE, only normal distributions are used.
panel_id	an optional variable name (string) or vector defining the panel structure (repeated measures). If provided, the standard errors and likelihood are estimated accounting for the panel structure.
weights	variable name to be used as frequency weights.
offset	variable name to be used as an offset.
method	optimization method (e.g., "BHHH", "BFGS", "NM").
max.iters	maximum number of iterations.
start.vals	optional vector of starting values.
verbose	logical.
engine	the engine to use for the integration if the Poisson Lognormal (PLN) distribution is used. Options include "halton" which uses Halton draws or "poilog" which uses the poilog package. The default is "halton".

Value

An object of class countreg which is a list with the following components:

• model: the fitted model object.

• data: the data frame used to fit the model.

• call: the matched call.

• formula: the formula used to fit the model.

Examples

# Load data
data("washington_roads")
washington_roads$AADT10kplus <- ifelse(washington_roads$AADT > 10000, 1, 0)

# 1. Basic Random Parameters Negative Binomial (NB2)
rp_nb2 <- countreg.rp(Total_crashes ~ lnaadt + lnlength,
                      rpar_formula = ~ -1 + speed50,
                      data = washington_roads,
                      family = "NB2",
                      rpardists = c(speed50 = "n"),
                      ndraws = 100,
                      method = "BHHH")
summary(rp_nb2)

# 2. Random Parameters with Panel Structure (if 'site_id' exists)
# rp_panel <- countreg.rp(Total_crashes ~ -1 + lnaadt,
#                        rpar_formula = ~ speed50,
#                        data = washington_roads,
#                        panel_id = "site_id",
#                        family = "NB2",
#                        ndraws = 100)

# 3. Generalized Random Parameters Model with Heterogeneity
rp_gen <- countreg.rp(Total_crashes ~ lnaadt,
                      rpar_formula = ~ -1 + speed50,
                      dis_param_formula_1 = ~ lnlength, 
                      het_mean_formula = ~ AADT10kplus,
                      data = washington_roads,
                      family = "NB2",
                      rpardists = c(speed50 = "n"),
                      ndraws = 100)
summary(rp_gen)

# 4. Random Parameters Poisson Model with panel specification
rp_poisson <- countreg.rp(Total_crashes ~ lnaadt,
                      rpar_formula = ~ -1 + speed50,
                      dis_param_formula_1 = ~ lnlength, 
                      het_mean_formula = ~ AADT10kplus,
                      data = washington_roads,
                      family = "POISSON",
                      rpardists = c(speed50 = "n"),
                      ndraws = 100,
                      panel_id = "ID")
summary(rp_poisson)

---

Cumulative Residuals (CURE) Plot for Count Models

Description

This function generates a Cumulative Residuals (CURE) plot for count models, including those with random parameters, estimated using the flexCountReg package.

Usage

cureplot(
  model,
  data = NULL,
  indvar = NULL,
  method = "Simulated",
  n_resamples = 0,
  ...
)

Arguments

model	A model object estimated using this R package.
data	Optional dataframe. If not provided, the data used to fit the model will be used.
indvar	Optional independent variable name (character string). This is the continuous independent variable to plot the cumulative residuals against. If not provided, the plot will be against the predicted values.
method	Optional parameter to pass to the predict function. This is only used for random parameters models (e.g., "Simulated" or "Individual"). For further details, see predict.flexCountReg.
n_resamples	Number of resamples for potential resampling in the CURE plot confidence bands. Default is 0 (no bands).
...	Additional arguments passed to cure_plot.

Value

A CURE plot generated with cureplots.

Examples

## Example using a Negative Binomial model
data("washington_roads")
washington_roads$AADTover10k <- ifelse(washington_roads$AADT>10000,1,0)

nb_model <- countreg(Total_crashes ~ lnaadt + lnlength + speed50 + 
                            ShouldWidth04 + AADTover10k,
                            data = washington_roads, family = 'nb2', 
                            method = 'NM', max.iters = 500)
                            
# 1. Plot against fitted values (default) with confidence bands
cureplot(nb_model, n_resamples = 20)

# 2. Plot against a specific covariate (e.g., lnlength)
cureplot(nb_model, indvar = "lnlength", n_resamples = 20)

---

One-Parameter Lindley Distribution

Description

Distribution function for the one-parameter Lindley distribution with parameter theta.

Usage

dlindley(x, theta = 1, log = FALSE)

plindley(q, theta = 1, lower.tail = TRUE, log.p = FALSE)

qlindley(p, theta = 1, lower.tail = TRUE, log.p = FALSE)

rlindley(n, theta = 1)

Arguments

x	a single value or vector of positive values.
theta	distribution parameter value. Default is 1.
log, log.p	logical; If TRUE, probabilities p are given as log(p). If FALSE, probabilities p are given directly. Default is FALSE.
q	a single value or vector of quantiles.
lower.tail	logical; If TRUE, (default), $P(X \leq x)$ are returned, otherwise $P(X > x)$ is returned. Default is TRUE.
p	a single value or vector of probabilities.
n	number of random values to generate.

Details

Probability density function (PDF)

$f(x\mid \theta )=\frac{\theta ^{2}}{(1+\theta )}(1+x)e^{-\theta x}$

Cumulative distribution function (CDF)

$F(x\mid \theta ) = 1 - \left(1+ \frac{\theta x}{1+\theta }\right)e^{-\theta x}$

Quantile function (Inverse CDF)

$Q(p\mid\theta) = -1 - \frac{1}{\theta} - \frac{1}{\theta} W_{-1}\!\left((1+\theta)(p-1)e^{-(1+\theta)}\right)$

where $W_{-1}()$ is the negative branch of the Lambert W function.

The moment generating function (MGF) is:

$M_X(t)=\frac{\theta^2(\theta-t+1)}{(\theta+1)(\theta-t)^2}$

The distribution mean and variance are:

$\mu=\frac{\theta+2}{\theta(1+\theta)}$

$\sigma^2=\frac{\mu}{\theta+2}\left(\frac{6}{\theta}-4\right)-\mu^2$

Value

dlindley gives the density, plindley gives the distribution function, qlindley gives the quantile function, and rlindley generates random deviates.

The length of the result is determined by n for rlindley, and is the maximum of the lengths of the numerical arguments for the other functions.

Examples

x <- seq(0, 5, by = 0.1)
p <- seq(0.1, 0.9, by = 0.1)
q <- c(0.2, 3, 0.2)
dlindley(x, theta = 1.5)
dlindley(x, theta=0.5, log=TRUE)
plindley(q, theta = 1.5)
plindley(q, theta = 0.5, lower.tail = FALSE)
qlindley(p, theta = 1.5)
qlindley(p, theta = 0.5)

set.seed(123154)
rlindley(5, theta = 1.5)
rlindley(5, theta = 0.5)

---

flexCountReg Class

Description

A class to represent various objects created by the flexCountReg package.

Slots

model

The fitted model object (can be any type).

data

The data used for fitting the model (data frame).

formula

The R formula used in the regression model (main formula, not including random parameters, etc.)

call

The matched call (language).

additional

Any additional information (list).

---

Generalized Waring Distribution

Description

These functions provide the probability mass function, distribution function, quantile function, and random number generation for the Generalized Waring distribution.

Usage

dgwar(y, mu, k, rho, log = FALSE)

pgwar(q, mu, k, rho, lower.tail = TRUE, log.p = FALSE)

qgwar(p, mu, k, rho)

rgwar(n, mu, k, rho)

Arguments

y	non-negative integer vector of count outcomes.
mu	positive numeric vector of distribution means.
k	positive numeric vector of shape parameters.
rho	numeric vector of shape parameters greater than 1. A finite variance requires $\rho > 2$.
log	logical; if TRUE, log probabilities are returned.
q	non-negative integer vector of quantiles.
lower.tail	logical; if TRUE, probabilities are $P[X \leq q]$; otherwise, $P[X > q]$.
log.p	logical; if TRUE, probabilities are returned as log probabilities.
p	numeric vector of probabilities.
n	integer number of random numbers to generate.

Details

The Generalized Waring distribution is a three-parameter count distribution used to model overdispersed count data.

dgwar computes the probability mass function (PMF) of the Generalized Waring distribution.

pgwar computes the CDF of the Generalized Waring distribution.

qgwar computes the quantile function of the Generalized Waring distribution.

rgwar generates random numbers from the Generalized Waring distribution.

The probability mass function for the Generalized Waring distribution is:

$f(y \mid a_x, k, \rho) = \frac{ \Gamma(a_x + \rho)\Gamma(k + \rho) (a_x)_y(k)_y }{ y!\Gamma(\rho)\Gamma(a_x + k + \rho) (a_x + k + \rho)_y }$

where $(\alpha)_r = \frac{\Gamma(\alpha+r)}{\Gamma(\alpha)}$, and $a_x > 0$, $k > 0$, and $\rho > 0$. Under the mean parameterization used by these functions, $\rho > 1$ is required.

When $\rho > 1$, the mean is:

$E[Y] = \frac{a_x k}{\rho - 1}.$

Therefore, the distribution can be parameterized in terms of its mean using:

$a_x = \frac{\mu(\rho - 1)}{k}.$

For a regression model:

$\mu = \exp(X\beta).$

When $\rho > 2$, the variance under this mean parameterization is:

$\mathrm{Var}(Y) = \frac{\mu(k+\mu)(k+\rho-1)}{k(\rho-2)}.$

Value

dgwar gives the probability mass function, pgwar gives the distribution function, qgwar gives the quantile function, and rgwar generates random deviates.

The length of the result is determined by n for rgwar and is the maximum of the lengths of the numerical arguments for the other functions.

Examples

dgwar(0, mu = 1, k = 2, rho = 3)
pgwar(c(0, 1, 2, 3), mu = 1, k = 2, rho = 3)
qgwar(0.8, mu = 1, k = 2, rho = 3)
rgwar(10, mu = 1, k = 2, rho = 3)

---

Generalized Poisson Version 1 Distribution

Description

These functions provide the density function, distribution function, quantile function, and random number generation for the Generalized Poisson Version 1 (GP-1) Distribution

Usage

dgp1(x, mu = 1, phi = 1, log = FALSE)

pgp1(q, mu = 1, phi = 1, lower.tail = TRUE, log.p = FALSE)

qgp1(p, mu = 1, phi = 1, lower.tail = TRUE, log.p = FALSE)

rgp1(n, mu = 1, phi = 1)

Arguments

x	numeric value or a vector of values.
mu	numeric value or vector of mean values for the distribution (the values have to be greater than 0).
phi	single value or vector of values for the scale parameter of the distribution (the values have to be greater than -1).
log	logical; if TRUE, probabilities p are given as log(p).
q	quantile or a vector of quantiles.
lower.tail	logical; if TRUE, probabilities p are $P[X\leq x]$ otherwise, $P[X>x]$.
log.p	logical; if TRUE, probabilities p are given as log(p).
p	probability or a vector of probabilities.
n	the number of random numbers to generate.

Details

dgp1 computes the density (PDF) of the Generalized Poisson Distribution.

pgp1 computes the CDF of the Generalized Poisson Distribution.

qgp1 computes the quantile function of the Generalized Poisson Distribution.

rgp1 generates random numbers from the Generalized Poisson Distribution.

The compound Probability Mass Function (PMF) for the Generalized Poisson distribution, version 1, (GP-1) is:

$f(y|\phi,\mu)=\frac{\mu(\mu+\phi y)^{y-1} exp\left(-\frac{\mu+\phi y} {1+\phi}\right)}{(1+\phi)^y y!}$

Where $\phi$ is a scale parameter with the restriction that $\mu>0$ is the mean value, and $y$ is a non-negative integer. This formulation uses the mean directly.

The variance of the GP-1 distribution is:

$\sigma^2=(1+\phi)^2 \mu$

If $\phi>0$, the distribution is overdispersed. If $\phi=0$, the distribution is equidispersed. If $\phi<0$, the distribution is underdispersed.

Furthermore, $\phi>-1$ is required for this distribution. When $\phi<0$, there is also a maximum value of support for the integar value $y$. This is $y_max=\left\lfloor -\frac{\mu}{\phi} \right\rfloor$.

Value

dgp1 gives the density, pgp1 gives the distribution function, qgp1 gives the quantile function, and rgp1 generates random deviates.

The length of the result is determined by n for rgp1, and is the maximum of the lengths of the numerical arguments for the other functions.

References

Consul, PoC, and Felix Famoye. "Generalized Poisson regression model." Communications in Statistics-Theory and Methods 21.1 (1992): 89-109.

Zamani, Hossein, and Noriszura Ismail. "Functional form for the generalized Poisson regression model." Communications in Statistics-Theory and Methods 41.20 (2012): 3666-3675.

Examples

dgp1(1, mu=0.75, phi=-0.1)
pgp1(c(0,1,2,3,5,7,9,10), mu=0.75, phi=3)
qgp1(c(0.1,0.3,0.5,0.9,0.95), mu=0.75, phi=0.5)
rgp1(30, mu=0.75, phi=1.5)

---

Generalized Poisson Version 2 Distribution

Description

These functions provide the density function, distribution function, quantile function, and random number generation for the Generalized Poisson Version 2 (GP-2) Distribution.

Usage

dgp2(x, mu = 1, alpha = 0, log = FALSE)

pgp2(q, mu = 1, alpha = 0, lower.tail = TRUE, log.p = FALSE)

qgp2(p, mu = 1, alpha = 0, lower.tail = TRUE, log.p = FALSE)

rgp2(n, mu = 1, alpha = 0)

Arguments

x	numeric value or a vector of values.
mu	numeric value or vector of mean values for the distribution (the values have to be greater than 0).
alpha	numeric value or vector of values for the dispersion parameter of the distribution. Values may be negative, zero, or positive, provided that $1 + \alpha \mu > 0$.
log	logical; if TRUE, probabilities are given as log(p).
q	quantile or a vector of quantiles.
lower.tail	logical; if TRUE, probabilities p are $P[X\leq x]$ otherwise, $P[X>x]$.
log.p	logical; if TRUE, probabilities p are given as log(p).
p	probability or a vector of probabilities.
n	the number of random numbers to generate.

Details

dgp2 computes the density (PMF) of the Generalized Poisson Version 2 Distribution.

pgp2 computes the CDF of the Generalized Poisson Version 2 Distribution.

qgp2 computes the quantile function of the Generalized Poisson Version 2 Distribution.

rgp2 generates random numbers from the Generalized Poisson Version 2 Distribution.

The probability mass function (PMF) for the Generalized Poisson Version 2 distribution (GP-2) is:

$f(y|\mu,\alpha)= \frac{\mu(\mu+\alpha\mu y)^{y-1} \exp\left(-\frac{\mu+\alpha\mu y}{1+\alpha\mu}\right)} {(1+\alpha\mu)^y y!}$

where $\mu>0$ is the mean and $y$ is a non-negative integer. This formulation uses the mean directly.

The variance of the GP-2 distribution is:

$\sigma^2=\mu(1+\alpha\mu)^2$

If $\alpha > 0$, the distribution is overdispersed. If $\alpha = 0$, the distribution is equidispersed and reduces to the ordinary Poisson distribution. If $\alpha < 0$ and $1 + \alpha\mu > 0$, the distribution is underdispersed and has finite support.

Under the underdispersed case, the largest possible value of $y$ is bounded above by $\left\lceil -1/\alpha \right\rceil - 1$, i.e., the largest integer satisfying $1 + \alpha y > 0$.

Value

dgp2 gives the density, pgp2 gives the distribution function, qgp2 gives the quantile function, and rgp2 generates random deviates.

The length of the result is determined by n for rgp2, and is the maximum of the lengths of the numerical arguments for the other functions.

References

Consul, PoC, and Felix Famoye. "Generalized Poisson regression model." Communications in Statistics-Theory and Methods 21.1 (1992): 89-109.

Wang, Weizhen, and Felix Famoye. "Modeling household fertility decisions with generalized Poisson regression." Journal of Population Economics 10.3 (1997): 273-283.

Yang, Zhao, James W. Hardin, and Cheryl L. Addy. "A score test for overdispersion in Poisson regression based on the generalized Poisson-2 model." Journal of Statistical Planning and Inference 139.4 (2009): 1514-1521.

Harris, Tammy, Zhao Yang, and James W. Hardin. "Modeling underdispersed count data with generalized Poisson regression." Stata Journal 12.4 (2012): 736-747.

Examples

dgp2(1, mu = 0.75, alpha = 0.3)
pgp2(c(0, 1, 2, 3, 5, 7, 9, 10), mu = 0.75, alpha = 0.25)
qgp2(c(0.1, 0.3, 0.5, 0.9, 0.95), mu = 0.75, alpha = -0.2)
rgp2(30, mu = 0.75, alpha = 0.2)

---

Generate pseudo-random draws from specified distributions using Halton draws

Description

Generate pseudo-random draws from specified distributions using Halton draws

Usage

halton_dists(dist, mean, sdev, hdraw = NULL, ndraws = 500)

Arguments

dist	The distribution type to use. The distribution options include normal ("n"), lognormal ("ln"), triangular ("t"), uniform ("u"), and gamma ("g").
mean	The mean value for the random draws.
sdev	The standard deviation value for the random draws.
hdraw	An optional vector of Halton draws to convert to the specified distribution. If not provided, the function will generate Halton draws.
ndraws	The number of random draws to generate. This is only used if hdraw is not provided.

Details

This function is used to convert Halton draws to the specified distribution. The function can be used to generate random draws for use in random parameter models, generating Halton-based pseudo-random draws for specified distributions, etc.

The distributions generated all use the mean (

$\mu$

) and sdev (

$\sigma$

) parameters to generate the random draws. The density functions for the distributions are as follows: The Normal distribution is: $f(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)$

The Lognormal distribution is: $f(x) = \frac{1}{x\sigma\sqrt{2\pi}} \exp\left(-\frac{(\log(x) - \mu)^2}{2\sigma^2}\right)$

The Triangular distribution is (note that this is a symmetrical triangular distribution where

$\mu$

is the median and

$\sigma$

is the half-width): $f(x) = \begin{cases} \frac{(x - \mu + \sigma)}{\sigma^2}, & \text{for } \mu - \sigma \leq x \leq \mu \\ \frac{(\mu + \sigma - x)}{\sigma^2}, & \text{for } \mu < x \leq \mu + \sigma \\0, & \text{otherwise}\end{cases}$

The Uniform distribution is (note that $\mu$ is the midpoint and $\sigma$ is the half-width): $f(x) = \frac{1}{(\beta_{\mu}+\beta_{\sigma}) - (\beta_{\mu}-\beta_{\sigma})}=\frac{1}{2\beta_{\sigma}}$

The Gamma distribution is based on

$\mu = \frac{\alpha}{\beta}$

and

$\sigma^2 = \frac{\alpha}{\beta^2}$

: $f(x) = \frac{\left(\frac{\mu}{\sigma^2}\right)^ {\frac{\mu^2}{\sigma^2}}}{\Gamma\left(\frac{\mu^2}{\sigma^2}\right)} x^{\frac{\mu^2}{\sigma^2} - 1} e^{-\frac{\mu}{\sigma^2} x}$

Value

A vector of psudo-random draws from the specified distribution, based on Halton draws.

Examples

# Generate 500 random draws from a normal distribution 
halton_dists(dist="n", mean=3, sdev=2, ndraws=500)

# Generate 500 random draws from a lognormal distribution 
halton_dists(dist="ln", mean=2, sdev=1.5, ndraws=500)

# Generate 500 random draws from a triangular distribution
halton_dists(dist="t", mean=1, sdev=0.5, ndraws=500)

# Generate 500 random draws from a uniform distribution
halton_dists(dist="u", mean=8, sdev=3, ndraws=500)

# Generate 500 random draws from a gamma distribution
halton_dists(dist="g", mean=0.5, sdev=1.5, ndraws=500)

---