From 40d76ba84f02fb1aab15a2eb06c53b053cc39758 Mon Sep 17 00:00:00 2001
From: CamDavidsonPilon <cam.davidson.pilon@gmail.com>
Date: Thu, 14 May 2020 11:25:48 -0400
Subject: [PATCH 01/10] adding npmle

---
 lifelines/fitters/npmle.py | 173 +++++++++++++++++++++++++++++++++++++
 1 file changed, 173 insertions(+)
 create mode 100644 lifelines/fitters/npmle.py

diff --git a/lifelines/fitters/npmle.py b/lifelines/fitters/npmle.py
new file mode 100644
index 000000000..07187657c
--- /dev/null
+++ b/lifelines/fitters/npmle.py
@@ -0,0 +1,173 @@
+# -*- coding: utf-8 -*-
+from collections import defaultdict, namedtuple
+import numpy as np
+from numpy.linalg import norm
+import pandas as pd
+
+interval = namedtuple("Interval", ["left", "right"])
+
+
+def E_step_M_step(observation_intervals, p_old, turnball_interval_lookup):
+
+    N = len(observation_intervals)
+    T = p_old.shape[0]
+
+    p_new = np.zeros_like(p_old)
+
+    for observation_interval in observation_intervals:
+        p_temp = np.zeros_like(p_old)
+
+        # find all turnball intervals, t, that are contained in (ol, or). Call this set T
+        # the denominator is sum of p_old[T] probabilities
+        # the numerator is p_old[t]
+        for ix_t in turnball_interval_lookup[observation_interval]:
+            t_p = p_old[ix_t]
+            p_temp[ix_t] = t_p
+
+        p_new = p_new + p_temp / p_temp.sum()
+
+    return p_new / N
+
+
+def create_turnball_intervals(left, right):
+    """
+    TIHI X 10000
+    """
+
+    left = [[l, 0, "1l"] for l in left]
+    right = [[r, 0, "0r"] for r in right]
+
+    for l, r in zip(left, right):
+        if l[0] == r[0]:
+            l[1] -= 0.01
+            r[1] += 0.01
+
+    import copy
+
+    union = sorted(list(left) + list(right))
+
+    """
+    # fix ties
+    for k in range(len(union)-1):
+        e_, e__ = union[k], union[k+1]
+        if e_[2] == "1l" and e__[2] == "0r" and e_[0] == e__[0]:
+            union_[k][1] += 0.01
+    """
+    intervals = []
+
+    for k in range(len(union) - 1):
+        e_, e__ = union[k], union[k + 1]
+        if e_[2] == "1l" and e__[2] == "0r":
+            intervals.append(interval(e_[0], e__[0]))
+
+    return sorted(set(intervals))
+
+
+def is_subset(query_interval, super_interval):
+    return super_interval.left <= query_interval.left and query_interval.right <= super_interval.right
+
+
+def create_turnball_lookup(turnball_intervals, observation_intervals):
+
+    turnball_lookup = defaultdict(set)
+
+    for i, turnball_interval in enumerate(turnball_intervals):
+        # ask: which observations is this t_interval part of?
+        for observation_interval in observation_intervals:
+            # since left and right are sorted by left, we can stop after left > turnball_interval[1] value
+            if observation_interval.left > turnball_interval.right:
+                break
+            if is_subset(turnball_interval, observation_interval):
+                turnball_lookup[observation_interval].add(i)
+
+    return turnball_lookup
+
+
+def check_convergence(p_new, p_old, tol, i, verbose=False):
+    if verbose:
+        print("Iteration %d: norm(p_new - p_old): %.6f" % (i, norm(p_new - p_old)))
+    if norm(p_new - p_old) < tol:
+        return True
+    return False
+
+
+def create_observation_intervals(left, right):
+    return [interval(l, r) for l, r in zip(left, right)]
+
+
+def npmle(left, right, tol=1e-5, verbose=False):
+
+    left, right = np.asarray(left, dtype=float), np.asarray(right, dtype=float)
+    assert left.shape == right.shape
+
+    # sort left, right arrays by (left, right).
+    ix = np.lexsort((right, left))
+    left = left[ix]
+    right = right[ix]
+
+    turnball_intervals = list(create_turnball_intervals(left, right))
+    observation_intervals = create_observation_intervals(left, right)
+    turnball_lookup = create_turnball_lookup(turnball_intervals, sorted(set(observation_intervals)))
+
+    T = len(turnball_intervals)
+
+    converged = False
+
+    # initialize to equal weight
+    p = 1 / T * np.ones(T)
+    i = 0
+    while not converged:
+        i += 1
+        p_new = E_step_M_step(observation_intervals, p, turnball_lookup)
+        converged = check_convergence(p_new, p, tol, i, verbose=verbose)
+        p = p_new
+
+    return p, turnball_intervals
+
+
+def reconstruct_survival_function(probabilities, turnball_intervals, timeline, label="NPMLE"):
+    index = [0.0]
+    values = [1.0]
+
+    for p, interval in zip(probabilities, turnball_intervals):
+        if interval.left != index[-1]:
+            index.append(interval.left)
+            values.append(values[-1])
+
+        if interval.left == interval.right:
+            values[-1] -= p
+        else:
+            index.append(interval.right)
+            values.append(values[-1] - p)
+
+    full_dataframe = pd.DataFrame(index=timeline, columns=[label + "_upper"])
+
+    turnball_dataframe = pd.DataFrame(values, index=index, columns=[label + "_upper"])
+
+    dataframe = full_dataframe.combine_first(turnball_dataframe).ffill()
+    dataframe[label + "_lower"] = dataframe[label + "_upper"].shift(1).fillna(1)
+    return dataframe
+
+
+def npmle_compute_confidence_intervals(left, right, mle_, alpha=0.05, samples=1000):
+    """
+    uses basic bootstrap
+    """
+    left, right = np.asarray(left, dtype=float), np.asarray(right, dtype=float)
+    all_times = np.unique(np.concatenate((left, right, [0])))
+
+    N = left.shape[0]
+
+    bootstrapped_samples = np.empty((all_times.shape[0], samples))
+
+    for i in range(samples):
+        ix = np.random.randint(low=0, high=N, size=N)
+        left_ = left[ix]
+        right_ = right[ix]
+
+        bootstrapped_samples[:, i] = reconstruct_survival_function(*npmle(left_, right_), all_times).values[:, 0]
+
+    return (
+        2 * mle_.squeeze() - pd.Series(np.percentile(bootstrapped_samples, (alpha / 2) * 100, axis=1), index=all_times),
+        2 * mle_.squeeze() - pd.Series(np.percentile(bootstrapped_samples, (1 - alpha / 2) * 100, axis=1), index=all_times),
+    )

From dd34b2749a4904b00eed8cb17cad4e25bab386e4 Mon Sep 17 00:00:00 2001
From: Cameron Davidson-Pilon <cam.davidson.pilon@gmail.com>
Date: Thu, 14 May 2020 11:34:05 -0400
Subject: [PATCH 02/10] Revert "Revert "Npmle""

---
 experiments/working_npmle.py             | 187 +++++++++++++++++++++++
 lifelines/fitters/kaplan_meier_fitter.py |  70 ++++++++-
 lifelines/plotting.py                    |   2 +-
 3 files changed, 257 insertions(+), 2 deletions(-)
 create mode 100644 experiments/working_npmle.py

diff --git a/experiments/working_npmle.py b/experiments/working_npmle.py
new file mode 100644
index 000000000..28cb32774
--- /dev/null
+++ b/experiments/working_npmle.py
@@ -0,0 +1,187 @@
+# -*- coding: utf-8 -*-
+from matplotlib import pyplot as plt
+from collections import defaultdict, namedtuple
+from numpy.linalg import norm
+
+interval = namedtuple("Interval", ["left", "right"])
+
+
+def E_step_M_step(observation_intervals, p_old, turnball_interval_lookup):
+
+    N = len(observation_intervals)
+    T = p_old.shape[0]
+
+    p_new = np.zeros_like(p_old)
+
+    for observation_interval in observation_intervals:
+        p_temp = np.zeros_like(p_old)
+
+        # find all turnball intervals, t, that are contained in (ol, or). Call this set T
+        # the denominator is sum of p_old[T] probabilities
+        # the numerator is p_old[t]
+        for ix_t in turnball_interval_lookup[observation_interval]:
+            t_p = p_old[ix_t]
+            p_temp[ix_t] = t_p
+
+        p_new = p_new + p_temp / p_temp.sum()
+
+    return p_new / N
+
+
+def create_turnball_intervals(left, right):
+    """
+    TIHI X 10000
+    """
+
+    left = [[l, 0, "1l"] for l in left]
+    right = [[r, 0, "0r"] for r in right]
+
+    for l, r in zip(left, right):
+        if l[0] == r[0]:
+            l[1] -= 0.01
+            r[1] += 0.01
+
+    import copy
+
+    union = sorted(list(left) + list(right))
+
+    """
+    # fix ties
+    for k in range(len(union)-1):
+        e_, e__ = union[k], union[k+1]
+        if e_[2] == "1l" and e__[2] == "0r" and e_[0] == e__[0]:
+            union_[k][1] += 0.01
+    """
+    intervals = []
+
+    for k in range(len(union) - 1):
+        e_, e__ = union[k], union[k + 1]
+        if e_[2] == "1l" and e__[2] == "0r":
+            intervals.append(interval(e_[0], e__[0]))
+
+    return sorted(set(intervals))
+
+
+def is_subset(query_interval, super_interval):
+    return super_interval.left <= query_interval.left and query_interval.right <= super_interval.right
+
+
+def create_turnball_lookup(turnball_intervals, observation_intervals):
+
+    turnball_lookup = defaultdict(set)
+
+    for i, turnball_interval in enumerate(turnball_intervals):
+        # ask: which observations is this t_interval part of?
+        for observation_interval in observation_intervals:
+            # since left and right are sorted by left, we can stop after left > turnball_interval[1] value
+            if observation_interval.left > turnball_interval.right:
+                break
+            if is_subset(turnball_interval, observation_interval):
+                turnball_lookup[observation_interval].add(i)
+
+    return turnball_lookup
+
+
+def check_convergence(p_new, p_old):
+    if norm(p_new - p_old) < 1e-4:
+        return True
+    return False
+
+
+def create_observation_intervals(left, right):
+    return [interval(l, r) for l, r in zip(left, right)]
+
+
+def npmle(left, right):
+
+    left, right = np.asarray(left, dtype=float), np.asarray(right, dtype=float)
+    assert left.shape == right.shape
+
+    # sort left, right arrays by (left, right).
+    ix = np.lexsort((right, left))
+    left = left[ix]
+    right = right[ix]
+
+    turnball_intervals = list(create_turnball_intervals(left, right))
+    observation_intervals = create_observation_intervals(left, right)
+    turnball_lookup = create_turnball_lookup(turnball_intervals, sorted(set(observation_intervals)))
+
+    T = len(turnball_intervals)
+
+    converged = False
+
+    # initialize to equal weight
+    p = 1 / T * np.ones(T)
+
+    while not converged:
+
+        p_new = E_step_M_step(observation_intervals, p, turnball_lookup)
+        converged = check_convergence(p_new, p)
+        p = p_new
+
+    return p, turnball_intervals
+
+
+def reconstruct_survival_function(probabilities, turnball_intervals, timeline):
+    index = [0.0]
+    values = [1.0]
+
+    for p, interval in zip(probabilities, turnball_intervals):
+        if interval.left != index[-1]:
+            index.append(interval.left)
+            values.append(values[-1])
+
+        if interval.left == interval.right:
+            values[-1] -= p
+        else:
+            index.append(interval.right)
+            values.append(values[-1] - p)
+
+    full_dataframe = pd.DataFrame(index=timeline, columns=["survival function"])
+
+    turnball_dataframe = pd.DataFrame(values, index=index, columns=["survival function"])
+
+    dataframe = full_dataframe.combine_first(turnball_dataframe).ffill()
+    return dataframe
+
+
+def compute_confidence_intervals(left, right, mle_, alpha=0.05, samples=1000):
+    """
+    uses basic bootstrap
+    """
+    left, right = np.asarray(left, dtype=float), np.asarray(right, dtype=float)
+    all_times = np.unique(np.concatenate((left, right, [np.inf, 0])))
+
+    N = left.shape[0]
+
+    bootstrapped_samples = np.empty((all_times.shape[0], samples))
+
+    for i in range(samples):
+        ix = np.random.randint(low=0, high=N, size=N)
+        left_ = left[ix]
+        right_ = right[ix]
+
+        bootstrapped_samples[:, i] = reconstruct_survival_function(*npmle(left_, right_), all_times).values[:, 0]
+
+    return (
+        2 * mle_.squeeze() - pd.Series(np.percentile(bootstrapped_samples, (alpha / 2) * 100, axis=1), index=all_times),
+        2 * mle_.squeeze() - pd.Series(np.percentile(bootstrapped_samples, (1 - alpha / 2) * 100, axis=1), index=all_times),
+    )
+
+
+from lifelines.datasets import load_diabetes
+
+data = load_diabetes()
+
+
+left = [1, 7, 8, 7, 7, 17, 37, 46, 46, 45]
+right = [7, 8, 10, 16, 14, 100, 44, 100, 100, 100]
+
+# left, right = data['left'], data['right']
+
+results = npmle(left, right)
+
+
+timeline = np.unique(np.concatenate((left, right, [np.inf, 0])))
+df = reconstruct_survival_function(*results, timeline)
+# CIs = compute_confidence_intervals(left, right, df)
diff --git a/lifelines/fitters/kaplan_meier_fitter.py b/lifelines/fitters/kaplan_meier_fitter.py
index facbd39fa..4b2e0d3b5 100644
--- a/lifelines/fitters/kaplan_meier_fitter.py
+++ b/lifelines/fitters/kaplan_meier_fitter.py
@@ -17,8 +17,11 @@
     StatisticalWarning,
     coalesce,
     CensoringType,
+    pass_for_numeric_dtypes_or_raise_array,
+    check_nans_or_infs,
 )
 from lifelines.plotting import loglogs_plot, _plot_estimate
+from lifelines.fitters.npmle import npmle, reconstruct_survival_function, npmle_compute_confidence_intervals
 
 
 class KaplanMeierFitter(UnivariateFitter):
@@ -112,6 +115,62 @@ def fit(
 
         return self._fit(durations, event_observed, timeline, entry, label, alpha, ci_labels, weights)
 
+    @CensoringType.interval_censoring
+    def fit_interval_censoring(
+        self,
+        lower_bound,
+        upper_bound,
+        event_observed=None,
+        timeline=None,
+        label=None,
+        alpha=None,
+        ci_labels=None,
+        show_progress=False,
+        entry=None,
+        weights=None,
+        tol=1e-5,
+    ) -> "KaplanMeierFitter":
+
+        if entry is not None or weights is not None:
+            raise NotImplementedError("entry / weights is not supported yet")
+
+        self.upper_bound = np.atleast_1d(pass_for_numeric_dtypes_or_raise_array(upper_bound))
+        self.lower_bound = np.atleast_1d(pass_for_numeric_dtypes_or_raise_array(lower_bound))
+        check_nans_or_infs(self.lower_bound)
+
+        self.timeline = coalesce(timeline, np.unique(np.concatenate((self.upper_bound, self.lower_bound))))
+
+        if (self.upper_bound < self.lower_bound).any():
+            raise ValueError("All upper_bound times must be greater than or equal to lower_bound times.")
+
+        if event_observed is None:
+            event_observed = self.upper_bound == self.lower_bound
+
+        if ((self.lower_bound == self.upper_bound) != event_observed).any():
+            raise ValueError(
+                "For all rows, lower_bound == upper_bound if and only if event observed = 1 (uncensored). Likewise, lower_bound < upper_bound if and only if event observed = 0 (censored)"
+            )
+
+        self._label = coalesce(label, self._label, "KM_estimate")
+
+        probs, t_intervals = npmle(self.lower_bound, self.upper_bound, verbose=show_progress)
+        self.survival_function_ = reconstruct_survival_function(probs, t_intervals, self.timeline, label=self._label)
+        self.cumulative_density_ = 1 - self.survival_function_
+
+        self._median = median_survival_times(self.survival_function_)
+        self.percentile = functools.partial(qth_survival_time, model_or_survival_function=self.survival_function_)
+
+        """
+        self.confidence_interval_ = npmle_compute_confidence_intervals(self.lower_bound, self.upper_bound, self.survival_function_, self.alpha)
+        self.confidence_interval_survival_function_ = self.confidence_interval_
+        self.confidence_interval_cumulative_density_ = 1 - self.confidence_interval_
+        """
+        # estimation methods
+        self._estimation_method = "survival_function_"
+        self._estimate_name = "survival_function_"
+        self._update_docstrings()
+        return self
+
     @CensoringType.left_censoring
     def fit_left_censoring(
         self, durations, event_observed=None, timeline=None, entry=None, label=None, alpha=None, ci_labels=None, weights=None
@@ -147,6 +206,7 @@ def fit_left_censoring(
           self with new properties like ``survival_function_``, ``plot()``, ``median_survival_time_``
 
         """
+        # left censoring is then defined in CensoringType.is_left_censoring(self)
         return self._fit(durations, event_observed, timeline, entry, label, alpha, ci_labels, weights)
 
     def _fit(
@@ -293,9 +353,17 @@ def cumulative_density_at_times(self, times, label=None) -> pd.Series:
         label = coalesce(label, self._label)
         return pd.Series(1 - self.predict(times), index=_to_1d_array(times), name=label)
 
+    def plot(self, **kwargs):
+        return self.plot_survival_function(**kwargs)
+
     def plot_survival_function(self, **kwargs):
         """Alias of ``plot``"""
-        return _plot_estimate(self, estimate="survival_function_", **kwargs)
+        if not CensoringType.is_interval_censoring(self):
+            return _plot_estimate(self, estimate="survival_function_", **kwargs)
+        else:
+            # hack for now.
+            color = coalesce(kwargs.get("c"), kwargs.get("color"), "k")
+            self.survival_function_.plot(drawstyle="steps", color=color, **kwargs)
 
     def plot_cumulative_density(self, **kwargs):
         """
diff --git a/lifelines/plotting.py b/lifelines/plotting.py
index ec44f1518..720d9dbdc 100644
--- a/lifelines/plotting.py
+++ b/lifelines/plotting.py
@@ -732,7 +732,7 @@ def _plot_estimate(
 ):
 
     """
-    Plots a pretty figure of {0}.{1}
+    Plots a pretty figure of estimates
 
     Matplotlib plot arguments can be passed in inside the kwargs, plus
 

From 56c140ba98f2dc57d1fc18bc54abfb91a36b1a69 Mon Sep 17 00:00:00 2001
From: CamDavidsonPilon <cam.davidson.pilon@gmail.com>
Date: Thu, 14 May 2020 11:35:11 -0400
Subject: [PATCH 03/10] adding AIC

---
 lifelines/fitters/__init__.py                | 13 ++++++
 lifelines/fitters/cox_time_varying_fitter.py | 37 +++++----------
 lifelines/fitters/coxph_fitter.py            |  4 +-
 lifelines/utils/printer.py                   | 47 +++++++++++++++-----
 4 files changed, 62 insertions(+), 39 deletions(-)

diff --git a/lifelines/fitters/__init__.py b/lifelines/fitters/__init__.py
index c17ef1539..ec49248a3 100644
--- a/lifelines/fitters/__init__.py
+++ b/lifelines/fitters/__init__.py
@@ -1214,6 +1214,15 @@ def predict_expectation(self, df, conditional_after=None) -> pd.Series:
         raise NotImplementedError()
 
 
+class SemiParametricRegressionFittter(RegressionFitter):
+    @property
+    def AIC_partial_(self) -> float:
+        """
+        "partial" because the log-likelihood is partial
+        """
+        return -2 * self.log_likelihood_ + 2 * self.params_.shape[0]
+
+
 class ParametricRegressionFitter(RegressionFitter):
 
     _scipy_fit_method = "BFGS"
@@ -2482,6 +2491,10 @@ def concordance_index_(self) -> float:
             return self.concordance_index_
         return self._concordance_index_
 
+    @property
+    def AIC_(self) -> float:
+        return -2 * self.log_likelihood_ + 2 * self.params_.shape[0]
+
 
 class ParametericAFTRegressionFitter(ParametricRegressionFitter):
 
diff --git a/lifelines/fitters/cox_time_varying_fitter.py b/lifelines/fitters/cox_time_varying_fitter.py
index b27fe5657..a54151839 100644
--- a/lifelines/fitters/cox_time_varying_fitter.py
+++ b/lifelines/fitters/cox_time_varying_fitter.py
@@ -16,7 +16,7 @@
 from autograd import numpy as anp
 
 
-from lifelines.fitters import RegressionFitter
+from lifelines.fitters import SemiParametricRegressionFittter
 from lifelines.fitters.mixins import ProportionalHazardMixin
 from lifelines.utils.printer import Printer
 from lifelines.statistics import _chisq_test_p_value, StatisticalResult
@@ -47,7 +47,7 @@
 matrix_axis_0_sum_to_1d_array = lambda m: np.sum(m, 0)
 
 
-class CoxTimeVaryingFitter(RegressionFitter, ProportionalHazardMixin):
+class CoxTimeVaryingFitter(SemiParametricRegressionFittter, ProportionalHazardMixin):
     r"""
     This class implements fitting Cox's time-varying proportional hazard model:
 
@@ -173,9 +173,7 @@ def fit(
 
         if weights_col is None:
             self.weights_col = None
-            assert (
-                "__weights" not in df.columns
-            ), "__weights is an internal lifelines column, please rename your column first."
+            assert "__weights" not in df.columns, "__weights is an internal lifelines column, please rename your column first."
             df["__weights"] = 1.0
         else:
             self.weights_col = weights_col
@@ -217,9 +215,7 @@ def fit(
         )
 
         self.params_ = pd.Series(params_, index=df.columns, name="coef") / self._norm_std
-        self.variance_matrix_ = pd.DataFrame(
-            -inv(self._hessian_) / np.outer(self._norm_std, self._norm_std), index=df.columns
-        )
+        self.variance_matrix_ = pd.DataFrame(-inv(self._hessian_) / np.outer(self._norm_std, self._norm_std), index=df.columns)
         self.standard_errors_ = self._compute_standard_errors(
             normalize(df, self._norm_mean, self._norm_std), events, start, stop, weights
         )
@@ -260,10 +256,9 @@ def _partition_by_strata(self, X, events, start, stop, weights):
             ), stratum
 
     def _partition_by_strata_and_apply(self, X, events, start, stop, weights, function, *args):
-        for (
-            (stratified_X, stratified_events, stratified_start, stratified_stop, stratified_W),
-            _,
-        ) in self._partition_by_strata(X, events, start, stop, weights):
+        for ((stratified_X, stratified_events, stratified_start, stratified_stop, stratified_W), _) in self._partition_by_strata(
+            X, events, start, stop, weights
+        ):
             yield function(stratified_X, stratified_events, stratified_start, stratified_stop, stratified_W, *args)
 
     def _compute_z_values(self):
@@ -375,9 +370,7 @@ def _newton_rhaphson(
             i += 1
 
             if self.strata is None:
-                h, g, ll = self._get_gradients(
-                    df.values, events.values, start.values, stop.values, weights.values, beta
-                )
+                h, g, ll = self._get_gradients(df.values, events.values, start.values, stop.values, weights.values, beta)
             else:
                 g = np.zeros_like(beta)
                 h = np.zeros((d, d))
@@ -714,11 +707,7 @@ def log_likelihood_ratio_test(self):
         degrees_freedom = self.params_.shape[0]
         p_value = _chisq_test_p_value(test_stat, degrees_freedom=degrees_freedom)
         return StatisticalResult(
-            p_value,
-            test_stat,
-            name="log-likelihood ratio test",
-            degrees_freedom=degrees_freedom,
-            null_distribution="chi squared",
+            p_value, test_stat, name="log-likelihood ratio test", degrees_freedom=degrees_freedom, null_distribution="chi squared"
         )
 
     def plot(self, columns=None, ax=None, **errorbar_kwargs):
@@ -777,18 +766,14 @@ def plot(self, columns=None, ax=None, **errorbar_kwargs):
 
         return ax
 
-    def _compute_cumulative_baseline_hazard(
-        self, tv_data, events, start, stop, weights
-    ):  # pylint: disable=too-many-locals
+    def _compute_cumulative_baseline_hazard(self, tv_data, events, start, stop, weights):  # pylint: disable=too-many-locals
 
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             hazards = self.predict_partial_hazard(tv_data).values
 
         unique_death_times = np.unique(stop[events.values])
-        baseline_hazard_ = pd.DataFrame(
-            np.zeros_like(unique_death_times), index=unique_death_times, columns=["baseline hazard"]
-        )
+        baseline_hazard_ = pd.DataFrame(np.zeros_like(unique_death_times), index=unique_death_times, columns=["baseline hazard"])
 
         for t in unique_death_times:
             ix = (start.values < t) & (t <= stop.values)
diff --git a/lifelines/fitters/coxph_fitter.py b/lifelines/fitters/coxph_fitter.py
index 90c0af63b..93fcbfdab 100644
--- a/lifelines/fitters/coxph_fitter.py
+++ b/lifelines/fitters/coxph_fitter.py
@@ -15,7 +15,7 @@
 from autograd import numpy as anp
 from autograd import elementwise_grad
 
-from lifelines.fitters import RegressionFitter, ParametricRegressionFitter
+from lifelines.fitters import RegressionFitter, SemiParametricRegressionFittter, ParametricRegressionFitter
 from lifelines.fitters.mixins import SplineFitterMixin, ProportionalHazardMixin
 from lifelines.plotting import set_kwargs_drawstyle
 from lifelines.statistics import _chisq_test_p_value, StatisticalResult
@@ -125,7 +125,7 @@ def decide(self, batch_mode: Optional[bool], n_unique: int, n_total: int, n_vars
         return self.SINGLE
 
 
-class CoxPHFitter(RegressionFitter, ProportionalHazardMixin):
+class CoxPHFitter(SemiParametricRegressionFittter, ProportionalHazardMixin):
     r"""
     This class implements fitting Cox's proportional hazard model using Efron's method for ties.
 
diff --git a/lifelines/utils/printer.py b/lifelines/utils/printer.py
index 5206165d1..7363acb7d 100644
--- a/lifelines/utils/printer.py
+++ b/lifelines/utils/printer.py
@@ -6,11 +6,20 @@
 
 
 class Printer:
-    def __init__(self, model, headers: List[Tuple[str, Any]], justify: Callable, decimals: int, header_kwargs: Dict):
+    def __init__(
+        self,
+        model,
+        headers: List[Tuple[str, Any]],
+        footers: List[Tuple[str, Any]],
+        justify: Callable,
+        decimals: int,
+        header_kwargs: Dict,
+    ):
         self.headers = headers
         self.model = model
         self.decimals = decimals
         self.justify = justify
+        self.footers = footers  # TODO: use this variable
 
         for tuple_ in header_kwargs.items():
             self.add_to_headers(tuple_)
@@ -85,18 +94,25 @@ def to_html(self):
             except AttributeError:
                 pass
 
+            try:
+                footers.append(("AIC", "{:.{prec}f}".format(self.model.AIC_, prec=decimals)))
+            except AttributeError:
+                pass
+
+            try:
+                footers.append(("Partial AIC", "{:.{prec}f}".format(self.model.AIC_partial_, prec=decimals)))
+            except AttributeError:
+                pass
+
             try:
                 sr = self.model.log_likelihood_ratio_test()
-                footers.extend(
-                    [
-                        ("Concordance", "{:.{prec}f}".format(self.model.concordance_index_, prec=decimals)),
-                        (
-                            "Log-likelihood ratio test",
-                            "{:.{prec}f} on {} df".format(sr.test_statistic, sr.degrees_freedom, prec=decimals),
-                        ),
-                        ("-log2(p) of ll-ratio test", "{:.{prec}f}".format(-np.log2(sr.p_value), prec=decimals)),
-                    ]
+                footers.append(
+                    (
+                        "Log-likelihood ratio test",
+                        "{:.{prec}f} on {} df".format(sr.test_statistic, sr.degrees_freedom, prec=decimals),
+                    )
                 )
+                footers.append(("-log2(p) of ll-ratio test", "{:.{prec}f}".format(-np.log2(sr.p_value), prec=decimals)))
             except AttributeError:
                 pass
 
@@ -177,12 +193,21 @@ def ascii_print(self):
 
         with np.errstate(invalid="ignore", divide="ignore"):
 
+            print("---")
             try:
-                print("---")
                 if utils.CensoringType.is_right_censoring(self.model) and self.model._KNOWN_MODEL:
                     print("Concordance = {:.{prec}f}".format(self.model.concordance_index_, prec=decimals))
             except AttributeError:
                 pass
+            try:
+                print("AIC = {:.{prec}f}".format(self.model.AIC_, prec=decimals))
+            except AttributeError:
+                pass
+
+            try:
+                print("Partial AIC = {:.{prec}f}".format(self.model.AIC_partial_, prec=decimals))
+            except AttributeError:
+                pass
 
             try:
                 sr = self.model.log_likelihood_ratio_test()

From 06871e92de9d144302aed9b13919be8497e8098e Mon Sep 17 00:00:00 2001
From: CamDavidsonPilon <cam.davidson.pilon@gmail.com>
Date: Thu, 14 May 2020 12:22:11 -0400
Subject: [PATCH 04/10] adding support for weighted penalty

---
 lifelines/fitters/__init__.py                 | 19 ++++---
 lifelines/fitters/cox_time_varying_fitter.py  | 14 +++--
 lifelines/fitters/coxph_fitter.py             | 20 ++++---
 .../generalized_gamma_regression_fitter.py    |  9 ++--
 lifelines/fitters/log_logistic_aft_fitter.py  |  5 +-
 lifelines/fitters/log_normal_aft_fitter.py    |  7 +--
 ...piecewise_exponential_regression_fitter.py |  2 +-
 lifelines/fitters/weibull_aft_fitter.py       |  4 +-
 lifelines/tests/test_estimation.py            | 53 +++++++++++++++++++
 9 files changed, 100 insertions(+), 33 deletions(-)

diff --git a/lifelines/fitters/__init__.py b/lifelines/fitters/__init__.py
index ec49248a3..b103266e8 100644
--- a/lifelines/fitters/__init__.py
+++ b/lifelines/fitters/__init__.py
@@ -1228,7 +1228,7 @@ class ParametricRegressionFitter(RegressionFitter):
     _scipy_fit_method = "BFGS"
     _scipy_fit_options: Dict[str, Any] = dict()
 
-    def __init__(self, alpha=0.05, penalizer=0.0, l1_ratio=0.0):
+    def __init__(self, alpha: float = 0.05, penalizer: Union[float, np.array] = 0.0, l1_ratio: float = 0.0):
         super(ParametricRegressionFitter, self).__init__(alpha=alpha)
         self.penalizer = penalizer
         self.l1_ratio = l1_ratio
@@ -1744,13 +1744,18 @@ def _add_penalty(self, params: Dict, neg_ll: float):
         params_array, _ = flatten(params)
         # remove intercepts from being penalized
         params_array = params_array[~self._constant_cols]
-        if self.penalizer > 0 and self.l1_ratio > 0:
-            penalty = self.l1_ratio * anp.abs(params_array).sum() + 0.5 * (1.0 - self.l1_ratio) * (params_array ** 2).sum()
-        elif self.penalizer > 0 and self.l1_ratio <= 0:
-            penalty = 0.5 * (params_array ** 2).sum()
+        if (isinstance(self.penalizer, np.ndarray) or self.penalizer > 0) and self.l1_ratio > 0:
+            penalty = (
+                self.l1_ratio * (self.penalizer * anp.abs(params_array)).sum()
+                + 0.5 * (1.0 - self.l1_ratio) * (self.penalizer * (params_array) ** 2).sum()
+            )
+
+        elif (isinstance(self.penalizer, np.ndarray) or self.penalizer > 0) and self.l1_ratio <= 0:
+            penalty = 0.5 * (self.penalizer * (params_array) ** 2).sum()
+
         else:
             penalty = 0
-        return neg_ll + self.penalizer * penalty
+        return neg_ll + penalty
 
     def _create_neg_likelihood_with_penalty_function(
         self, params_array, Ts, E, weights, entries, Xs, likelihood=None, penalty=None
@@ -2096,7 +2101,7 @@ def print_summary(self, decimals=2, style=None, **kwargs):
             headers.append(("weights col", "'%s'" % self.weights_col))
         if self.entry_col:
             headers.append(("entry col", "'%s'" % self.entry_col))
-        if self.penalizer > 0:
+        if isinstance(self.penalizer, np.ndarray) or self.penalizer > 0:
             headers.append(("penalizer", self.penalizer))
         if self.robust:
             headers.append(("robust variance", True))
diff --git a/lifelines/fitters/cox_time_varying_fitter.py b/lifelines/fitters/cox_time_varying_fitter.py
index a54151839..cf104b2af 100644
--- a/lifelines/fitters/cox_time_varying_fitter.py
+++ b/lifelines/fitters/cox_time_varying_fitter.py
@@ -86,9 +86,6 @@ class CoxTimeVaryingFitter(SemiParametricRegressionFittter, ProportionalHazardMi
 
     def __init__(self, alpha=0.05, penalizer=0.0, l1_ratio: float = 0.0, strata=None):
         super(CoxTimeVaryingFitter, self).__init__(alpha=alpha)
-        if penalizer < 0:
-            raise ValueError("penalizer parameter must be >= 0.")
-
         self.alpha = alpha
         self.penalizer = penalizer
         self.strata = strata
@@ -338,14 +335,15 @@ def _newton_rhaphson(
         """
         assert precision <= 1.0, "precision must be less than or equal to 1."
 
-        # soft penalizer functions, from https://www.cs.ubc.ca/cgi-bin/tr/2009/TR-2009-19.pdf
         # soft penalizer functions, from https://www.cs.ubc.ca/cgi-bin/tr/2009/TR-2009-19.pdf
         soft_abs = lambda x, a: 1 / a * (anp.logaddexp(0, -a * x) + anp.logaddexp(0, a * x))
         penalizer = (
             lambda beta, a: n
             * 0.5
-            * self.penalizer
-            * (self.l1_ratio * soft_abs(beta, a).sum() + (1 - self.l1_ratio) * ((beta) ** 2).sum())
+            * (
+                self.l1_ratio * (self.penalizer * soft_abs(beta, a)).sum()
+                + (1 - self.l1_ratio) * (self.penalizer * beta ** 2).sum()
+            )
         )
         d_penalizer = elementwise_grad(penalizer)
         dd_penalizer = elementwise_grad(d_penalizer)
@@ -388,7 +386,7 @@ def _newton_rhaphson(
                 # if the user supplied a non-trivial initial point, we need to delay this.
                 self._log_likelihood_null = ll
 
-            if self.penalizer > 0:
+            if isinstance(self.penalizer, np.ndarray) or self.penalizer > 0:
                 ll -= penalizer(beta, 1.5 ** i)
                 g -= d_penalizer(beta, 1.5 ** i)
                 h[np.diag_indices(d)] -= dd_penalizer(beta, 1.5 ** i)
@@ -653,7 +651,7 @@ def print_summary(self, decimals=2, style=None, **kwargs):
             headers.append(("event col", "'%s'" % self.event_col))
         if self.weights_col:
             headers.append(("weights col", "'%s'" % self.weights_col))
-        if self.penalizer > 0:
+        if isinstance(self.penalizer, np.ndarray) or self.penalizer > 0:
             headers.append(("penalizer", self.penalizer))
         if self.strata:
             headers.append(("strata", self.strata))
diff --git a/lifelines/fitters/coxph_fitter.py b/lifelines/fitters/coxph_fitter.py
index 93fcbfdab..05113e99b 100644
--- a/lifelines/fitters/coxph_fitter.py
+++ b/lifelines/fitters/coxph_fitter.py
@@ -143,12 +143,15 @@ class CoxPHFitter(SemiParametricRegressionFittter, ProportionalHazardMixin):
       baseline_estimation_method: string, optional
         specify how the fitter should estimate the baseline. ``"breslow"`` or ``"spline"``
 
-      penalizer: float, optional (default=0.0)
+      penalizer: float or array, optional (default=0.0)
         Attach a penalty to the size of the coefficients during regression. This improves
         stability of the estimates and controls for high correlation between covariates.
         For example, this shrinks the magnitude value of :math:`\beta_i`. See ``l1_ratio`` below.
         The penalty term is :math:`\frac{1}{2} \text{penalizer} \left( (1-\text{l1_ratio}) ||\beta||_2^2 + \text{l1_ratio}||\beta||_1\right)`.
 
+        Alternatively, penalizer is an array equal in size to the number of parameters, with penalty coefficients for specific variables. For
+        example, `penalizer=0.01 * np.ones(p)` is the same as `penalizer=0.01`
+
       l1_ratio: float, optional (default=0.0)
         Specify what ratio to assign to a L1 vs L2 penalty. Same as scikit-learn. See ``penalizer`` above.
 
@@ -202,7 +205,7 @@ class CoxPHFitter(SemiParametricRegressionFittter, ProportionalHazardMixin):
     def __init__(
         self,
         baseline_estimation_method: str = "breslow",
-        penalizer: float = 0.0,
+        penalizer: Union[float, np.ndarray] = 0.0,
         strata: Optional[Union[List[str], str]] = None,
         l1_ratio: float = 0.0,
         n_baseline_knots: int = 1,
@@ -210,8 +213,7 @@ def __init__(
     ) -> None:
 
         super(CoxPHFitter, self).__init__(**kwargs)
-        if penalizer < 0:
-            raise ValueError("penalizer parameter must be >= 0.")
+
         if l1_ratio < 0 or l1_ratio > 1:
             raise ValueError("l1_ratio parameter must in [0, 1].")
 
@@ -555,8 +557,10 @@ def _newton_rhapson_for_efron_model(
         penalizer = (
             lambda beta, a: n
             * 0.5
-            * self.penalizer
-            * (self.l1_ratio * soft_abs(beta, a).sum() + (1 - self.l1_ratio) * ((beta) ** 2).sum())
+            * (
+                self.l1_ratio * (self.penalizer * soft_abs(beta, a)).sum()
+                + (1 - self.l1_ratio) * (self.penalizer * beta ** 2).sum()
+            )
         )
         d_penalizer = elementwise_grad(penalizer)
         dd_penalizer = elementwise_grad(d_penalizer)
@@ -605,7 +609,7 @@ def _newton_rhapson_for_efron_model(
                 # if the user supplied a non-trivial initial point, we need to delay this.
                 self._ll_null_ = ll_
 
-            if self.penalizer > 0:
+            if isinstance(self.penalizer, np.ndarray) or self.penalizer > 0:
                 ll_ -= penalizer(beta, 1.3 ** i)
                 g -= d_penalizer(beta, 1.3 ** i)
                 h[np.diag_indices(d)] -= dd_penalizer(beta, 1.3 ** i)
@@ -1323,7 +1327,7 @@ def print_summary(self, decimals: int = 2, style: Optional[str] = None, **kwargs
             headers.append(("weights col", "'%s'" % self.weights_col))
         if self.cluster_col:
             headers.append(("cluster col", "'%s'" % self.cluster_col))
-        if self.penalizer > 0:
+        if isinstance(self.penalizer, np.ndarray) or self.penalizer > 0:
             headers.append(("penalizer", self.penalizer))
             headers.append(("l1 ratio", self.l1_ratio))
         if self.robust or self.cluster_col:
diff --git a/lifelines/fitters/generalized_gamma_regression_fitter.py b/lifelines/fitters/generalized_gamma_regression_fitter.py
index a11f84e4d..f1af8e425 100644
--- a/lifelines/fitters/generalized_gamma_regression_fitter.py
+++ b/lifelines/fitters/generalized_gamma_regression_fitter.py
@@ -58,7 +58,10 @@ class GeneralizedGammaRegressionFitter(ParametricRegressionFitter):
     -----------
     alpha: float, optional (default=0.05)
         the level in the confidence intervals.
-
+    penalizer: float or array, optional (default=0.0)
+        the penalizer coefficient to the size of the coefficients. See `l1_ratio`. Must be equal to or greater than 0.
+        Alternatively, penalizer is an array equal in size to the number of parameters, with penalty coefficients for specific variables. For
+        example, `penalizer=0.01 * np.ones(p)` is the same as `penalizer=0.01`
 
     Examples
     --------
@@ -126,9 +129,7 @@ def _create_initial_point(self, Ts, E, entries, weights, Xs):
             # we may use this later in print_summary
             self._ll_null_ = uni_model.log_likelihood_
 
-            default_point = super(GeneralizedGammaRegressionFitter, self)._create_initial_point(
-                Ts, E, entries, weights, Xs
-            )
+            default_point = super(GeneralizedGammaRegressionFitter, self)._create_initial_point(Ts, E, entries, weights, Xs)
             nested_point = {}
 
             nested_point["mu_"] = np.array([0.0] * (len(Xs.mappings["mu_"])))
diff --git a/lifelines/fitters/log_logistic_aft_fitter.py b/lifelines/fitters/log_logistic_aft_fitter.py
index 1502d8dc1..0d99d6bee 100644
--- a/lifelines/fitters/log_logistic_aft_fitter.py
+++ b/lifelines/fitters/log_logistic_aft_fitter.py
@@ -35,8 +35,11 @@ class LogLogisticAFTFitter(ParametericAFTRegressionFitter):
     fit_intercept: boolean, optional (default=True)
         Allow lifelines to add an intercept column of 1s to df, and ancillary_df if applicable.
 
-    penalizer: float, optional (default=0.0)
+    penalizer: float or array, optional (default=0.0)
         the penalizer coefficient to the size of the coefficients. See `l1_ratio`. Must be equal to or greater than 0.
+        Alternatively, penalizer is an array equal in size to the number of parameters, with penalty coefficients for specific variables. For
+        example, `penalizer=0.01 * np.ones(p)` is the same as `penalizer=0.01`
+
 
     l1_ratio: float, optional (default=0.0)
         how much of the penalizer should be attributed to an l1 penalty (otherwise an l2 penalty). The penalty function looks like
diff --git a/lifelines/fitters/log_normal_aft_fitter.py b/lifelines/fitters/log_normal_aft_fitter.py
index 75d7d06e2..d8d4e3dca 100644
--- a/lifelines/fitters/log_normal_aft_fitter.py
+++ b/lifelines/fitters/log_normal_aft_fitter.py
@@ -37,8 +37,10 @@ class LogNormalAFTFitter(ParametericAFTRegressionFitter):
     fit_intercept: bool, optional (default=True)
         Allow lifelines to add an intercept column of 1s to df, and ancillary_df if applicable.
 
-    penalizer: float, optional (default=0.0)
+    penalizer: float or array, optional (default=0.0)
         the penalizer coefficient to the size of the coefficients. See `l1_ratio`. Must be equal to or greater than 0.
+        Alternatively, penalizer is an array equal in size to the number of parameters, with penalty coefficients for specific variables. For
+        example, `penalizer=0.01 * np.ones(p)` is the same as `penalizer=0.01`
 
     l1_ratio: float, optional (default=0.0)
         how much of the penalizer should be attributed to an l1 penalty (otherwise an l2 penalty). The penalty function looks like
@@ -158,8 +160,7 @@ def predict_percentile(
             S = norm.sf(Z)
 
             return pd.Series(
-                exp_mu_ * np.exp(np.sqrt(2) * sigma_ * erfinv(2 * (1 - p * S) - 1)) - conditional_after,
-                index=_get_index(df),
+                exp_mu_ * np.exp(np.sqrt(2) * sigma_ * erfinv(2 * (1 - p * S) - 1)) - conditional_after, index=_get_index(df)
             )
 
     def predict_expectation(self, df: pd.DataFrame, ancillary_df: Optional[pd.DataFrame] = None) -> pd.Series:
diff --git a/lifelines/fitters/piecewise_exponential_regression_fitter.py b/lifelines/fitters/piecewise_exponential_regression_fitter.py
index e96c3743c..5efca2634 100644
--- a/lifelines/fitters/piecewise_exponential_regression_fitter.py
+++ b/lifelines/fitters/piecewise_exponential_regression_fitter.py
@@ -53,13 +53,13 @@ def __init__(self, breakpoints, alpha=0.05, penalizer=0.0):
         self.breakpoints = breakpoints
         self.n_breakpoints = len(self.breakpoints)
 
+        assert isinstance(self.penalizer, float), "penalizer must be a float"
         self.penalizer = penalizer
         self._fitted_parameter_names = ["lambda_%d_" % i for i in range(self.n_breakpoints + 1)]
 
     def _add_penalty(self, params, neg_ll):
         params_stacked = np.stack(params.values())
         coef_penalty = 0
-
         if self.penalizer > 0:
             for i in range(params_stacked.shape[1]):
                 if not self._constant_cols[i]:
diff --git a/lifelines/fitters/weibull_aft_fitter.py b/lifelines/fitters/weibull_aft_fitter.py
index d9370fe46..ad627ba2d 100644
--- a/lifelines/fitters/weibull_aft_fitter.py
+++ b/lifelines/fitters/weibull_aft_fitter.py
@@ -43,8 +43,10 @@ class WeibullAFTFitter(ParametericAFTRegressionFitter, ProportionalHazardMixin):
     fit_intercept: boolean, optional (default=True)
         Allow lifelines to add an intercept column of 1s to df, and ancillary_df if applicable.
 
-    penalizer: float, optional (default=0.0)
+    penalizer: float or array, optional (default=0.0)
         the penalizer coefficient to the size of the coefficients. See `l1_ratio`. Must be equal to or greater than 0.
+        Alternatively, penalizer is an array equal in size to the number of parameters, with penalty coefficients for specific variables. For
+        example, `penalizer=0.01 * np.ones(p)` is the same as `penalizer=0.01`
 
     l1_ratio: float, optional (default=0.0)
         how much of the penalizer should be attributed to an l1 penalty (otherwise an l2 penalty). The penalty function looks like
diff --git a/lifelines/tests/test_estimation.py b/lifelines/tests/test_estimation.py
index 602dd85f8..b09239bfd 100644
--- a/lifelines/tests/test_estimation.py
+++ b/lifelines/tests/test_estimation.py
@@ -1481,6 +1481,21 @@ def rossi(self):
         rossi["_int"] = 1.0
         return rossi
 
+    def test_penalizer_can_be_an_array(self, rossi):
+
+        wf_array = WeibullAFTFitter(penalizer=0.01 * np.ones(7)).fit(rossi, "week", "arrest")
+        wf_float = WeibullAFTFitter(penalizer=0.01).fit(rossi, "week", "arrest")
+
+        assert_frame_equal(wf_array.summary, wf_float.summary)
+
+    def test_penalizer_can_be_an_array_and_check_it_behaves_as_expected(self, rossi):
+
+        penalty = np.array([0, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01])
+        wf_array = WeibullAFTFitter(penalizer=penalty).fit(rossi, "week", "arrest")
+        wf_float = WeibullAFTFitter(penalizer=0.01).fit(rossi, "week", "arrest")
+
+        assert abs(wf_array.summary.loc[("lambda_", "fin"), "coef"]) > abs(wf_float.summary.loc[("lambda_", "fin"), "coef"])
+
     def test_custom_weibull_model_gives_the_same_data_as_implemented_weibull_model(self, rossi):
         class CustomWeibull(ParametricRegressionFitter):
             _scipy_fit_method = "SLSQP"
@@ -2534,6 +2549,21 @@ def cph(self):
     def cph_spline(self):
         return CoxPHFitter(baseline_estimation_method="spline")
 
+    def test_penalizer_can_be_an_array(self, rossi):
+
+        cph_array = CoxPHFitter(penalizer=0.01 * np.ones(7)).fit(rossi, "week", "arrest")
+        cph_float = CoxPHFitter(penalizer=0.01).fit(rossi, "week", "arrest")
+
+        assert_frame_equal(cph_array.summary, cph_float.summary)
+
+    def test_penalizer_can_be_an_array_and_check_it_behaves_as_expected(self, rossi):
+
+        penalty = np.array([0, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01])
+        cph_array = CoxPHFitter(penalizer=penalty).fit(rossi, "week", "arrest")
+        cph_float = CoxPHFitter(penalizer=0.01).fit(rossi, "week", "arrest")
+
+        assert abs(cph_array.summary.loc["fin", "coef"]) > abs(cph_float.summary.loc["fin", "coef"])
+
     def test_compute_followup_hazard_ratios(self, cph, cph_spline, rossi):
         cph.fit(rossi, "week", "arrest")
         cph.compute_followup_hazard_ratios(rossi, [15, 25, 35, 45])
@@ -4353,6 +4383,29 @@ def dfcv(self):
     def heart(self):
         return load_stanford_heart_transplants()
 
+    def test_penalizer_can_be_an_array(self, dfcv):
+
+        cph_array = CoxTimeVaryingFitter(penalizer=0.01 * np.ones(2)).fit(
+            dfcv, id_col="id", start_col="start", stop_col="stop", event_col="event"
+        )
+        cph_float = CoxTimeVaryingFitter(penalizer=0.01).fit(
+            dfcv, id_col="id", start_col="start", stop_col="stop", event_col="event"
+        )
+
+        assert_frame_equal(cph_array.summary, cph_float.summary)
+
+    def test_penalizer_can_be_an_array_and_check_it_behaves_as_expected(self, dfcv):
+
+        penalty = np.array([0, 0.01])
+        cph_array = CoxTimeVaryingFitter(penalizer=penalty).fit(
+            dfcv, id_col="id", start_col="start", stop_col="stop", event_col="event"
+        )
+        cph_float = CoxTimeVaryingFitter(penalizer=0.01).fit(
+            dfcv, id_col="id", start_col="start", stop_col="stop", event_col="event"
+        )
+
+        assert abs(cph_array.summary.loc["z", "coef"]) > abs(cph_float.summary.loc["z", "coef"])
+
     def test_model_can_accept_null_covariates(self, ctv, dfcv):
         ctv.fit(dfcv[["id", "start", "stop", "event"]], id_col="id", start_col="start", stop_col="stop", event_col="event")
 

From 348281a21bb77b7d0bdd5de73d7731c0fbd1243a Mon Sep 17 00:00:00 2001
From: CamDavidsonPilon <cam.davidson.pilon@gmail.com>
Date: Thu, 14 May 2020 13:42:07 -0400
Subject: [PATCH 05/10] refactor footer stuff in Printer

---
 docs/Survival Regression.rst                 | 57 ++++++++++++------
 lifelines/fitters/__init__.py                | 16 ++++-
 lifelines/fitters/cox_time_varying_fitter.py | 15 ++++-
 lifelines/fitters/coxph_fitter.py            | 15 ++++-
 lifelines/utils/printer.py                   | 63 ++------------------
 5 files changed, 88 insertions(+), 78 deletions(-)

diff --git a/docs/Survival Regression.rst b/docs/Survival Regression.rst
index 3471d0fe6..480dd78ac 100644
--- a/docs/Survival Regression.rst	
+++ b/docs/Survival Regression.rst	
@@ -103,7 +103,9 @@ The implementation of the Cox model in *lifelines* is under :class:`~lifelines.f
     prio  3.19 <0.005       9.48
     ---
     Concordance = 0.64
-    Log-likelihood ratio test = 33.27 on 7 df, -log2(p)=15.37
+    Partial AIC = 1331.50
+    log-likelihood ratio test = 33.27 on 7 df
+    -log2(p) of ll-ratio test = 15.37
     """
 
 To access the coefficients and the baseline hazard directly, you can use :attr:`~lifelines.fitters.coxph_fitter.CoxPHFitter.params_` and :attr:`~lifelines.fitters.coxph_fitter.CoxPHFitter.baseline_hazard_` respectively. Taking a look at these coefficients for a moment, ``prio`` (the number of prior arrests) has a coefficient of about 0.09. Thus, a one unit increase in ``prio`` means the the baseline hazard will increase by a factor of :math:`\exp{(0.09)} = 1.10` - about a 10% increase. Recall, in the Cox proportional hazard model, a higher hazard means more at risk of the event occurring. The value :math:`\exp{(0.09)}` is called the *hazard ratio*, a name that will be clear with another example.
@@ -181,7 +183,7 @@ Back to our original problem of predicting the event time of censored individual
 Penalties and sparse regression
 -----------------------------------------------
 
-It's possible to add a penalizer term to the Cox regression as well. One can use these to i) stabilize the coefficients, ii) shrink the estimates to 0, iii) encourages a Bayesian interpretation, and iv) create sparse coefficients. Regression models, including the Cox model, include both an L1 and L2 penalty:
+It's possible to add a penalizer term to the Cox regression as well. One can use these to i) stabilize the coefficients, ii) shrink the estimates to 0, iii) encourages a Bayesian viewpoint, and iv) create sparse coefficients. Regression models, including the Cox model, include both an L1 and L2 penalty:
 
 .. math:: \frac{1}{2} \text{penalizer} \left((1-\text{l1\_ratio}) \cdot ||\beta||_2^2 + \text{l1\_ratio} \cdot ||\beta||_1\right)
 
@@ -204,6 +206,29 @@ To use this in *lifelines*, both the ``penalizer`` and ``l1_ratio`` can be speci
     cph.print_summary()
 
 
+Instead of a float, an *array* can be provided that is the same size as the number of estimated parameters. The values in the array
+are specific penalty coefficients for each covariate. This is useful for more complicated covariate structure. Some examples:
+
+i) you have lots of confounders you wish to penalizer, but not the main treatment(s).
+
+.. code:: python
+
+    from lifelines import CoxPHFitter
+    from lifelines.datasets import load_rossi
+
+    rossi = load_rossi()
+
+    # variable `fin` is the treatment of interest so don't penalize it at all
+    penalty = np.array([0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5])
+
+    cph = CoxPHFitter(penalizer=penalty)
+    cph.fit(rossi, 'week', 'arrest')
+    cph.print_summary()
+
+ii) you have to `fuse categories together <https://stats.stackexchange.com/questions/146907/principled-way-of-collapsing-categorical-variables-with-many-levels>`_.
+
+
+
 Plotting the coefficients
 ------------------------------
 
@@ -300,15 +325,7 @@ To specify variables to be used in stratification, we define them in the call to
 
 .. code:: python
 
-    from lifelines.datasets import load_rossi
-    from lifelines import CoxPHFitter
-
-    rossi_dataset = load_rossi()
-    cph = CoxPHFitter()
-    cph.fit(rossi_dataset, 'week', event_col='arrest', strata=['race'])
-
-    cph.print_summary()  # access the results using cph.summary
-
+p
     """
     <lifelines.CoxPHFitter: fitted with 432 observations, 318 censored>
           duration col = 'week'
@@ -328,8 +345,10 @@ To specify variables to be used in stratification, we define them in the call to
     paro -0.09       0.92      0.20 -0.44   0.66      0.60       -0.47        0.30
     prio  0.09       1.10      0.03  3.21 <0.005      9.56        0.04        0.15
     ---
-    Concordance = 0.64
-    Likelihood ratio test = 109.63 on 6 df, -log2(p)=68.48
+    Concordance = 0.63
+    Partial AIC = 1253.13
+    log-likelihood ratio test = 32.73 on 6 df
+    -log2(p) of ll-ratio test = 16.37
     """
 
     cph.baseline_cumulative_hazard_.shape
@@ -492,7 +511,9 @@ The Weibull AFT model is implemented under :class:`~lifelines.fitters.weibull_af
     rho_    _intercept  0.339      1.404     0.089  3.809 <0.0005    12.808       0.165       0.514
     ---
     Concordance = 0.640
-    Log-likelihood ratio test = 33.416 on 7 df, -log2(p)=15.462
+    AIC = 1377.833
+    log-likelihood ratio test = 33.416 on 7 df
+    -log2(p) of ll-ratio test = 15.462
     """
 
 From above, we can see that ``prio``, which is the number of previous incarcerations, has a large negative coefficient. This means that each addition incarcerations changes a subject's mean/median survival time by :math:`\exp(-0.066) = 0.936`, approximately a 7% decrease in mean/median survival time. What is the mean/median survival time?
@@ -646,7 +667,7 @@ When predicting time remaining for uncensored individuals, you can use the `cond
     aft.predict_percentile(X, p=0.9, ancillary_df=censored_X, conditional_after=censored_subjects_last_obs)
 
 
-There are two hyper-parameters that can be used to to achieve a better test score. These are ``penalizer`` and ``l1_ratio`` in the call to :class:`~lifelines.fitters.weibull_aft_fitter.WeibullAFTFitter`. The penalizer is similar to scikit-learn's ``ElasticNet`` model, see their `docs <https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html>`_.
+There are two hyper-parameters that can be used to to achieve a better test score. These are ``penalizer`` and ``l1_ratio`` in the call to :class:`~lifelines.fitters.weibull_aft_fitter.WeibullAFTFitter`. The penalizer is similar to scikit-learn's ``ElasticNet`` model, see their `docs <https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html>`_. (However, *lifelines* will also accept an array for custom penalizer per variable, see `Cox docs above <https://lifelines.readthedocs.io/en/latest/Survival%20Regression.html#penalties-and-sparse-regression>`_)
 
 .. code:: python
 
@@ -680,8 +701,10 @@ There are two hyper-parameters that can be used to to achieve a better test scor
             _intercept  0.00       1.00      0.19  0.00   1.00      0.00       -0.38        0.38
     rho_    _intercept -0.00       1.00       nan   nan    nan       nan         nan         nan
     ---
-    Concordance = 0.60
-    Log-likelihood ratio test = -4028.65 on 7 df, -log2(p)=-0.00
+    Concordance = 0.64
+    AIC = 1377.91
+    log-likelihood ratio test = 33.34 on 7 df
+    -log2(p) of ll-ratio test = 15.42
     """
 
 
diff --git a/lifelines/fitters/__init__.py b/lifelines/fitters/__init__.py
index b103266e8..4da7d9a8b 100644
--- a/lifelines/fitters/__init__.py
+++ b/lifelines/fitters/__init__.py
@@ -2115,7 +2115,21 @@ def print_summary(self, decimals=2, style=None, **kwargs):
             ]
         )
 
-        p = Printer(self, headers, justify, decimals, kwargs)
+        sr = self.log_likelihood_ratio_test()
+        footers = []
+        footers.extend(
+            [
+                ("Concordance", "{:.{prec}f}".format(self.concordance_index_, prec=decimals)),
+                ("AIC", "{:.{prec}f}".format(self.AIC_, prec=decimals)),
+                (
+                    "log-likelihood ratio test",
+                    "{:.{prec}f} on {} df".format(sr.test_statistic, sr.degrees_freedom, prec=decimals),
+                ),
+                ("-log2(p) of ll-ratio test", "{:.{prec}f}".format(-np.log2(sr.p_value), prec=decimals)),
+            ]
+        )
+
+        p = Printer(self, headers, footers, justify, decimals, kwargs)
 
         p.print(style=style)
 
diff --git a/lifelines/fitters/cox_time_varying_fitter.py b/lifelines/fitters/cox_time_varying_fitter.py
index cf104b2af..71f568795 100644
--- a/lifelines/fitters/cox_time_varying_fitter.py
+++ b/lifelines/fitters/cox_time_varying_fitter.py
@@ -666,8 +666,21 @@ def print_summary(self, decimals=2, style=None, **kwargs):
             ]
         )
 
-        p = Printer(self, headers, justify, decimals, kwargs)
+        sr = self.log_likelihood_ratio_test()
+        footers = []
+        footers.extend(
+            [
+                ("Concordance", "{:.{prec}f}".format(self.concordance_index_, prec=decimals)),
+                ("Partial AIC", "{:.{prec}f}".format(self.AIC_partial_, prec=decimals)),
+                (
+                    "log-likelihood ratio test",
+                    "{:.{prec}f} on {} df".format(sr.test_statistic, sr.degrees_freedom, prec=decimals),
+                ),
+                ("-log2(p) of ll-ratio test", "{:.{prec}f}".format(-np.log2(sr.p_value), prec=decimals)),
+            ]
+        )
 
+        p = Printer(self, headers, footers, justify, decimals, kwargs)
         p.print(style=style)
 
     def log_likelihood_ratio_test(self):
diff --git a/lifelines/fitters/coxph_fitter.py b/lifelines/fitters/coxph_fitter.py
index 05113e99b..fdcbcd75d 100644
--- a/lifelines/fitters/coxph_fitter.py
+++ b/lifelines/fitters/coxph_fitter.py
@@ -1350,8 +1350,21 @@ def print_summary(self, decimals: int = 2, style: Optional[str] = None, **kwargs
             ]
         )
 
-        p = Printer(self, headers, justify, decimals, kwargs)
+        sr = self.log_likelihood_ratio_test()
+        footers = []
+        footers.extend(
+            [
+                ("Concordance", "{:.{prec}f}".format(self.concordance_index_, prec=decimals)),
+                ("Partial AIC", "{:.{prec}f}".format(self.AIC_partial_, prec=decimals)),
+                (
+                    "log-likelihood ratio test",
+                    "{:.{prec}f} on {} df".format(sr.test_statistic, sr.degrees_freedom, prec=decimals),
+                ),
+                ("-log2(p) of ll-ratio test", "{:.{prec}f}".format(-np.log2(sr.p_value), prec=decimals)),
+            ]
+        )
 
+        p = Printer(self, headers, footers, justify, decimals, kwargs)
         p.print(style=style)
 
     def _trivial_log_likelihood(self):
diff --git a/lifelines/utils/printer.py b/lifelines/utils/printer.py
index 7363acb7d..5c80f7efc 100644
--- a/lifelines/utils/printer.py
+++ b/lifelines/utils/printer.py
@@ -19,7 +19,7 @@ def __init__(
         self.model = model
         self.decimals = decimals
         self.justify = justify
-        self.footers = footers  # TODO: use this variable
+        self.footers = footers
 
         for tuple_ in header_kwargs.items():
             self.add_to_headers(tuple_)
@@ -85,39 +85,8 @@ def to_html(self):
             },
         )
 
-        footers = []
-        with np.errstate(invalid="ignore", divide="ignore"):
-
-            try:
-                if utils.CensoringType.is_right_censoring(self.model) and self.model._KNOWN_MODEL:
-                    footers.append(("Concordance", "{:.{prec}f}".format(self.model.score_, prec=decimals)))
-            except AttributeError:
-                pass
-
-            try:
-                footers.append(("AIC", "{:.{prec}f}".format(self.model.AIC_, prec=decimals)))
-            except AttributeError:
-                pass
-
-            try:
-                footers.append(("Partial AIC", "{:.{prec}f}".format(self.model.AIC_partial_, prec=decimals)))
-            except AttributeError:
-                pass
-
-            try:
-                sr = self.model.log_likelihood_ratio_test()
-                footers.append(
-                    (
-                        "Log-likelihood ratio test",
-                        "{:.{prec}f} on {} df".format(sr.test_statistic, sr.degrees_freedom, prec=decimals),
-                    )
-                )
-                footers.append(("-log2(p) of ll-ratio test", "{:.{prec}f}".format(-np.log2(sr.p_value), prec=decimals)))
-            except AttributeError:
-                pass
-
-        if footers:
-            footer_df = pd.DataFrame.from_records(footers).set_index(0)
+        if self.footers:
+            footer_df = pd.DataFrame.from_records(self.footers).set_index(0)
             footer_html = footer_df.to_html(header=False, notebook=True, index_names=False)
         else:
             footer_html = ""
@@ -194,28 +163,6 @@ def ascii_print(self):
         with np.errstate(invalid="ignore", divide="ignore"):
 
             print("---")
-            try:
-                if utils.CensoringType.is_right_censoring(self.model) and self.model._KNOWN_MODEL:
-                    print("Concordance = {:.{prec}f}".format(self.model.concordance_index_, prec=decimals))
-            except AttributeError:
-                pass
-            try:
-                print("AIC = {:.{prec}f}".format(self.model.AIC_, prec=decimals))
-            except AttributeError:
-                pass
-
-            try:
-                print("Partial AIC = {:.{prec}f}".format(self.model.AIC_partial_, prec=decimals))
-            except AttributeError:
-                pass
-
-            try:
-                sr = self.model.log_likelihood_ratio_test()
-                print(
-                    "Log-likelihood ratio test = {:.{prec}f} on {} df, -log2(p)={:.{prec}f}".format(
-                        sr.test_statistic, sr.degrees_freedom, -np.log2(sr.p_value), prec=decimals
-                    )
-                )
-            except AttributeError:
-                pass
+            for string, value in self.footers:
+                print("{} = {}".format(string, value))
         print()

From 9660f6c42187d75b05bb9f228b3c319fdc44bc84 Mon Sep 17 00:00:00 2001
From: CamDavidsonPilon <cam.davidson.pilon@gmail.com>
Date: Fri, 15 May 2020 09:40:20 -0400
Subject: [PATCH 06/10] adding new functionality to find_best..

---
 lifelines/tests/utils/test_utils.py | 21 ++++++++
 lifelines/utils/__init__.py         | 83 ++++++++++++++++++++++++-----
 2 files changed, 92 insertions(+), 12 deletions(-)

diff --git a/lifelines/tests/utils/test_utils.py b/lifelines/tests/utils/test_utils.py
index 93c4c5724..e616ab3da 100644
--- a/lifelines/tests/utils/test_utils.py
+++ b/lifelines/tests/utils/test_utils.py
@@ -1065,6 +1065,27 @@ def test_find_best_parametric_model_with_BIC():
     assert True
 
 
+def test_find_best_parametric_model_works_for_left_censoring():
+    T = np.random.exponential(2, 100)
+    model, score = utils.find_best_parametric_model(T, censoring_type="left", show_progress=True)
+    assert True
+
+
+def test_find_best_parametric_model_works_for_interval_censoring():
+    T_1 = np.random.exponential(2, 100)
+    T_2 = T_1 + 1
+    model, score = utils.find_best_parametric_model((T_1, T_2), censoring_type="interval", show_progress=True)
+    assert True
+
+
+def test_find_best_parametric_model_works_with_weights_and_entry():
+    T = np.random.exponential(2, 100)
+    W = np.random.randint(1, 5, size=100)
+    entry = np.random.exponential(0.01, 100)
+    model, score = utils.find_best_parametric_model(T, weights=W, entry=entry, show_progress=True)
+    assert True
+
+
 def test_safe_exp():
     from lifelines.utils.safe_exp import MAX
 
diff --git a/lifelines/utils/__init__.py b/lifelines/utils/__init__.py
index 724d08aa7..f5cf5e47b 100644
--- a/lifelines/utils/__init__.py
+++ b/lifelines/utils/__init__.py
@@ -42,6 +42,10 @@ class CensoringType:
     INTERVAL = 2
     RIGHT = 3
 
+    MAP = {"right": RIGHT, "left": LEFT, "interval": INTERVAL}
+
+    HUMAN_MAP = {LEFT: "left", RIGHT: "right", INTERVAL: "interval"}
+
     @classmethod
     def right_censoring(cls, function: Callable) -> Callable:
         @wraps(function)
@@ -84,11 +88,13 @@ def is_interval_censoring(cls, model) -> bool:
     @classmethod
     def get_human_readable_censoring_type(cls, model) -> str:
         if cls.is_interval_censoring(model):
-            return "interval"
+            return cls.HUMAN_MAP[cls.INTERVAL]
         elif cls.is_right_censoring(model):
-            return "right"
+            return cls.HUMAN_MAP[cls.RIGHT]
+        elif cls.is_left_censoring(model):
+            return cls.HUMAN_MAP[cls.LEFT]
         else:
-            return "left"
+            return
 
 
 class StatError(Exception):
@@ -1675,7 +1681,19 @@ def iterdicts(self):
             yield DataframeSliceDict(x.to_frame().T, self.mappings)
 
 
-def find_best_parametric_model(event_times, event_observed=None, scoring_method: str = "AIC", additional_models=None):
+def find_best_parametric_model(
+    event_times,
+    event_observed=None,
+    scoring_method: str = "AIC",
+    additional_models=None,
+    censoring_type="right",
+    timeline=None,
+    alpha=None,
+    ci_labels=None,
+    entry=None,
+    weights=None,
+    show_progress=False,
+):
     """
     To quickly determine the best¹ univariate model, this function will iterate through each
     parametric model available in lifelines and select the one that minimizes a particular measure of fit.
@@ -1685,13 +1703,27 @@ def find_best_parametric_model(event_times, event_observed=None, scoring_method:
     Parameters
     -------------
     event_times: list, np.array, pd.Series
-        a (n,) array of observed survival times.
+        a (n,) array of observed survival times. If interval censoring, a tuple of (lower_bound, upper_bound).
     event_observed: list, np.array, pd.Series
         a (n,) array of censored flags, 1 if observed,  0 if not. Default None assumes all observed.
     scoring_method: string
         one of {"AIC", "BIC"}
     additional_models: list
         list of other parametric models that implement the lifelines API.
+    censoring_type: str
+        {"right", "left", "interval"}
+    timeline: list, optional
+        return the model at the values in timeline (positively increasing)
+    alpha: float, optional
+        the alpha value in the confidence intervals. Overrides the initializing
+       alpha for this call to fit only.
+    ci_labels: list, optional
+        add custom column names to the generated confidence intervals as a length-2 list: [<lower-bound name>, <upper-bound name>]. Default: <label>_lower_<alpha>
+    entry: an array, or pd.Series, of length n
+        relative time when a subject entered the study. This is useful for left-truncated (not left-censored) observations. If None, all members of the population
+        entered study when they were "born": time zero.
+    weights: an array, or pd.Series, of length n
+        integer weights per observation
 
     Returns
     ----------
@@ -1711,20 +1743,37 @@ def find_best_parametric_model(event_times, event_observed=None, scoring_method:
     if additional_models is None:
         additional_models = []
 
+    censoring_type = CensoringType.MAP[censoring_type]
+
     evaluation_lookup = {
         "AIC": lambda model: 2 * len(model._fitted_parameter_names) - 2 * model.log_likelihood_,
-        "BIC": lambda model: 2 * len(model._fitted_parameter_names) - 2 * model.log_likelihood_ * np.log(event_times.shape[0]),
+        "BIC": lambda model: 2 * len(model._fitted_parameter_names)
+        - 2 * model.log_likelihood_ * np.log(model.event_observed.shape[0]),
     }
 
     eval = evaluation_lookup[scoring_method]
 
+    if censoring_type != CensoringType.INTERVAL:
+        event_times = (event_times,)
+
+    n = event_times[0].shape
+
     if event_observed is None:
-        event_observed = np.ones_like(event_times, dtype=bool)
+        if censoring_type == CensoringType.INTERVAL:
+            event_observed = event_times[0] == event_times[1]
+        else:
+            event_observed = np.ones(n, dtype=bool)
 
-    observed_T = event_times[event_observed.astype(bool)]
-    knots1 = np.percentile(observed_T, 100 * np.linspace(0.05, 0.95, 3))
-    knots2 = np.percentile(observed_T, 100 * np.linspace(0.05, 0.95, 4))
-    knots3 = np.percentile(observed_T, 100 * np.linspace(0.05, 0.95, 5))
+    try:
+        observed_T = event_times[0][event_observed.astype(bool)]
+        knots1 = np.percentile(observed_T, 100 * np.linspace(0.05, 0.95, 3))
+        knots2 = np.percentile(observed_T, 100 * np.linspace(0.05, 0.95, 4))
+        knots3 = np.percentile(observed_T, 100 * np.linspace(0.05, 0.95, 5))
+    except:
+        observed_T = event_times[0]
+        knots1 = np.percentile(observed_T, 100 * np.linspace(0.05, 0.95, 3))
+        knots2 = np.percentile(observed_T, 100 * np.linspace(0.05, 0.95, 4))
+        knots3 = np.percentile(observed_T, 100 * np.linspace(0.05, 0.95, 5))
 
     best_model = None
     best_score = np.inf
@@ -1745,9 +1794,19 @@ def find_best_parametric_model(event_times, event_observed=None, scoring_method:
         try:
             with warnings.catch_warnings():
                 warnings.simplefilter("ignore")
-                model.fit(event_times, event_observed)
+                getattr(model, "fit_" + CensoringType.HUMAN_MAP[censoring_type] + "_censoring")(
+                    *event_times,
+                    event_observed=event_observed,
+                    weights=weights,
+                    entry=entry,
+                    alpha=alpha,
+                    ci_labels=ci_labels,
+                    timeline=timeline
+                )
             score_ = eval(model)
 
+            if show_progress:
+                print(model._label, ", Score: %.2f" % score_)
             if score_ < best_score:
                 best_score = score_
                 best_model = model

From 39f9d192e4cdca8784496109ecdf2bf4228f37c1 Mon Sep 17 00:00:00 2001
From: CamDavidsonPilon <cam.davidson.pilon@gmail.com>
Date: Fri, 15 May 2020 11:55:57 -0400
Subject: [PATCH 07/10] adding new functionality to find_best..

---
 lifelines/fitters/__init__.py                |  27 ++--
 lifelines/fitters/aalen_additive_fitter.py   |  17 +--
 lifelines/fitters/cox_time_varying_fitter.py |   1 -
 lifelines/tests/test_estimation.py           | 132 +++++++++++++++++++
 lifelines/tests/utils/test_utils.py          |   2 +-
 5 files changed, 158 insertions(+), 21 deletions(-)

diff --git a/lifelines/fitters/__init__.py b/lifelines/fitters/__init__.py
index 4da7d9a8b..343c29bea 100644
--- a/lifelines/fitters/__init__.py
+++ b/lifelines/fitters/__init__.py
@@ -46,6 +46,7 @@ def __init__(self, alpha: float = 0.05, label: str = None):
         self.alpha = alpha
         self._class_name = self.__class__.__name__
         self._label = label
+        self._censoring_type = None
 
     def __repr__(self) -> str:
         classname = self._class_name
@@ -487,10 +488,12 @@ def _fit_model(self, Ts, E, entry, weights, show_progress=True):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
 
-            minimizing_results, minimizing_ll = None, np.inf
-            for method in [self._scipy_fit_method, "Nelder-Mead"]:
+            minimizing_results, previous_results, minimizing_ll = None, None, np.inf
+            for method, option in zip(
+                ["Nelder-Mead", self._scipy_fit_method], [{"maxiter": 25}, {**{"disp": show_progress}, **self._scipy_fit_options}]
+            ):
 
-                initial_value = self._initial_values if minimizing_results is None else utils._to_1d_array(minimizing_results.x)
+                initial_value = self._initial_values if previous_results is None else utils._to_1d_array(previous_results.x)
 
                 results = minimize(
                     value_and_grad(negative_log_likelihood),  # pylint: disable=no-value-for-parameter
@@ -499,15 +502,18 @@ def _fit_model(self, Ts, E, entry, weights, show_progress=True):
                     method=method,
                     args=(Ts, E, entry, weights),
                     bounds=self._bounds,
-                    options={**{"disp": show_progress}, **self._scipy_fit_options},
+                    options=option,
                 )
+                previous_results = results
 
-                if results.success and (results.fun < minimizing_ll):
+                if results.success and ~np.isnan(results.x).any() and (results.fun < minimizing_ll):
                     minimizing_ll = results.fun
                     minimizing_results = results
 
             # convergence successful.
-            if minimizing_results and minimizing_results.success:
+            # I still need to check for ~np.isnan(minimizing_results.x).any() since minimize will happily
+            # return nans even when criteria is satisified.
+            if minimizing_results and minimizing_results.success and ~np.isnan(minimizing_results.x).any():
                 sol = utils._to_1d_array(minimizing_results.x)
                 # pylint: disable=no-value-for-parameter
                 hessian_ = hessian(negative_log_likelihood)(sol, Ts, E, entry, weights)
@@ -516,7 +522,8 @@ def _fit_model(self, Ts, E, entry, weights, show_progress=True):
                 return sol, -minimizing_results.fun * weights.sum(), hessian_ * weights.sum()
 
             # convergence failed.
-            print(minimizing_results)
+            if show_progress:
+                print(minimizing_results)
             if self._KNOWN_MODEL:
                 raise utils.ConvergenceError(
                     dedent(
@@ -631,6 +638,7 @@ def print_summary(self, decimals=2, style=None, **kwargs):
                     ),
                 ),
             ],
+            [],
             justify,
             decimals,
             kwargs,
@@ -2117,9 +2125,12 @@ def print_summary(self, decimals=2, style=None, **kwargs):
 
         sr = self.log_likelihood_ratio_test()
         footers = []
+
+        if CensoringType.is_right_censoring(self):
+            footers.apoend(("Concordance", "{:.{prec}f}".format(self.concordance_index_, prec=decimals)))
+
         footers.extend(
             [
-                ("Concordance", "{:.{prec}f}".format(self.concordance_index_, prec=decimals)),
                 ("AIC", "{:.{prec}f}".format(self.AIC_, prec=decimals)),
                 (
                     "log-likelihood ratio test",
diff --git a/lifelines/fitters/aalen_additive_fitter.py b/lifelines/fitters/aalen_additive_fitter.py
index d01178a7f..cfe2b0f43 100644
--- a/lifelines/fitters/aalen_additive_fitter.py
+++ b/lifelines/fitters/aalen_additive_fitter.py
@@ -201,9 +201,7 @@ def _fit_model(self, X, T, E, weights, show_progress):
 
         hazards = pd.DataFrame(hazards_, columns=columns, index=index).iloc[:stop]
         cumulative_hazards_ = hazards.cumsum()
-        cumulative_variance_hazards_ = (
-            pd.DataFrame(variance_hazards_, columns=columns, index=index).iloc[:stop].cumsum()
-        )
+        cumulative_variance_hazards_ = pd.DataFrame(variance_hazards_, columns=columns, index=index).iloc[:stop].cumsum()
 
         return hazards, cumulative_hazards_, cumulative_variance_hazards_
 
@@ -299,9 +297,7 @@ def _preprocess_dataframe(self, df):
         self._check_values(df, T, E)
 
         if self.fit_intercept:
-            assert (
-                "_intercept" not in df.columns
-            ), "_intercept is an internal lifelines column, please rename your column first."
+            assert "_intercept" not in df.columns, "_intercept is an internal lifelines column, please rename your column first."
             X["_intercept"] = 1.0
 
         return X, T, E, W
@@ -333,9 +329,7 @@ def predict_cumulative_hazard(self, X):
         X_ = X_ if not self.fit_intercept else np.c_[X_, np.ones((n, 1))]
 
         timeline = self._index
-        individual_cumulative_hazards_ = pd.DataFrame(
-            np.dot(self.cumulative_hazards_, X_.T), index=timeline, columns=cols
-        )
+        individual_cumulative_hazards_ = pd.DataFrame(np.dot(self.cumulative_hazards_, X_.T), index=timeline, columns=cols)
 
         return individual_cumulative_hazards_
 
@@ -490,7 +484,7 @@ def smoothed_hazards_(self, bandwidth=1):
         )
 
     @property
-    def score_(self):
+    def concordance_index_(self):
         """
         The concordance score (also known as the c-index) of the fit.  The c-index is a generalization of the ROC AUC
         to survival data, including censorships.
@@ -572,7 +566,8 @@ def print_summary(self, decimals=2, style=None, **kwargs):
             ]
         )
 
-        p = Printer(self, headers, justify, decimals, kwargs)
+        footers = [("Concordance", "{:.{prec}f}".format(self.concordance_index_, prec=decimals))]
+        p = Printer(self, headers, footers, justify, decimals, kwargs)
 
         p.print(style=style)
 
diff --git a/lifelines/fitters/cox_time_varying_fitter.py b/lifelines/fitters/cox_time_varying_fitter.py
index 71f568795..3f5fbf57c 100644
--- a/lifelines/fitters/cox_time_varying_fitter.py
+++ b/lifelines/fitters/cox_time_varying_fitter.py
@@ -670,7 +670,6 @@ def print_summary(self, decimals=2, style=None, **kwargs):
         footers = []
         footers.extend(
             [
-                ("Concordance", "{:.{prec}f}".format(self.concordance_index_, prec=decimals)),
                 ("Partial AIC", "{:.{prec}f}".format(self.AIC_partial_, prec=decimals)),
                 (
                     "log-likelihood ratio test",
diff --git a/lifelines/tests/test_estimation.py b/lifelines/tests/test_estimation.py
index b09239bfd..ad15c1f63 100644
--- a/lifelines/tests/test_estimation.py
+++ b/lifelines/tests/test_estimation.py
@@ -1357,6 +1357,138 @@ def test_late_entry_with_tied_entry_and_death(self):
         npt.assert_allclose(rf["KM_lifelines_latest"].values, rf["KM_lateenterafter"].values, rtol=10e-2)
         npt.assert_allclose(rf["KM_lifelines_latest"].values, rf["KM_true"].values, rtol=10e-2)
 
+    def test_interval_censoring_to_r_1(self):
+        kmf = KaplanMeierFitter()
+        left = [1, 7, 8, 7, 7, 17, 37, 46, 46, 45]
+        right = [7, 8, 10, 16, 14, 100, 44, 100, 100, 100]
+
+        kmf.fit_interval_censoring(left, right)
+
+        npt.assert_allclose(kmf.survival_function_.loc[7.0].values, np.array([0.9, 1.0]))
+        npt.assert_allclose(kmf.survival_function_.loc[8.0].values, np.array([0.7, 0.9]))
+        npt.assert_allclose(kmf.survival_function_.loc[10.0].values, np.array([0.5, 0.7]))
+        npt.assert_allclose(kmf.survival_function_.loc[44.0].values, np.array([0.375, 0.5]), rtol=1e-3)
+        npt.assert_allclose(kmf.survival_function_.iloc[-1].values, np.array([0.0, 0.375]), atol=1e-3)
+
+    def test_interval_censoring_to_r_2_test_observed_event(self):
+        kmf = KaplanMeierFitter()
+        left = [1, 8, 8, 7, 7, 17, 37, 46, 46, 45]
+        right = [7, 8, 10, 16, 14, 100, 44, 100, 100, 100]
+
+        kmf.fit_interval_censoring(left, right)
+
+        npt.assert_allclose(kmf.survival_function_.loc[7.0].values, np.array([0.9, 1.0]))
+        npt.assert_allclose(kmf.survival_function_.loc[8.0].values, np.array([0.5, 0.9]))
+        npt.assert_allclose(kmf.survival_function_.loc[10.0].values, np.array([0.5, 0.5]))
+        npt.assert_allclose(kmf.survival_function_.loc[44.0].values, np.array([0.375, 0.5]), rtol=1e-3)
+        npt.assert_allclose(kmf.survival_function_.iloc[-1].values, np.array([0.0, 0.375]), atol=1e-3)
+
+    def test_interval_censoring_to_r_3_test_0_and_inf(self):
+        kmf = KaplanMeierFitter()
+        left = [0, 8, 8, 7, 7, 17, 37, 46, 46, 45]
+        right = [7, 8, 10, 16, 14, np.inf, 44, np.inf, np.inf, np.inf]
+
+        kmf.fit_interval_censoring(left, right)
+
+        npt.assert_allclose(kmf.survival_function_.loc[0.0].values, np.array([1.0, 1.0]))
+        npt.assert_allclose(kmf.survival_function_.loc[7.0].values, np.array([0.9, 1.0]))
+        npt.assert_allclose(kmf.survival_function_.loc[8.0].values, np.array([0.5, 0.9]))
+        npt.assert_allclose(kmf.survival_function_.loc[10.0].values, np.array([0.5, 0.5]))
+        npt.assert_allclose(kmf.survival_function_.loc[44.0].values, np.array([0.375, 0.5]), rtol=1e-3)
+        npt.assert_allclose(kmf.survival_function_.iloc[-1].values, np.array([0.0, 0.375]), atol=1e-3)
+
+    def test_interval_censoring_to_r_3_test_ties_and_overlapping_intervals(self):
+        kmf = KaplanMeierFitter()
+        left = [6, 7, 8, 7, 5]
+        right = [7, 8, 10, 16, 20]
+
+        kmf.fit_interval_censoring(left, right)
+
+        npt.assert_allclose(kmf.survival_function_.loc[5.0].values, np.array([1.0, 1.0]))
+        npt.assert_allclose(kmf.survival_function_.loc[6.0].values, np.array([1.0, 1.0]))
+        npt.assert_allclose(kmf.survival_function_.loc[7.0].values, np.array([0.75, 1.0]), rtol=1e-05)
+        npt.assert_allclose(kmf.survival_function_.loc[8.0].values, np.array([0.375, 0.75]), rtol=1e-05)
+        npt.assert_allclose(kmf.survival_function_.loc[10.0].values, np.array([0, 0.375]), rtol=1e-05)
+
+    def test_interval_censoring_with_custom_index(self):
+        kmf = KaplanMeierFitter()
+        left = [6, 7, 8, 7, 5]
+        right = [7, 8, 10, 16, 20]
+
+        kmf.fit_interval_censoring(left, right, timeline=np.arange(10))
+        npt.assert_allclose(kmf.survival_function_.index.values, np.arange(10))
+
+    def test_interval_censoring_fit_with_exact_observations_is_equal_to_km(self):
+        left = np.array([6, 7, 8])
+        right = np.array([6, 7, 8])
+
+        kmf_interval = KaplanMeierFitter()
+        kmf_interval.fit_interval_censoring(left, right)
+
+        kmf_right = KaplanMeierFitter().fit(left, left == right)
+
+        npt.assert_allclose(
+            kmf_interval.survival_function_["NPMLE_estimate_lower"],
+            kmf_right.survival_function_["KM_estimate"],
+            rtol=1e-3,
+            atol=1e-3,
+        )
+
+    def test_interval_censoring_fit_with_exact_observations_is_equal_to_km2(self):
+        left = np.array([6, 7, 8, 9])
+        right = np.array([6, 7, 8, np.inf])
+
+        kmf_interval = KaplanMeierFitter()
+        kmf_interval.fit_interval_censoring(left, right)
+
+        kmf_right = KaplanMeierFitter().fit(left, left == right)
+
+        npt.assert_allclose(
+            kmf_interval.survival_function_["NPMLE_estimate_lower"].iloc[:-1],
+            kmf_right.survival_function_["KM_estimate"],
+            rtol=1e-3,
+            atol=1e-3,
+        )
+
+    def test_interval_censoring_fit_with_exact_observations_is_equal_to_km3(self):
+        left = np.array([6, 7, 8, 7.5])
+        right = np.array([6, 7, 8, np.inf])
+
+        kmf_interval = KaplanMeierFitter()
+        kmf_interval.fit_interval_censoring(left, right)
+
+        kmf_right = KaplanMeierFitter().fit(left, left == right)
+
+        npt.assert_allclose(
+            kmf_interval.survival_function_["NPMLE_estimate_lower"].iloc[:-1],
+            kmf_right.survival_function_["KM_estimate"],
+            rtol=1e-3,
+            atol=1e-3,
+        )
+
+    def test_interval_censoring_fit_with_exact_observations_is_equal_to_km4(self):
+        left = np.array([6, 7, 8, 7.0])
+        right = np.array([6, 7, 8, np.inf])
+
+        kmf_interval = KaplanMeierFitter()
+        kmf_interval.fit_interval_censoring(left, right)
+
+        kmf_right = KaplanMeierFitter().fit(left, left == right)
+
+        npt.assert_allclose(
+            kmf_interval.survival_function_["NPMLE_estimate_lower"].iloc[:-1],
+            kmf_right.survival_function_["KM_estimate"],
+            rtol=1e-3,
+            atol=1e-3,
+        )
+
+    def test_interval_model_has_all_attributes_as_non_interval(self):
+        left = [6, 7, 8, 7, 5]
+        right = [7, 8, 10, 16, 20]
+        kmf_interval = KaplanMeierFitter().fit(left, right)
+
+        kmf_right = KaplanMeierFitter().fit(np.arange(10))
+
 
 class TestNelsonAalenFitter:
     def nelson_aalen(self, lifetimes, observed=None):
diff --git a/lifelines/tests/utils/test_utils.py b/lifelines/tests/utils/test_utils.py
index e616ab3da..5a98befa2 100644
--- a/lifelines/tests/utils/test_utils.py
+++ b/lifelines/tests/utils/test_utils.py
@@ -1079,7 +1079,7 @@ def test_find_best_parametric_model_works_for_interval_censoring():
 
 
 def test_find_best_parametric_model_works_with_weights_and_entry():
-    T = np.random.exponential(2, 100)
+    T = np.random.exponential(5, 100)
     W = np.random.randint(1, 5, size=100)
     entry = np.random.exponential(0.01, 100)
     model, score = utils.find_best_parametric_model(T, weights=W, entry=entry, show_progress=True)

From 79367410b43a54f3fc7674067fa962ed55b1d6bb Mon Sep 17 00:00:00 2001
From: CamDavidsonPilon <cam.davidson.pilon@gmail.com>
Date: Sat, 16 May 2020 10:11:07 -0400
Subject: [PATCH 08/10] prepare for release

---
 CHANGELOG.md                             |  14 +
 docs/Changelog.rst                       | 367 ++++++++++++-----------
 docs/Survival Regression.rst             |  32 +-
 lifelines/fitters/__init__.py            |   4 +
 lifelines/fitters/kaplan_meier_fitter.py |  57 +++-
 lifelines/fitters/npmle.py               |  66 +++-
 lifelines/plotting.py                    |  16 +-
 lifelines/statistics.py                  |   1 +
 lifelines/tests/test_plotting.py         |   9 +
 lifelines/utils/concordance.py           |  43 ++-
 lifelines/version.py                     |   2 +-
 11 files changed, 412 insertions(+), 199 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1ecfc2716..e16650b3f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,19 @@
 ## Changelog
 
+#### 0.24.7
+
+##### New features
+ - `find_best_parametric_model` can handle left and interval censoring. Also allows for more fitting options.
+ - `AIC_` is a property on parametric models, and `AIC_partial_` is a property on Cox models.
+ - `penalizer` in all regression models can now be an array instead of a float. This enables new functionality and better
+ control over penalization. This is similar (but not identical) to `penalty.factors` in glmnet in R.
+ - some convergence tweaks which should help recent performance regressions.
+
+##### Bug fixes
+ - fixed bug where `cdf_plot` and `qq_plot` were not factoring in the weights correctly.
+
+
+
 #### 0.24.6 - 2020-05-05
 
 ##### New features
diff --git a/docs/Changelog.rst b/docs/Changelog.rst
index 64ba664b6..0d4bbb929 100644
--- a/docs/Changelog.rst
+++ b/docs/Changelog.rst
@@ -1,9 +1,36 @@
 Changelog
 ---------
 
+0.24.7
+^^^^^^
+
+New features
+''''''''''''
+
+-  ``find_best_parametric_model`` can handle left and interval
+   censoring. Also allows for more fitting options.
+-  ``AIC_`` is a property on parametric models, and ``AIC_partial_`` is
+   a property on Cox models.
+-  ``penalizer`` in all regression models can now be an array instead of
+   a float. This enables new functionality and better control over
+   penalization. This is similar (but not identical) to
+   ``penalty.factors`` in glmnet in R.
+-  some convergence tweaks which should help recent performance
+   regressions.
+
+Bug fixes
+'''''''''
+
+-  fixed bug where ``cdf_plot`` and ``qq_plot`` were not factoring in
+   the weights correctly.
+
+.. _section-1:
+
 0.24.6 - 2020-05-05
 ^^^^^^^^^^^^^^^^^^^
 
+.. _new-features-1:
+
 New features
 ''''''''''''
 
@@ -12,25 +39,27 @@ New features
 -  New ``lifelines.plotting.plot_interval_censored_lifetimes`` for
    plotting interval censored data - thanks @sean-reed!
 
+.. _bug-fixes-1:
+
 Bug fixes
 '''''''''
 
 -  fixed bug where ``cdf_plot`` and ``qq_plot`` were not factoring in
    the weights correctly.
 
-.. _section-1:
+.. _section-2:
 
 0.24.5 - 2020-05-01
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-1:
+.. _new-features-2:
 
 New features
 ''''''''''''
 
 -  ``plot_lifetimes`` accepts pandas Series.
 
-.. _bug-fixes-1:
+.. _bug-fixes-2:
 
 Bug fixes
 '''''''''
@@ -40,12 +69,12 @@ Bug fixes
 -  Improved ``at_risk_counts`` for subplots.
 -  More data validation checks for ``CoxTimeVaryingFitter``
 
-.. _section-2:
+.. _section-3:
 
 0.24.4 - 2020-04-13
 ^^^^^^^^^^^^^^^^^^^
 
-.. _bug-fixes-2:
+.. _bug-fixes-3:
 
 Bug fixes
 '''''''''
@@ -54,12 +83,12 @@ Bug fixes
 -  setting a dataframe in ``ancillary_df`` works for interval censoring
 -  ``.score`` works for interval censored models
 
-.. _section-3:
+.. _section-4:
 
 0.24.3 - 2020-03-25
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-2:
+.. _new-features-3:
 
 New features
 ''''''''''''
@@ -69,7 +98,7 @@ New features
    the hazard ratio would be at previous times. This is useful because
    the final hazard ratio is some weighted average of these.
 
-.. _bug-fixes-3:
+.. _bug-fixes-4:
 
 Bug fixes
 '''''''''
@@ -77,12 +106,12 @@ Bug fixes
 -  Fixed error in HTML printer that was hiding concordance index
    information.
 
-.. _section-4:
+.. _section-5:
 
 0.24.2 - 2020-03-15
 ^^^^^^^^^^^^^^^^^^^
 
-.. _bug-fixes-4:
+.. _bug-fixes-5:
 
 Bug fixes
 '''''''''
@@ -94,12 +123,12 @@ Bug fixes
 -  Fixed a keyword bug in ``plot_covariate_groups`` for parametric
    models.
 
-.. _section-5:
+.. _section-6:
 
 0.24.1 - 2020-03-05
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-3:
+.. _new-features-4:
 
 New features
 ''''''''''''
@@ -107,14 +136,14 @@ New features
 -  Stability improvements for GeneralizedGammaRegressionFitter and
    CoxPHFitter with spline estimation.
 
-.. _bug-fixes-5:
+.. _bug-fixes-6:
 
 Bug fixes
 '''''''''
 
 -  Fixed bug with plotting hazards in NelsonAalenFitter.
 
-.. _section-6:
+.. _section-7:
 
 0.24.0 - 2020-02-20
 ^^^^^^^^^^^^^^^^^^^
@@ -123,7 +152,7 @@ This version and future versions of lifelines no longer support py35.
 Pandas 1.0 is fully supported, along with previous versions. Minimum
 Scipy has been bumped to 1.2.0.
 
-.. _new-features-4:
+.. _new-features-5:
 
 New features
 ''''''''''''
@@ -173,7 +202,7 @@ API Changes
    to ``scoring_method``.
 -  removed ``_score_`` and ``path`` from Cox model.
 
-.. _bug-fixes-6:
+.. _bug-fixes-7:
 
 Bug fixes
 '''''''''
@@ -186,12 +215,12 @@ Bug fixes
 -  Cox models now incorporate any penalizers in their
    ``log_likelihood_``
 
-.. _section-7:
+.. _section-8:
 
 0.23.9 - 2020-01-28
 ^^^^^^^^^^^^^^^^^^^
 
-.. _bug-fixes-7:
+.. _bug-fixes-8:
 
 Bug fixes
 '''''''''
@@ -202,12 +231,12 @@ Bug fixes
    of ``GeneralizedGammaRegressionFitter`` and any custom regression
    models should update their code as soon as possible.
 
-.. _section-8:
+.. _section-9:
 
 0.23.8 - 2020-01-21
 ^^^^^^^^^^^^^^^^^^^
 
-.. _bug-fixes-8:
+.. _bug-fixes-9:
 
 Bug fixes
 '''''''''
@@ -218,19 +247,19 @@ Bug fixes
    ``GeneralizedGammaRegressionFitter`` and any custom regression models
    should update their code as soon as possible.
 
-.. _section-9:
+.. _section-10:
 
 0.23.7 - 2020-01-14
 ^^^^^^^^^^^^^^^^^^^
 
 Bug fixes for py3.5.
 
-.. _section-10:
+.. _section-11:
 
 0.23.6 - 2020-01-07
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-5:
+.. _new-features-6:
 
 New features
 ''''''''''''
@@ -244,12 +273,12 @@ New features
 -  custom parametric regression models can now do left and interval
    censoring.
 
-.. _section-11:
+.. _section-12:
 
 0.23.5 - 2020-01-05
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-6:
+.. _new-features-7:
 
 New features
 ''''''''''''
@@ -258,7 +287,7 @@ New features
 -  New lymph node cancer dataset, originally from *H.F. for the German
    Breast Cancer Study Group (GBSG) (1994)*
 
-.. _bug-fixes-9:
+.. _bug-fixes-10:
 
 Bug fixes
 '''''''''
@@ -268,26 +297,26 @@ Bug fixes
 -  fixed bug where large exponential numbers in ``print_summary`` were
    not being suppressed correctly.
 
-.. _section-12:
+.. _section-13:
 
 0.23.4 - 2019-12-15
 ^^^^^^^^^^^^^^^^^^^
 
 -  Bug fix for PyPI
 
-.. _section-13:
+.. _section-14:
 
 0.23.3 - 2019-12-11
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-7:
+.. _new-features-8:
 
 New features
 ''''''''''''
 
 -  ``StatisticalResult.print_summary`` supports html output.
 
-.. _bug-fixes-10:
+.. _bug-fixes-11:
 
 Bug fixes
 '''''''''
@@ -295,12 +324,12 @@ Bug fixes
 -  fix import in ``printer.py``
 -  fix html printing with Univariate models.
 
-.. _section-14:
+.. _section-15:
 
 0.23.2 - 2019-12-07
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-8:
+.. _new-features-9:
 
 New features
 ''''''''''''
@@ -312,7 +341,7 @@ New features
 -  performance improvements on regression models’ preprocessing. Should
    make datasets with high number of columns more performant.
 
-.. _bug-fixes-11:
+.. _bug-fixes-12:
 
 Bug fixes
 '''''''''
@@ -321,12 +350,12 @@ Bug fixes
 -  fixed repr for ``sklearn_adapter`` classes.
 -  fixed ``conditional_after`` in Cox model with strata was used.
 
-.. _section-15:
+.. _section-16:
 
 0.23.1 - 2019-11-27
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-9:
+.. _new-features-10:
 
 New features
 ''''''''''''
@@ -336,7 +365,7 @@ New features
 -  performance improvements for ``CoxPHFitter`` - up to 30% performance
    improvements for some datasets.
 
-.. _bug-fixes-12:
+.. _bug-fixes-13:
 
 Bug fixes
 '''''''''
@@ -348,12 +377,12 @@ Bug fixes
 -  fixed bug when using ``print_summary`` with left censored models.
 -  lots of minor bug fixes.
 
-.. _section-16:
+.. _section-17:
 
 0.23.0 - 2019-11-17
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-10:
+.. _new-features-11:
 
 New features
 ''''''''''''
@@ -362,7 +391,7 @@ New features
    Jupyter notebooks!
 -  silenced some warnings.
 
-.. _bug-fixes-13:
+.. _bug-fixes-14:
 
 Bug fixes
 '''''''''
@@ -384,7 +413,7 @@ API Changes
 -  ``left_censorship`` in ``fit`` has been removed in favour of
    ``fit_left_censoring``.
 
-.. _section-17:
+.. _section-18:
 
 0.22.10 - 2019-11-08
 ^^^^^^^^^^^^^^^^^^^^
@@ -392,7 +421,7 @@ API Changes
 The tests were re-factored to be shipped with the package. Let me know
 if this causes problems.
 
-.. _bug-fixes-14:
+.. _bug-fixes-15:
 
 Bug fixes
 '''''''''
@@ -402,12 +431,12 @@ Bug fixes
 -  fixed bug in plot_covariate_groups for AFT models when >1d arrays
    were used for values arg.
 
-.. _section-18:
+.. _section-19:
 
 0.22.9 - 2019-10-30
 ^^^^^^^^^^^^^^^^^^^
 
-.. _bug-fixes-15:
+.. _bug-fixes-16:
 
 Bug fixes
 '''''''''
@@ -419,12 +448,12 @@ Bug fixes
 -  ``CoxPHFitter`` now displays correct columns values when changing
    alpha param.
 
-.. _section-19:
+.. _section-20:
 
 0.22.8 - 2019-10-06
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-11:
+.. _new-features-12:
 
 New features
 ''''''''''''
@@ -434,19 +463,19 @@ New features
 -  ``conditional_after`` now available in ``CoxPHFitter.predict_median``
 -  Suppressed some unimportant warnings.
 
-.. _bug-fixes-16:
+.. _bug-fixes-17:
 
 Bug fixes
 '''''''''
 
 -  fixed initial_point being ignored in AFT models.
 
-.. _section-20:
+.. _section-21:
 
 0.22.7 - 2019-09-29
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-12:
+.. _new-features-13:
 
 New features
 ''''''''''''
@@ -454,7 +483,7 @@ New features
 -  new ``ApproximationWarning`` to tell you if the package is making an
    potentially mislead approximation.
 
-.. _bug-fixes-17:
+.. _bug-fixes-18:
 
 Bug fixes
 '''''''''
@@ -473,19 +502,19 @@ API Changes
 -  Some previous ``StatisticalWarnings`` have been replaced by
    ``ApproximationWarning``
 
-.. _section-21:
+.. _section-22:
 
 0.22.6 - 2019-09-25
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-13:
+.. _new-features-14:
 
 New features
 ''''''''''''
 
 -  ``conditional_after`` works for ``CoxPHFitter`` prediction models 😅
 
-.. _bug-fixes-18:
+.. _bug-fixes-19:
 
 Bug fixes
 '''''''''
@@ -501,12 +530,12 @@ API Changes
 -  ``utils.dataframe_interpolate_at_times`` renamed to
    ``utils.interpolate_at_times_and_return_pandas``.
 
-.. _section-22:
+.. _section-23:
 
 0.22.5 - 2019-09-20
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-14:
+.. _new-features-15:
 
 New features
 ''''''''''''
@@ -515,7 +544,7 @@ New features
    weights.
 -  Better support for predicting on Pandas Series
 
-.. _bug-fixes-19:
+.. _bug-fixes-20:
 
 Bug fixes
 '''''''''
@@ -532,12 +561,12 @@ API Changes
 -  ``_get_initial_value`` in parametric univariate models is renamed
    ``_create_initial_point``
 
-.. _section-23:
+.. _section-24:
 
 0.22.4 - 2019-09-04
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-15:
+.. _new-features-16:
 
 New features
 ''''''''''''
@@ -556,7 +585,7 @@ API changes
 -  ``KaplanMeierFitter.survival_function_``\ ‘s’ index is no longer
    given the name “timeline”.
 
-.. _bug-fixes-20:
+.. _bug-fixes-21:
 
 Bug fixes
 '''''''''
@@ -564,12 +593,12 @@ Bug fixes
 -  Fixed issue where ``concordance_index`` would never exit if NaNs in
    dataset.
 
-.. _section-24:
+.. _section-25:
 
 0.22.3 - 2019-08-08
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-16:
+.. _new-features-17:
 
 New features
 ''''''''''''
@@ -593,7 +622,7 @@ API changes
    gains only in Cox models, and only a small fraction of the API was
    being used.
 
-.. _bug-fixes-21:
+.. _bug-fixes-22:
 
 Bug fixes
 '''''''''
@@ -605,19 +634,19 @@ Bug fixes
 -  Fixed an error in the ``predict_percentile`` of
    ``LogLogisticAFTFitter``. New tests have been added around this.
 
-.. _section-25:
+.. _section-26:
 
 0.22.2 - 2019-07-25
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-17:
+.. _new-features-18:
 
 New features
 ''''''''''''
 
 -  lifelines is now compatible with scipy>=1.3.0
 
-.. _bug-fixes-22:
+.. _bug-fixes-23:
 
 Bug fixes
 '''''''''
@@ -628,12 +657,12 @@ Bug fixes
    errors when using the library. The correctly numpy has been pinned
    (to 1.14.0+)
 
-.. _section-26:
+.. _section-27:
 
 0.22.1 - 2019-07-14
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-18:
+.. _new-features-19:
 
 New features
 ''''''''''''
@@ -661,7 +690,7 @@ API changes
    ``.print_summary`` includes confidence intervals for the exponential
    of the value.
 
-.. _bug-fixes-23:
+.. _bug-fixes-24:
 
 Bug fixes
 '''''''''
@@ -671,12 +700,12 @@ Bug fixes
 -  fixed an overflow bug in ``KaplanMeierFitter`` confidence intervals
 -  improvements in data validation for ``CoxTimeVaryingFitter``
 
-.. _section-27:
+.. _section-28:
 
 0.22.0 - 2019-07-03
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-19:
+.. _new-features-20:
 
 New features
 ''''''''''''
@@ -710,7 +739,7 @@ API changes
    could set ``fit_intercept`` to False and not have to set
    ``ancillary_df`` - now one must specify a DataFrame.
 
-.. _bug-fixes-24:
+.. _bug-fixes-25:
 
 Bug fixes
 '''''''''
@@ -719,21 +748,21 @@ Bug fixes
    is now exact instead of an approximation.
 -  fixed a name error bug in ``CoxTimeVaryingFitter.plot``
 
-.. _section-28:
+.. _section-29:
 
 0.21.5 - 2019-06-22
 ^^^^^^^^^^^^^^^^^^^
 
 I’m skipping 0.21.4 version because of deployment issues.
 
-.. _new-features-20:
+.. _new-features-21:
 
 New features
 ''''''''''''
 
 -  ``scoring_method`` now a kwarg on ``sklearn_adapter``
 
-.. _bug-fixes-25:
+.. _bug-fixes-26:
 
 Bug fixes
 '''''''''
@@ -743,12 +772,12 @@ Bug fixes
 -  fixed visual bug that misaligned x-axis ticks and at-risk counts.
    Thanks @christopherahern!
 
-.. _section-29:
+.. _section-30:
 
 0.21.3 - 2019-06-04
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-21:
+.. _new-features-22:
 
 New features
 ''''''''''''
@@ -762,19 +791,19 @@ New features
 -  ``CoxPHFitter.check_assumptions`` now accepts a ``columns`` parameter
    to specify only checking a subset of columns.
 
-.. _bug-fixes-26:
+.. _bug-fixes-27:
 
 Bug fixes
 '''''''''
 
 -  ``covariates_from_event_matrix`` handle nulls better
 
-.. _section-30:
+.. _section-31:
 
 0.21.2 - 2019-05-16
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-22:
+.. _new-features-23:
 
 New features
 ''''''''''''
@@ -798,17 +827,17 @@ API changes
 -  removing ``_compute_likelihood_ratio_test`` on regression models. Use
    ``log_likelihood_ratio_test`` now.
 
-.. _bug-fixes-27:
+.. _bug-fixes-28:
 
 Bug fixes
 '''''''''
 
-.. _section-31:
+.. _section-32:
 
 0.21.1 - 2019-04-26
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-23:
+.. _new-features-24:
 
 New features
 ''''''''''''
@@ -825,19 +854,19 @@ API changes
 -  output of ``survival_table_from_events`` when collapsing rows to
    intervals now removes the “aggregate” column multi-index.
 
-.. _bug-fixes-28:
+.. _bug-fixes-29:
 
 Bug fixes
 '''''''''
 
 -  fixed bug in CoxTimeVaryingFitter when ax is provided, thanks @j-i-l!
 
-.. _section-32:
+.. _section-33:
 
 0.21.0 - 2019-04-12
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-24:
+.. _new-features-25:
 
 New features
 ''''''''''''
@@ -862,7 +891,7 @@ API changes
 -  ``entries`` property in multivariate parametric models has a new
    Series name: ``entry``
 
-.. _bug-fixes-29:
+.. _bug-fixes-30:
 
 Bug fixes
 '''''''''
@@ -872,12 +901,12 @@ Bug fixes
 -  Fixed an error that didn’t let users use Numpy arrays in prediction
    for AFT models
 
-.. _section-33:
+.. _section-34:
 
 0.20.5 - 2019-04-08
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-25:
+.. _new-features-26:
 
 New features
 ''''''''''''
@@ -894,7 +923,7 @@ API changes
 -  in ``AalenJohansenFitter``, the ``variance`` parameter is renamed to
    ``variance_`` to align with the usual lifelines convention.
 
-.. _bug-fixes-30:
+.. _bug-fixes-31:
 
 Bug fixes
 '''''''''
@@ -903,12 +932,12 @@ Bug fixes
    test when using strata.
 -  Fixed some plotting bugs with ``AalenJohansenFitter``
 
-.. _section-34:
+.. _section-35:
 
 0.20.4 - 2019-03-27
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-26:
+.. _new-features-27:
 
 New features
 ''''''''''''
@@ -927,7 +956,7 @@ API changes
 -  Pandas is now correctly pinned to >= 0.23.0. This was always the
    case, but not specified in setup.py correctly.
 
-.. _bug-fixes-31:
+.. _bug-fixes-32:
 
 Bug fixes
 '''''''''
@@ -936,12 +965,12 @@ Bug fixes
 -  ``PiecewiseExponentialFitter`` is available with
    ``from lifelines import *``.
 
-.. _section-35:
+.. _section-36:
 
 0.20.3 - 2019-03-23
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-27:
+.. _new-features-28:
 
 New features
 ''''''''''''
@@ -954,12 +983,12 @@ New features
    ``plot_survival_function`` and
    ``confidence_interval_survival_function_``.
 
-.. _section-36:
+.. _section-37:
 
 0.20.2 - 2019-03-21
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-28:
+.. _new-features-29:
 
 New features
 ''''''''''''
@@ -983,7 +1012,7 @@ API changes
    @vpolimenov!
 -  The ``C`` column in ``load_lcd`` dataset is renamed to ``E``.
 
-.. _bug-fixes-32:
+.. _bug-fixes-33:
 
 Bug fixes
 '''''''''
@@ -999,7 +1028,7 @@ Bug fixes
    the q parameter was below the truncation limit. This should have been
    ``-np.inf``
 
-.. _section-37:
+.. _section-38:
 
 0.20.1 - 2019-03-16
 ^^^^^^^^^^^^^^^^^^^
@@ -1023,7 +1052,7 @@ API changes
    This is no longer the case. A 0 will still be added if there is a
    duration (observed or not) at 0 occurs however.
 
-.. _section-38:
+.. _section-39:
 
 0.20.0 - 2019-03-05
 ^^^^^^^^^^^^^^^^^^^
@@ -1032,7 +1061,7 @@ API changes
    recent installs where Py3.
 -  Updated minimum dependencies, specifically Matplotlib and Pandas.
 
-.. _new-features-29:
+.. _new-features-30:
 
 New features
 ''''''''''''
@@ -1052,19 +1081,19 @@ API changes
    transposed now (previous parameters where columns, now parameters are
    rows).
 
-.. _bug-fixes-33:
+.. _bug-fixes-34:
 
 Bug fixes
 '''''''''
 
 -  Fixed a bug with plotting and ``check_assumptions``.
 
-.. _section-39:
+.. _section-40:
 
 0.19.5 - 2019-02-26
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-30:
+.. _new-features-31:
 
 New features
 ''''''''''''
@@ -1074,24 +1103,24 @@ New features
    features or categorical variables.
 -  Convergence improvements for AFT models.
 
-.. _section-40:
+.. _section-41:
 
 0.19.4 - 2019-02-25
 ^^^^^^^^^^^^^^^^^^^
 
-.. _bug-fixes-34:
+.. _bug-fixes-35:
 
 Bug fixes
 '''''''''
 
 -  remove some bad print statements in ``CoxPHFitter``.
 
-.. _section-41:
+.. _section-42:
 
 0.19.3 - 2019-02-25
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-31:
+.. _new-features-32:
 
 New features
 ''''''''''''
@@ -1103,12 +1132,12 @@ New features
 -  Performance increase to ``print_summary`` in the ``CoxPHFitter`` and
    ``CoxTimeVaryingFitter`` model.
 
-.. _section-42:
+.. _section-43:
 
 0.19.2 - 2019-02-22
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-32:
+.. _new-features-33:
 
 New features
 ''''''''''''
@@ -1116,7 +1145,7 @@ New features
 -  ``ParametricUnivariateFitters``, like ``WeibullFitter``, have
    smoothed plots when plotting (vs stepped plots)
 
-.. _bug-fixes-35:
+.. _bug-fixes-36:
 
 Bug fixes
 '''''''''
@@ -1126,12 +1155,12 @@ Bug fixes
 -  Univariate fitters are more flexiable and can allow 2-d and
    DataFrames as inputs.
 
-.. _section-43:
+.. _section-44:
 
 0.19.1 - 2019-02-21
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-33:
+.. _new-features-34:
 
 New features
 ''''''''''''
@@ -1148,12 +1177,12 @@ API changes
    ``PiecewiseExponential`` to the same as ``ExponentialFitter`` (from
    ``\lambda * t`` to ``t / \lambda``).
 
-.. _section-44:
+.. _section-45:
 
 0.19.0 - 2019-02-20
 ^^^^^^^^^^^^^^^^^^^
 
-.. _new-features-34:
+.. _new-features-35:
 
 New features
 ''''''''''''
@@ -1192,7 +1221,7 @@ API changes
    means that the *default* for alpha is set to 0.05 in the latest
    lifelines, instead of 0.95 in previous versions.
 
-.. _bug-fixes-36:
+.. _bug-fixes-37:
 
 Bug Fixes
 '''''''''
@@ -1209,7 +1238,7 @@ Bug Fixes
    models. Thanks @airanmehr!
 -  Fixed some Pandas <0.24 bugs.
 
-.. _section-45:
+.. _section-46:
 
 0.18.6 - 2019-02-13
 ^^^^^^^^^^^^^^^^^^^
@@ -1219,7 +1248,7 @@ Bug Fixes
    ``rank`` and ``km`` p-values now.
 -  some performance improvements to ``qth_survival_time``.
 
-.. _section-46:
+.. _section-47:
 
 0.18.5 - 2019-02-11
 ^^^^^^^^^^^^^^^^^^^
@@ -1240,7 +1269,7 @@ Bug Fixes
    that can be used to turn off variance calculations since this can
    take a long time for large datasets. Thanks @pzivich!
 
-.. _section-47:
+.. _section-48:
 
 0.18.4 - 2019-02-10
 ^^^^^^^^^^^^^^^^^^^
@@ -1250,7 +1279,7 @@ Bug Fixes
 -  adding left-truncation support to parametric univarite models with
    the ``entry`` kwarg in ``.fit``
 
-.. _section-48:
+.. _section-49:
 
 0.18.3 - 2019-02-07
 ^^^^^^^^^^^^^^^^^^^
@@ -1260,7 +1289,7 @@ Bug Fixes
    warnings are more noticeable.
 -  Improved some warning and error messages.
 
-.. _section-49:
+.. _section-50:
 
 0.18.2 - 2019-02-05
 ^^^^^^^^^^^^^^^^^^^
@@ -1276,7 +1305,7 @@ Bug Fixes
    Moved them all (most) to use ``autograd``.
 -  ``LogNormalFitter`` no longer models ``log_sigma``.
 
-.. _section-50:
+.. _section-51:
 
 0.18.1 - 2019-02-02
 ^^^^^^^^^^^^^^^^^^^
@@ -1287,7 +1316,7 @@ Bug Fixes
 -  use the ``autograd`` lib to help with gradients.
 -  New ``LogLogisticFitter`` univariate fitter available.
 
-.. _section-51:
+.. _section-52:
 
 0.18.0 - 2019-01-31
 ^^^^^^^^^^^^^^^^^^^
@@ -1324,7 +1353,7 @@ Bug Fixes
    ``LinAlgError: Matrix is singular.`` and report back to the user
    advice.
 
-.. _section-52:
+.. _section-53:
 
 0.17.5 - 2019-01-25
 ^^^^^^^^^^^^^^^^^^^
@@ -1332,7 +1361,7 @@ Bug Fixes
 -  more bugs in ``plot_covariate_groups`` fixed when using non-numeric
    strata.
 
-.. _section-53:
+.. _section-54:
 
 0.17.4 -2019-01-25
 ^^^^^^^^^^^^^^^^^^
@@ -1344,7 +1373,7 @@ Bug Fixes
 -  ``groups`` is now called ``values`` in
    ``CoxPHFitter.plot_covariate_groups``
 
-.. _section-54:
+.. _section-55:
 
 0.17.3 - 2019-01-24
 ^^^^^^^^^^^^^^^^^^^
@@ -1352,7 +1381,7 @@ Bug Fixes
 -  Fix in ``compute_residuals`` when using ``schoenfeld`` and the
    minumum duration has only censored subjects.
 
-.. _section-55:
+.. _section-56:
 
 0.17.2 2019-01-22
 ^^^^^^^^^^^^^^^^^
@@ -1363,7 +1392,7 @@ Bug Fixes
    ``for`` loop. The downside is the code is more esoteric now. I’ve
    added comments as necessary though 🤞
 
-.. _section-56:
+.. _section-57:
 
 0.17.1 - 2019-01-20
 ^^^^^^^^^^^^^^^^^^^
@@ -1380,7 +1409,7 @@ Bug Fixes
 -  Fixes a Pandas performance warning in ``CoxTimeVaryingFitter``.
 -  Performances improvements to ``CoxTimeVaryingFitter``.
 
-.. _section-57:
+.. _section-58:
 
 0.17.0 - 2019-01-11
 ^^^^^^^^^^^^^^^^^^^
@@ -1401,7 +1430,7 @@ Bug Fixes
 
 -  some plotting improvemnts to ``plotting.plot_lifetimes``
 
-.. _section-58:
+.. _section-59:
 
 0.16.3 - 2019-01-03
 ^^^^^^^^^^^^^^^^^^^
@@ -1409,7 +1438,7 @@ Bug Fixes
 -  More ``CoxPHFitter`` performance improvements. Up to a 40% reduction
    vs 0.16.2 for some datasets.
 
-.. _section-59:
+.. _section-60:
 
 0.16.2 - 2019-01-02
 ^^^^^^^^^^^^^^^^^^^
@@ -1420,14 +1449,14 @@ Bug Fixes
    has lots of duplicate times. See
    https://github.com/CamDavidsonPilon/lifelines/issues/591
 
-.. _section-60:
+.. _section-61:
 
 0.16.1 - 2019-01-01
 ^^^^^^^^^^^^^^^^^^^
 
 -  Fixed py2 division error in ``concordance`` method.
 
-.. _section-61:
+.. _section-62:
 
 0.16.0 - 2019-01-01
 ^^^^^^^^^^^^^^^^^^^
@@ -1463,7 +1492,7 @@ Bug Fixes
    ``lifelines.utils.to_episodic_format``.
 -  ``CoxTimeVaryingFitter`` now accepts ``strata``.
 
-.. _section-62:
+.. _section-63:
 
 0.15.4
 ^^^^^^
@@ -1471,14 +1500,14 @@ Bug Fixes
 -  bug fix for the Cox model likelihood ratio test when using
    non-trivial weights.
 
-.. _section-63:
+.. _section-64:
 
 0.15.3 - 2018-12-18
 ^^^^^^^^^^^^^^^^^^^
 
 -  Only allow matplotlib less than 3.0.
 
-.. _section-64:
+.. _section-65:
 
 0.15.2 - 2018-11-23
 ^^^^^^^^^^^^^^^^^^^
@@ -1489,7 +1518,7 @@ Bug Fixes
 -  removed ``entry`` from ``ExponentialFitter`` and ``WeibullFitter`` as
    it was doing nothing.
 
-.. _section-65:
+.. _section-66:
 
 0.15.1 - 2018-11-23
 ^^^^^^^^^^^^^^^^^^^
@@ -1498,7 +1527,7 @@ Bug Fixes
 -  Raise NotImplementedError if the ``robust`` flag is used in
    ``CoxTimeVaryingFitter`` - that’s not ready yet.
 
-.. _section-66:
+.. _section-67:
 
 0.15.0 - 2018-11-22
 ^^^^^^^^^^^^^^^^^^^
@@ -1569,7 +1598,7 @@ Bug Fixes
    When Estimating Risks in Pharmacoepidemiology” for a nice overview of
    the model.
 
-.. _section-67:
+.. _section-68:
 
 0.14.6 - 2018-07-02
 ^^^^^^^^^^^^^^^^^^^
@@ -1577,7 +1606,7 @@ Bug Fixes
 -  fix for n > 2 groups in ``multivariate_logrank_test`` (again).
 -  fix bug for when ``event_observed`` column was not boolean.
 
-.. _section-68:
+.. _section-69:
 
 0.14.5 - 2018-06-29
 ^^^^^^^^^^^^^^^^^^^
@@ -1585,7 +1614,7 @@ Bug Fixes
 -  fix for n > 2 groups in ``multivariate_logrank_test``
 -  fix weights in KaplanMeierFitter when using a pandas Series.
 
-.. _section-69:
+.. _section-70:
 
 0.14.4 - 2018-06-14
 ^^^^^^^^^^^^^^^^^^^
@@ -1602,7 +1631,7 @@ Bug Fixes
 -  New ``delay`` parameter in ``add_covariate_to_timeline``
 -  removed ``two_sided_z_test`` from ``statistics``
 
-.. _section-70:
+.. _section-71:
 
 0.14.3 - 2018-05-24
 ^^^^^^^^^^^^^^^^^^^
@@ -1614,7 +1643,7 @@ Bug Fixes
 -  adds a ``column`` argument to ``CoxTimeVaryingFitter`` and
    ``CoxPHFitter`` ``plot`` method to plot only a subset of columns.
 
-.. _section-71:
+.. _section-72:
 
 0.14.2 - 2018-05-18
 ^^^^^^^^^^^^^^^^^^^
@@ -1622,7 +1651,7 @@ Bug Fixes
 -  some quality of life improvements for working with
    ``CoxTimeVaryingFitter`` including new ``predict_`` methods.
 
-.. _section-72:
+.. _section-73:
 
 0.14.1 - 2018-04-01
 ^^^^^^^^^^^^^^^^^^^
@@ -1640,7 +1669,7 @@ Bug Fixes
    faster completion of ``fit`` for large dataframes, and up to 10%
    faster for small dataframes.
 
-.. _section-73:
+.. _section-74:
 
 0.14.0 - 2018-03-03
 ^^^^^^^^^^^^^^^^^^^
@@ -1662,7 +1691,7 @@ Bug Fixes
    of a ``RuntimeWarning``
 -  New checks for complete separation in the dataset for regressions.
 
-.. _section-74:
+.. _section-75:
 
 0.13.0 - 2017-12-22
 ^^^^^^^^^^^^^^^^^^^
@@ -1691,7 +1720,7 @@ Bug Fixes
    group the same subjects together and give that observation a weight
    equal to the count. Altogether, this means a much faster regression.
 
-.. _section-75:
+.. _section-76:
 
 0.12.0
 ^^^^^^
@@ -1708,7 +1737,7 @@ Bug Fixes
 -  Additional functionality to ``utils.survival_table_from_events`` to
    bin the index to make the resulting table more readable.
 
-.. _section-76:
+.. _section-77:
 
 0.11.3
 ^^^^^^
@@ -1720,7 +1749,7 @@ Bug Fixes
    observation or censorship.
 -  More accurate prediction methods parametrics univariate models.
 
-.. _section-77:
+.. _section-78:
 
 0.11.2
 ^^^^^^
@@ -1728,14 +1757,14 @@ Bug Fixes
 -  Changing liscense to valilla MIT.
 -  Speed up ``NelsonAalenFitter.fit`` considerably.
 
-.. _section-78:
+.. _section-79:
 
 0.11.1 - 2017-06-22
 ^^^^^^^^^^^^^^^^^^^
 
 -  Python3 fix for ``CoxPHFitter.plot``.
 
-.. _section-79:
+.. _section-80:
 
 0.11.0 - 2017-06-21
 ^^^^^^^^^^^^^^^^^^^
@@ -1749,14 +1778,14 @@ Bug Fixes
    of a new ``loc`` kwarg. This is to align with Pandas deprecating
    ``ix``
 
-.. _section-80:
+.. _section-81:
 
 0.10.1 - 2017-06-05
 ^^^^^^^^^^^^^^^^^^^
 
 -  fix in internal normalization for ``CoxPHFitter`` predict methods.
 
-.. _section-81:
+.. _section-82:
 
 0.10.0
 ^^^^^^
@@ -1771,7 +1800,7 @@ Bug Fixes
    mimic R’s ``basehaz`` API.
 -  new ``predict_log_partial_hazards`` to ``CoxPHFitter``
 
-.. _section-82:
+.. _section-83:
 
 0.9.4
 ^^^^^
@@ -1794,7 +1823,7 @@ Bug Fixes
 -  performance improvements in ``CoxPHFitter`` - should see at least a
    10% speed improvement in ``fit``.
 
-.. _section-83:
+.. _section-84:
 
 0.9.2
 ^^^^^
@@ -1803,7 +1832,7 @@ Bug Fixes
 -  throw an error if no admissable pairs in the c-index calculation.
    Previously a NaN was returned.
 
-.. _section-84:
+.. _section-85:
 
 0.9.1
 ^^^^^
@@ -1811,7 +1840,7 @@ Bug Fixes
 -  add two summary functions to Weibull and Exponential fitter, solves
    #224
 
-.. _section-85:
+.. _section-86:
 
 0.9.0
 ^^^^^
@@ -1827,7 +1856,7 @@ Bug Fixes
 -  Default predict method in ``k_fold_cross_validation`` is now
    ``predict_expectation``
 
-.. _section-86:
+.. _section-87:
 
 0.8.1 - 2015-08-01
 ^^^^^^^^^^^^^^^^^^
@@ -1844,7 +1873,7 @@ Bug Fixes
    -  scaling of smooth hazards in NelsonAalenFitter was off by a factor
       of 0.5.
 
-.. _section-87:
+.. _section-88:
 
 0.8.0
 ^^^^^
@@ -1863,7 +1892,7 @@ Bug Fixes
    ``lifelines.statistics. power_under_cph``.
 -  fixed a bug when using KaplanMeierFitter for left-censored data.
 
-.. _section-88:
+.. _section-89:
 
 0.7.1
 ^^^^^
@@ -1882,7 +1911,7 @@ Bug Fixes
 -  refactor each fitter into it’s own submodule. For now, the tests are
    still in the same file. This will also *not* break the API.
 
-.. _section-89:
+.. _section-90:
 
 0.7.0 - 2015-03-01
 ^^^^^^^^^^^^^^^^^^
@@ -1901,7 +1930,7 @@ Bug Fixes
    duration remaining until the death event, given survival up until
    time t.
 
-.. _section-90:
+.. _section-91:
 
 0.6.1
 ^^^^^
@@ -1913,7 +1942,7 @@ Bug Fixes
    your work is to sum up the survival function (for expected values or
    something similar), it’s more difficult to make a mistake.
 
-.. _section-91:
+.. _section-92:
 
 0.6.0 - 2015-02-04
 ^^^^^^^^^^^^^^^^^^
@@ -1936,7 +1965,7 @@ Bug Fixes
 -  In ``KaplanMeierFitter``, ``epsilon`` has been renamed to
    ``precision``.
 
-.. _section-92:
+.. _section-93:
 
 0.5.1 - 2014-12-24
 ^^^^^^^^^^^^^^^^^^
@@ -1957,7 +1986,7 @@ Bug Fixes
    ``lifelines.plotting.add_at_risk_counts``.
 -  Fix bug Epanechnikov kernel.
 
-.. _section-93:
+.. _section-94:
 
 0.5.0 - 2014-12-07
 ^^^^^^^^^^^^^^^^^^
@@ -1970,7 +1999,7 @@ Bug Fixes
 -  add test for summary()
 -  Alternate metrics can be used for ``k_fold_cross_validation``.
 
-.. _section-94:
+.. _section-95:
 
 0.4.4 - 2014-11-27
 ^^^^^^^^^^^^^^^^^^
@@ -1982,7 +2011,7 @@ Bug Fixes
 -  Fixes bug in 1-d input not returning in CoxPHFitter
 -  Lots of new tests.
 
-.. _section-95:
+.. _section-96:
 
 0.4.3 - 2014-07-23
 ^^^^^^^^^^^^^^^^^^
@@ -2003,7 +2032,7 @@ Bug Fixes
 -  Adds option ``include_likelihood`` to CoxPHFitter fit method to save
    the final log-likelihood value.
 
-.. _section-96:
+.. _section-97:
 
 0.4.2 - 2014-06-19
 ^^^^^^^^^^^^^^^^^^
@@ -2023,7 +2052,7 @@ Bug Fixes
    from failing so often (this a stop-gap)
 -  pep8 everything
 
-.. _section-97:
+.. _section-98:
 
 0.4.1.1
 ^^^^^^^
@@ -2036,7 +2065,7 @@ Bug Fixes
 -  Adding more robust cross validation scheme based on issue #67.
 -  fixing ``regression_dataset`` in ``datasets``.
 
-.. _section-98:
+.. _section-99:
 
 0.4.1 - 2014-06-11
 ^^^^^^^^^^^^^^^^^^
@@ -2055,7 +2084,7 @@ Bug Fixes
 -  Adding a Changelog.
 -  more sanitizing for the statistical tests =)
 
-.. _section-99:
+.. _section-100:
 
 0.4.0 - 2014-06-08
 ^^^^^^^^^^^^^^^^^^
diff --git a/docs/Survival Regression.rst b/docs/Survival Regression.rst
index 480dd78ac..e419ef367 100644
--- a/docs/Survival Regression.rst	
+++ b/docs/Survival Regression.rst	
@@ -1035,13 +1035,17 @@ The sections `Testing the Proportional Hazard Assumptions`_ and `Assessing Cox m
 .. note:: Work is being done to extend residual methods to all regression models. Stay tuned.
 
 
-Model selection based on predictive power
------------------------------------------------
+Model selection based on predictive power and fit
+---------------------------------------------------
 
 If censoring is present, it's not appropriate to use a loss function like mean-squared-error or
 mean-absolute-loss. This is because the difference between a censored value and the predicted value could be
 due to poor prediction *or* due to censoring. Below we introduce alternative ways to measure prediction performance.
 
+Out-of-sample validation
+****************************
+
+
 In this author's opinion, the best way to measure predictive performance is evaluating the log-likelihood on out-of-sample data. The log-likelihood correctly handles any type of censoring, and is precisely what we are maximizing in the model training. The in-sample log-likelihood is available under ``log_likelihood_`` of any regression model. For out-of-sample data, the  :meth:`~lifelines.fitters.cox_ph_fitter.CoxPHFitter.score` method (available on all regression models) can be used. This returns the *average evaluation of the out-of-sample log-likelihood*. We want to maximize this.
 
 .. code:: python
@@ -1060,6 +1064,27 @@ In this author's opinion, the best way to measure predictive performance is eval
     print(cph_l1.score(test_rossi)) # better model
 
 
+Within-sample validation
+****************************
+
+For within-sample validation, the AIC is a great metric for comparing models as it relies on the log-likelihood. It's available under ``AIC_`` for parametric models, and ``AIC_partial_`` for Cox models (because the Cox model maximizes a *partial* log-likelihood, it can't be reliably compared to parametric model's AIC.)
+
+
+.. code:: python
+
+    from lifelines import CoxPHFitter
+    from lifelines.datasets import load_rossi
+
+    rossi = load_rossi()
+
+    cph_l2 = CoxPHFitter(penalizer=0.1, l1_ratio=0.).fit(rossi, 'week', 'arrest')
+    cph_l1 = CoxPHFitter(penalizer=0.1, l1_ratio=1.).fit(rossi, 'week', 'arrest')
+
+    print(cph_l2.AIC_partial_) # lower is better
+    print(cph_l1.AIC_partial_)
+
+
+
 Another censoring-sensitive measure is the concordance-index, also known as the c-index. This measure
 evaluates the accuracy of the *ranking* of predicted time. It is in fact a generalization
 of AUC, another common loss function, and is interpreted similarly:
@@ -1097,6 +1122,9 @@ Fitted survival models typically have a concordance index between 0.55 and 0.75
 .. note:: Remember, the concordance score evaluates the relative rankings of subject's event times. Thus, it is scale and shift invariant (i.e. you can multiple by a positive constant, or add a constant, and the rankings won't change). A model maximized for concordance-index does not necessarily give good predicted *times*, but will give good predicted *rankings*.
 
 
+Cross validation
+****************************
+
 *lifelines* has an implementation of k-fold cross validation under :func:`lifelines.utils.k_fold_cross_validation`. This function accepts an instance of a regression fitter (either :class:`~lifelines.fitters.coxph_fitter.CoxPHFitter` of :class:`~lifelines.fitters.aalen_additive_fitter.AalenAdditiveFitter`), a dataset, plus ``k`` (the number of folds to perform, default 5). On each fold, it splits the data
 into a training set and a testing set fits itself on the training set and evaluates itself on the testing set (using the concordance measure by default).
 
diff --git a/lifelines/fitters/__init__.py b/lifelines/fitters/__init__.py
index 343c29bea..87db76201 100644
--- a/lifelines/fitters/__init__.py
+++ b/lifelines/fitters/__init__.py
@@ -293,6 +293,10 @@ def __init__(self, *args, **kwargs):
         if "alpha" in self._fitted_parameter_names:
             raise NameError("'alpha' in _fitted_parameter_names is a lifelines reserved word. Try 'alpha_' instead.")
 
+    @property
+    def AIC_(self) -> float:
+        return -2 * self.log_likelihood_ + 2 * self.params_.shape[0]
+
     def _check_cumulative_hazard_is_monotone_and_positive(self, durations, values):
         class_name = self._class_name
 
diff --git a/lifelines/fitters/kaplan_meier_fitter.py b/lifelines/fitters/kaplan_meier_fitter.py
index 4b2e0d3b5..7fff8c176 100644
--- a/lifelines/fitters/kaplan_meier_fitter.py
+++ b/lifelines/fitters/kaplan_meier_fitter.py
@@ -128,16 +128,58 @@ def fit_interval_censoring(
         show_progress=False,
         entry=None,
         weights=None,
-        tol=1e-5,
+        tol=1e-7,
     ) -> "KaplanMeierFitter":
+        """
+        Fit the model to a interval-censored dataset using non-parametric MLE. This estimator is
+        also called the Turball Estimator.
+
+        Note
+        ------
+        This is new and experimental, and many feature are missing.
+
+        Parameters
+        ----------
+          lower_bound: an array, list, pd.DataFrame or pd.Series
+            length n -- lower bound of observations
+          upper_bound: an array, list, pd.DataFrame or pd.Series
+            length n -- upper bound of observations
+          event_observed: an array, list, pd.DataFrame, or pd.Series, optional
+             True if the the death was observed, False if the event was lost (right-censored). This can be computed from
+             the lower_bound and upper_bound, and can be left blank.
+          timeline: an array, list, pd.DataFrame, or pd.Series, optional
+            return the best estimate at the values in timelines (positively increasing)
+          entry: an array, list, pd.DataFrame, or pd.Series, optional
+             relative time when a subject entered the study. This is useful for left-truncated (not left-censored) observations. If None, all members of the population
+             entered study when they were "born".
+          label: string, optional
+            a string to name the column of the estimate.
+          alpha: float, optional
+            the alpha value in the confidence intervals. Overrides the initializing alpha for this call to fit only.
+          ci_labels: tuple, optional
+                add custom column names to the generated confidence intervals as a length-2 list: [<lower-bound name>, <upper-bound name>]. Default: <label>_lower_<1-alpha/2>
+          weights: an array, list, pd.DataFrame, or pd.Series, optional
+              if providing a weighted dataset. For example, instead
+              of providing every subject as a single element of `durations` and `event_observed`, one could
+              weigh subject differently.
+
+        Returns
+        -------
+        self: KaplanMeierFitter
+          self with new properties like ``survival_function_``, ``plot()``, ``median_survival_time_``
+        """
+        warnings.warn("This is new and experimental, many feature are missing and accuracy is not reliable", UserWarning)
 
         if entry is not None or weights is not None:
             raise NotImplementedError("entry / weights is not supported yet")
+        self.weights = np.ones_like(upper_bound)
 
         self.upper_bound = np.atleast_1d(pass_for_numeric_dtypes_or_raise_array(upper_bound))
         self.lower_bound = np.atleast_1d(pass_for_numeric_dtypes_or_raise_array(lower_bound))
         check_nans_or_infs(self.lower_bound)
 
+        self.event_observed = self.lower_bound == self.upper_bound
+
         self.timeline = coalesce(timeline, np.unique(np.concatenate((self.upper_bound, self.lower_bound))))
 
         if (self.upper_bound < self.lower_bound).any():
@@ -151,10 +193,12 @@ def fit_interval_censoring(
                 "For all rows, lower_bound == upper_bound if and only if event observed = 1 (uncensored). Likewise, lower_bound < upper_bound if and only if event observed = 0 (censored)"
             )
 
-        self._label = coalesce(label, self._label, "KM_estimate")
+        self._label = coalesce(label, self._label, "NPMLE_estimate")
 
         probs, t_intervals = npmle(self.lower_bound, self.upper_bound, verbose=show_progress)
-        self.survival_function_ = reconstruct_survival_function(probs, t_intervals, self.timeline, label=self._label)
+        self.survival_function_ = reconstruct_survival_function(probs, t_intervals, self.timeline, label=self._label).loc[
+            self.timeline
+        ]
         self.cumulative_density_ = 1 - self.survival_function_
 
         self._median = median_survival_times(self.survival_function_)
@@ -405,7 +449,12 @@ def plot_cumulative_density(self, **kwargs):
         ax:
             a pyplot axis object
         """
-        return _plot_estimate(self, estimate="cumulative_density_", **kwargs)
+        if not CensoringType.is_interval_censoring(self):
+            return _plot_estimate(self, estimate="cumulative_density_", **kwargs)
+        else:
+            # hack for now.
+            color = coalesce(kwargs.get("c"), kwargs.get("color"), "k")
+            self.cumulative_density_.plot(drawstyle="steps", color=color, **kwargs)
 
     def _bounds(self, cumulative_sq_, alpha, ci_labels):
         # This method calculates confidence intervals using the exponential Greenwood formula.
diff --git a/lifelines/fitters/npmle.py b/lifelines/fitters/npmle.py
index 07187657c..3d28c78d1 100644
--- a/lifelines/fitters/npmle.py
+++ b/lifelines/fitters/npmle.py
@@ -1,4 +1,12 @@
 # -*- coding: utf-8 -*-
+"""
+This code isn't to be called directly, but is the core logic of the KaplanMeierFitter.fit_interval_censoring
+
+References
+https://upcommons.upc.edu/bitstream/handle/2117/93831/01Rop01de01.pdf
+
+
+"""
 from collections import defaultdict, namedtuple
 import numpy as np
 from numpy.linalg import norm
@@ -20,9 +28,12 @@ def E_step_M_step(observation_intervals, p_old, turnball_interval_lookup):
         # find all turnball intervals, t, that are contained in (ol, or). Call this set T
         # the denominator is sum of p_old[T] probabilities
         # the numerator is p_old[t]
+
+        # TODO: I think I can remove the bottom for loop:
+        # > ix = turnball_interval_lookup[observation_interval]
+        # > p_temp[ix] = p_old[ix]
         for ix_t in turnball_interval_lookup[observation_interval]:
-            t_p = p_old[ix_t]
-            p_temp[ix_t] = t_p
+            p_temp[ix_t] = p_old[ix_t]
 
         p_new = p_new + p_temp / p_temp.sum()
 
@@ -108,6 +119,17 @@ def npmle(left, right, tol=1e-5, verbose=False):
     turnball_intervals = list(create_turnball_intervals(left, right))
     observation_intervals = create_observation_intervals(left, right)
     turnball_lookup = create_turnball_lookup(turnball_intervals, sorted(set(observation_intervals)))
+    print(turnball_lookup)
+
+    """ is correct
+    turnball_lookup = {
+        interval(6.0, 6.0): {0},
+        interval(7.0, 7.0): {1},
+        interval(7.0, np.inf): {2},  # {1, 2}
+        interval(8., 8.): {2},
+    }
+    """
+    print(turnball_lookup)
 
     T = len(turnball_intervals)
 
@@ -126,10 +148,21 @@ def npmle(left, right, tol=1e-5, verbose=False):
 
 
 def reconstruct_survival_function(probabilities, turnball_intervals, timeline, label="NPMLE"):
-    index = [0.0]
-    values = [1.0]
+    """
+    TIHI
+
+    """
+    index = []
+    values = []
+
+    for i, (p, interval) in enumerate(zip(probabilities, turnball_intervals)):
+        if i == 0:
+            index.append(interval.left)
+            index.append(interval.right)
+            values.append(1.0)
+            values.append(1 - p)
+            continue
 
-    for p, interval in zip(probabilities, turnball_intervals):
         if interval.left != index[-1]:
             index.append(interval.left)
             values.append(values[-1])
@@ -140,12 +173,12 @@ def reconstruct_survival_function(probabilities, turnball_intervals, timeline, l
             index.append(interval.right)
             values.append(values[-1] - p)
 
-    full_dataframe = pd.DataFrame(index=timeline, columns=[label + "_upper"])
+    full_dataframe = pd.DataFrame(index=timeline, columns=[label + "_lower"])
 
-    turnball_dataframe = pd.DataFrame(values, index=index, columns=[label + "_upper"])
+    turnball_dataframe = pd.DataFrame(values, index=index, columns=[label + "_lower"])
 
-    dataframe = full_dataframe.combine_first(turnball_dataframe).ffill()
-    dataframe[label + "_lower"] = dataframe[label + "_upper"].shift(1).fillna(1)
+    dataframe = full_dataframe.combine_first(turnball_dataframe).ffill().fillna(1)
+    dataframe[label + "_upper"] = dataframe[label + "_lower"].shift(1).fillna(1)
     return dataframe
 
 
@@ -171,3 +204,18 @@ def npmle_compute_confidence_intervals(left, right, mle_, alpha=0.05, samples=10
         2 * mle_.squeeze() - pd.Series(np.percentile(bootstrapped_samples, (alpha / 2) * 100, axis=1), index=all_times),
         2 * mle_.squeeze() - pd.Series(np.percentile(bootstrapped_samples, (1 - alpha / 2) * 100, axis=1), index=all_times),
     )
+
+
+"""
+Test cases
+"""
+
+
+ti = create_turnball_intervals(*zip(*[(0, 1), (4, 6), (2, 6), (0, 3), (2, 4), (5, 7)]))
+
+assert ti == [interval(0, 1), interval(2, 3), interval(5, 6)]
+
+
+ti = create_turnball_intervals(*zip(*[(4, 7), (3, 5), (0, 2), (1, 4), (6, 9), (8, 10)]))
+
+assert ti == [interval(1, 2), interval(3, 4), interval(4, 5), interval(6, 7), interval(8, 9)]
diff --git a/lifelines/plotting.py b/lifelines/plotting.py
index 720d9dbdc..a15a1a468 100644
--- a/lifelines/plotting.py
+++ b/lifelines/plotting.py
@@ -95,7 +95,9 @@ def cdf_plot(model, timeline=None, ax=None, **plot_kwargs):
             model.durations, model.event_observed, label=COL_EMP, timeline=timeline, weights=model.weights, entry=model.entry
         )
     elif CensoringType.is_interval_censoring(model):
-        raise NotImplementedError("lifelines does not have a non-parametric interval model yet.")
+        empirical_kmf = KaplanMeierFitter().fit_interval_censoring(
+            model.lower_bound, model.upper_bound, label=COL_EMP, timeline=timeline, weights=model.weights, entry=model.entry
+        )
 
     empirical_kmf.plot_cumulative_density(ax=ax, **plot_kwargs)
 
@@ -258,16 +260,22 @@ def qq_plot(model, ax=None, **plot_kwargs):
         kmf = KaplanMeierFitter().fit_left_censoring(
             model.durations, model.event_observed, label=COL_EMP, weights=model.weights, entry=model.entry
         )
+        sf, cdf = kmf.survival_function_[COL_EMP], kmf.cumulative_density_[COL_EMP]
     elif CensoringType.is_right_censoring(model):
         kmf = KaplanMeierFitter().fit_right_censoring(
             model.durations, model.event_observed, label=COL_EMP, weights=model.weights, entry=model.entry
         )
+        sf, cdf = kmf.survival_function_[COL_EMP], kmf.cumulative_density_[COL_EMP]
+
     elif CensoringType.is_interval_censoring(model):
-        raise NotImplementedError("lifelines does not have a non-parametric interval model yet.")
+        kmf = KaplanMeierFitter().fit_interval_censoring(
+            model.lower_bound, model.upper_bound, label=COL_EMP, weights=model.weights, entry=model.entry
+        )
+        sf, cdf = kmf.survival_function_[COL_EMP + "_lower"], kmf.cumulative_density_[COL_EMP + "_lower"]
 
-    q = np.unique(kmf.cumulative_density_.values[:, 0])
+    q = np.unique(cdf.values)
 
-    quantiles = qth_survival_times(1 - q, kmf.survival_function_)
+    quantiles = qth_survival_times(1 - q, sf)
     quantiles[COL_THEO] = dist_object.ppf(q)
     quantiles = quantiles.replace([-np.inf, 0, np.inf], np.nan).dropna()
 
diff --git a/lifelines/statistics.py b/lifelines/statistics.py
index 1bc7c7043..6ee75e5f3 100644
--- a/lifelines/statistics.py
+++ b/lifelines/statistics.py
@@ -29,6 +29,7 @@
     "proportional_hazard_test",
     "power_under_cph",
     "sample_size_necessary_under_cph",
+    "somers_d",
 ]
 
 
diff --git a/lifelines/tests/test_plotting.py b/lifelines/tests/test_plotting.py
index 7f5b16cef..b099b842d 100644
--- a/lifelines/tests/test_plotting.py
+++ b/lifelines/tests/test_plotting.py
@@ -32,6 +32,7 @@
     load_rossi,
     load_multicenter_aids_cohort_study,
     load_nh4,
+    load_diabetes,
 )
 from lifelines.generate_datasets import cumulative_integral
 
@@ -215,6 +216,14 @@ def test_aalen_additive_plot(self, block):
         self.plt.show(block=block)
         return
 
+    def test_kmf_with_interval_censoring_plotting(self, block):
+        kmf = KaplanMeierFitter()
+        left, right = load_diabetes()[["left", "right"]]
+        kmf.fit(left, right)
+        kmf.plot(color="r")
+        self.plt.show(block=block)
+        return
+
     def test_aalen_additive_smoothed_plot(self, block):
         # this is a visual test of the fitting the cumulative
         # hazards.
diff --git a/lifelines/utils/concordance.py b/lifelines/utils/concordance.py
index 231442049..0d113e981 100644
--- a/lifelines/utils/concordance.py
+++ b/lifelines/utils/concordance.py
@@ -4,6 +4,35 @@
 from lifelines.utils.btree import _BTree
 
 
+def somers_d(event_times, x, event_observed=None) -> float:
+    """
+    A measure of rank association between [-1, 1] between a censored variable, event_times,
+    and another (uncensored) variable, x. -1 is strong anti-correlation, 1 is strong correlation.
+
+
+    event_times: iterable
+         a length-n iterable of observed survival times.
+    x: iterable
+        a length-n iterable to compare against
+    event_observed: iterable, optional
+        a length-n iterable censoring flags, 1 if observed, 0 if not. Default None assumes all observed.
+
+
+    Examples
+    --------
+    .. code:: python
+        from lifelines.datasets import load_rossi
+        from lifelines.utils
+
+        T, E = df['week'], df['arrest']
+        x = df['age']
+        somers_d(T, x, E)
+
+
+    """
+    return 2 * concordance_index(event_times, x, event_observed) - 1
+
+
 def concordance_index(event_times, predicted_scores, event_observed=None) -> float:
     """
     Calculates the concordance index (C-index) between two series
@@ -29,7 +58,7 @@ def concordance_index(event_times, predicted_scores, event_observed=None) -> flo
     predicted_scores: iterable
         a length-n iterable of predicted scores - these could be survival times, or hazards, etc. See https://stats.stackexchange.com/questions/352183/use-median-survival-time-to-calculate-cph-c-statistic/352435#352435
     event_observed: iterable, optional
-        a length-n iterable censorship flags, 1 if observed, 0 if not. Default None assumes all observed.
+        a length-n iterable censoring flags, 1 if observed, 0 if not. Default None assumes all observed.
 
     Returns
     -------
@@ -51,9 +80,7 @@ def concordance_index(event_times, predicted_scores, event_observed=None) -> flo
         concordance_index(df['T'], -cph.predict_partial_hazard(df), df['E'])
 
     """
-    event_times, predicted_scores, event_observed = _preprocess_scoring_data(
-        event_times, predicted_scores, event_observed
-    )
+    event_times, predicted_scores, event_observed = _preprocess_scoring_data(event_times, predicted_scores, event_observed)
     num_correct, num_tied, num_pairs = _concordance_summary_statistics(event_times, predicted_scores, event_observed)
 
     return _concordance_ratio(num_correct, num_tied, num_pairs)
@@ -65,9 +92,7 @@ def _concordance_ratio(num_correct: int, num_tied: int, num_pairs: int) -> float
     return (num_correct + num_tied / 2) / num_pairs
 
 
-def _concordance_summary_statistics(
-    event_times, predicted_event_times, event_observed
-):  # pylint: disable=too-many-locals
+def _concordance_summary_statistics(event_times, predicted_event_times, event_observed):  # pylint: disable=too-many-locals
     """Find the concordance index in n * log(n) time.
 
     Assumes the data has been verified by lifelines.utils.concordance_index first.
@@ -234,9 +259,7 @@ def naive_concordance_index(event_times, predicted_event_times, event_observed=N
     event_times, predicted_event_times, event_observed = _preprocess_scoring_data(
         event_times, predicted_event_times, event_observed
     )
-    return _concordance_ratio(
-        *_naive_concordance_summary_statistics(event_times, predicted_event_times, event_observed)
-    )
+    return _concordance_ratio(*_naive_concordance_summary_statistics(event_times, predicted_event_times, event_observed))
 
 
 def _preprocess_scoring_data(event_times, predicted_scores, event_observed):
diff --git a/lifelines/version.py b/lifelines/version.py
index 0663adbfd..168a9842d 100644
--- a/lifelines/version.py
+++ b/lifelines/version.py
@@ -1,4 +1,4 @@
 # -*- coding: utf-8 -*-
 from __future__ import unicode_literals
 
-__version__ = "0.24.6"
+__version__ = "0.24.7"

From 8578b564f5970ceb2eeb489946edd705faff906b Mon Sep 17 00:00:00 2001
From: CamDavidsonPilon <cam.davidson.pilon@gmail.com>
Date: Sat, 16 May 2020 10:30:20 -0400
Subject: [PATCH 09/10] fix tests

---
 lifelines/fitters/__init__.py      |  8 ++++----
 lifelines/tests/test_estimation.py | 24 ++++++++++++++++++++----
 2 files changed, 24 insertions(+), 8 deletions(-)

diff --git a/lifelines/fitters/__init__.py b/lifelines/fitters/__init__.py
index 87db76201..d05f8ba92 100644
--- a/lifelines/fitters/__init__.py
+++ b/lifelines/fitters/__init__.py
@@ -295,7 +295,7 @@ def __init__(self, *args, **kwargs):
 
     @property
     def AIC_(self) -> float:
-        return -2 * self.log_likelihood_ + 2 * self.params_.shape[0]
+        return -2 * self.log_likelihood_ + 2 * self._fitted_parameters_.shape[0]
 
     def _check_cumulative_hazard_is_monotone_and_positive(self, durations, values):
         class_name = self._class_name
@@ -1940,7 +1940,7 @@ def _compute_variance_matrix(self) -> np.array:
                 Some ways to possible ways fix this:
 
                 0. Are there any lifelines warnings outputted during the `fit`?
-                1. Inspect your DataFrame: does everything look as expected?
+                1. Inspect your DataFrame: does everything look as expected? Do you need to add/drop a constant (intercept) column?
                 2. Is there high-collinearity in the dataset? Try using the variance inflation factor (VIF) to find redundant variables.
                 3. Trying adding a small penalizer (or changing it, if already present). Example: `%s(penalizer=0.01).fit(...)`.
                 4. Are there any extreme outliers? Try modeling them or dropping them to see if it helps convergence.
@@ -2130,8 +2130,8 @@ def print_summary(self, decimals=2, style=None, **kwargs):
         sr = self.log_likelihood_ratio_test()
         footers = []
 
-        if CensoringType.is_right_censoring(self):
-            footers.apoend(("Concordance", "{:.{prec}f}".format(self.concordance_index_, prec=decimals)))
+        if utils.CensoringType.is_right_censoring(self):
+            footers.append(("Concordance", "{:.{prec}f}".format(self.concordance_index_, prec=decimals)))
 
         footers.extend(
             [
diff --git a/lifelines/tests/test_estimation.py b/lifelines/tests/test_estimation.py
index ad15c1f63..2c63a4386 100644
--- a/lifelines/tests/test_estimation.py
+++ b/lifelines/tests/test_estimation.py
@@ -228,6 +228,14 @@ def _cumulative_hazard(self, params, times):
         assert coef - std > lower
         assert coef + std < upper
 
+    def test_AIC_on_models(self, known_parametric_univariate_fitters):
+        T = np.random.exponential(1, size=1000)
+
+        for fitter in known_parametric_univariate_fitters:
+            f = fitter().fit(T)
+            assert f.AIC_ > 0
+            npt.assert_allclose(f.AIC_, -2 * f.log_likelihood_ + 2 * f.summary.shape[0])
+
     def test_models_can_handle_really_large_duration_values(self, known_parametric_univariate_fitters):
         T1 = np.random.exponential(1e12, size=1000)
         T2 = np.random.exponential(1e12, size=1000)
@@ -1613,18 +1621,22 @@ def rossi(self):
         rossi["_int"] = 1.0
         return rossi
 
+    def test_AIC_on_models(self, rossi):
+        model = WeibullAFTFitter(fit_intercept=False).fit(rossi, "week", "arrest")
+        npt.assert_allclose(model.AIC_, -2 * model.log_likelihood_ + 2 * model.summary.shape[0])
+
     def test_penalizer_can_be_an_array(self, rossi):
 
-        wf_array = WeibullAFTFitter(penalizer=0.01 * np.ones(7)).fit(rossi, "week", "arrest")
-        wf_float = WeibullAFTFitter(penalizer=0.01).fit(rossi, "week", "arrest")
+        wf_array = WeibullAFTFitter(penalizer=0.01 * np.ones(7), fit_intercept=False).fit(rossi, "week", "arrest")
+        wf_float = WeibullAFTFitter(penalizer=0.01, fit_intercept=False).fit(rossi, "week", "arrest")
 
         assert_frame_equal(wf_array.summary, wf_float.summary)
 
     def test_penalizer_can_be_an_array_and_check_it_behaves_as_expected(self, rossi):
 
         penalty = np.array([0, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01])
-        wf_array = WeibullAFTFitter(penalizer=penalty).fit(rossi, "week", "arrest")
-        wf_float = WeibullAFTFitter(penalizer=0.01).fit(rossi, "week", "arrest")
+        wf_array = WeibullAFTFitter(penalizer=penalty, fit_intercept=False).fit(rossi, "week", "arrest")
+        wf_float = WeibullAFTFitter(penalizer=0.01, fit_intercept=False).fit(rossi, "week", "arrest")
 
         assert abs(wf_array.summary.loc[("lambda_", "fin"), "coef"]) > abs(wf_float.summary.loc[("lambda_", "fin"), "coef"])
 
@@ -2696,6 +2708,10 @@ def test_penalizer_can_be_an_array_and_check_it_behaves_as_expected(self, rossi)
 
         assert abs(cph_array.summary.loc["fin", "coef"]) > abs(cph_float.summary.loc["fin", "coef"])
 
+    def test_AIC_partial_(self, cph, rossi):
+        cph.fit(rossi, "week", "arrest")
+        npt.assert_allclose(cph.AIC_partial_, -2 * cph.log_likelihood_ + 2 * cph.summary.shape[0])
+
     def test_compute_followup_hazard_ratios(self, cph, cph_spline, rossi):
         cph.fit(rossi, "week", "arrest")
         cph.compute_followup_hazard_ratios(rossi, [15, 25, 35, 45])

From 94b2193958184e0870af7d190b601ff06e88fad2 Mon Sep 17 00:00:00 2001
From: CamDavidsonPilon <cam.davidson.pilon@gmail.com>
Date: Sat, 16 May 2020 13:10:43 -0400
Subject: [PATCH 10/10] fix test

---
 lifelines/fitters/__init__.py                 |  3 +-
 lifelines/fitters/generalized_gamma_fitter.py |  2 +-
 lifelines/tests/test_estimation.py            | 80 ++++++++++++-------
 lifelines/tests/test_plotting.py              |  4 +-
 4 files changed, 57 insertions(+), 32 deletions(-)

diff --git a/lifelines/fitters/__init__.py b/lifelines/fitters/__init__.py
index d05f8ba92..8457fb3be 100644
--- a/lifelines/fitters/__init__.py
+++ b/lifelines/fitters/__init__.py
@@ -494,7 +494,8 @@ def _fit_model(self, Ts, E, entry, weights, show_progress=True):
 
             minimizing_results, previous_results, minimizing_ll = None, None, np.inf
             for method, option in zip(
-                ["Nelder-Mead", self._scipy_fit_method], [{"maxiter": 25}, {**{"disp": show_progress}, **self._scipy_fit_options}]
+                ["Nelder-Mead", self._scipy_fit_method],
+                [{"maxiter": 100}, {**{"disp": show_progress}, **self._scipy_fit_options}],
             ):
 
                 initial_value = self._initial_values if previous_results is None else utils._to_1d_array(previous_results.x)
diff --git a/lifelines/fitters/generalized_gamma_fitter.py b/lifelines/fitters/generalized_gamma_fitter.py
index c41fbf06e..9ec638f3e 100644
--- a/lifelines/fitters/generalized_gamma_fitter.py
+++ b/lifelines/fitters/generalized_gamma_fitter.py
@@ -115,7 +115,7 @@ def _create_initial_point(self, Ts, E, *args):
             log_data = log(Ts[1])
         elif CensoringType.is_interval_censoring(self):
             # this fails if Ts[1] == Ts[0], so we add a some fudge factors.
-            log_data = log(Ts[1] - Ts[0] + 0.01)
+            log_data = log(Ts[1] - Ts[0] + 0.1)
         return np.array([log_data.mean(), log(log_data.std() + 0.01), 0.1])
 
     def _cumulative_hazard(self, params, times):
diff --git a/lifelines/tests/test_estimation.py b/lifelines/tests/test_estimation.py
index 2c63a4386..372957383 100644
--- a/lifelines/tests/test_estimation.py
+++ b/lifelines/tests/test_estimation.py
@@ -260,6 +260,9 @@ def test_models_can_handle_really_small_duration_values_for_left_censorship(self
         T = np.maximum(T1, T2)
 
         for fitter in known_parametric_univariate_fitters:
+            if isinstance(fitter(), PiecewiseExponentialFitterTesting):
+                # not a good model since the "break" is at 5.
+                continue
             fitter().fit_left_censoring(T, E)
 
     def test_parametric_univariate_fitters_can_print_summary(
@@ -383,7 +386,8 @@ def test_default_alpha_is_005(self, univariate_fitters):
             assert f().alpha == 0.05
 
     def test_univariate_fitters_accept_late_entries(self, positive_sample_lifetimes, univariate_fitters):
-        entries = 0.1 * positive_sample_lifetimes[0]
+        positive_sample_lifetimes = positive_sample_lifetimes
+        entries = positive_sample_lifetimes[0] - 3
         for fitter in univariate_fitters:
             f = fitter().fit(positive_sample_lifetimes[0], entry=entries)
             assert f.entry is not None
@@ -3124,26 +3128,38 @@ def test_print_summary(self, rossi, cph):
                     repr(cph)
                     + "\n"
                     + """
-      duration col = week
-         event col = arrest
-number of subjects = 432
-  number of events = 114
-    log-likelihood = -658.748
-  time fit was run = 2018-10-23 02:40:45 UTC
+<lifelines.CoxPHFitter: fitted with 432 total observations, 318 right-censored observations>
+             duration col = 'week'
+                event col = 'arrest'
+      baseline estimation = breslow
+   number of observations = 432
+number of events observed = 114
+   partial log-likelihood = -658.75
+         time fit was run = 2018-10-23 02:40:45 UTC
 
 ---
-        coef  exp(coef)  se(coef)       z      p  coef lower 95%  coef upper 95%
-fin  -0.3794     0.6843    0.1914 -1.9826 0.0474     -0.7545     -0.0043
-age  -0.0574     0.9442    0.0220 -2.6109 0.0090     -0.1006     -0.0143
-race  0.3139     1.3688    0.3080  1.0192 0.3081     -0.2898      0.9176
-wexp -0.1498     0.8609    0.2122 -0.7058 0.4803     -0.5657      0.2662
-mar  -0.4337     0.6481    0.3819 -1.1358 0.2561     -1.1821      0.3147
-paro -0.0849     0.9186    0.1958 -0.4336 0.6646     -0.4685      0.2988
-prio  0.0915     1.0958    0.0286  3.1939 0.0014      0.0353      0.1476
+       coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%
+fin   -0.38       0.68       0.19            -0.75            -0.00                 0.47                 1.00
+age   -0.06       0.94       0.02            -0.10            -0.01                 0.90                 0.99
+race   0.31       1.37       0.31            -0.29             0.92                 0.75                 2.50
+wexp  -0.15       0.86       0.21            -0.57             0.27                 0.57                 1.30
+mar   -0.43       0.65       0.38            -1.18             0.31                 0.31                 1.37
+paro  -0.08       0.92       0.20            -0.47             0.30                 0.63                 1.35
+prio   0.09       1.10       0.03             0.04             0.15                 1.04                 1.16
+
+         z      p   -log2(p)
+fin  -1.98   0.05       4.40
+age  -2.61   0.01       6.79
+race  1.02   0.31       1.70
+wexp -0.71   0.48       1.06
+mar  -1.14   0.26       1.97
+paro -0.43   0.66       0.59
+prio  3.19 <0.005       9.48
 ---
-
-Concordance = 0.640
-Log-likelihood ratio test = 33.27 on 7 df, -log2(p)=15.37
+Concordance = 0.64
+Partial AIC = 1331.50
+log-likelihood ratio test = 33.27 on 7 df
+-log2(p) of ll-ratio test = 15.37
 """
                 )
                 .strip()
@@ -4949,22 +4965,30 @@ def test_print_summary(self, ctv, heart):
                     repr(ctv)
                     + "\n"
                     + """
-         event col = event
+<lifelines.CoxTimeVaryingFitter: fitted with 172 periods, 103 subjects, 75 events>
+         event col = 'event'
 number of subjects = 103
  number of periods = 172
   number of events = 75
-    log-likelihood = -290.566
+partial log-likelihood = -290.57
   time fit was run = 2018-10-23 02:41:45 UTC
 
 ---
-              coef  exp(coef)  se(coef)       z      p  coef lower 95%  coef upper 95%
-age         0.0272     1.0275    0.0137  1.9809 0.0476      0.0003      0.0540
-year       -0.1463     0.8639    0.0705 -2.0768 0.0378     -0.2845     -0.0082
-surgery    -0.6372     0.5288    0.3672 -1.7352 0.0827     -1.3570      0.0825
-transplant -0.0103     0.9898    0.3138 -0.0327 0.9739     -0.6252      0.6047
+             coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%
+age          0.03       1.03       0.01             0.00             0.05                 1.00                 1.06
+year        -0.15       0.86       0.07            -0.28            -0.01                 0.75                 0.99
+surgery     -0.64       0.53       0.37            -1.36             0.08                 0.26                 1.09
+transplant  -0.01       0.99       0.31            -0.63             0.60                 0.54                 1.83
+
+               z    p   -log2(p)
+age         1.98 0.05       4.39
+year       -2.08 0.04       4.72
+surgery    -1.74 0.08       3.60
+transplant -0.03 0.97       0.04
 ---
-
-Likelihood ratio test = 15.11 on 4 df, -log2(p)=7.80
+Partial AIC = 589.13
+log-likelihood ratio test = 15.11 on 4 df
+-log2(p) of ll-ratio test = 7.80
 """
                 )
                 .strip()
@@ -5169,4 +5193,4 @@ def _create_initial_point(self, Ts, E, entry, weights):
         wmc.fit(T, E)
         mcfitter.fit(T, E)
 
-        assert abs(wmc.c_ - mcfitter.cured_fraction_) < 0.001
+        assert_frame_equal(wmc.summary.reset_index(drop=True), mcfitter.summary.reset_index(drop=True), check_less_precise=0)
diff --git a/lifelines/tests/test_plotting.py b/lifelines/tests/test_plotting.py
index b099b842d..79c3c671e 100644
--- a/lifelines/tests/test_plotting.py
+++ b/lifelines/tests/test_plotting.py
@@ -218,8 +218,8 @@ def test_aalen_additive_plot(self, block):
 
     def test_kmf_with_interval_censoring_plotting(self, block):
         kmf = KaplanMeierFitter()
-        left, right = load_diabetes()[["left", "right"]]
-        kmf.fit(left, right)
+        left, right = load_diabetes()["left"], load_diabetes()["right"]
+        kmf.fit_interval_censoring(left, right)
         kmf.plot(color="r")
         self.plt.show(block=block)
         return