1library(survival)
2library(R2nparray)
3
4ixd = list(c(20,1), c(50,1), c(50,2), c(100,5), c(1000,10))
5
6res = list()
7
8for (ix in ixd) {
9    fname = sprintf("results/survival_data_%d_%d.csv", ix[1], ix[2])
10    data = read.table(fname)
11
12    time = data[,1]
13    status = data[,2]
14    entry = data[,3]
15    exog = data[,4:dim(data)[2]]
16    exog = as.matrix(exog)
17    n = dim(exog)[1]
18    p = dim(exog)[2]
19
20    # Needs to match the kronecker statement in test_phreg.py
21    strata = kronecker(seq(5), array(1, n/5))
22
23    for (ties in c("breslow", "efron")) {
24
25        ti = substr(ties, 1, 3)
26
27        # Base model
28        surv = Surv(time, status)
29        md = coxph(surv ~ exog, ties=ties)
30        tag = sprintf("%d_%d_%s", n, p, ti)
31        res[[sprintf("coef_%s", tag)]] = md$coef
32        res[[sprintf("se_%s", tag)]] = sqrt(diag(md$var))
33        #bhaz = basehaz(md)
34        bhaz = survfit(md, type="aalen")
35        #bhaz = survfit(md, type="efron")
36        res[[sprintf("time_%s", tag)]] = bhaz$time
37        res[[sprintf("hazard_%s", tag)]] = -log(bhaz$surv)
38
39        # With entry time
40        surv = Surv(entry, time, status)
41        md = coxph(surv ~ exog, ties=ties)
42        tag = sprintf("%d_%d_et_%s", n, p, ti)
43        res[[sprintf("coef_%s", tag)]] = md$coef
44        res[[sprintf("se_%s", tag)]] = sqrt(diag(md$var))
45        res[[sprintf("time_%s", tag)]] = c(0)
46        res[[sprintf("hazard_%s", tag)]] = c(0)
47
48        # With strata
49        surv = Surv(time, status)
50        md = coxph(surv ~ exog + strata(strata), ties=ties)
51        tag = sprintf("%d_%d_st_%s", n, p, ti)
52        res[[sprintf("coef_%s", tag)]] = md$coef
53        res[[sprintf("se_%s", tag)]] = sqrt(diag(md$var))
54        res[[sprintf("time_%s", tag)]] = c(0)
55        res[[sprintf("hazard_%s", tag)]] = c(0)
56
57        # With entry time and strata
58        surv = Surv(entry, time, status)
59        md = coxph(surv ~ exog + strata(strata), ties=ties)
60        tag = sprintf("%d_%d_et_st_%s", n, p, ti)
61        res[[sprintf("coef_%s", tag)]] = md$coef
62        res[[sprintf("se_%s", tag)]] = sqrt(diag(md$var))
63        res[[sprintf("time_%s", tag)]] = c(0)
64        res[[sprintf("hazard_%s", tag)]] = c(0)
65    }
66}
67
68R2nparray(res, fname="survival_r_results.py")
69