I don't really have another place to write up, so here is a computational problem I keep coming back to: getting individual probs of survival from an array of times and an array of probs

Let me sharpen the problem, I have 2 arrays for probability of survival and the corresponding time. Everything following will be written in Python
>>> pr = [0.9, 0.8, 0.7, 0.5, 0.2, 0.1]
>>> t = [1, 2, 3, 5, 6, 9]
So S(t=1)=0.9, S(t=2)=0.8, and so on.

Show thread

Given an array of individual times
>>> ti = [0.9, 5, 7.5, 6, 9, 10]
we want to get the probability of survival at that time. So our output should look like
>>> [1.0, 0.5, 0.2, 0.2, 0.1, 0.1]
The first element is 1 because S(0)=1 by definition. So, how can we get this array?

Show thread

My first thought is that we can write a loop. We look over each element in `ti`, then use that element to find the last valid index in `t` (for ti=7.5 that would be 6), then use that index to look up the element in `pr`.

Show thread

It might be faster to convert `pr` and `t` to a dictionary (in Python) but this is still expensive as we need to loop. If there are lots of elements in `ti`, this can become far too slow. So the idea is to vectorize this procedure. The following is the best solution I’ve come up with at this point

Show thread

So first, I am going to convert `t`, `ti`, `pr` to NumPy arrays so we can use NumPy to help speed things along. After that
>>> t = np.insert(t, 0, 0)
>>> pr = np.insert(pr, 0, 1)
>>> shift_t = np.append(t, np.max(ti)+1)[1:]
>>> upper = (ti >= t[:, None]).astype(int)
>>> lower = (ti < shift_t[:, None]).astype(int)
>>> t_ind = upper * lower
>>> np.sum(t_ind * pr[:, None], axis=0)

Show thread

So what does this do. Setup is done on lines 1-3. Lines 1-2 tack on S(0)=1 to the start of the `t` and `pr` arrays. Line 3 adds the maximum time we saw in `ti` then drops the first element (zero)

Show thread

The clever part (in my opinion) is the remainder. Line 4 creates a matrix of indicators where the columns are the individuals and rows are whether their value in `ti` was >= the corresponding `t`. I do a similar thing in line 5 but now < the shifted times. When these are multiplied together in 6, we get a matrix where the only non-zero value in a column corresponds to the final time the person was observed

Then in the final step, we multiply that by `pr`. Therefore, we can sum over the rows, as the only non-zero row in this matrix will be the probability at their corresponding individual time. There is probably a better way to do this, but I think it’s a clever vectorization of the problem

Show thread

The censoring events before the first event time in `t`, and censoring times greater than the last event time in `t` were the annoying parts to deal with. The first 3 lines were primarily to setup the input arrays to deal with those issues

Show thread
Sign in to participate in the conversation
Qoto Mastodon

QOTO: Question Others to Teach Ourselves
An inclusive, Academic Freedom, instance
All cultures welcome.
Hate speech and harassment strictly forbidden.