Source code for naginterfaces.library.examples.surviv.coxmodel_ex

#!/usr/bin/env python3
"``naginterface.library.surviv.coxmodel`` Python Example."

# NAG Copyright 2018-2019.

# pylint: disable=invalid-name,too-many-locals

from math import log

import numpy as np

from naginterfaces.library import correg, surviv

[docs]def main(): """ Example for :func:`naginterfaces.library.surviv.coxmodel`. Fits Cox's proportional hazard model. >>> main() naginterfaces.library.surviv.coxmodel Python Example Results. Parameter estimates, survival function for a group of leukaemia patients. Parameter Estimate Standard Error 1 -1.5091 0.4096 Deviance = 1.7276e+02 Time Survivor Function 1 0.9640 2 0.9264 3 0.9065 4 0.8661 5 0.8235 6 0.7566 7 0.7343 8 0.6506 10 0.6241 11 0.5724 12 0.5135 13 0.4784 15 0.4447 16 0.4078 17 0.3727 22 0.2859 23 0.1908 """ print('naginterfaces.library.surviv.coxmodel Python Example Results.') print( 'Parameter estimates, survival function for a group ' 'of leukaemia patients.' ) # The number of strata: ns = 0 # Iterations limit: maxit = 12 # The data: t, z, ic, isz = get_data() # The tolerance to use: tol = 5.e-5 # Get initial estimates by fitting an exponential model. y = 1. - np.array(ic, dtype=float) # The number of parameters in the model: ip = isz.count(1) # The offset values (we are fitting a mean in the exponential model): v = np.empty((len(y), ip+1+7)) v[:, 6] = [log(tt) for tt in t] fit = correg.glm_poisson( z, isz, y, v=v, tol=tol, maxit=maxit, ) # Perform the Cox fit. # Stratum indicators, not required: isi = [0] # Move all parameter estimates down one so as to drop the parameter # estimate for the mean: cox_fit = surviv.coxmodel( ns, z, isz, t, ic, isi, fit.b[1:], len(y), tol=tol, maxit=maxit, ) print('Parameter Estimate Standard Error') for i, b_i in enumerate(cox_fit.b): print( '{:4d}'.format(i+1) + ' '*10 + '{:8.4f}'.format(b_i) + ' '*10 + '{:8.4f}'.format(cox_fit.se[i]) ) print('Deviance = {:13.4e}'.format(cox_fit.dev)) print(' Time Survivor Function') for i, tp_i in enumerate(cox_fit.tp[:cox_fit.nd]): print( '{:7.0f}'.format(tp_i) + ' '*6 + ''.join(['{:8.4f}']*max(ns, 1)).format(*cox_fit.sur[i, :]) )
def get_data(): """ The data for this problem. """ # The failure censoring times: t = [ 1., 1., 2., 2., 3., 4., 4., 5., 5., 8., 8., 8., 8., 11., 11., 12., 12., 15., 17., 22., 23., 6., 6., 6., 7., 10., 13., 16., 22., 23., 6., 9., 10., 11., 17., 19., 20., 25., 32., 32., 34., 35., ] # The covariances: z = [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., ] z = np.reshape(z, (len(z), 1)) # The status indicators: ic = [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ] # The covariate inclusions: isz = [1] return t, z, ic, isz if __name__ == '__main__': import doctest import sys sys.exit( doctest.testmod( None, verbose=True, report=False, optionflags=doctest.REPORT_NDIFF, ).failed )