R Under development (unstable) (2024-03-19 r86153 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 
   2.51    0.25    2.75