ash_pois mode estimation fails on near-Poisson data

[aksarkar]

Test data is available here https://users.rcc.uchicago.edu/~aksarkar/pois-ex2.Rds

> dat <- readRDS("pois-ex2.Rds")
> ashr::ash_pois(dat$x, dat$s, mixcompdist='halfuniform')$loglik
[1] -223.0521
Warning messages:
1: In estimate_mixprop(data, g, prior, optmethod = optmethod, control = control,  :
  Optimization failed to converge. Results may be unreliable. Try increasing maxiter and rerunning.
2: In estimate_mixprop(data, g, prior, optmethod = optmethod, control = control,  :
  Optimization failed to converge. Results may be unreliable. Try increasing maxiter and rerunning.
3: In estimate_mixprop(data, g, prior, optmethod = optmethod, control = control,  :
  Optimization failed to converge. Results may be unreliable. Try increasing maxiter and rerunning.
4: In estimate_mixprop(data, g, prior, optmethod = optmethod, control = control,  :
  Optimization failed to converge. Results may be unreliable. Try increasing maxiter and rerunning.
5: In estimate_mixprop(data, g, prior, optmethod = optmethod, control = control,  :
  Optimization failed to converge. Results may be unreliable. Try increasing maxiter and rerunning.
> sum(dpois(dat$x, dat$s * sum(dat$x) / sum(dat$s), log=TRUE))
[1] -142.5174
> sessionInfo()
sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /scratch/midway2/aksarkar/miniconda3/envs/scmodes/lib/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] MASS_7.3-50       compiler_3.5.1    Matrix_1.2-14     parallel_3.5.1   
 [5] mixsqp_0.1-97     Rcpp_0.12.18      codetools_0.2-15  SQUAREM_2017.10-1
 [9] pscl_1.5.2        doParallel_1.0.11 grid_3.5.1        iterators_1.0.10 
[13] foreach_1.4.4     truncnorm_1.0-8   ashr_2.2-39       lattice_0.20-35  
> 

[aksarkar]

Estimating the mode assuming g is a point mass, and fixing it in ash does better:

> ashr::ash_pois(dat$x, dat$s, mixcompdist="halfuniform", mode=sum(dat$x) / sum(dat$s))$loglik
[1] -137.1668

[stephens999]

i tried this - i get different results. I will update ashr in a minute to make sure I have latest version, but want to record this before I do.

> library(ashr)
> dat <- readRDS("~/Downloads/pois-ex2.Rds")
> ashr::ash_pois(dat$x, dat$s, mixcompdist='halfuniform')$loglik
[1] -136.6626

The sessionInfo says I have: ashr_2.2-39

[stephens999]

ok, when i went to update it says i have the latest version. What version of ashr do you have? what version of mixsqp?

[stephens999]

oh, sorry i just noticed your sessionInfo.

[stephens999]

can you update your mixsqp? I have mixsqp_0.2-2

[aksarkar]

OK, updating mixsqp fixed this example.