<s> '''<STR_LIT>''' <EOL> from __future__ import print_function <EOL> from statsmodels . compat . python import iterkeys , range <EOL> import shutil <EOL> import tempfile <EOL> import numpy as np <EOL> from statsmodels . base import model <EOL> from . import survival <EOL> class Discrete ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , x , w = None ) : <EOL> self . x = np . squeeze ( x ) <EOL> if self . x . shape == ( ) : <EOL> self . x = np . array ( [ self . x ] ) <EOL> self . n = self . x . shape [ <NUM_LIT:0> ] <EOL> if w is None : <EOL> w = np . ones ( self . n , np . float64 ) <EOL> else : <EOL> if w . shape [ <NUM_LIT:0> ] != self . n : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> if np . any ( np . less ( w , <NUM_LIT:0> ) ) : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> self . w = w * <NUM_LIT:1.0> / w . sum ( ) <EOL> def mean ( self , f = None ) : <EOL> if f is None : <EOL> fx = self . x <EOL> else : <EOL> fx = f ( self . x ) <EOL> return ( fx * self . w ) . sum ( ) <EOL> def cov ( self ) : <EOL> mu = self . mean ( ) <EOL> dx = self . x - mu <EOL> return np . dot ( dx , np . transpose ( dx ) ) <EOL> class Observation ( survival . RightCensored ) : <EOL> def __getitem__ ( self , item ) : <EOL> if self . namespace is not None : <EOL> return self . namespace [ item ] <EOL> else : <EOL> return getattr ( self , item ) <EOL> def __init__ ( self , time , delta , namespace = None ) : <EOL> self . namespace = namespace <EOL> survival . RightCensored . __init__ ( self , time , delta ) <EOL> def __call__ ( self , formula , time = None , ** extra ) : <EOL> return formula ( namespace = self , time = time , ** extra ) <EOL> class CoxPH ( model . LikelihoodModel ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , subjects , formula , time_dependent = False ) : <EOL> self . subjects , self . formula = subjects , formula <EOL> self . time_dependent = time_dependent <EOL> self . initialize ( self . subjects ) <EOL> def initialize ( self , subjects ) : <EOL> print ( '<STR_LIT>' ) <EOL> self . failures = { } <EOL> for i in range ( len ( subjects ) ) : <EOL> s = subjects [ i ] <EOL> if s . delta : <EOL> if s . time not in self . failures : <EOL> self . failures [ s . time ] = [ i ] <EOL> else : <EOL> self . failures [ s . time ] . append ( i ) <EOL> self . failure_times = list ( iterkeys ( self . failures ) ) <EOL> self . failure_times . sort ( ) <EOL> def cache ( self ) : <EOL> if self . time_dependent : <EOL> self . cachedir = tempfile . mkdtemp ( ) <EOL> self . design = { } <EOL> self . risk = { } <EOL> first = True <EOL> for t in iterkeys ( self . failures ) : <EOL> if self . time_dependent : <EOL> d = np . array ( [ s ( self . formula , time = t ) <EOL> for s in self . subjects ] ) . astype ( float ) [ : , None ] <EOL> dshape = d . shape <EOL> dfile = file ( tempfile . mkstemp ( dir = self . cachedir ) [ <NUM_LIT:1> ] , '<STR_LIT:w>' ) <EOL> d . tofile ( dfile ) <EOL> dfile . close ( ) <EOL> del ( d ) <EOL> self . design [ t ] = np . memmap ( dfile . name , <EOL> dtype = np . dtype ( float ) , <EOL> shape = dshape ) <EOL> elif first : <EOL> d = np . array ( [ s ( self . formula , time = t ) <EOL> for s in self . subjects ] ) . astype ( np . float64 ) <EOL> self . design [ t ] = d <EOL> else : <EOL> self . design [ t ] = d <EOL> self . risk [ t ] = np . compress ( [ s . atrisk ( t ) for s in self . subjects ] , <EOL> np . arange ( self . design [ t ] . shape [ <NUM_LIT:0> ] ) , axis = - <NUM_LIT:1> ) <EOL> def __del__ ( self ) : <EOL> try : <EOL> shutil . rmtree ( self . cachedir , ignore_errors = True ) <EOL> except AttributeError : <EOL> print ( "<STR_LIT>" ) <EOL> pass <EOL> def loglike ( self , b , ties = '<STR_LIT>' ) : <EOL> logL = <NUM_LIT:0> <EOL> for t in iterkeys ( self . failures ) : <EOL> fail = self . failures [ t ] <EOL> d = len ( fail ) <EOL> risk = self . risk [ t ] <EOL> Zb = np . dot ( self . design [ t ] , b ) <EOL> logL += Zb [ fail ] . sum ( ) <EOL> if ties == '<STR_LIT>' : <EOL> s = np . exp ( Zb [ risk ] ) . sum ( ) <EOL> logL -= np . log ( np . exp ( Zb [ risk ] ) . sum ( ) ) * d <EOL> elif ties == '<STR_LIT>' : <EOL> s = np . exp ( Zb [ risk ] ) . sum ( ) <EOL> r = np . exp ( Zb [ fail ] ) . sum ( ) <EOL> for j in range ( d ) : <EOL> logL -= np . log ( s - j * r / d ) <EOL> elif ties == '<STR_LIT>' : <EOL> raise NotImplementedError ( '<STR_LIT>' ) <EOL> else : <EOL> raise NotImplementedError ( '<STR_LIT>' ) <EOL> return logL <EOL> def score ( self , b , ties = '<STR_LIT>' ) : <EOL> score = <NUM_LIT:0> <EOL> for t in iterkeys ( self . failures ) : <EOL> fail = self . failures [ t ] <EOL> d = len ( fail ) <EOL> risk = self . risk [ t ] <EOL> Z = self . design [ t ] <EOL> score += Z [ fail ] . sum ( ) <EOL> if ties == '<STR_LIT>' : <EOL> w = np . exp ( np . dot ( Z , b ) ) <EOL> rv = Discrete ( Z [ risk ] , w = w [ risk ] ) <EOL> score -= rv . mean ( ) * d <EOL> elif ties == '<STR_LIT>' : <EOL> w = np . exp ( np . dot ( Z , b ) ) <EOL> score += Z [ fail ] . sum ( ) <EOL> for j in range ( d ) : <EOL> efron_w = w <EOL> efron_w [ fail ] -= i * w [ fail ] / float ( d ) <EOL> rv = Discrete ( Z [ risk ] , w = efron_w [ risk ] ) <EOL> score -= rv . mean ( ) <EOL> elif ties == '<STR_LIT>' : <EOL> raise NotImplementedError ( '<STR_LIT>' ) <EOL> else : <EOL> raise NotImplementedError ( '<STR_LIT>' ) <EOL> return np . array ( [ score ] ) <EOL> def information ( self , b , ties = '<STR_LIT>' ) : <EOL> info = <NUM_LIT:0> <EOL> score = <NUM_LIT:0> <EOL> for t in iterkeys ( self . failures ) : <EOL> fail = self . failures [ t ] <EOL> d = len ( fail ) <EOL> risk = self . risk [ t ] <EOL> Z = self . design [ t ] <EOL> if ties == '<STR_LIT>' : <EOL> w = np . exp ( np . dot ( Z , b ) ) <EOL> rv = Discrete ( Z [ risk ] , w = w [ risk ] ) <EOL> info += rv . cov ( ) <EOL> elif ties == '<STR_LIT>' : <EOL> w = np . exp ( np . dot ( Z , b ) ) <EOL> score += Z [ fail ] . sum ( ) <EOL> for j in range ( d ) : <EOL> efron_w = w <EOL> efron_w [ fail ] -= i * w [ fail ] / d <EOL> rv = Discrete ( Z [ risk ] , w = efron_w [ risk ] ) <EOL> info += rv . cov ( ) <EOL> elif ties == '<STR_LIT>' : <EOL> raise NotImplementedError ( '<STR_LIT>' ) <EOL> else : <EOL> raise NotImplementedError ( '<STR_LIT>' ) <EOL> return score <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> import numpy . random as R <EOL> n = <NUM_LIT:100> <EOL> X = np . array ( [ <NUM_LIT:0> ] * n + [ <NUM_LIT:1> ] * n ) <EOL> b = <NUM_LIT> <EOL> lin = <NUM_LIT:1> + b * X <EOL> Y = R . standard_exponential ( ( <NUM_LIT:2> * n , ) ) / lin <EOL> delta = R . binomial ( <NUM_LIT:1> , <NUM_LIT> , size = ( <NUM_LIT:2> * n , ) ) <EOL> subjects = [ Observation ( Y [ i ] , delta [ i ] ) for i in range ( <NUM_LIT:2> * n ) ] <EOL> for i in range ( <NUM_LIT:2> * n ) : <EOL> subjects [ i ] . X = X [ i ] <EOL> import statsmodels . sandbox . formula as F <EOL> x = F . Quantitative ( '<STR_LIT:X>' ) <EOL> f = F . Formula ( x ) <EOL> c = CoxPH ( subjects , f ) <EOL> c = CoxPH ( subjects , f , time_dependent = True ) <EOL> c . cache ( ) <EOL> res = c . fit ( [ <NUM_LIT> ] , method = "<STR_LIT>" ) <EOL> print ( res . params ) <EOL> print ( dir ( c ) ) <EOL> '''<STR_LIT>''' </s>
