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

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

# 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
)
```