190 lines
6.7 KiB
Plaintext
190 lines
6.7 KiB
Plaintext
|
|
||
|
|
R Under development (unstable) (2024-04-17 r86441) -- "Unsuffered Consequences"
|
||
|
|
Copyright (C) 2024 The R Foundation for Statistical Computing
|
||
|
|
Platform: aarch64-unknown-linux-gnu
|
||
|
|
|
||
|
|
R is free software and comes with ABSOLUTELY NO WARRANTY.
|
||
|
|
You are welcome to redistribute it under certain conditions.
|
||
|
|
Type 'license()' or 'licence()' for distribution details.
|
||
|
|
|
||
|
|
R is a collaborative project with many contributors.
|
||
|
|
Type 'contributors()' for more information and
|
||
|
|
'citation()' on how to cite R or R packages in publications.
|
||
|
|
|
||
|
|
Type 'demo()' for some demos, 'help()' for on-line help, or
|
||
|
|
'help.start()' for an HTML browser interface to help.
|
||
|
|
Type 'q()' to quit R.
|
||
|
|
|
||
|
|
> library(survival)
|
||
|
|
> options(na.action=na.exclude) # preserve missings
|
||
|
|
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
|
||
|
|
>
|
||
|
|
> #
|
||
|
|
> # Simple tests of concordance
|
||
|
|
> #
|
||
|
|
> aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
|
||
|
|
>
|
||
|
|
> grank <- function(x, time, grp, wt)
|
||
|
|
+ unlist(tapply(x, grp, rank))
|
||
|
|
> grank2 <- function(x, time, grp, wt) { #for case weights
|
||
|
|
+ if (length(wt)==0) wt <- rep(1, length(x))
|
||
|
|
+ z <- double(length(x))
|
||
|
|
+ for (i in unique(grp)) {
|
||
|
|
+ indx <- which(grp==i)
|
||
|
|
+ temp <- tapply(wt[indx], x[indx], sum)
|
||
|
|
+ temp <- temp/2 + c(0, cumsum(temp)[-length(temp)])
|
||
|
|
+ z[indx] <- temp[match(x[indx], names(temp))]
|
||
|
|
+ }
|
||
|
|
+ z
|
||
|
|
+ }
|
||
|
|
>
|
||
|
|
> # Concordance by brute force. O(n^2) algorithm, but ok for n<500 or so
|
||
|
|
> allpair <- function(x, time, status, wt, all=FALSE) {
|
||
|
|
+ if (missing(wt)) wt <- rep(1, length(x))
|
||
|
|
+ count <- sapply(which(status==1), function(i) {
|
||
|
|
+ atrisk <- (time > time[i]) | (time==time[i] & status==0)
|
||
|
|
+ temp <- tapply(wt[atrisk], factor(sign(x[i] -x[atrisk]), c(1, -1, 0)),
|
||
|
|
+ sum)
|
||
|
|
+ wt[i]* c(ifelse(is.na(temp), 0, temp),
|
||
|
|
+ (sum(wt[time==time[i] & status==1]) - wt[i])/2)
|
||
|
|
+ })
|
||
|
|
+ rownames(count) <- c("concordant", "discordant", "tied.x", "tied.y")
|
||
|
|
+ if (all) {
|
||
|
|
+ colnames(count) <- time[status==1]
|
||
|
|
+ t(count)
|
||
|
|
+ }
|
||
|
|
+ else rowSums(count)
|
||
|
|
+ }
|
||
|
|
>
|
||
|
|
>
|
||
|
|
> # The std of C = std(numerator)/(number of comparable pairs)
|
||
|
|
> # The information matrix of a Cox model is = to the var(C-D)
|
||
|
|
> cfun <- function(fit) fit$cvar * sum(fit$count[1:3])^2
|
||
|
|
>
|
||
|
|
> tdata <- aml[aml$x=='Maintained', c("time", "status")]
|
||
|
|
> tdata$x <- c(1,6,2,7,3,7,3,8,4,4,5)
|
||
|
|
> tdata$wt <- c(1,2,3,2,1,2,3,4,3,2,1)
|
||
|
|
> fit <- concordance(Surv(time, status) ~x, tdata)
|
||
|
|
>
|
||
|
|
> aeq(fit$count[1:4], c(24,14,2,0))
|
||
|
|
[1] TRUE
|
||
|
|
> cfit <- coxph(Surv(time, status) ~ tt(x), tdata, tt=grank, method='breslow',
|
||
|
|
+ iter=0, x=T)
|
||
|
|
> cdt <- coxph.detail(cfit)
|
||
|
|
> aeq(sum(cdt$imat), cfun(fit))
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(sum(2*cdt$score), diff(fit$count[1:2]))
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(with(tdata, allpair(x, time, status)), c(14,24,2,0))
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> # Lots of ties
|
||
|
|
> tempy <- Surv(c(1,2,2,2,3,4,4,4,5,2), c(1,0,1,0,1,0,1,1,0,1))
|
||
|
|
> tempx <- c(5,5,4,4,3,3,7,6,5,4)
|
||
|
|
> fit2 <- concordance(tempy ~ tempx)
|
||
|
|
> addxy <- function(x) c(x[1:3], sum(x[4:5]))
|
||
|
|
> aeq(addxy(fit2$count), allpair(tempx, tempy[,1], tempy[,2]))
|
||
|
|
[1] TRUE
|
||
|
|
> cfit2 <- coxph(tempy ~ tt(tempx), tt=grank, method='breslow', iter=0)
|
||
|
|
> aeq(cfit2$var, 1/cfun(fit2))
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> # Direct call
|
||
|
|
> fit2b <- concordancefit(tempy, tempx)
|
||
|
|
> fit2c <- concordancefit(tempy, tempx, std.err=FALSE)
|
||
|
|
> all.equal(fit2[1:5], fit2b)
|
||
|
|
[1] TRUE
|
||
|
|
> all.equal(fit2b[1:3], fit2c)
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> # Bigger data
|
||
|
|
> fit3 <- concordance(Surv(time, status) ~ age, lung, reverse=TRUE)
|
||
|
|
> aeq(addxy(fit3$count), allpair(lung$age, lung$time, lung$status-1))
|
||
|
|
[1] TRUE
|
||
|
|
> cfit3 <- coxph(Surv(time, status) ~ tt(age), lung,
|
||
|
|
+ iter=0, method='breslow', tt=grank, x=T)
|
||
|
|
> cdt <- coxph.detail(cfit3)
|
||
|
|
> aeq(sum(cdt$imat), cfun(fit3))
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(2*sum(cdt$score), diff(fit3$count[2:1]))
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
>
|
||
|
|
> # More ties
|
||
|
|
> fit4 <- concordance(Surv(time, status) ~ ph.ecog, lung, reverse=TRUE)
|
||
|
|
> aeq(addxy(fit4$count), allpair(lung$ph.ecog, lung$time, lung$status-1))
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(fit4$count[1:5], c(8392, 4258, 7137, 21, 7))
|
||
|
|
[1] TRUE
|
||
|
|
> cfit4 <- coxph(Surv(time, status) ~ tt(ph.ecog), lung,
|
||
|
|
+ iter=0, method='breslow', tt=grank)
|
||
|
|
> aeq(1/cfit4$var, cfun(fit4))
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> # Case weights
|
||
|
|
> fit5 <- concordance(Surv(time, status) ~ x, tdata, weights=wt, reverse=TRUE)
|
||
|
|
> fit6 <- concordance(Surv(time, status) ~x, tdata[rep(1:11,tdata$wt),])
|
||
|
|
> aeq(addxy(fit5$count), with(tdata, allpair(x, time, status, wt)))
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(fit5$count[1:4], c(70, 91, 7, 0)) # checked by hand
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(fit5$count[1:3], fit6$count[c(2,1,3)]) #spurious "tied on time" values, ignore
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(fit5$std, fit6$std)
|
||
|
|
[1] TRUE
|
||
|
|
> cfit5 <- coxph(Surv(time, status) ~ tt(x), tdata, weights=wt,
|
||
|
|
+ iter=0, method='breslow', tt=grank2)
|
||
|
|
> cfit6 <- coxph(Surv(time, status) ~ tt(x), tdata[rep(1:11,tdata$wt),],
|
||
|
|
+ iter=0, method='breslow', tt=grank)
|
||
|
|
> aeq(1/cfit6$var, cfun(fit6))
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(cfit5$var, cfit6$var)
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> # Start, stop simplest cases
|
||
|
|
> fit7 <- concordance(Surv(rep(0,11), time, status) ~ x, tdata)
|
||
|
|
> aeq(fit7$count, fit$count)
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(fit7$std.err, fit$std.err)
|
||
|
|
[1] TRUE
|
||
|
|
> fit7 <- concordance(Surv(rep(0,11), time, status) ~ x, tdata, weights=wt)
|
||
|
|
> aeq(fit5$count, fit7$count[c(2,1,3:5)]) #one reversed, one not
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> # Multiple intervals for some, but same risk sets as tdata
|
||
|
|
> tdata2 <- data.frame(time1=c(0,3, 5, 6,7, 0, 4,17, 7, 0,16, 2, 0,
|
||
|
|
+ 0,9, 5),
|
||
|
|
+ time2=c(3,9, 13, 7,13, 18, 17,23, 28, 16,31, 34, 45,
|
||
|
|
+ 9,48, 60),
|
||
|
|
+ status=c(0,1, 1, 0,0, 1, 0,1, 0, 0,1, 1, 0, 0,1, 0),
|
||
|
|
+ x = c(1,1, 6, 2,2, 7, 3,3, 7, 3,3, 8, 4, 4,4, 5),
|
||
|
|
+ wt= c(1,1, 2, 3,3, 2, 1,1, 2, 3,3, 4, 3, 2,2, 1))
|
||
|
|
> fit8 <- concordance(Surv(time1, time2, status) ~x, tdata2, weights=wt,
|
||
|
|
+ reverse=TRUE)
|
||
|
|
> aeq(fit5$count, fit8$count)
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(fit5$std.err, fit8$std.err)
|
||
|
|
[1] TRUE
|
||
|
|
> cfit8 <- coxph(Surv(time1, time2, status) ~ tt(x), tdata2, weights=wt,
|
||
|
|
+ iter=0, method='breslow', tt=grank2)
|
||
|
|
> aeq(1/cfit8$var, cfun(fit8))
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> # Stratified
|
||
|
|
> tdata3 <- data.frame(time1=c(tdata2$time1, rep(0, nrow(lung))),
|
||
|
|
+ time2=c(tdata2$time2, lung$time),
|
||
|
|
+ status = c(tdata2$status, lung$status -1),
|
||
|
|
+ x = c(tdata2$x, lung$ph.ecog),
|
||
|
|
+ wt= c(tdata2$wt, rep(1, nrow(lung))),
|
||
|
|
+ grp=rep(1:2, c(nrow(tdata2), nrow(lung))))
|
||
|
|
> fit9 <- concordance(Surv(time1, time2, status) ~x + strata(grp),
|
||
|
|
+ data=tdata3, weights=wt, reverse=TRUE)
|
||
|
|
> aeq(fit9$count[1,], fit5$count)
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(fit9$count[2,], fit4$count)
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> proc.time()
|
||
|
|
user system elapsed
|
||
|
|
0.506 0.016 0.519
|