Kolmogorov-Smirnov and AIC giving opposite goodness-of-fit results. Is this possible/surprising/normal?

Question

I have some data on the duration of several activities (rounded to the nearest half hour). I'm trying to add up these random variables (one per activity) so that I can calculate the total duration of a project, as well as extract some summary statistics from it.

In order to do that, I'm trying to determine which distribution fits it reasonably well. This is for two reasons:

• if I can model these RVs by using a "named" distribution, it becomes simpler to combine them
• by fitting a "named" distibution, I assume I'll be able to infer the generator behind the data. For instance, if a lognormal is a good fit for the data, then it might be generated by some sort of random walk process.

I've selected five candidate distributions: lognorm, exponweib, norm, t and dweibull. Three of them are unbounded; I've chosen them just to see if they fit my data reasonably well, even though time durations can't possibly be negative.

I've also selected two criteria with which to judge goodness-of-fit, K-S and AIC. I wanted something that would apply to all of the distributions I selected, in an automated manner. I also wanted a criterion that would penalize more parameters.

However, something apparently strange happened. Here are the distributions' CDFs and the duration data empirical CDF (dashed line):

Judging by this plot, the norm, dweibull and t are all reasonably good fits, which is confimed by their K-S scores.

Now here are the same distributions in PDF form (as well as the data histogram):

Judging by this plot, the lognorm and exponweib are the clear winners, with dweibull a distant third. This is also confirmed by their AIC scores.

Here is the code in Python that calculates the AIC:

def aic(dist, dados, second_order = True):
    fit = dist.fit(dados)
    k = len(fit)
    lnL = dist(*fit).logpdf(dados).sum()
    aic = 2 * k - 2 * lnL
    if second_order:
        n = len(dados)
        aicc_p = 2 * ((k ** 2 + k) / (n - k - 1))
        aic += aicc_p

    return aic

In light of this, I have a few questions:

1. Is this a good strategy (fitting a distribution to the data) for my problem (adding RVs up and extracting summary statistics)?
2. Is it possible/surprising/normal that the K-S and AIC statistics give opposing results when used as goodness-of-fit criteria?
3. Which distribution should I pick? Am I misusing these statistics? Is there some other consideration I should be making before choosing?

Answer 1

1. The AIC is normally a criterion that is used to compare models with different numbers of parameters to fit. It is not really clear exactly how many parameters you have fitted here (chances are 2-3 for most if not all of these, but for example it is not clear where you got the degrees of freedom of the t from; also it seems you have shifted or cut off distributions so that they start at 1/2 where the data start, rather than at zero (like, e.g., lognormal) or at $-\infty$ (like t or normal), so you used these in a nonstandard way and it may even debatable whether the shift to 1/2 should count as an additional parameter in the sense of AIC or not. In fact AIC differences are so striking that they cannot be explained in my opinion by different numbers of parameters. I wonder even whether the AIC was computed correctly (or whether the images reflect precisely what you did), because I haven't seen any such calculations for a truncated normal or t-distribution yet (no idea what dist.fit and len(fit) exactly do in your code, as I don't use python).

2. In any case, even assuming that AIC values are correct, I wouldn't put much trust in them. This is because (as was stated already in the other answer) the data are highly discrete, and density (pdf) is a somewhat unstable concept in the sense that discretisation and changes of low amounts of probability can change densities a lot. The good results of lognorm and exponweib in terms of the density seem to be entirely due to the fact that they correctly capture the peak at 1/2; the quality of their fit at higher values (which is clearly worse than Weibull and somewhat worse than the other two) is in the density-based AIC-computation totally dominated by what goes on at 1/2. And this only works for those two distributions because you let them start at 1/2 rather than zero, which they theoretically do. (Of course this is not wrong in the sense that shifting them to 1/2 still leaves them as well defined distributions, but it would be hard to argue that this captures a real process). In fact it may be that the underlying more continuous real data that probably start at zero do not have a density peak at zero but somewhat higher, which would make a Weibull shape clearly more appropriate compared to lognormal or exponweib, but as the rounded data don't show anything below 1/2, the visible density peak is at 1/2, favouring in somewhat artificial manner the distributions that peak at their origin and then are shifted to 1/2.

3. I recommend to look at K-S p-values, because these give you an indication if and whether these distributions are compatible with the data, which seems to be what you are interested in. This may rule out some (or all) of these distributions, or it may tell you that all are compatible, or it may tell you that some work and some don't, so it wouldn't necessarily pick a clear winner, but that will just be an honest result; it may be that the data cannot distinguish between all or some of these distributions apart from random variation. If you need to pick one that is best, you can still go for the one with lowest K-S distance. Note though that a proper computation of the K-S p-value needs to take into account the number of fitted parameters, and not all existing software does this. This of course brings in the difficulty of parameter counting mentioned for the AIC above (and in fact honest results here should depend on the number of degrees of freedom used in fitting), but I'd think even ignoring this, K-S p-values at least give some orientation, interpreted in a cautious "exploratory" way.

4. An alternative approach would be to not shift or truncate distributions to 1/2, but rather to use them to generate continuous data, and then model applying the rounding/binning procedure as is done for the real data. This will give you probability mass functions starting at 1/2, which you can compare with the one from your real data in a more appropriate way than comparing the data pmf to continuous densities. (Also K-S distances can be computed in this way, and may then give a result more compatible to comparing based on densities/pmfs.)

Answer 2

This seems to be an example of where binning a continuous variable is leading to problems, as it often does. With half-hour bins, over 3/4 of your observations are restricted to the first 3 bins. It's going to be hard to gauge any model fit very well, whether visually or by a statistic related to a fit, with data values restricted this way.

I'm also a bit worried that there seem to be no observed values in the time bin between 0 and 30 minutes in these histograms. If there were any actual times less than 15 minutes, rounding "to the nearest half hour" should lead to times of 0 for those cases.

And if you did round to the nearest half hour, it's not clear whether what you show in the bin from 1 to 1.5 hours, for example, represents actual task durations between 0.75 and 1.25 hours or between 1.25 and 1.75 hours. At least, for visual comparisons against parametric fits, the histogram bars should be shifted so that they are centered on the correct "nearest half hour" values rather than starting or ending at them.

It would be best to get actual durations instead of these highly binned data. If that's not possible, as this is a time-to-event (end-of-task event) analysis, you might consider parametric survival models with the times properly treated as interval-censored. If you do use an interval-censoring approach, make sure to use the correct boundary times for intervals if the current data display is for the "closest" half hour. For example, cases in the bin shown here from 1 to 1.5 hours would be coded to represent durations either between 0.75 and 1.25 hours or 1.25 to 1.75 hours, depending on how you have done the time-rounding.

In a situation like this I would tend to trust the integration provided by the cumulative distribution over whatever you would get from the highly discretized "pdf" (more like a probability mass function here). From that cumulative plot I suspect that the Weibull, with its relationship to extreme-value distributions, will work pretty well as @BruceET suggests in a comment.