A bug in vcmp - Githubissues

andrewraim commented 3 years ago

An issue was reported with the newly-exposed vcmp function in version v0.7.1 of COMPoissonReg.

library(COMPoissonReg)

## parameter 1
lam = 0.6785115
nu = 0.007096562

# Compare vmp to calculation by truncated sums
vcmp(lam, nu)
var(rcmp(100000, lam, nu))

## parameter 2
lam = 1.199938
nu = 0.4124449
vcmp(lam, nu)

# Compare vmp to calculation by truncated sums
vcmp(lam, nu)
var(rcmp(100000, lam, nu))

The vcmp results for the two given cases are -985.2381 and 912.5801, which are way off the empirical variances (and obviously a negative variance result isn't right). The underlying issue appears to be in the use of the grad.fwd function. It is not producing accurate second derivatives in these cases.

Here is a workaround that uses the numDeriv package, until we can fix the underlying problem in vcmp.

vcmp_alt = function(lambda, nu) 
{
    # Only handle univariate inputs for now
    stopifnot(length(lambda) == 1)
    stopifnot(length(nu) == 1)

    # Make sure these options are set
    hybrid_tol = getOption("COMPoissonReg.hybrid_tol")
    truncate_tol = getOption("COMPoissonReg.truncate_tol")
    if (is.null(hybrid_tol)) {
        options(COMPoissonReg.hybrid_tol = 1e-2)
    }
    if (is.null(truncate_tol)) {
        options(COMPoissonReg.truncate_tol = 1e-6)
    }

    # Variance calculation using derivatives
    ev = ecmp(lambda, nu)
    dd = lambda^2 * numDeriv::hessian(ncmp, lambda, nu = nu, log = TRUE)
    as.numeric(dd + ev)
}

Here is a quick test of vcmp_alt on the inputs above.

library(numDeriv)

lam = 0.6785115
nu = 0.007096562
vcmp_alt(lam, nu)
var(rcmp(100000, lam, nu))

lam = 1.199938
nu = 0.4124449
vcmp_alt(lam, nu)
var(rcmp(100000, lam, nu))

andrewraim commented 3 years ago

The other function that would be affected by this issue is vzicmp. Here is a workaround version of that:

vzicmp_alt = function(lambda, nu, p)
{
    ee  = ecmp(lambda, nu)
    vv  = vcmp_alt(lambda, nu)
    (1-p) * (p*ee^2 + vv)
}

andrewraim commented 3 years ago

This is corrected in the vignette branch and should be included in the next release.

lotze / COMPoissonReg

A bug in vcmp #12