lotze / COMPoissonReg

COMPoissonReg R package
GNU General Public License v2.0
9 stars 4 forks source link

Simple fitting? #13

Closed jimmymathews closed 2 years ago

jimmymathews commented 2 years ago

Can this package be used to do basic parameter estimation for a Conway-Maxwell-Poisson distribution? I see lots of packages are doing something quite a bit fancier than this. Something like the following function signature:

x = c(1, 2, 3)
y = c(0.9, 0.05, 0.05)
result = fit(x, y)

result["lambda"]
result["nu"]
...
andrewraim commented 2 years ago

Dear @jimmymathews,

Yes - this package can be used to fit a regression using the Conway-Maxwell-Poisson distribution. See the function glm.cmp. Some examples are given in the package manual, which is available at https://CRAN.R-project.org/package=COMPoissonReg.

jimmymathews commented 2 years ago

Sorry, I looked at glm.cmp before. The main readme refers to arguments:

none of which I have. I'm not trying to do regression; there is no response variable. I am trying to do parameter estimation.

andrewraim commented 2 years ago

I think I understand your question better now. It's possible to fit an independent and identically distributed CMP model (with no predictors) to count observations using glm.cmp. To do this, specify only an intercept as the predictor as follows.

library(COMPoissonReg)
set.seed(1234)

# Generate counts from CMP, for demonstration
n = 1000
lambda_true = 1.2
nu_true = 0.4
y = rcmp(n, lambda_true, nu_true)

# Fit an intercept-only model
fit_out = glm.cmp(y ~ 1)

Here are the results.

R> print(fit_out)
CMP coefficients
              Estimate     SE z.value   p.value
X:(Intercept)   0.1847 0.0485  3.8069 0.0001407
S:(Intercept)  -0.9553 0.0959 -9.9595 2.293e-23
--
Transformed intercept-only parameters
       Estimate     SE
lambda   1.2028 0.0584
nu       0.3847 0.0369
--
Chi-squared test for equidispersion
X^2 = 204.1333, df = 1, p-value = 2.6175e-46
--
Elapsed Sec: 0.26   Sample size: 1000   SEs via Hessian
LogLik: -2014.6274   AIC: 4033.2549   BIC: 4043.0704   
Optimization Method: L-BFGS-B   Converged status: 0
Message: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH

The results for "CMP coefficients" are on the coefficient scale, which you may not be interested in. The "transformed" parameters are shown when the regression formula contains only an intercept. These are maximum likelihood estimates for the parameters lambda and nu, which sound like what you are looking for.

dgbsg commented 1 year ago

How to obtain the maximum likelihood estimation of parameter zero-inflated p , if I generate data in the following way

y = rzicmp(n, lambda_true, nu_true, p_true)
andrewraim commented 1 year ago

Dear @dgbsg,

Is your question specifically about fitting intercept-only ZICMP models? Here are two ways to accomplish that.

# Generate data
set.seed(1234)
n = 400
y = rzicmp(n, lambda = 5, nu = 0.5, p = 0.10)

# Formula interface
out1 = glm.cmp(y ~ 1, formula.nu = ~1, formula.p = ~1)

# Raw interface
X = S = W = matrix(1, n, 1)
out2 = glm.zicmp.raw(y, X, S, W)

Both results look something like this.

R> print(out2)
ZICMP coefficients
    Estimate     SE  z-value   p-value
X:1   1.5188 0.1194  12.7225 4.432e-37
S:1  -0.7556 0.0774  -9.7579 1.706e-22
W:1  -2.0655 0.1582 -13.0529 6.119e-39
--
Transformed intercept-only parameters
       Estimate     SE
lambda   4.5667 0.5452
nu       0.4697 0.0364
p        0.1125 0.0158
--
Chi-squared test for equidispersion
X^2 = 126.3064, df = 1, p-value = 2.6350e-29
--
Elapsed: 0.52 sec   Sample size: 400   raw interface
LogLik: -1349.7863   AIC: 2705.5725   BIC: 2717.5469   
Optimization Method: L-BFGS-B   Converged status: 0
Message: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH