Semi-parametric estimation of interval duration distributions

[JoKra1]

Hey @GWeindel ,

This PR outlines how to perform semi-parametric estimation of the probability mass functions for the peak-2-peak interval durations.

Motivation:

• We want zero durations for the first and last interval, but only the Weibull handles this gracefully and you mentioned that it is not really flexible in simulations.
• We want flexibility in the distributions. Not always can we zoom in on a suitable candidate in advance/is the same distribution going to work for all intervals.
• Related to 2: Picking a distribution can be a daunting task for the applied user who just wants to run this but might hesitate because the choice can sometimes feel arbitrary.

Proposed Solution:
So, if parametric distributions don't do it for us, we let the estimation routine figure out the appropriate distribution itself. That is, rather than estimating the scale (or shape and scale) for a given parametric distribution, we estimate the probabilities $p(\tau_j=d)$ instead, where $\tau_j$ is the duration of interval $j$. The downside is that we have to estimate up to max_dur parameters per interval/group combination, but we can apply regularization to still make this work.

Comparisons to existing approaches:
This is essentially the most complex/flexible approach that we can take. So far, we have the most restrictive implemented (just fit scale). #267 allows to estimate both scale and shape. This PR estimates the entire distribution. Funny: while this is the most complex approach, the math is actually simpler and borrows from the M step in #267.

Current issues:

• I probably messed up some things with the paramter initialization but it works fine for the tutorial 1 example. We will need to clean this up. I updated it to show you the new functionality.
• We need to think about how to check for convergence appropriately. I implemented a rough first criterion, but it's not ideal.
• Again this is much slower than the current approach, but we can parallelize the computations in estim_d_probs.
• I have not updated the doc strings yet. Will do that if you think this is something worthwhile to add.

Math:
Like I said, the following borrows a lot from the M step discussed in #267. Let $p(\tau_{s,j} = d | \mathbf{C}, \boldsymbol{\theta})$ again denote the posterior probabilities over the latent variables; the duration $\tau$ of peak-2-peak interval $j$ on trial $s$ given data $\mathbf{C}$ and current parameters $\boldsymbol{\theta}$. I outline how to compute them in #267 (as mentioned, the notation incorrectly works with stages but the math stays the same). Now, we are interested in some estimate of $p(\tau_j=d | \mathbf{C}, \boldsymbol{\theta}^k)$, the probability of the duration $\tau$ of peak-2-peak interval on any given trial given the data and current ($k$ is the iteration of the EM algorithm) parameter estimates. We can then use $p(\tau_j=d | \mathbf{C}, \boldsymbol{\theta}^{k-1})$ as the pmf for interval durations in the E step of the EM algorithm to estimate $p(\tau_j=d | \mathbf{C}, \boldsymbol{\theta}^k)$; so the probabilities themselves - from the previous iteration - are part of $\boldsymbol{\theta}$.

$p(\tau_j=d | \mathbf{C}, \boldsymbol{\theta})$ is sort of the optimal distribution we can think of. It is technically a Bayesian posterior distribution, but I will use a Frequentist approach to derive an estimate for the probabilities (sue me). From that perspective, $p(\tau_j=d)$, and I have dropped the conditioning to kind of emphasize the Frequentist take, is just the limiting behavior of the proportion of independent trials on which interval $j$ has a duration $\tau_j=d$ if we were to sample an infinite number of trials. That immediately suggests the following estimate

$$\hat{p}_{\boldsymbol{\theta}^{k}}(\tau_j=d) = \frac{\sum_{s=1}^Sp(\tau_{s,j}=d|\mathbf{C},\boldsymbol{\theta}^{k-1})}{\sum_{d'=1}^D\sum_{s=1}^Sp(\tau_{s,j}=d'|\mathbf{C},\boldsymbol{\theta}^{k-1})}$$

where $S$ is the total number of series. First note, that the denominator simplifies to $S$, since $\sum_{d'=1}^Dp(\tau_{s,j}=d'|\mathbf{C},\boldsymbol{\theta}^{k-1})=1\ \forall\ s$. Then note, that the numerator is essentially an estimate of the expected number of trials on which $\tau_j=d$ so that the ratio becomes a valid estimate of the proportion. Note also, that because trials are assumed to be mutually independent, the numerator converges to the expected number of trials on which $\tau_j=d$ as $S$ increases (theoreticlally to infinity). This estimate is also essentially the same one used for more general HSMMs (see for example section 2.3.1 of Yu, 2011).

$\hat{p}_{\boldsymbol{\theta}^{k}}(\tau_j=d)$ is a great estimate, but it can be quite noisy and overfits easily, especially when we don't have a lot of data. We can address this problem by introducing bias. Specifically, we can impose that $\hat{p}_{\boldsymbol{\theta}^{k}}(\tau_j=d)$ should only vary smoothly over $d$. In practice, this is enforced by fitting a smoothing spline with smoothing penalty $\lambda$ to the estimates of $\hat{p}_{\boldsymbol{\theta}^{k}}(\tau_j=d)\ \forall\ d \in {0,D}$ where $D$ is the max_dur (followed by re-normalizing). $\lambda$ is a hyper parameter, but we can just set this to a quite high value by default, since we would a priori expect that the pmfs are rather smooth and do not have a lot of sudden dips/bumps. If we have more data, we can lower it.

Current implementation

• Setting model.semi_parametric=True after init triggers semi-parametric estimation of the pmfs. I propose to implement it in a way that the scale (and optionally shape) parameters are still estimated alongside the probabilities. They will not be MLE estimates, but together with the specified distribution we can think of them as the best parametric description of the semi-parametric distribtution estimates.
• model.smooth_spmfs = True and model.smooth_spmfs_lam = 100 control the optional smoothing step, which is enabled and quite strong by default. These parameters can also be changed after init.
• model.spmfs holds the semi-parametric estimates of the distributions after calling fit. See the updated tutorial for a plot - if you like this addition I will clean that up.

Here an image of the true (solid) vs. recovered (dashed) distributions from the tutorial 1 example (for $\lambda=200$):

Happy to discuss this in the meeting with @jelmerborst and @lvmaanen as well!

Cheers,
Josh

[GWeindel]

Hey Joshua,

Looks very cool! Let's discuss that soon indeed.
Can you try to make it run for the data of the 3rd tutorial?

[JoKra1]

Done! I updated tutorials 1-3 and tried to set all models to semi_parametric. For tutorials 1-2 differences are completely negligible. For tutorial 3 things get more interesting. The semi_parametric model finds an additional event when using cumulative for the accuracy condition.

Two possible explanations:

• greater risk for over-fitting with more parameters
• the fast-forward mechanism of the cumulative method actually acts like a regulariser, by constraining the explored solution space, and thus likely already counter-acts this a little bit (sadly to an unknown extent). However, I did not immediately see how to generalise the implemented logic to the semi-parametric distribution estimates. So for now it just falls back on the parametric estimates obtained alongside the semi-parametric distributions. The semi-parametric distributions do however differ quite a lot from the best parametric estimates, so this might be far from optimal.

I also fit eliminative models to the accuracy data and both the original/semi-parametric find essentially an identical 4 intervals + probable duplicate solutions for 5 and 6 intervals. Generally, the likelihood is consistently higher for the estimates proposed here. But also, the likelihood just increases with events, so I was really hoping to do some cross-validation - but I realised that is not actually implemented at the moment. We should really add that. Neither the cumulative nor the eliminative method provide an estimate of the expected loss/how well the estimated number of events generalises - so across experiments/simulations we are guaranteed to overfit eventually. I personally would thus only ever trust event number estimates based on an estimate of the expected loss rather than the empirical loss.

Anyway, that is a different problem. Generally the proposal here works really well. On the real data we see the benefit of not specifying a parametric family - the estimated durations are really quite different between intervals and also not necessarily close to a Gamma. Especially for the last interval we also see the bias introduced by the Gamma with a minimum duration of 1: the semi-parametric estimate attributes a significant mass to zero durations for that interval. The Gamma cannot. The approach might just require cross-validation even more than the current implementation because of the increased number of parameters/risk of over-fitting.

[GWeindel]

I really like the principle, it's also elegant! Shouldn't we think of reducing $\lambda$ with the iterations? Logic being that the smooting should be inversely related to the precision of the estimates based on the magnitudes alone.

[JoKra1]

Hmm, I am not so sure about that. The point of the smoothing (bias) is to ensure better generalisation of the estimated distributions: the un-smoothed estimate can really result in pmfs that we would immediately consider an over-fitted estimate. So it can easily miss the true distributions. See the plot below for the example from tutorial 1:

To be clear: the estimate is not bad! It's actually much better for the data at hand (llk ~ -138 vs. llk ~ -170 for the smoothed one). It achieves this by really clinging to the most plausible interval durations for this data. But with little data, we cannot really be sure that these interval durations are really reflective of what we would find in a different/larger dataset (this should improve with more data). The smoothed estimate does not suffer that problem (as much), and that is precisely the point: if we get $\lambda$ just right we could expect a lower expected loss with the smoothed estimate because of the bias. In other words, with the right degree of smoothing we should be more likely to recover the true distributions (as is evident in the plot I included in my original post).

That said, there is definitely a point in trying to estimate $\lambda$. This is essentially a smoothing model problem and thus at the heart of my own research - so I can tell you from experience that this'll cause you headaches :-). The big problem is which criterion to optimise for $\lambda$. Established criteria are for example generalised cross-validation (gcv) or Bayesian marginal likelihood (REML in the frequentist world). gcv is built into make_smoothing_spline from scipy, so you can actually do this already when setting model.smooth_spmfs_lam=None (Note that the implementation assumes a gaussian likelihood/undertakes least squares fit for the proportions). However, this criterion is well known to under-smooth in the finite data case (anything practical) and for tutorial 1 actually produces the plot I included above.

The more stable route would be REML, which requires some model of the proportions $\hat{p}_{\boldsymbol{\theta}^{k}}(\tau_j=d)$ and for REML the Gaussian/least-squares route is tricky to justify. A suitable option would be to assume that $\hat{p}_{\boldsymbol{\theta}^{k}}(\tau_j=d) \sim B(\alpha_d,\beta)$ where $B$ denotes the beta distribution and $\alpha_d = f(d)$ where $f$ is again a smooth function to be estimated alongside smoothing penalty $\lambda$. Essentially, we assume that the estimated proportions are noisy realisations of a beta distribution with shared $\beta$ (this could be relaxed of course) across durations $d$ and alpha parameter $\alpha_d$ varying smoothly as a function of $d$. The estimated mean of the fitted beta distributions then becomes the smoothed estimate of the proportions.

The resulting smooth functions (most importantly their $\lambda$) can be estimated with mssm. The resulting HMP model has the lowest log-likelihood (~-206) but undeniably achieves the best fit out of all options - in the sense that the estimate is most likely to generalise/closest to the actual distributions:

mssm is GPL-licensed, so I will not commit this version, but we could expose an attribute of the EventModel that accepts any function, passes along the proportions, and accepts smoothed proportions in return. If this argument is not specified, the instance of EventModel simply falls back to the scipy implementation. Then I could make the code available in a standalone repo that users can just download if they want (and anyone could implement their own smoother).