CamDavidsonPilon/lifelines — pull request #57 by spacecowboy (2014-05-31 15:39:51+00:00)
Title: Add concordance index function

This commit includes a function for calculating Harrell's concordance
index, which can be calculated in R using 'hmisc'. The function is
implemented in Fortran, with a small Python wrapper. The reason for this
is that since calculating the C-index is an O(n^2) process, it quickly
becomes unacceptably slow with pure Python. Comparing a pure Python
implementation with the Fortran version, on arrays with length 1000,
Python required 434 ms while Fortran did it in 4.73 ms. So almost a
factor of 100 difference.

As a consequence of the addition of the Fortran module, the setup
script now utilizes numpy's setup function which will handle the
compilation of the native code.

In addition, a small unit test has been added. To be able to run the
unit tests, it is likely necessary to compile the native code first with:

python setup.py build_ext --inplace

I'm not sure how you want to organize the source code, so I opted for
naming the file "_statistics.f90" which compiles to a module "_statistics".
The function inside is then imported and wrapped in "statistics.py". My
thinking is that any "module.py" might have some native code related
to it in a "_module.f90" or "_module.c" file.

As a reference, here is a pure python version of the function:

```python
def concordance_index(event_times, predicted_event_times, event_observed=None):
"""
Calculates the concordance index (C-index) between two series
of event times. The first is the real survival times from
the experimental data, and the other is the predicted survival
times from a model of some kind.

The concordance index is a value between 0 and 1 where,
0.5 is the expected result from random predictions,
1.0 is perfect concordance and,
0.0 is perfect anti-concordance (multiply predictions with -1 to get 1.0)

Parameters:
event_times: a (nx1) array of observed survival times.
predicted_event_times: a (nx1) array of predicted survival times.
event_observed: a (nx1) array of censorship flags, 1 if observed,
0 if not. Default assumes all observed.

Returns:
c-index: a value between 0 and 1.
"""
event_times = np.array(event_times, dtype=float)
predicted_event_times = np.array(predicted_event_times, dtype=float)

if event_observed is None:
event_observed = np.ones(event_times.shape[0], dtype=float)

if event_times.shape != predicted_event_times.shape:
raise ValueError("Event times arrays must have the same shape!")

def valid_comparison(time_a, time_b, event_a, event_b):
"""True if times can be compared."""
if event_a and event_b:
return True
elif event_a and time_a < time_b:
return True
elif event_b and time_b < time_a:
return True
else:
return False

def concordance_value(time_a, time_b, pred_a, pred_b):
if pred_a == pred_b:
# Same as random
return 0.5
elif time_a < time_b and pred_a < pred_b:
return 1.0
elif time_b < time_a and pred_b < pred_a:
return 1.0
else:
return 0.0

paircount = 0.0
csum = 0.0

for a, (time_a, pred_a, event_a) in enumerate(zip(event_times,
predicted_event_times,
event_observed)):
# Don't want to double count
for b in range(a + 1, len(event_times)):
time_b = event_times[b]
pred_b = predicted_event_times[b]
event_b = event_observed[b]

if valid_comparison(time_a, time_b, event_a, event_b):
paircount += 1.0
csum += concordance_value(time_a, time_b, pred_a, pred_b)

return csum / paircount
```