from dsm import datasets, DeepSurvivalMachines
import numpy as np
x, t, e = datasets.load_dataset('SUPPORT')
times = np.quantile(t[e==1], [0.25, 0.5, 0.75]).tolist()
cv_folds = 5
folds = list(range(cv_folds))*10000
folds = np.array(folds[:len(x)])
from sksurv.metrics import concordance_index_ipcw, brier_score
cis = []
brs = []
for fold in range(cv_folds):
print ("On Fold:", fold)
x_train, t_train, e_train = x[folds!=fold], t[folds!=fold], e[folds!=fold]
x_test, t_test, e_test = x[folds==fold], t[folds==fold], e[folds==fold]
print (x_train.shape)
model = DeepSurvivalMachines(distribution='Weibull', layers=[100])
model.fit(x_train, t_train, e_train, iters=100, learning_rate=1e-3,batch_size=10)
et_train = np.array([(e_train[i], t_train[i]) for i in range(len(e_train))],
dtype=[('e', bool), ('t', int)])
et_test = np.array([(e_test[i], t_test[i]) for i in range(len(e_test))],
dtype=[('e', bool), ('t', int)])
out_risk = model.predict_risk(x_test, times)
out_survival = model.predict_survival(x_test, times)
cis_ = []
for i in range(len(times)):
cis_.append(concordance_index_ipcw(et_train, et_test, out_risk[:,i], times[i])[0])
cis.append(cis_)
brs.append(brier_score(et_train, et_test, out_survival, times )[1])
1%| | 54/10000 [00:00<00:18, 536.02it/s]
On Fold: 0 (7284, 44) Pretraining the Underlying Distributions... torch.Size([6192]) torch.Size([6192]) torch.Size([6192]) torch.Size([6192])
13%|█▎ | 1312/10000 [00:02<00:13, 628.77it/s] 0%| | 0/100 [00:00<?, ?it/s]
-0.7721933727751713 -6.453758404040128
4%|▍ | 4/100 [00:07<02:56, 1.84s/it] 1%| | 60/10000 [00:00<00:16, 594.79it/s]
On Fold: 1 (7284, 44) Pretraining the Underlying Distributions... torch.Size([6192]) torch.Size([6192]) torch.Size([6192]) torch.Size([6192])
14%|█▍ | 1438/10000 [00:02<00:13, 640.39it/s] 0%| | 0/100 [00:00<?, ?it/s]
-0.7678721461202748 -6.509151589080709
7%|▋ | 7/100 [00:11<02:29, 1.61s/it] 1%| | 63/10000 [00:00<00:15, 629.21it/s]
On Fold: 2 (7284, 44) Pretraining the Underlying Distributions... torch.Size([6192]) torch.Size([6192]) torch.Size([6192]) torch.Size([6192])
16%|█▌ | 1582/10000 [00:02<00:13, 641.78it/s] 0%| | 0/100 [00:00<?, ?it/s]
-0.7727229786575149 -6.511681205334374
6%|▌ | 6/100 [00:10<02:38, 1.68s/it] 1%| | 53/10000 [00:00<00:19, 519.58it/s]
On Fold: 3 (7284, 44) Pretraining the Underlying Distributions... torch.Size([6192]) torch.Size([6192]) torch.Size([6192]) torch.Size([6192])
14%|█▍ | 1385/10000 [00:02<00:16, 515.82it/s] 0%| | 0/100 [00:00<?, ?it/s]
-0.7719732355296647 -6.47692660327643
5%|▌ | 5/100 [00:14<04:36, 2.91s/it] 0%| | 40/10000 [00:00<00:25, 391.81it/s]
On Fold: 4 (7284, 44) Pretraining the Underlying Distributions... torch.Size([6192]) torch.Size([6192]) torch.Size([6192]) torch.Size([6192])
18%|█▊ | 1799/10000 [00:04<00:18, 433.18it/s] 0%| | 0/100 [00:00<?, ?it/s]
-0.7750461172086791 -6.534539612585216
5%|▌ | 5/100 [00:11<03:30, 2.22s/it]
print ("Concordance Index:", np.mean(cis,axis=0))
print ("Brier Score:", np.mean(brs,axis=0))
Concordance Index: [0.74546862 0.706156 0.67491647] Brier Score: [0.12746457 0.19834306 0.21601803]
from sksurv.linear_model import CoxPHSurvivalAnalysis
cis = []
for fold in range(cv_folds):
print ("On Fold:", fold)
x_train, t_train, e_train = x[folds!=fold], t[folds!=fold], e[folds!=fold]
x_test, t_test, e_test = x[folds==fold], t[folds==fold], e[folds==fold]
et_train = np.array([(e_train[i], t_train[i]) for i in range(len(e_train))],
dtype=[('e', bool), ('t', int)])
et_test = np.array([(e_test[i], t_test[i]) for i in range(len(e_test))],
dtype=[('e', bool), ('t', int)])
model = CoxPHSurvivalAnalysis(alpha=1e-3)
model.fit(x_test, et_test)
out_risk = model.predict_survival_function(x_test)
cis_ = []
for i in range(len(times)):
cis_.append(concordance_index_ipcw(et_train, et_test, out_risk, times[i])[0])
cis.append(cis_)
On Fold: 0
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-34-105a4002e245> in <module> 22 23 for i in range(len(times)): ---> 24 cis_.append(concordance_index_ipcw(et_train, et_test, out_risk, times[i])[0]) 25 26 cis.append(cis_) ~/anaconda3/lib/python3.8/site-packages/sksurv/metrics.py in concordance_index_ipcw(survival_train, survival_test, estimate, tau, tied_tol) 300 survival_test = survival_test[mask] 301 --> 302 estimate = _check_estimate(estimate, test_time) 303 304 cens = CensoringDistributionEstimator() ~/anaconda3/lib/python3.8/site-packages/sksurv/metrics.py in _check_estimate(estimate, test_time) 29 30 def _check_estimate(estimate, test_time): ---> 31 estimate = check_array(estimate, ensure_2d=False) 32 if estimate.ndim != 1: 33 raise ValueError( ~/anaconda3/lib/python3.8/site-packages/sklearn/utils/validation.py in inner_f(*args, **kwargs) 71 FutureWarning) 72 kwargs.update({k: arg for k, arg in zip(sig.parameters, args)}) ---> 73 return f(**kwargs) 74 return inner_f 75 ~/anaconda3/lib/python3.8/site-packages/sklearn/utils/validation.py in check_array(array, accept_sparse, accept_large_sparse, dtype, order, copy, force_all_finite, ensure_2d, allow_nd, ensure_min_samples, ensure_min_features, estimator) 597 array = array.astype(dtype, casting="unsafe", copy=False) 598 else: --> 599 array = np.asarray(array, order=order, dtype=dtype) 600 except ComplexWarning: 601 raise ValueError("Complex data not supported\n" ~/anaconda3/lib/python3.8/site-packages/numpy/core/_asarray.py in asarray(a, dtype, order) 83 84 """ ---> 85 return array(a, dtype, copy=False, order=order) 86 87 TypeError: float() argument must be a string or a number, not 'StepFunction'
time = 6
int(np.where(out_risk[0].x == time)[0])
3
out_risk[0].x
array([ 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 74, 75, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 88, 90, 91, 92, 93, 94, 95, 96, 97, 98, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 114, 116, 117, 118, 119, 120, 121, 122, 124, 126, 127, 128, 129, 130, 132, 133, 134, 136, 137, 139, 142, 143, 145, 146, 147, 148, 149, 151, 152, 153, 156, 157, 160, 162, 163, 164, 165, 166, 167, 168, 170, 171, 172, 173, 174, 176, 180, 181, 182, 183, 185, 186, 187, 189, 191, 193, 194, 195, 197, 198, 199, 200, 201, 202, 203, 204, 205, 207, 208, 212, 213, 214, 215, 217, 218, 220, 223, 224, 225, 227, 229, 230, 231, 233, 234, 235, 236, 237, 240, 242, 244, 247, 248, 251, 252, 253, 254, 258, 259, 260, 263, 264, 265, 266, 268, 269, 273, 274, 276, 277, 279, 281, 283, 287, 288, 290, 291, 292, 294, 295, 297, 299, 300, 303, 309, 310, 311, 312, 313, 314, 316, 318, 319, 320, 321, 322, 323, 324, 326, 328, 330, 335, 338, 339, 340, 343, 344, 346, 347, 348, 350, 351, 352, 353, 355, 356, 359, 360, 361, 363, 365, 366, 368, 370, 372, 377, 379, 380, 381, 382, 384, 385, 386, 387, 389, 392, 393, 394, 395, 396, 397, 399, 400, 401, 403, 404, 405, 406, 407, 408, 409, 410, 411, 413, 415, 417, 418, 420, 421, 422, 423, 425, 428, 430, 432, 433, 434, 435, 436, 440, 442, 444, 446, 447, 448, 449, 450, 451, 453, 455, 459, 460, 461, 463, 464, 465, 467, 468, 469, 470, 472, 473, 474, 477, 479, 482, 484, 485, 486, 487, 489, 491, 492, 493, 494, 496, 497, 499, 500, 501, 503, 504, 507, 509, 511, 513, 515, 517, 518, 521, 523, 524, 526, 527, 528, 529, 531, 533, 534, 536, 541, 542, 546, 548, 551, 552, 553, 554, 555, 557, 558, 560, 562, 563, 564, 566, 567, 573, 575, 576, 577, 578, 582, 584, 585, 586, 587, 588, 589, 591, 595, 597, 599, 603, 604, 605, 608, 609, 610, 613, 615, 616, 617, 618, 619, 620, 621, 623, 624, 626, 627, 628, 629, 631, 632, 633, 634, 636, 637, 641, 643, 644, 648, 649, 650, 652, 653, 655, 656, 657, 658, 659, 661, 662, 664, 665, 666, 667, 668, 669, 670, 671, 674, 675, 677, 679, 680, 682, 685, 686, 690, 692, 695, 702, 703, 705, 706, 707, 708, 709, 710, 712, 714, 716, 717, 719, 720, 721, 724, 726, 727, 734, 738, 741, 744, 745, 746, 747, 751, 756, 757, 760, 761, 763, 765, 766, 768, 770, 772, 773, 774, 776, 777, 779, 781, 783, 784, 786, 789, 790, 795, 797, 798, 799, 800, 803, 804, 807, 808, 809, 811, 812, 814, 815, 816, 817, 818, 819, 820, 821, 823, 824, 825, 827, 829, 830, 831, 833, 835, 839, 842, 844, 845, 847, 849, 851, 852, 853, 855, 857, 858, 861, 864, 867, 868, 869, 872, 873, 875, 877, 878, 879, 883, 885, 887, 889, 890, 891, 892, 897, 898, 904, 910, 914, 917, 918, 919, 923, 926, 928, 929, 930, 934, 936, 937, 940, 941, 944, 946, 950, 951, 954, 958, 965, 969, 970, 971, 972, 973, 977, 978, 982, 984, 985, 986, 987, 988, 989, 992, 996, 998, 999, 1000, 1006, 1009, 1011, 1012, 1017, 1018, 1021, 1022, 1023, 1029, 1034, 1036, 1037, 1043, 1045, 1046, 1047, 1049, 1050, 1051, 1055, 1059, 1060, 1063, 1064, 1068, 1070, 1072, 1073, 1074, 1075, 1078, 1079, 1082, 1087, 1088, 1099, 1109, 1116, 1126, 1134, 1138, 1142, 1162, 1164, 1172, 1174, 1177, 1185, 1201, 1212, 1213, 1224, 1227, 1232, 1238, 1250, 1253, 1265, 1269, 1289, 1301, 1304, 1307, 1310, 1312, 1320, 1321, 1326, 1327, 1328, 1342, 1344, 1345, 1346, 1347, 1349, 1352, 1355, 1356, 1360, 1363, 1365, 1369, 1371, 1373, 1377, 1379, 1380, 1382, 1384, 1385, 1388, 1391, 1392, 1396, 1398, 1401, 1406, 1409, 1410, 1411, 1416, 1418, 1421, 1422, 1427, 1439, 1441, 1442, 1444, 1449, 1452, 1455, 1458, 1466, 1467, 1474, 1475, 1484, 1485, 1486, 1487, 1489, 1492, 1495, 1497, 1503, 1510, 1512, 1514, 1517, 1518, 1519, 1521, 1530, 1531, 1534, 1539, 1542, 1543, 1547, 1551, 1552, 1558, 1560, 1563, 1566, 1567, 1568, 1572, 1573, 1578, 1579, 1593, 1596, 1599, 1600, 1605, 1610, 1613, 1614, 1618, 1622, 1623, 1629, 1636, 1642, 1647, 1648, 1654, 1655, 1657, 1659, 1665, 1670, 1671, 1676, 1677, 1681, 1683, 1686, 1688, 1689, 1691, 1697, 1699, 1701, 1705, 1712, 1717, 1718, 1719, 1722, 1723, 1728, 1732, 1733, 1734, 1739, 1740, 1742, 1745, 1747, 1748, 1761, 1763, 1767, 1769, 1772, 1778, 1782, 1783, 1785, 1788, 1790, 1792, 1795, 1798, 1801, 1807, 1812, 1814, 1819, 1820, 1823, 1825, 1826, 1830, 1845, 1853, 1857, 1863, 1866, 1867, 1882, 1885, 1886, 1887, 1892, 1910, 1911, 1915, 1916, 1918, 1921, 1928, 1938, 1940, 1944, 1945, 1948, 1949, 1951, 1952, 1963, 1971, 1976, 1978, 1979, 1980, 1984, 1990, 1992, 1995, 1998, 1999, 2001, 2007, 2009, 2010, 2012, 2014, 2016, 2019, 2022, 2024, 2026, 2027, 2028, 2029])
model = CoxPHSurvivalAnalysis(alpha=1e-3)
model.fit(x_test, et_test)
np.mean(cis,axis=0)
array([0.74335312, 0.7045087 , 0.68096073])
out_risk = model.predict_risk(x, times)
model.torch_model.eval()
DeepSurvivalMachinesTorch( (act): SELU() (embedding): Sequential( (0): Linear(in_features=44, out_features=100, bias=False) (1): ReLU6() (2): Linear(in_features=100, out_features=100, bias=False) (3): ReLU6() ) (gate): Sequential( (0): Linear(in_features=100, out_features=3, bias=False) ) (scaleg): Sequential( (0): Linear(in_features=100, out_features=3, bias=True) ) (shapeg): Sequential( (0): Linear(in_features=100, out_features=3, bias=True) ) )
out_survival = model.predict_survival(x, times)
from matplotlib import pyplot as plt
from sksurv.metrics import brier_score, concordance_index_ipcw
import numpy as np
et = np.array([(e[i], t[i]) for i in range(len(e))],
dtype=[('e', bool), ('t', int)])
brier_score(et, et, out_survival, times )
array([0.13039755, 0.20234974, 0.21643684])
for i in range(len(times)):
print(concordance_index_ipcw(et, et, out_risk[:,i], times[i])[0])
0.7519513749695589 0.7074775823879251 0.678728630898966
from sksurv.linear_model import CoxPHSurvivalAnalysis
estimator = CoxPHSurvivalAnalysis(alpha=1e-3).fit(x, et,)
surv_funcs = estimator.predict(x)
surv_funcs
array([ 0.86249313, 0.16849345, -0.45380257, ..., -0.14997697, 0.35619347, -0.12209867])
for i in range(len(times)):
print(concordance_index_ipcw(et, et, surv_funcs, times[i])[0])
0.6924659134706312 0.6741630293711603 0.6724802772351569