Supplementary materials for "A Comparative Study of Two-Sample Tests for Interval-Censored Data" R code for the analysis of the breast cosmesis data


library ( survival ) library ( icenReg ) library ( flexsurv ) library ( interval ) library ( glrt ) library ( FHtest ) # # Analysis of the brest cosmesis data data ( cosmesis ) mydata = cosmesis mydata [ , 1 ][ mydata [ , 1 ]== 0 ] < -0 . 0 1 mydata [ , 2 ][ mydata [ , 2 ]== 0 ] < -0 . 0 1 mydata = data . frame ( mydata ) names ( mydata ) < -c ( " myl " , " myu " , " tr " ) # ## For a plot pdf ( " cosmesis . pdf " ) np _ fit = ic _ np ( cbind ( myl , myu )~tr , data = mydata ) plot ( np _ fit , col = c ( " Black " , " Blue " ) , xlim = c ( 0 , 1 0 0 ) , xlab = " Time " , ylab = " Estimated survival " , lwd = 2 ) dev . off () # #### preparation status = rep ( 3 , 9 4 ) newa = Surv ( mydata [ , 1 ] , mydata [ , 2 ] , status , type = " interval " ) newa 1 = Surv ( mydata [ 1 : 4 6 ,1 ] , mydata [ 1 : 4 6 ,2 ] , status [ 1 : 4 6 ] , type = " interval " ) newa 2 = Surv ( mydata [ 4 7 : 9 4 ,1 ] , mydata [ 4 7 : 9 4 ,2 ] , status [ 4 7 : 9 4 ] 
, type = " interval " ) # ### # Parametric methods # ### Exponential model fitting modelc <-survreg ( newa ~1 , dist = " exponential " ) model 1 <-survreg ( newa 1 ~1 , dist = " exponential " ) model 2 <-survreg ( newa 2 ~1 , dist = " exponential " ) myts = 2 *( model 1 $ loglik [ 1 ]+ model 2 $ loglik [ 1 ] -modelc $ loglik [ 1 ]) pvale = 1 -pchisq ( myts , 1 ) cat ( " pvale = " , pvale , " \ n " ) # ### Weibull model fitting modelc <-survreg ( newa ~1 , dist = " weibull " ) model 1 <-survreg ( newa 1 ~1 , dist = " weibull " ) model 2 <-survreg ( newa 2 ~1 , dist = " weibull " ) myts = 2 *( model 1 $ loglik [ 1 ]+ model 2 $ loglik [ 1 ] -modelc $ loglik [ 1 ]) pvalw = 1 -pchisq ( myts , 2 ) cat ( " pvalw = " , pvalw , " \ n " ) # ### Loglogistic model fitting modelc <-survreg ( newa ~1 , dist = " loglogistic " ) model 1 <-survreg ( newa 1 ~1 , dist = " loglogistic " ) model 2 <-survreg ( newa 2 ~1 , dist = " loglogistic " ) myts = 2 *( model 1 $ loglik [ 1 ]+ model 2 $ loglik [ 1 ] -modelc $ loglik [ 1 ]) pvalll = 1 -pchisq ( myts , 2 ) cat ( " pvalll = " , pvalll , " \ n " ) # ### Lognormal model fitting modelc <-survreg ( newa ~1 , dist = " lognormal " ) model 1 <-survreg ( newa 1 ~1 , dist = " lognormal " ) model 2 <-survreg ( newa 2 ~1 , dist = " lognormal " ) myts = 2 *( model 1 $ loglik [ 1 ]+ model 2 $ loglik [ 1 ] -modelc $ loglik [ 1 ]
) pvalln = 1 -pchisq ( myts , 2 ) cat ( " pvalln = " , pvalln , " \ n " ) # ### Generalized F modelc <-flexsurvreg ( newa ~1 , dist = " genf " , control = list ( fnscale = 2 5 0 )) model 1 <-flexsurvreg ( newa 1 ~1 , dist = " genf " , control = list ( fnscale = 2 5 0 )) model 2 <-flexsurvreg ( newa 2 ~1 , dist = " genf " , control = list ( fnscale = 2 5 0 )) myts = 2 *( model 1 $ loglik + model 2 $ loglik -modelc $ loglik ) pvalgf = 1 -pchisq ( myts , 4 ) cat ( " pvalgf = " , pvalgf , " \ n " ) generalization Wilcoxon -type ( Fay , 1 9 9 6 
# ####### # ####### # ####### Nonparametic methods # # Zhao and Sun mat <-matrix ( c ( cosmesis [ , 1 ] , cosmesis [ , 2 ] , cosmesis [ , 3 ]) , nrow = length ( cosmesis [ , 1 ])) p _ tem 1 = gLRT 1 ( mat ) p _ value 1 = as . vector ( p _ tem 1 $ p ) # # Sun df = data . frame ( " left " = cosmesis [ , 1 ] , " right " = cosmesis [ , 2 ] , " Treat " = cosmesis [ , 3 ]) icout = ictest ( Surv ( left , right , type = ' interval 2 ' )~Treat , data = df ) fit 3 = ictest ( Surv ( left , right , type = ' interval 2 ' )~Treat , data = df , scores = ' logrank 1 ') p _ value 2 = fit 3 $ p . value # # Finkelstein ' s fit 4 = ictest ( Surv ( left , right , type = ' interval 2 ' )~Treat , data = df , initfit = icout $ fit , scores = ' logrank 2 ') # Finkelstein ' s p _ value 3 = fit 4 $ p . value # # Peto & Peto
*( d 1 + d 2 )/( y 1 + y 2 ) V _ tem 1 [ is . nan ( V _ tem 1 )] <-0 V _ tem 2 [ is . nan ( V _ tem 2 )] <-0 out = list ( sum ( G _ tem ) , sum ( V _ tem 1 )+ sum ( V _ tem 2 )) names ( out ) < -c ( " G " , " V " ) return ( out ) } # # sfunction < -function ( left ,
_ dj ))) new _ s = c ( 1 , cumprod ( 1 -turn _ dj / turn _ yj )) new _ s [ which ( is . na ( new _ s )== TRUE )]= 0 turn _ err = sum ( abs (( new _ s [ -len _ tt ] -old _ s [ -len _ tt ]))) } new _ s [ len _ tt ]= 0 . 5 * new _ s [ len _ tt -1 ] return ( new _ s ) } # # M = 1 0 # the number of imputations rho 0 = 0 rho 1 = 1 rho 5 = 0 . 5 cosmesis 0 <-cosmesis [ cosmesis [ , 3 ]== 0 , 1 : 2 ] cosmesis 1 <-cosmesis [ cosmesis [ , 3 ]== 1 , 1 : 2 ] li 1 = cosmesis 0 [ , 1 ] ri 1 = cosmesis 0 [ , 2 ] li 2 = cosmesis 1 [ , 1 ] ri 2 = cosmesis 1 [ , 2 ] Li = c ( li 1 , li 2 ) Ri = c ( ri 1 , ri 2 ) nh 1 = nrow ( cosmesis 0 ) nh 2 = nrow ( cosmesis 1 ) # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # G 0 = V 0 = G 1 = V 1 = G 5 = V 5 = rep ( 0 ,M ) rnoncensored 1 = which ( ri 1 <1 0 0 ) rnoncensored 2 = which ( ri 2 <1 0 0 ) nrnonc 1 = length ( rnoncensored 1 ) nrnonc 2 = length ( rnoncensored 2 ) for ( m in 1 : M ){ # bootstrap from each (U , V ) df 1 = data . frame ( li 1 , ri 1 ) df 1 = df 1 [ sample (( 1 : nrow ( df 1 )) ,( nh 1 ) , replace = TRUE ) , ] df 2 = data . frame ( li 2 , ri 2 ) df 2 = df 2 [ sample (( 1 : nrow ( df 2 )) ,( nh 2 ) , replace = TRUE ) , ] # # newli 1 = df 1 [ , " li 1 " ] newri 1 = df 1 [ , " ri 1 " ] newli 2 = df 2 [ , " li 2 " ] newri 2 = df 2 [ , " ri 2 " ] # mydata = data . frame ( newli 1 , newri 1 ) # s 1 = sfunction ( newli 1 , newri 1 ) s 2 = sfunction ( newli 2 , newri 2 ) # L 1 = sort ( unique ( c ( newli 1 , newri 1 ))) L 2 = sort ( unique ( c ( newli 2 , newri 2 ))) # x 1 = li 1 x 2 = li 2 statu 1 = rep ( 0 ,( nh 1 )) statu 1 [ rnoncensored 1 ]= 1 # statu 2 = rep ( 0 ,( nh 2 )) statu 2 [ rnoncensored 2 ]= 1 # # if ( min ( L 1 ) > 0 ){ ax = c ( 0 ,L 1 



  ) L = with ( df , left ) R = with ( df , right ) trt = with ( df , Treat ) fit 5 = ictest (L ,R , trt , initfit = icout $ fit , scores = ' wmw ') p _ value 4 = fit 5 $ p . value # # Huang et al fit 6 = ictest ( Surv ( left , right , type = ' interval 2 ' )~Treat , data = df , initfit = icout $ fit , method = ' wsr . HLY ' , mcontrol = mControl ( nwsr = 9 9 ) , scores = ' logrank 1 ') p _ value 5 = fit 6 $ p . value # # Oller and Gomez fit 7 = FHtesticp ( Surv ( left , right , type = ' interval 2 ' )~Treat , data = df , rho = 0 , lambda = 0 ) p _ value 6 = fit 7 $ pvalue # fit 8 = FHtesticp ( Surv ( left , right , type = ' interval 2 ' )~Treat , data = df , rho = 1 , lambda = 0 ) # early diff p _ value 7 = fit 8 $ pvalue # fit 9 = FHtesticp ( Surv ( left , right , type = ' interval 2 ' )~Treat , data = df , rho = 0 , lambda = 1 ) # middle or late diff p _ value 8 = fit 9 $ pvalue # fit 1 0 = FHtesticp ( Surv ( left , right , type = ' interval 2 ' )~Treat , data = df , rho = 1 , lambda = 1 ) # middle or late diff p _ value 9 = fit 1 0 $ pvalue # # Pan ' s method # The following two functions are needed for Pan ' s approach pan = function ( rho ){ G _ tem = su ^rho *( y 1 * y 2 /( y 1 + y 2 ))*( d 1 / y 1 -d 2 / y 2 ) G _ tem [ is . nan ( G _ tem )] 




  = sort ( unique ( c ( r 6 _ left , r 6 _ right ))) len _ tt = length ( turn _ tau ) survobj = survfit ( Surv ( turn _ tau [ 1 :( len _ tt -1 )] , rep ( 1 , len _ tt -1 ))~1 ) init _ s = c ( 1 , survobj $ surv ) turn _ p = -diff ( init _ s ) # turn _ tau 1 2 = cbind ( turn _ tau [ -len _ tt ] , turn _ tau [ -1 ]) alpha _ fn = function ( xx , infi , supr ) ifelse ( xx [ 1 ] >= infi & xx [ 2 ] <= supr , 1 , 0 ) alpha _ mat = apply ( turn _ tau 1 2 ,1 , alpha _ fn , infi = left , supr = right 1 ) # id _ zero = which ( apply ( alpha _ mat == 0 ,1 , all )) if ( length ( id _ zero ) > 0 ) alpha _ mat = alpha _ mat [ -id _ zero ,= new _ sturn _ p = -diff ( old _ s ) dj _ temp = t ( t ( alpha _ mat )* turn _ p )/ as . vector ( alpha _ mat \%*\% turn _ p ) turn _ dj = apply ( dj _ temp , 2 , sum ) turn _ yj = rev ( cumsum ( rev ( turn 




  ); ay = c ( 1 ,s 1 )} else { ax = L 1 ; ay = s 1 } # lpoint = findInterval ( li 1 [ rnoncensored 1 ] , ax ); rpoint = findInterval ( ri 1 [ rnoncensored 1 ] , ax ) deno = ay [ lpoint ] -ay [ rpoint + 1 ] p = runif ( nrnonc 1 , 0 , 1 ) ry = ay [ lpoint ] -deno * p myindex 1 = findInterval ( -ry , -ay ) myslope =( ax [ myindex 1 ] -ax [ myindex 1 + 1 ])/( ay [ myindex 1 ] -ay [ myindex 1 + 1 ]) myintercept = ax [ myindex 1 ] -ay [ myindex 1 ]* myslope x 1 [ rnoncensored 1 ]= myslope * ry + myintercept # # For the second group if ( min ( L 2 ) > 0 ){ ax = c ( 0 ,L 2 ); ay = c ( 1 ,s 2 )} else { ax = L 2 ; ay = s 2 } # lpoint = findInterval ( li 2 [ rnoncensored 2 ] , ax ); rpoint = findInterval ( ri 2 [ rnoncensored 2 ] , ax ) deno = ay [ lpoint ] -ay [ rpoint + 1 ] p = runif ( nrnonc 2 , 0 , 1 ) ry = ay [ lpoint ] -deno * p myindex 2 = findInterval ( -ry , -ay ) myslope =( ax [ myindex 2 ] -ax [ myindex 2 + 1 ])/( ay [ myindex 2 ] -ay [ myindex 2 + 1 ]) myintercept = ax [ myindex 2 ] -ay [ myindex 2 ]* myslope x 2 [ rnoncensored 2 ]= myslope * ry + myintercept # # x 1 = round ( x 1 ,5 ) x 2 = round ( x 2 ,5 ) out 1 = survfit ( Surv ( x 1 , statu 1 )~1 ) out 2 = survfit ( Surv ( x 2 , statu 2 )~1 ) outcomb = survfit ( Surv ( c ( x 1 , x 2 ) , c ( statu 1 , statu 2 ))~1 ) y 1 = d 1 = rep ( 0 , length ( outcomb $ time )) y 2 = d 2 = rep ( 0 , length ( outcomb $ time )) d 1 [ match ( out 1 $ time , outcomb $ time )]= out 1 $ n . event d 2 [ match ( out 2 $ time , outcomb $ time )]= out 2 $ n . event y 1 [ match ( out 1 $ time , outcomb $ time )]= out 1 $ n . risk y 2 [ match ( out 2 $ time , outcomb $ time )]= out 2 $ n . risk ind 1 = which ( y 1 == 0 ) ind 2 = which ( y 2 == 0 ) # ny 1 = y 1 ny 1 [ ind 1 ]= sapply ( 1 : length ( ind 1 ) , function( k ) { temp 1 = y 1 [ -c ( 1 : ind 1 [ k ])]; temp 1 1 = temp 1 [ temp 1 >0 ]; if ( length ( temp 1 1 ) > 0 ) { y 1 [ ind 1 [ k ]] < -max ( temp 1 1 )} else { y 1 [ ind 1 [ k ]]= ind 2 ]= sapply ( 1 : length ( ind 2 ) , function ( k ) { temp 2 = y 2 [ -c ( 1 : ind 2 [ k ])]; temp 2 1 = temp 2 [ temp 2 >0 ]; if ( length ( temp 2 1 ) > 0 ) { y 2 [ ind 2 [ k ]] < -max ( temp 2 1 )} else { y 2 [ ind 2 [ k ]]= 0 } }) y 2 = unlist ( ny 2 ) # su = c ( 1 , outcomb $ surv ) su = su [ -length ( su )] # out _ pan = pan ( 0 ) G 0 [ m ]= out _ pan $ G V 0 [ m ]= out _ pan $ sum ( G 0 )/ M finalv 0 = sum ( V 0 )/ M +( 1 + 1 / M )* var ( G 0 ) finalG 1 = sum ( G 1 )/ M finalv 1 = sum ( V 1 )/ M +( 1 + 1 / M )* var ( G 1 ) finalG 5 = sum ( G 5 )/ M finalv 5 = sum ( V 5 )/ M +( 1 + 1 / M )* var ( G 5 ) # p _ value 1 0 = 2 *pnorm ( abs ( finalG 0 ) , mean = 0 , sd = sqrt ( finalv 0 ) , lower . tail = FALSE ) # Pan ( 0 ) p _ value 1 1 = 2 * pnorm ( abs ( finalG 5 ) , mean = 0 , sd = sqrt ( finalv 5 ) , lower . tail = FALSE ) # Pan ( 0 . 5 ) p _ value 1 2 = 2 * pnorm ( abs ( finalG 1 ) , mean = 0 , sd = sqrt ( finalv 1 ) , lower . tail = FALSE ) # Pan ( 1 ) # 


For the scenario where assessment time distributions were different for the two groups, the following code was used to get the intervals (L, R] for group 1 and group 2:
Supplementary materials for "A Comparative Study of Two-Sample Tests for Interval-Censored Data" R code for the analysis of the breast cosmesis data