#!/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
)
