bug in `samep.test()`

[RobinHankin]

reproducible, and not the same thing as the "wmmin not finite" bug:

library("hyper2")
#> Loading required package: magrittr
#> Loading required package: cubature
samep.test(skating,c("slutskaya","kwan","hughes"))
#> Error in constrOptim(theta = startq, f = objective, grad = NULL, ui = UI, : initial value is not in the interior of the feasible region

[RobinHankin]

Sometimes it returns a value and sometimes an error is thrown. I hesitate to say "works" because:

samep.test(skating,c("slutskaya","kwan","hughes"))
## 
##  Constrained support maximization
## 
## data:  skating
## null hypothesis: slutskaya = kwan = hughes
## null estimate:
##   babiakova  butyrskaya       cohen     fontana     giunchi   gusmeroli 
##  3.4434e-06  4.9650e-03  8.6052e-02  3.2119e-04  8.3833e-06  6.4434e-05 
##       hegel      hubert      hughes    kettunen       kopac        kwan 
##  6.9880e-06  2.2718e-04  2.9328e-01  1.1148e-03  2.2230e-06  2.9328e-01 
##   liashenko        luca maniachenko       meier        onda    robinson 
##  4.5101e-04  5.0343e-07  7.4024e-04  2.5567e-04  4.2841e-04  5.0280e-03 
##   sebestyen   slutskaya   soldatova      suguri   volchkova 
##  2.4650e-03  2.9328e-01  8.9931e-06  1.7130e-02  8.9147e-04 
## (argmax, constrained optimization)
## Support for null:  -227.89 + K
## 
## alternative hypothesis:  sum p_i=1 
## alternative estimate:
##   babiakova  butyrskaya       cohen     fontana     giunchi   gusmeroli 
##  4.3112e-06  4.8978e-03  8.6031e-02  3.3198e-04  9.1760e-06  7.1191e-05 
##       hegel      hubert      hughes    kettunen       kopac        kwan 
##  8.2830e-06  2.4259e-04  3.0212e-01  1.2563e-03  2.6104e-06  2.6156e-01 
##   liashenko        luca maniachenko       meier        onda    robinson 
##  4.7050e-04  1.0000e-06  7.3469e-04  2.6326e-04  4.4246e-04  4.8840e-03 
##   sebestyen   slutskaya   soldatova      suguri   volchkova 
##  2.3920e-03  3.1616e-01  1.0698e-05  1.7073e-02  1.0381e-03 
## (argmax, free optimization)
## Support for alternative:  -229.22 + K
## 
## degrees of freedom: 2
## support difference = -1.3235
## p-value: 1

Above we see an impossible situation in which the constrained likelihood is higher than the unconstrained likelihood.

[RobinHankin]

It's possible to ameliorate the problem a tiny bit by using the (null) evaluate as a starting point:

library("hyper2")
#> Loading required package: cubature
samep.test(skating, c("slutskaya","kwan","hughes"),startp=skating_maxp) 
#> 
#>  Constrained support maximization
#> 
#> data:  skating
#> null hypothesis: slutskaya = kwan = hughes
#> null estimate:
#>    babiakova   butyrskaya        cohen      fontana      giunchi    gusmeroli 
#> 5.047230e-06 5.304464e-03 9.024140e-02 3.683279e-04 1.137324e-05 8.316906e-05 
#>        hegel       hubert       hughes     kettunen        kopac         kwan 
#> 9.785622e-06 2.689828e-04 2.909049e-01 1.376411e-03 2.996664e-06 2.909049e-01 
#>    liashenko         luca  maniachenko        meier         onda     robinson 
#> 5.188620e-04 1.000001e-06 8.068784e-04 2.928086e-04 4.880185e-04 5.308256e-03 
#>    sebestyen    slutskaya    soldatova       suguri    volchkova 
#> 2.588300e-03 2.909049e-01 1.269274e-05 1.845713e-02 1.139378e-03 
#> (argmax, constrained optimization)
#> Support for null:  -229.341 + K
#> 
#> alternative hypothesis:  sum p_i=1 
#> alternative estimate:
#>    babiakova   butyrskaya        cohen      fontana      giunchi    gusmeroli 
#> 5.047235e-06 5.304470e-03 9.024149e-02 3.683282e-04 1.137326e-05 8.316914e-05 
#>        hegel       hubert       hughes     kettunen        kopac         kwan 
#> 9.785631e-06 2.689831e-04 2.983166e-01 1.376412e-03 2.996667e-06 2.655059e-01 
#>    liashenko         luca  maniachenko        meier         onda     robinson 
#> 5.188626e-04 1.000000e-06 8.068792e-04 2.928089e-04 4.880190e-04 5.308262e-03 
#>    sebestyen    slutskaya    soldatova       suguri    volchkova 
#> 2.588302e-03 3.088930e-01 1.269275e-05 1.845715e-02 1.138379e-03 
#> (argmax, free optimization)
#> Support for alternative:  -229.2836 + K
#> 
#> degrees of freedom: 2
#> support difference = 0.05737755
#> p-value: 0.9442375

^{Created on 2024-03-07 with reprex v2.1.0}

We see that the unconstrained optimization is unable to find a better estimate, but at least it's not worse.