R Under development (unstable) (2024-03-08 r86070 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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(survey) Loading required package: grid Loading required package: Matrix Loading required package: survival Attaching package: 'survey' The following object is masked from 'package:graphics': dotchart > example(svycoxph, ask=FALSE) svycxp> ## Somewhat unrealistic example of nonresponse bias. svycxp> data(pbc, package="survival") svycxp> pbc$randomized<-with(pbc, !is.na(trt) & trt>0) svycxp> biasmodel<-glm(randomized~age*edema,data=pbc,family=binomial) svycxp> pbc$randprob<-fitted(biasmodel) svycxp> if (is.null(pbc$albumin)) pbc$albumin<-pbc$alb ##pre2.9.0 svycxp> dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized)) svycxp> rpbc<-as.svrepdesign(dpbc) svycxp> (model<-svycoxph(Surv(time,status>0)~log(bili)+protime+albumin,design=dpbc)) Call: svycoxph(formula = Surv(time, status > 0) ~ log(bili) + protime + albumin, design = dpbc) coef exp(coef) se(coef) robust se z p log(bili) 0.88592 2.42522 0.09140 0.09048 9.791 < 2e-16 protime 0.24487 1.27745 0.07825 0.08122 3.015 0.00257 albumin -1.04298 0.35240 0.21211 0.20454 -5.099 3.41e-07 Likelihood ratio test= on 3 df, p= n= 312, number of events= 144 svycxp> svycoxph(Surv(time,status>0)~log(bili)+protime+albumin,design=rpbc) Call: svycoxph.svyrep.design(formula = Surv(time, status > 0) ~ log(bili) + protime + albumin, design = rpbc) coef exp(coef) se(coef) z p log(bili) 0.88592 2.42522 0.09838 9.005 < 2e-16 protime 0.24487 1.27745 0.09373 2.612 0.00899 albumin -1.04298 0.35240 0.21966 -4.748 2.05e-06 Likelihood ratio test=NA on 3 df, p=NA n= 312, number of events= 144 svycxp> s<-predict(model,se=TRUE, type="curve", svycxp+ newdata=data.frame(bili=c(3,9), protime=c(10,10), albumin=c(3.5,3.5))) svycxp> plot(s[[1]],ci=TRUE,col="sienna") svycxp> lines(s[[2]], ci=TRUE,col="royalblue") svycxp> quantile(s[[1]], ci=TRUE) 0.75 0.5 0.25 1435 2503 3574 attr(,"ci") 0.025 0.975 0.75 1217 1786 0.5 2256 3170 0.25 3222 Inf svycxp> confint(s[[2]], parm=365*(1:5)) 0.025 0.975 365 0.8375139 0.9453781 730 0.7382750 0.8999016 1095 0.4784105 0.7478460 1460 0.3192009 0.6206764 1825 0.2149475 0.5292978 > m<-update(model, .~.+I(protime^2)) > a<-anova(m,model) > b<-anova(m, model,force=TRUE) > stopifnot(isTRUE(all.equal(b[2:6],a[c(3,4,6,7,8)]))) > > proc.time() user system elapsed 1.42 0.21 1.64