Differing model fits pre/post mrds 3.0.0 to binned data

[erex]

Distance issue DistanceDevelopment/Distance#191 submitted by a user suggests "new" Distance package produces different model fit than "old".

The AIC score from the "new" package seems to be larger than the score from previous package.

The data are binned into four distance categories. Could that be the cause?

[fpetersma]

I am not sure, but one thing I am noticing is the detection function seemingly increasing between 170 and 190, thus not being strictly monotonic. Whether the optimiser and check.mono() function pick this up depends on their settings (e.g., which points to evaluate and what the tolerance is) which have been modified in the latest release. It could have something to do with that, maybe? Would be interesting to see what happens when monotonicity="weak".

[fpetersma]

Another thing that the user could try is setting mono_method="solnp" to check whether it has anything to do with the solnp solver vs slsqp.

[erex]

Thanks @fpetersma. User is asking permission to pass the data on to us. Did the test suite contain binned data, particularly with only 4 points of support? Could that induce strange behaviour?

I'll keep you posted.

[Globicephala]

@erex @fpetersma I'm pasting here my response on the original issue, so that there is visibility in both:

@erex Thank you for the help. I've been following the other issue that you raised. Here are the results.

mono_method:

With this change in the argument, the results approach what was done with the previous versions of the packages.

Monotonicity:

As you can see, if you compare with the previous functions I've pasted, it gives a different result from either.

Regarding sharing data, Tiago told me to wait a bit untill he returns from the US, so he can also take a look into this properly. Now, considering all these results, what would you recommend me to o for now? Obviously this will most likely be easier if/when you have the dataset, but I'm all ears.

[fpetersma]

Sorry for the delay on my end. I had another look at it hoping that I would find a small bug somewhere that explained the difference. Sadly, I did not. A bit of a technical answer, but hopefully it helps.

I ran the following line to get a better understanding of why slsqp did not find the best solution, with debug_level = 3.

fit <- ds(obs_data, nadj = 1, mono_method = "slsqp", optimizer = "R", debug_level = 3)

What this showed is that slsqp very quickly finds the solution that it eventually returns (nll = 1036.83) and after that finds many potential solutions that are almost perfect, except for the fact that the solutions just do not meet the monotonicity constraint. I am not sure what solnp() does better.

If you supply slsqp with the solution that solnp finds for nadj = 1, it finds pretty much the same solution as solnp. So, for some reason, it struggles to get to the point that solnp seems to be able to find, but once it is there it does recognise it as the best solution to the problem.

The most likely explanation I can think off for this is simply that these two search algorithms are different, and that this is one of those cases in which one perform significantly better than the other. This is a problem when solving difficult functions in general, which is why it can be best practice to test a solution with more than one solver if multiple solvers are available (which is the case here).

Another explanation (less likely but potentially worth exploring at some point in the future) is that there is a small bug hidden somewhere in mrds related to the implementation of the slsqp algorithm.

Hope this helps a bit, and happy holidays!

[lenthomas]

Hi @fpetersma, This is really helpful, thanks for looking into it. I also had a quick look, and one thing I struggled to understand was what the g(x) line is reporting in the debug output (example below). I assume it's the log-transformed value of g(x) at the 20 constraint points, but I was not sure of the ordering of the points. Perhaps we can have a quick chat in the new year about all this? I'll drop you an email...

DEBUG: initial values = 4.049167 0 
par =  4.049167 1e-16 
nll =  1480.917 
par =  4.049167 1e-16 
nll =  1480.917 
iteration: 1
	x = (4.049167, 0.000000)
par =  4.049167 1e-16 
nll =  1480.917 
	f(x) = 1480.917423
	g(x) = (-0.013589, -0.090090, -0.205183, -0.274548, -0.235771, -0.128758, -0.042902, -0.008253, -0.000861, -0.000045, -0.986411, -0.896321, -0.691138, -0.416591, -0.180820, -0.052062, -0.009160, -0.000907, -0.000047, -0.000001)
iteration: 2
	x = (5.271586, 1.000000)
par =  5.271586 1 
nll =  986.6635 
	f(x) = 986.663470
	g(x) = (-0.011012, -0.074612, -0.181863, -0.277767, -0.284777, -0.160522, 0.037605, 0.166768, 0.124441, -0.033059, -0.988988, -0.914377, -0.732514, -0.454746, -0.169969, -0.009447, -0.047052, -0.213820, -0.338261, -0.305202)
iteration: 3
	x = (5.413112, 0.127111)
par =  5.413112 0.1271109 
nll =  1036.833 
	f(x) = 1036.832940
	g(x) = (-0.003111, -0.021239, -0.053146, -0.086626, -0.103568, -0.091618, -0.059074, -0.035643, -0.048005, -0.089042, -0.996889, -0.975650, -0.922504, -0.835878, -0.732310, -0.640692, -0.581618, -0.545975, -0.497969, -0.408928)
iteration: 4
	x = (5.335655, 0.604844)
par =  5.335655 0.604844 
nll =  956.4749 
	f(x) = 956.474888
	g(x) = (-0.008451, -0.057350, -0.140535, -0.217431, -0.229707, -0.142667, 0.005499, 0.108305, 0.080293, -0.045930, -0.991549, -0.934199, -0.793664, -0.576233, -0.346526, -0.203859, -0.209358, -0.317663, -0.397956, -0.352026)

[Globicephala]

Hi, all.

Thank you for looking into this issue during this time of the year. I'll be showing your inputs to Tiago after New Year's, and share his input here. In the meantime, if there is any other analysis I can try to run with the data I provided to you let me know.

Again, thanks @erex @fpetersma and @lenthomas for your availability!