Poisson mode estimation incorrectly handles scaling factor

[aksarkar]

By default, the mode estimation procedure searches over the range of the observations instead of over the range of observation-specific MLEs for the Poisson rate. This leads to a much worse fit on real data.

Test case:

test_that("Mode estimation for pois_lik finds an acceptable solution", {
    set.seed(1)
    # Load example 10X Genomics data
    dat = readRDS("test_pois_data.Rds")
    m0 = ashr::ash(rep(0, nrow(dat)), 1, lik=ashr::lik_pois(dat$x, scale=dat$scale, link="identity"), mode="estimate")
    lam = dat$x / dat$scale
    m1 = ashr::ash(rep(0, nrow(dat)), 1, lik=ashr::lik_pois(dat$x, scale=dat$scale, link="identity"), mode=c(min(lam), max(lam)))
    expect_equal(m0$loglik, m1$loglik, tolerance=1, scale=1)
})

Test data: https://users.rcc.uchicago.edu/~aksarkar/singlecell-ideas/test_pois_data.Rds

Test result:

test_pois.R:85: failure: Mode estimation for pois_lik finds an acceptable solution
m0$loglik not equal to m1$loglik.
1/1 mismatches
[1] -80770 - -43477 == -37293