Skip to content

Instantly share code, notes, and snippets.

@CSJCampbell

CSJCampbell/survfit.coxph.patch

Created April 30, 2021 23:07
Show Gist options
  • Save CSJCampbell/09dfc0cee3365b3b38fccf0b108b209f to your computer and use it in GitHub Desktop.
Save CSJCampbell/09dfc0cee3365b3b38fccf0b108b209f to your computer and use it in GitHub Desktop.
Download ZIP
#130 - Survival curves of stratified coxph objects
Raw
survfit.coxph.patch
diff --git a/R/survfit.coxph.R b/R/survfit.coxph.R
index 43b10c4..25872e5 100644
--- a/R/survfit.coxph.R
+++ b/R/survfit.coxph.R
@@ -110,15 +110,6 @@ survfit.coxph <-
}
}
- }
- if (FALSE) {
- if (!is.null(mf)){
- y2 <- object[['y']]
- if (!is.null(y2)) {
- if (ncol(y2) != ncol(Y) || length(y2) != length(Y))
- stop("Could not reconstruct the y vector")
- }
- }
}
type <- attr(Y, 'type')
if (!type %in% c("right", "counting", "mright", "mcounting"))
@@ -247,6 +238,7 @@ survfit.coxph <-
else strata2 <- factor(rep(0, nrow(mf2)))
if (!robust) cluster <- NULL
+ x2 <- y2 <- NULL
if (individual) {
if (missing(newdata))
stop("The newdata argument must be present when individual=TRUE")
@@ -301,34 +293,16 @@ survfit.coxph <-
y2, x2, risk2, strata2, id2)
}
else {
- result <- coxsurv.fit(ctype, stype, se.fit, varmat, cluster,
- Y, X, weights, risk, position, strata, oldid,
- y2, x2, risk2)
- if (has.strata && found.strata) {
- if (is.matrix(result$surv)) {
- nr <- nrow(result$surv) #a vector if newdata had only 1 row
- indx1 <- split(1:nr, rep(1:length(result$strata), result$strata))
- rows <- indx1[as.numeric(strata2)] #the rows for each curve
-
- indx2 <- unlist(rows) #index for time, n.risk, n.event, n.censor
- indx3 <- as.integer(strata2) #index for n and strata
-
- for(i in 2:length(rows)) rows[[i]] <- rows[[i]]+ (i-1)*nr #linear subscript
- indx4 <- unlist(rows) #index for surv and std.err
- temp <- result$strata[indx3]
- names(temp) <- row.names(mf2)
- new <- list(n = result$n[indx3],
- time= result$time[indx2],
- n.risk= result$n.risk[indx2],
- n.event=result$n.event[indx2],
- n.censor=result$n.censor[indx2],
- strata = temp,
- surv= result$surv[indx4],
- cumhaz = result$cumhaz[indx4])
- if (se.fit) new$std.err <- result$std.err[indx4]
- result <- new
- }
- }
+ if (has.strata && found.strata) {
+ Y <- Y[strata %in% strata2]
+ X <- X[strata %in% strata2, , drop = FALSE]
+ if (!is.null(weights)) { weights <- weights[strata %in% strata2] }
+ risk <- risk[strata %in% strata2]
+ strata <- factor(strata[strata %in% strata2])
+ }
+ result <- coxsurv.fit(ctype, stype, se.fit, varmat, cluster,
+ Y, X, weights, risk, position, strata, oldid,
+ y2, x2, risk2)
}
if (!censor) {
kfun <- function(x, keep){ if (is.matrix(x)) x[keep,,drop=F]
diff --git a/tests/survfit3.R b/tests/survfit3.R
new file mode 100644
index 0000000..899c8ce
--- /dev/null
+++ b/tests/survfit3.R
@@ -0,0 +1,249 @@
+library(survival)
+
+set.seed(345726)
+df <- data.frame(X1 = runif(200),
+ X2 = sample(c("A", "B"), 200, replace = TRUE),
+ Ev = sample(c(0,1), 200, replace = TRUE),
+ Time = rexp(200))
+
+testfit0 <- coxph(formula = Surv(Time, Ev) ~ X1, data = df)
+
+out0 <- survfit(testfit0, newdata = data.frame(X1 = 0.6))
+
+out_compare0 <- c(
+ records = 200,
+ n.max = 200,
+ n.start = 200,
+ events = 90,
+ median = 1.20017089876802,
+ `0.95LCL` = 0.953603094792439,
+ `0.95UCL` = 1.80559369245789)
+
+all.equal(out_compare0, survival:::survmean(out0, rmean = "none")$matrix)
+
+testfit1 <- coxph(formula = Surv(Time, Ev) ~ X1 + strata(X2), data = df)
+
+out1 <- survfit(testfit1, newdata = data.frame(X1 = 0.6, X2 = "A"))
+
+all.equal(sum(df$X2 == "A"), length(out1$surv))
+
+out_compare1 <- structure(
+ list(
+ n = 96L,
+ time = c(0.00598399853333831, 0.0267057245032998,
+ 0.050189660402173, 0.0644403744917489, 0.0701612638463647, 0.0724950589515161,
+ 0.0793271128505524, 0.093478377442807, 0.0947427659133707, 0.0982021834805559,
+ 0.129375483375043, 0.137765689287335, 0.154622706118971, 0.17002234675396,
+ 0.170190538301736, 0.186400834576651, 0.198276239447296, 0.20554883684963,
+ 0.209514521677677, 0.216919615399092, 0.226764952763915, 0.231043137754327,
+ 0.231377071334626, 0.285285478699579, 0.302137257799159, 0.312510287730083,
+ 0.320241346650024, 0.330649028066546, 0.338164258282632, 0.349865917112634,
+ 0.351110461633652, 0.356002631597221, 0.361543103586882, 0.387070584345154,
+ 0.390072865877301, 0.390173776540905, 0.39377862541005, 0.41481919772923,
+ 0.417203719727695, 0.448883372151561, 0.470826670527458, 0.473127941135317,
+ 0.500355682335794, 0.512959829065949, 0.543956416254501, 0.561118300072849,
+ 0.577319613657892, 0.590505231171846, 0.592247009277344, 0.602149887941778,
+ 0.603863363154233, 0.623027286957949, 0.688235442154109, 0.740739317568851,
+ 0.743342567193222, 0.746460048439813, 0.752596555235905, 0.792471590409887,
+ 0.794642869341088, 0.854227680653655, 0.875820309803661, 0.877021801231731,
+ 0.916470456826401, 0.941462716977012, 0.967143015330303, 1.00954419370265,
+ 1.0159810702373, 1.09290593672068, 1.11139687053384, 1.12096728982327,
+ 1.12445468544484, 1.12924358359679, 1.19186105388672, 1.30088131684334,
+ 1.3415423631248, 1.34419240891458, 1.44386141227057, 1.49093266381862,
+ 1.50939832102898, 1.57309243633392, 1.5840862540438, 1.64667304783108,
+ 1.74344241051915, 1.79478753722253, 1.92308481960537, 1.94059755458477,
+ 2.16018752554867, 2.97577757863932, 3.15856469071252, 3.17990924401169,
+ 3.50181933307077, 3.92315771656383, 4.47221537281699, 4.91030936099985,
+ 5.48181478836038, 5.67623017694721),
+ n.risk = c(96, 95, 94, 93,
+ 92, 91, 90, 89, 88, 87, 86, 85, 84, 83, 82, 81, 80, 79, 78, 77,
+ 76, 75, 74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 63, 62, 61,
+ 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45,
+ 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29,
+ 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13,
+ 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1),
+ n.event = c(1, 1, 0,
+ 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0,
+ 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1,
+ 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1,
+ 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0,
+ 0, 1, 1, 0, 0, 0, 0, 0, 1),
+ n.censor = c(0, 0, 1, 1, 0, 1, 0,
+ 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1,
+ 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0,
+ 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1,
+ 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1,
+ 1, 1, 1, 1, 0),
+ surv = c(0.989263204056545,
+ 0.978516304323273, 0.978516304323273, 0.978516304323273, 0.967527826035276,
+ 0.967527826035276, 0.956441726180785, 0.945348803653435, 0.934269566440011,
+ 0.934269566440011, 0.934269566440011, 0.92293425642073, 0.92293425642073,
+ 0.92293425642073, 0.92293425642073, 0.911167854431887, 0.911167854431887,
+ 0.899235443524779, 0.899235443524779, 0.899235443524779, 0.887019479777956,
+ 0.887019479777956, 0.874648697610133, 0.874648697610133, 0.862128637077813,
+ 0.849607778633846, 0.849607778633846, 0.849607778633846, 0.849607778633846,
+ 0.836570482616144, 0.836570482616144, 0.836570482616144, 0.836570482616144,
+ 0.836570482616144, 0.836570482616144, 0.836570482616144, 0.822166857718013,
+ 0.822166857718013, 0.807531588472591, 0.807531588472591, 0.79264203093055,
+ 0.777756282931739, 0.777756282931739, 0.762554609926638, 0.747371105288107,
+ 0.747371105288107, 0.747371105288107, 0.747371105288107, 0.731240535268028,
+ 0.715076965149398, 0.698902039504983, 0.68272489836221, 0.666577761088128,
+ 0.666577761088128, 0.650092208444908, 0.650092208444908, 0.650092208444908,
+ 0.632814187347251, 0.632814187347251, 0.632814187347251, 0.632814187347251,
+ 0.632814187347251, 0.632814187347251, 0.632814187347251, 0.632814187347251,
+ 0.611751463718706, 0.611751463718706, 0.611751463718706, 0.589204290906436,
+ 0.589204290906436, 0.565866423216841, 0.565866423216841, 0.541551101599392,
+ 0.517345256562477, 0.493207441124124, 0.493207441124124, 0.493207441124124,
+ 0.466462953525818, 0.439841802292848, 0.439841802292848, 0.439841802292848,
+ 0.410130748728784, 0.410130748728784, 0.410130748728784, 0.410130748728784,
+ 0.410130748728784, 0.410130748728784, 0.410130748728784, 0.36116684107885,
+ 0.312557263029036, 0.312557263029036, 0.312557263029036, 0.312557263029036,
+ 0.312557263029036, 0.312557263029036, 0.114470875549317),
+ cumhaz = c(0.0107948512630526,
+ 0.0217178297163296, 0.0217178297163296, 0.0217178297163296, 0.0330110937318073,
+ 0.0330110937318073, 0.044535415971588, 0.0562013151954944, 0.0679902673624694,
+ 0.0679902673624694, 0.0679902673624694, 0.0801972751626847, 0.0801972751626847,
+ 0.0801972751626847, 0.0801972751626847, 0.0930281457507738, 0.0930281457507738,
+ 0.106210383888236, 0.106210383888236, 0.106210383888236, 0.119888335495492,
+ 0.119888335495492, 0.13393296168815, 0.13393296168815, 0.148350788503535,
+ 0.162980472893901, 0.162980472893901, 0.162980472893901, 0.162980472893901,
+ 0.178444503149502, 0.178444503149502, 0.178444503149502, 0.178444503149502,
+ 0.178444503149502, 0.178444503149502, 0.178444503149502, 0.195811914601034,
+ 0.195811914601034, 0.213773105798942, 0.213773105798942, 0.232383570448668,
+ 0.251342064893415, 0.251342064893415, 0.271081153489211, 0.291193422982997,
+ 0.291193422982997, 0.291193422982997, 0.291193422982997, 0.313012823743045,
+ 0.33536509851569, 0.358244690339077, 0.381663284793446, 0.40559847536899,
+ 0.40559847536899, 0.430641067007499, 0.430641067007499, 0.430641067007499,
+ 0.457578442825448, 0.457578442825448, 0.457578442825448, 0.457578442825448,
+ 0.457578442825448, 0.457578442825448, 0.457578442825448, 0.457578442825448,
+ 0.491429183987586, 0.491429183987586, 0.491429183987586, 0.528982311822698,
+ 0.528982311822698, 0.569397230015335, 0.569397230015335, 0.613317846542379,
+ 0.659044819693656, 0.706825420368365, 0.706825420368365, 0.706825420368365,
+ 0.762576675512011, 0.821340156963448, 0.821340156963448, 0.821340156963448,
+ 0.891279270782872, 0.891279270782872, 0.891279270782872, 0.891279270782872,
+ 0.891279270782872, 0.891279270782872, 0.891279270782872, 1.01841526376066,
+ 1.16296758489947, 1.16296758489947, 1.16296758489947, 1.16296758489947,
+ 1.16296758489947, 1.16296758489947, 2.16743485035805),
+ std.err = c(0.0108125001325002,
+ 0.0154093640388869, 0.0154093640388869, 0.0154093640388869, 0.0191634362018011,
+ 0.0191634362018011, 0.0224313882575495, 0.0253671755680212, 0.0280648453735425,
+ 0.0280648453735425, 0.0280648453735425, 0.0307092324729056, 0.0307092324729056,
+ 0.0307092324729056, 0.0307092324729056, 0.0334061135217622, 0.0334061135217622,
+ 0.0360579337420049, 0.0360579337420049, 0.0360579337420049, 0.03871774105768,
+ 0.03871774105768, 0.0413506345532574, 0.0413506345532574, 0.0439628640337065,
+ 0.0465166967316136, 0.0465166967316136, 0.0465166967316136, 0.0465166967316136,
+ 0.0492068443844735, 0.0492068443844735, 0.0492068443844735, 0.0492068443844735,
+ 0.0492068443844735, 0.0492068443844735, 0.0492068443844735, 0.0524362425759656,
+ 0.0524362425759656, 0.0556987065758033, 0.0556987065758033, 0.0590197120954823,
+ 0.0623036830207577, 0.0623036830207577, 0.0657223491140423, 0.0691145142278423,
+ 0.0691145142278423, 0.0691145142278423, 0.0691145142278423, 0.0729203435886759,
+ 0.0767735746392451, 0.0806657607710129, 0.0846006061612719, 0.088550720081238,
+ 0.088550720081238, 0.0926827020934768, 0.0926827020934768, 0.0926827020934768,
+ 0.097225265717323, 0.097225265717323, 0.097225265717323, 0.097225265717323,
+ 0.097225265717323, 0.097225265717323, 0.097225265717323, 0.097225265717323,
+ 0.103825136842638, 0.103825136842638, 0.103825136842638, 0.111431387155048,
+ 0.111431387155048, 0.119630906618085, 0.119630906618085, 0.12872496782783,
+ 0.137903720147344, 0.147302655117893, 0.147302655117893, 0.147302655117893,
+ 0.159351333752637, 0.171769961340616, 0.171769961340616, 0.171769961340616,
+ 0.187404206167724, 0.187404206167724, 0.187404206167724, 0.187404206167724,
+ 0.187404206167724, 0.187404206167724, 0.187404206167724, 0.227492322295445,
+ 0.270207434823876, 0.270207434823876, 0.270207434823876, 0.270207434823876,
+ 0.270207434823876, 0.270207434823876, 1.04075302326175),
+ logse = TRUE,
+ lower = c(0.96851920944807, 0.949405169124627, 0.949405169124627,
+ 0.949405169124627, 0.93186181868013, 0.93186181868013, 0.915302992589028,
+ 0.899496545830592, 0.884266887010891, 0.884266887010891,
+ 0.884266887010891, 0.869022497757744, 0.869022497757744,
+ 0.869022497757744, 0.869022497757744, 0.853420467432664,
+ 0.853420467432664, 0.837878115999675, 0.837878115999675,
+ 0.837878115999675, 0.822198265971395, 0.822198265971395,
+ 0.80655860666671, 0.80655860666671, 0.79095324738692, 0.77557427567173,
+ 0.77557427567173, 0.77557427567173, 0.77557427567173, 0.759657089205706,
+ 0.759657089205706, 0.759657089205706, 0.759657089205706,
+ 0.759657089205706, 0.759657089205706, 0.759657089205706,
+ 0.741867174805844, 0.741867174805844, 0.724016885269468,
+ 0.724016885269468, 0.706056454691397, 0.688351934097492,
+ 0.688351934097492, 0.670390707705081, 0.652688441870397,
+ 0.652688441870397, 0.652688441870397, 0.652688441870397,
+ 0.633855625241047, 0.615181126579596, 0.596696496318155,
+ 0.578407046226118, 0.560371850125976, 0.560371850125976,
+ 0.542104861366039, 0.542104861366039, 0.542104861366039,
+ 0.523019535930763, 0.523019535930763, 0.523019535930763,
+ 0.523019535930763, 0.523019535930763, 0.523019535930763,
+ 0.523019535930763, 0.523019535930763, 0.499113021497944,
+ 0.499113021497944, 0.499113021497944, 0.473603967846141,
+ 0.473603967846141, 0.447593645969083, 0.447593645969083,
+ 0.420793044301925, 0.394817698059267, 0.369526316493336,
+ 0.369526316493336, 0.369526316493336, 0.341332040192808,
+ 0.314112808511793, 0.314112808511793, 0.314112808511793,
+ 0.284055747683803, 0.284055747683803, 0.284055747683803,
+ 0.284055747683803, 0.284055747683803, 0.284055747683803,
+ 0.284055747683803, 0.231241647486522, 0.18404694362627, 0.18404694362627,
+ 0.18404694362627, 0.18404694362627, 0.18404694362627, 0.18404694362627,
+ 0.0148869052793071),
+ upper = c(1, 1, 1, 1, 1, 1, 0.999429459956347,
+ 0.993538401799817, 0.987099749631646, 0.987099749631646,
+ 0.987099749631646, 0.980190551881825, 0.980190551881825,
+ 0.980190551881825, 0.980190551881825, 0.97282276513425, 0.97282276513425,
+ 0.965085932488443, 0.965085932488443, 0.965085932488443,
+ 0.956951127324476, 0.956951127324476, 0.94848698892779, 0.94848698892779,
+ 0.939708875746047, 0.930708250850577, 0.930708250850577,
+ 0.930708250850577, 0.930708250850577, 0.921271166068321,
+ 0.921271166068321, 0.921271166068321, 0.921271166068321,
+ 0.921271166068321, 0.921271166068321, 0.921271166068321,
+ 0.911158176134182, 0.911158176134182, 0.900679638346226,
+ 0.900679638346226, 0.889845826099439, 0.878772624403089,
+ 0.878772624403089, 0.867389011269196, 0.855788969418397,
+ 0.855788969418397, 0.855788969418397, 0.855788969418397,
+ 0.843587560204624, 0.831194332846806, 0.818613924898563,
+ 0.805856861331304, 0.792912619499665, 0.792912619499665,
+ 0.77959064675389, 0.77959064675389, 0.77959064675389, 0.765657433799898,
+ 0.765657433799898, 0.765657433799898, 0.765657433799898,
+ 0.765657433799898, 0.765657433799898, 0.765657433799898,
+ 0.765657433799898, 0.7498098371363, 0.7498098371363, 0.7498098371363,
+ 0.733021089332045, 0.733021089332045, 0.715391766187715,
+ 0.715391766187715, 0.696963981736076, 0.677897966082358,
+ 0.658284861247747, 0.658284861247747, 0.658284861247747,
+ 0.637466341832782, 0.615895964130853, 0.615895964130853,
+ 0.615895964130853, 0.592162744195105, 0.592162744195105,
+ 0.592162744195105, 0.592162744195105, 0.592162744195105,
+ 0.592162744195105, 0.592162744195105, 0.564091669959574,
+ 0.530799592470156, 0.530799592470156, 0.530799592470156,
+ 0.530799592470156, 0.530799592470156, 0.530799592470156,
+ 0.880208552629219),
+ conf.type = "log",
+ conf.int = 0.95,
+ call = as.call(survfit(formula = testfit1,
+ newdata = data.frame(X1 = 0.6, X2 = "A")))
+ ),
+ class = c("survfitcox", "survfit")
+)
+
+all.equal(`[.AsIs`(out_compare1, -14), `[.AsIs`(out1, -14))
+
+set.seed(35342)
+df2 <- data.frame(X1 = rpois(n = 300, lambda = 2),
+ X2 = c(rep("A", 10), sample(c("A", "B"), 290, replace = TRUE)),
+ X3 = sample(c("M", "F"), 300, replace = TRUE),
+ Ev = c(rep(1, 10), sample(c(0,1), 290, replace = TRUE)),
+ Time = 1:300)
+
+testfit2 <- coxph( Surv(Time, Ev) ~ X1 + strata(X2) + strata(X3), df2)
+
+out2 <- survfit(testfit2, newdata = data.frame(X1 = 2, X2 = "B", X3 = c("M", "F")))
+
+out_compare2 <- matrix(
+ c(53, 85, 53, 85, 53, 85, 53, 85, 53, 85, 53, 85, 25,
+ 56, 25, 56, 223, 226, 223, 226, 191, 194, 191, 194, NA, 257,
+ NA, 257),
+ nrow = 4,
+ ncol = 7,
+ dimnames = list(
+ c("B, F, 1", "B, M, 1",
+ "B, F, 2", "B, M, 2"),
+ c("records", "n.max", "n.start", "events", "median", "0.95LCL", "0.95UCL")
+ )
+)
+all.equal(out_compare2, survival:::survmean(out2, rmean = "none")$matrix)
+
+
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment