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