Python estimation.CoxPHFitter类代码示例

    def test_predict_log_hazard_relative_to_mean_with_normalization(self, rossi):
        cox = CoxPHFitter(normalize=True)
        cox.fit(rossi, 'week', 'arrest')

        # they are equal because the data is normalized, so the mean of the covarites is all 0,
        # thus exp(beta * 0) == 1, so exp(beta * X)/exp(beta * 0) = exp(beta * X)
        assert_frame_equal(cox.predict_log_hazard_relative_to_mean(rossi), np.log(cox.predict_partial_hazard(rossi)))

 def test_coxph_plotting_normalized(self, block):
     df = load_regression_dataset()
     cp = CoxPHFitter()
     cp.fit(df, "T", "E")
     cp.plot(True)
     self.plt.title('test_coxph_plotting')
     self.plt.show(block=block)

    def test_print_summary(self, rossi):

        import sys
        saved_stdout = sys.stdout
        try:
            out = StringIO()
            sys.stdout = out

            cp = CoxPHFitter()
            cp.fit(rossi, duration_col='week', event_col='arrest')
            cp.print_summary()
            output = out.getvalue().strip().split()
            expected = """n=432, number of events=114

           coef  exp(coef)  se(coef)          z         p  lower 0.95  upper 0.95
fin  -1.897e-01  8.272e-01 9.579e-02 -1.981e+00 4.763e-02  -3.775e-01  -1.938e-03   *
age  -3.500e-01  7.047e-01 1.344e-01 -2.604e+00 9.210e-03  -6.134e-01  -8.651e-02  **
race  1.032e-01  1.109e+00 1.012e-01  1.020e+00 3.078e-01  -9.516e-02   3.015e-01
wexp -7.486e-02  9.279e-01 1.051e-01 -7.124e-01 4.762e-01  -2.809e-01   1.311e-01
mar  -1.421e-01  8.675e-01 1.254e-01 -1.134e+00 2.570e-01  -3.880e-01   1.037e-01
paro -4.134e-02  9.595e-01 9.522e-02 -4.341e-01 6.642e-01  -2.280e-01   1.453e-01
prio  2.639e-01  1.302e+00 8.291e-02  3.182e+00 1.460e-03   1.013e-01   4.264e-01  **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Concordance = 0.640""".strip().split()
            for i in [0, 1, 2, -2, -1]:
                assert output[i] == expected[i]
        finally:
            sys.stdout = saved_stdout

 def test_strata_works_if_only_a_single_element_is_in_the_strata(self):
     df = load_holly_molly_polly()
     del df['Start(days)']
     del df['Stop(days)']
     del df['ID']
     cp = CoxPHFitter()
     cp.fit(df, 'T', 'Status', strata=['Stratum'])
     assert True

    def test_se_against_Survival_Analysis_by_John_Klein_and_Melvin_Moeschberger(self):
        # see table 8.1 in Survival Analysis by John P. Klein and Melvin L. Moeschberger, Second Edition
        df = load_larynx()
        cf = CoxPHFitter(normalize=False)
        cf.fit(df, duration_col='time', event_col='death')

        # standard errors
        actual_se = cf._compute_standard_errors().values
        expected_se = np.array([[0.0143,  0.4623,  0.3561,  0.4222]])
        npt.assert_array_almost_equal(actual_se, expected_se, decimal=2)

    def test_p_value_against_Survival_Analysis_by_John_Klein_and_Melvin_Moeschberger(self):
        # see table 8.1 in Survival Analysis by John P. Klein and Melvin L. Moeschberger, Second Edition
        df = load_larynx()
        cf = CoxPHFitter()
        cf.fit(df, duration_col='time', event_col='death')

        # p-values
        actual_p = cf._compute_p_values()
        expected_p = np.array([0.1847, 0.7644,  0.0730, 0.00])
        npt.assert_array_almost_equal(actual_p, expected_p, decimal=2)

    def test_fit_methods_require_duration_col(self):
        X = load_regression_dataset()

        aaf = AalenAdditiveFitter()
        cph = CoxPHFitter()

        with pytest.raises(TypeError):
            aaf.fit(X)
        with pytest.raises(TypeError):
            cph.fit(X)

    def test_predict_methods_in_regression_return_same_types(self):
        X = load_regression_dataset()

        aaf = AalenAdditiveFitter()
        cph = CoxPHFitter()

        aaf.fit(X, duration_col='T', event_col='E')
        cph.fit(X, duration_col='T', event_col='E')

        for fit_method in ['predict_percentile', 'predict_median', 'predict_expectation', 'predict_survival_function', 'predict_cumulative_hazard']:
            assert isinstance(getattr(aaf, fit_method)(X), type(getattr(cph, fit_method)(X)))

 def test_penalized_output_against_R(self, rossi):
     # R code:
     #
     # rossi <- read.csv('.../lifelines/datasets/rossi.csv')
     # mod.allison <- coxph(Surv(week, arrest) ~ ridge(fin, age, race, wexp, mar, paro, prio,
     #                                                 theta=1.0, scale=FALSE), data=rossi)
     # cat(round(mod.allison$coefficients, 4), sep=", ")
     expected = np.array([[-0.3641, -0.0580, 0.2894, -0.1496, -0.3837, -0.0822, 0.0913]])
     cf = CoxPHFitter(normalize=False, penalizer=1.0)
     cf.fit(rossi, duration_col='week', event_col='arrest')
     npt.assert_array_almost_equal(cf.hazards_.values, expected, decimal=3)

    def test_coef_output_against_Survival_Analysis_by_John_Klein_and_Melvin_Moeschberger(self):
        # see example 8.3 in Survival Analysis by John P. Klein and Melvin L. Moeschberger, Second Edition
        df = load_kidney_transplant(usecols=['time', 'death',
                                             'black_male', 'white_male',
                                             'black_female'])
        cf = CoxPHFitter(normalize=False)
        cf.fit(df, duration_col='time', event_col='death')

        # coefs
        actual_coefs = cf.hazards_.values
        expected_coefs = np.array([[0.1596, 0.2484, 0.6567]])
        npt.assert_array_almost_equal(actual_coefs, expected_coefs, decimal=4)

 def test_output_against_R(self, rossi):
     # from http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-cox-regression.pdf
     # Link is now broken, but this is the code:
     #
     # rossi <- read.csv('.../lifelines/datasets/rossi.csv')
     # mod.allison <- coxph(Surv(week, arrest) ~ fin + age + race + wexp + mar + paro + prio,
     #     data=rossi)
     # cat(round(mod.allison$coefficients, 4), sep=", ")
     expected = np.array([[-0.3794, -0.0574, 0.3139, -0.1498, -0.4337, -0.0849,  0.0915]])
     cf = CoxPHFitter(normalize=False)
     cf.fit(rossi, duration_col='week', event_col='arrest')
     npt.assert_array_almost_equal(cf.hazards_.values, expected, decimal=3)

    def test_summary(self, rossi):

        cp = CoxPHFitter()
        cp.fit(rossi, duration_col='week', event_col='arrest')
        summDf = cp.summary
        expectedColumns = ['coef',
                           'exp(coef)',
                           'se(coef)',
                           'z',
                           'p',
                           'lower 0.95',
                           'upper 0.95']
        assert all([col in summDf.columns for col in expectedColumns])

    def test_strata_against_r_output(self, rossi):
        """
        > r = coxph(formula = Surv(week, arrest) ~ fin + age + strata(race,
            paro, mar, wexp) + prio, data = rossi)
        > r
        > r$loglik
        """

        cp = CoxPHFitter(normalize=False)
        cp.fit(rossi, 'week', 'arrest', strata=['race', 'paro', 'mar', 'wexp'], include_likelihood=True)

        npt.assert_almost_equal(cp.summary['coef'].values, [-0.335, -0.059, 0.100], decimal=3)
        assert abs(cp._log_likelihood - -436.9339) / 436.9339 < 0.01

    def test_fit_methods_can_accept_optional_event_col_param(self):
        X = load_regression_dataset()

        aaf = AalenAdditiveFitter()
        aaf.fit(X, 'T', event_col='E')
        assert_series_equal(aaf.event_observed.sort_index(), X['E'].astype(bool), check_names=False)

        aaf.fit(X, 'T')
        npt.assert_array_equal(aaf.event_observed.values, np.ones(X.shape[0]))

        cph = CoxPHFitter()
        cph.fit(X, 'T', event_col='E')
        assert_series_equal(cph.event_observed.sort_index(), X['E'].astype(bool), check_names=False)

        cph.fit(X, 'T')
        npt.assert_array_equal(cph.event_observed.values, np.ones(X.shape[0]))

    def test_data_normalization(self, data_pred2):
        # During fit, CoxPH copies the training data and normalizes it.
        # Future calls should be normalized in the same way and
        # internal training set should not be saved in a normalized state.

        cf = CoxPHFitter(normalize=True)
        cf.fit(data_pred2, duration_col='t', event_col='E')

        # Internal training set
        ci_trn = concordance_index(cf.durations,
                                   -cf.predict_partial_hazard(cf.data).values,
                                   cf.event_observed)
        # New data should normalize in the exact same way
        ci_org = concordance_index(data_pred2['t'],
                                   -cf.predict_partial_hazard(data_pred2[['x1', 'x2']]).values,
                                   data_pred2['E'])

        assert ci_org == ci_trn