<s> from __future__ import print_function <EOL> import numpy as np <EOL> import pandas as pd <EOL> from numpy import dot , exp <EOL> from numpy . linalg import solve , norm , inv <EOL> from scipy . integrate import trapz <EOL> import scipy . stats as stats <EOL> from lifelines . fitters import BaseFitter <EOL> from lifelines . utils import survival_table_from_events , inv_normal_cdf , normalize , significance_code , concordance_index , _get_index , qth_survival_times <EOL> class CoxPHFitter ( BaseFitter ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , alpha = <NUM_LIT> , tie_method = '<STR_LIT>' , normalize = True , penalizer = <NUM_LIT:0.0> ) : <EOL> if not ( <NUM_LIT:0> < alpha <= <NUM_LIT:1.> ) : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> if penalizer < <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if tie_method != '<STR_LIT>' : <EOL> raise NotImplementedError ( "<STR_LIT>" ) <EOL> self . alpha = alpha <EOL> self . normalize = normalize <EOL> self . tie_method = tie_method <EOL> self . penalizer = penalizer <EOL> self . strata = None <EOL> def _get_efron_values ( self , X , beta , T , E , include_likelihood = False ) : <EOL> """<STR_LIT>""" <EOL> n , d = X . shape <EOL> hessian = np . zeros ( ( d , d ) ) <EOL> gradient = np . zeros ( ( <NUM_LIT:1> , d ) ) <EOL> log_lik = <NUM_LIT:0> <EOL> x_tie_sum = np . zeros ( ( <NUM_LIT:1> , d ) ) <EOL> risk_phi , tie_phi = <NUM_LIT:0> , <NUM_LIT:0> <EOL> risk_phi_x , tie_phi_x = np . zeros ( ( <NUM_LIT:1> , d ) ) , np . zeros ( ( <NUM_LIT:1> , d ) ) <EOL> risk_phi_x_x , tie_phi_x_x = np . zeros ( ( d , d ) ) , np . zeros ( ( d , d ) ) <EOL> tie_count = <NUM_LIT:0> <EOL> for i , ( ti , ei ) in reversed ( list ( enumerate ( zip ( T , E ) ) ) ) : <EOL> xi = X [ i : i + <NUM_LIT:1> ] <EOL> phi_i = exp ( dot ( xi , beta ) ) <EOL> phi_x_i = dot ( phi_i , xi ) <EOL> phi_x_x_i = dot ( xi . T , xi ) * phi_i <EOL> risk_phi += phi_i <EOL> risk_phi_x += phi_x_i <EOL> risk_phi_x_x += phi_x_x_i <EOL> if ei : <EOL> x_tie_sum += xi <EOL> tie_phi += phi_i <EOL> tie_phi_x += phi_x_i <EOL> tie_phi_x_x += phi_x_x_i <EOL> tie_count += <NUM_LIT:1> <EOL> if i > <NUM_LIT:0> and T [ i - <NUM_LIT:1> ] == ti : <EOL> continue <EOL> elif tie_count == <NUM_LIT:0> : <EOL> continue <EOL> partial_gradient = np . zeros ( ( <NUM_LIT:1> , d ) ) <EOL> for l in range ( tie_count ) : <EOL> c = l / tie_count <EOL> denom = ( risk_phi - c * tie_phi ) <EOL> z = ( risk_phi_x - c * tie_phi_x ) <EOL> if denom == <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> partial_gradient += z / denom <EOL> a1 = ( risk_phi_x_x - c * tie_phi_x_x ) / denom <EOL> a2 = dot ( z . T / denom , z / denom ) <EOL> hessian -= ( a1 - a2 ) <EOL> if include_likelihood : <EOL> log_lik -= np . log ( denom ) . ravel ( ) [ <NUM_LIT:0> ] <EOL> gradient += x_tie_sum - partial_gradient <EOL> if include_likelihood : <EOL> log_lik += dot ( x_tie_sum , beta ) . ravel ( ) [ <NUM_LIT:0> ] <EOL> tie_count = <NUM_LIT:0> <EOL> x_tie_sum = np . zeros ( ( <NUM_LIT:1> , d ) ) <EOL> tie_phi = <NUM_LIT:0> <EOL> tie_phi_x = np . zeros ( ( <NUM_LIT:1> , d ) ) <EOL> tie_phi_x_x = np . zeros ( ( d , d ) ) <EOL> if include_likelihood : <EOL> return hessian , gradient , log_lik <EOL> else : <EOL> return hessian , gradient <EOL> def _newton_rhaphson ( self , X , T , E , initial_beta = None , step_size = <NUM_LIT:1.> , <EOL> precision = <NUM_LIT> , show_progress = True , include_likelihood = False ) : <EOL> """<STR_LIT>""" <EOL> assert precision <= <NUM_LIT:1.> , "<STR_LIT>" <EOL> n , d = X . shape <EOL> E = E . astype ( bool ) <EOL> if initial_beta is not None : <EOL> assert initial_beta . shape == ( d , <NUM_LIT:1> ) <EOL> beta = initial_beta <EOL> else : <EOL> beta = np . zeros ( ( d , <NUM_LIT:1> ) ) <EOL> if self . tie_method == '<STR_LIT>' : <EOL> get_gradients = self . _get_efron_values <EOL> else : <EOL> raise NotImplementedError ( "<STR_LIT>" ) <EOL> i = <NUM_LIT:1> <EOL> converging = True <EOL> while converging and i < <NUM_LIT:50> and step_size > <NUM_LIT> : <EOL> if self . strata is None : <EOL> output = get_gradients ( X . values , beta , T . values , E . values , include_likelihood = include_likelihood ) <EOL> h , g = output [ : <NUM_LIT:2> ] <EOL> else : <EOL> g = np . zeros_like ( beta ) . T <EOL> h = np . zeros ( ( beta . shape [ <NUM_LIT:0> ] , beta . shape [ <NUM_LIT:0> ] ) ) <EOL> ll = <NUM_LIT:0> <EOL> for strata in np . unique ( X . index ) : <EOL> stratified_X , stratified_T , stratified_E = X . loc [ [ strata ] ] , T . loc [ [ strata ] ] , E . loc [ [ strata ] ] <EOL> output = get_gradients ( stratified_X . values , beta , stratified_T . values , stratified_E . values , include_likelihood = include_likelihood ) <EOL> _h , _g = output [ : <NUM_LIT:2> ] <EOL> g += _g <EOL> h += _h <EOL> ll += output [ <NUM_LIT:2> ] if include_likelihood else <NUM_LIT:0> <EOL> if self . penalizer > <NUM_LIT:0> : <EOL> g -= self . penalizer * beta . T <EOL> h . flat [ : : d + <NUM_LIT:1> ] -= self . penalizer <EOL> delta = solve ( - h , step_size * g . T ) <EOL> if np . any ( np . isnan ( delta ) ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> hessian , gradient = h , g <EOL> if norm ( delta ) < precision : <EOL> converging = False <EOL> if norm ( delta ) > <NUM_LIT:10> : <EOL> step_size *= <NUM_LIT:0.5> <EOL> continue <EOL> beta += delta <EOL> if ( ( i % <NUM_LIT:10> ) == <NUM_LIT:0> ) and show_progress : <EOL> print ( "<STR_LIT>" % ( i , norm ( delta ) ) ) <EOL> i += <NUM_LIT:1> <EOL> self . _hessian_ = hessian <EOL> self . _score_ = gradient <EOL> if include_likelihood : <EOL> self . _log_likelihood = output [ - <NUM_LIT:1> ] if self . strata is None else ll <EOL> if show_progress : <EOL> print ( "<STR_LIT>" % ( i ) ) <EOL> return beta <EOL> def fit ( self , df , duration_col , event_col = None , <EOL> show_progress = False , initial_beta = None , include_likelihood = False , <EOL> strata = None ) : <EOL> """<STR_LIT>""" <EOL> df = df . copy ( ) <EOL> df . sort_values ( by = duration_col , inplace = True ) <EOL> self . strata = strata <EOL> if strata is not None : <EOL> df = df . set_index ( strata ) <EOL> T = df [ duration_col ] <EOL> del df [ duration_col ] <EOL> if event_col is None : <EOL> E = pd . Series ( np . ones ( df . shape [ <NUM_LIT:0> ] ) , index = df . index ) <EOL> else : <EOL> E = df [ event_col ] <EOL> del df [ event_col ] <EOL> self . data = df if self . strata is None else df . reset_index ( ) <EOL> if self . normalize : <EOL> self . _norm_mean = df . mean ( <NUM_LIT:0> ) <EOL> self . _norm_std = df . std ( <NUM_LIT:0> ) <EOL> df = normalize ( df ) <EOL> E = E . astype ( bool ) <EOL> self . _check_values ( df ) <EOL> hazards_ = self . _newton_rhaphson ( df , T , E , initial_beta = initial_beta , <EOL> show_progress = show_progress , <EOL> include_likelihood = include_likelihood ) <EOL> self . hazards_ = pd . DataFrame ( hazards_ . T , columns = df . columns , <EOL> index = [ '<STR_LIT>' ] ) <EOL> self . confidence_intervals_ = self . _compute_confidence_intervals ( ) <EOL> self . durations = T <EOL> self . event_observed = E <EOL> self . baseline_hazard_ = self . _compute_baseline_hazard ( ) <EOL> self . baseline_cumulative_hazard_ = self . baseline_hazard_ . cumsum ( ) <EOL> self . baseline_survival_ = exp ( - self . baseline_cumulative_hazard_ ) <EOL> return self <EOL> def _check_values ( self , X ) : <EOL> low_var = ( X . var ( <NUM_LIT:0> ) < <NUM_LIT> ) <EOL> if low_var . any ( ) : <EOL> cols = str ( list ( X . columns [ low_var ] ) ) <EOL> print ( "<STR_LIT>" % cols ) <EOL> def _compute_confidence_intervals ( self ) : <EOL> alpha2 = inv_normal_cdf ( ( <NUM_LIT:1.> + self . alpha ) / <NUM_LIT> ) <EOL> se = self . _compute_standard_errors ( ) <EOL> hazards = self . hazards_ . values <EOL> return pd . DataFrame ( np . r_ [ hazards - alpha2 * se , <EOL> hazards + alpha2 * se ] , <EOL> index = [ '<STR_LIT>' , '<STR_LIT>' ] , <EOL> columns = self . hazards_ . columns ) <EOL> def _compute_standard_errors ( self ) : <EOL> se = np . sqrt ( inv ( - self . _hessian_ ) . diagonal ( ) ) <EOL> return pd . DataFrame ( se [ None , : ] , <EOL> index = [ '<STR_LIT>' ] , columns = self . hazards_ . columns ) <EOL> def _compute_z_values ( self ) : <EOL> return ( self . hazards_ . ix [ '<STR_LIT>' ] / <EOL> self . _compute_standard_errors ( ) . ix [ '<STR_LIT>' ] ) <EOL> def _compute_p_values ( self ) : <EOL> U = self . _compute_z_values ( ) ** <NUM_LIT:2> <EOL> return stats . chi2 . sf ( U , <NUM_LIT:1> ) <EOL> @ property <EOL> def summary ( self ) : <EOL> """<STR_LIT>""" <EOL> df = pd . DataFrame ( index = self . hazards_ . columns ) <EOL> df [ '<STR_LIT>' ] = self . hazards_ . ix [ '<STR_LIT>' ] . values <EOL> df [ '<STR_LIT>' ] = exp ( self . hazards_ . ix [ '<STR_LIT>' ] . values ) <EOL> df [ '<STR_LIT>' ] = self . _compute_standard_errors ( ) . ix [ '<STR_LIT>' ] . values <EOL> df [ '<STR_LIT:z>' ] = self . _compute_z_values ( ) <EOL> df [ '<STR_LIT:p>' ] = self . _compute_p_values ( ) <EOL> df [ '<STR_LIT>' % self . alpha ] = self . confidence_intervals_ . ix [ '<STR_LIT>' ] . values <EOL> df [ '<STR_LIT>' % self . alpha ] = self . confidence_intervals_ . ix [ '<STR_LIT>' ] . values <EOL> return df <EOL> def print_summary ( self ) : <EOL> """<STR_LIT>""" <EOL> df = self . summary <EOL> df [ '<STR_LIT>' ] = [ significance_code ( p ) for p in df [ '<STR_LIT:p>' ] ] <EOL> print ( '<STR_LIT>' . format ( self . data . shape [ <NUM_LIT:0> ] , <EOL> np . where ( self . event_observed ) [ <NUM_LIT:0> ] . shape [ <NUM_LIT:0> ] ) , <EOL> end = '<STR_LIT>' ) <EOL> print ( df . to_string ( float_format = lambda f : '<STR_LIT>' . format ( f ) ) ) <EOL> print ( '<STR_LIT>' ) <EOL> print ( "<STR_LIT>" , <EOL> end = '<STR_LIT>' ) <EOL> print ( "<STR_LIT>" <EOL> . format ( concordance_index ( self . durations , <EOL> - self . predict_partial_hazard ( self . data ) . values . ravel ( ) , <EOL> self . event_observed ) ) ) <EOL> return <EOL> def predict_partial_hazard ( self , X ) : <EOL> """<STR_LIT>""" <EOL> index = _get_index ( X ) <EOL> if isinstance ( X , pd . DataFrame ) : <EOL> order = self . hazards_ . columns <EOL> X = X [ order ] <EOL> if self . normalize : <EOL> X = normalize ( X , self . _norm_mean . values , self . _norm_std . values ) <EOL> return pd . DataFrame ( exp ( np . dot ( X , self . hazards_ . T ) ) , index = index ) <EOL> def predict_log_hazard_relative_to_mean ( self , X ) : <EOL> """<STR_LIT>""" <EOL> mean_covariates = self . data . mean ( <NUM_LIT:0> ) . to_frame ( ) . T <EOL> return np . log ( self . predict_partial_hazard ( X ) / self . predict_partial_hazard ( mean_covariates ) . squeeze ( ) ) <EOL> def predict_cumulative_hazard ( self , X ) : <EOL> """<STR_LIT>""" <EOL> v = self . predict_partial_hazard ( X ) <EOL> s_0 = self . baseline_survival_ <EOL> col = _get_index ( X ) <EOL> return pd . DataFrame ( - np . dot ( np . log ( s_0 ) , v . T ) , index = self . baseline_survival_ . index , columns = col ) <EOL> def predict_survival_function ( self , X ) : <EOL> """<STR_LIT>""" <EOL> return exp ( - self . predict_cumulative_hazard ( X ) ) <EOL> def predict_percentile ( self , X , p = <NUM_LIT:0.5> ) : <EOL> """<STR_LIT>""" <EOL> index = _get_index ( X ) <EOL> return qth_survival_times ( p , self . predict_survival_function ( X ) [ index ] ) <EOL> def predict_median ( self , X ) : <EOL> """<STR_LIT>""" <EOL> return self . predict_percentile ( X , <NUM_LIT:0.5> ) <EOL> def predict_expectation ( self , X ) : <EOL> """<STR_LIT>""" <EOL> index = _get_index ( X ) <EOL> v = self . predict_survival_function ( X ) [ index ] <EOL> return pd . DataFrame ( trapz ( v . values . T , v . index ) , index = index ) <EOL> def _compute_baseline_hazard ( self ) : <EOL> ind_hazards = self . predict_partial_hazard ( self . data ) . values <EOL> event_table = survival_table_from_events ( self . durations . values , <EOL> self . event_observed . values ) <EOL> baseline_hazard_ = pd . DataFrame ( np . zeros ( ( event_table . shape [ <NUM_LIT:0> ] , <NUM_LIT:1> ) ) , <EOL> index = event_table . index , <EOL> columns = [ '<STR_LIT>' ] ) <EOL> for t , s in event_table . iterrows ( ) : <EOL> less = np . array ( self . durations >= t ) <EOL> if ind_hazards [ less ] . sum ( ) == <NUM_LIT:0> : <EOL> v = <NUM_LIT:0> <EOL> else : <EOL> v = ( s [ '<STR_LIT>' ] / ind_hazards [ less ] . sum ( ) ) <EOL> baseline_hazard_ . ix [ t ] = v <EOL> return baseline_hazard_ </s>
