Compare survival distributions

Computes a weighted version of the logrank test for comparing two or more survival distributions. The generalization to complex samples is based on the characterization of the logrank test as the score test in a Cox model. Under simple random sampling with replacement, this function with rho=0 and gamma=0 is almost identical to the robust score test in the survival package. The rho=0 and gamma=0 version was proposed by Rader (2014).

svylogrank(formula, design, rho=0,gamma=0,method=c("small","large","score"), ...)

Arguments

formula

Model formula with a single predictor. The predictor must be a factor if it has more than two levels.

design

A survey design object

rho,gamma

Coefficients for the Harrington/Fleming G-rho-gamma tests. The default is the logrank test, rho=1 gives a generalised Wilcoxon test

method

"small" works faster when a matrix with dimension number of events by number of people fits easily in memory; "large" works faster for large data sets; "score" works by brute-force construction of an expanded data set, and is for debugging

...

for future expansion.

Value

A vector containing the z-statistic for comparing each level of the variable to the lowest, the chisquared statistic for the logrank test, and the p-value.

References

Rader, Kevin Andrew. 2014. Methods for Analyzing Survival and Binary Data in Complex Surveys. Doctoral dissertation, Harvard University.http://nrs.harvard.edu/urn-3:HUL.InstRepos:12274283

See also

svykm, svycoxph.

Examples

library("survival")
data(nwtco)
## stratified on case status
dcchs<-twophase(id=list(~seqno,~seqno), strata=list(NULL,~rel),
         subset=~I(in.subcohort | rel), data=nwtco, method="simple")
svylogrank(Surv(edrel,rel)~factor(stage),design=dcchs)
#> [[1]]
#>         score       se        z            p
#> [1,] 20.36150 13.86013 1.469070 1.418138e-01
#> [2,] 37.69020 13.92394 2.706863 6.792220e-03
#> [3,] 56.47217 11.15524 5.062388 4.140369e-07
#> 
#> [[2]]
#>        Chisq            p 
#> 7.008976e+01 4.083552e-15 
#> 
#> attr(,"class")
#> [1] "svylogrank"

data(pbc, package="survival")
pbc$randomized <- with(pbc, !is.na(trt) & trt>0)
biasmodel<-glm(randomized~age*edema,data=pbc)
pbc$randprob<-fitted(biasmodel)
dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized))

svylogrank(Surv(time,status==2)~trt,design=dpbc)
#> [[1]]
#>          score                              
#> [1,] -1.442224 7.500245 -0.1922902 0.8475149
#> 
#> [[2]]
#>      Chisq          p 
#> 0.03697554 0.84751486 
#> 
#> attr(,"class")
#> [1] "svylogrank"

svylogrank(Surv(time,status==2)~trt,design=dpbc,rho=1)
#> [[1]]
#>            score                               
#> [1,] -0.08777618 5.908643 -0.01485556 0.9881474
#> 
#> [[2]]
#>        Chisq            p 
#> 0.0002206876 0.9881474163 
#> 
#> attr(,"class")
#> [1] "svylogrank"

rpbc<-as.svrepdesign(dpbc)
svylogrank(Surv(time,status==2)~trt,design=rpbc)
#> [[1]]
#>          score       se          z         p
#> [1,] -1.442224 7.500245 -0.1922902 0.8475149
#> 
#> [[2]]
#>      Chisq          p 
#> 0.03697554 0.84751486 
#> 
#> attr(,"class")
#> [1] "svylogrank"