Summarizing the result of the product limit method in life-table format. Calculates the number of subjects at risk and counts events and censored observations at specified times or in specified time intervals.
# S3 method for class 'prodlim'
summary(
object,
times,
newdata,
max.tables = 20,
surv = TRUE,
cause,
intervals = FALSE,
percent = FALSE,
format = "df",
...
)An object with class `prodlim' derived with
prodlim
Vector of times at which to return the estimated probabilities.
A data frame with the same variable names as those
that appear on the right hand side of the 'prodlim' formula.
Defaults to object$X.
Integer. If newdata is not given the value
of max.tables decides about the maximal number of tables to
be shown. Defaults to 20.
Logical. If FALSE report event probabilities instead of
survival probabilities. Only available for
object$model=="survival".
For competing risk models. The event of interest for which predictions of the absolute risks are obtained by evaluating the cause-specific cumulative
incidence functions at times.
Logical. If TRUE count events and censored in
intervals between the values of times.
Logical. If TRUE all estimated values are multiplied by 100 and thus interpretable on a percent scale.
Control format of output. Since May 2021,
the result is a data.table and data.frame with attributes. When there are multiple
covariate strata or competing risks, these are indicated by columns.
Set format to "list" to get the old behaviour.
Further arguments that are passed to the print function.
A data.frame with the relevant information.
For cluster-correlated data the number of clusters at-risk are are also given. Confidence intervals are displayed when they are part of the fitted object.
library(lava)
set.seed(17)
m <- survModel()
distribution(m,~age) <- uniform.lvm(30,80)
distribution(m,~sex) <- binomial.lvm()
m <- categorical(m,~z,K=3)
regression(m,eventtime~age) <- 0.01
regression(m,eventtime~sex) <- -0.4
d <- sim(m,50)
d$sex <- factor(d$sex,levels=c(0,1),labels=c("female","male"))
d$Z <- factor(d$z,levels=c(1,0,2),labels=c("B","A","C"))
# Univariate Kaplan-Meier
# -----------------------------------------------------------------------------------------
fit0 <- prodlim(Hist(time,event)~1,data=d)
summary(fit0)
#> time n.risk n.event n.lost surv se.surv lower upper
#> 1 1.25 50 1 0 0.9800 0.0198 0.94119 1.000
#> 2 1.28 49 0 1 0.9800 0.0198 0.94119 1.000
#> 3 1.84 48 1 0 0.9596 0.0280 0.90470 1.000
#> 4 1.85 47 1 0 0.9392 0.0340 0.87244 1.000
#> 5 1.86 46 0 1 0.9392 0.0340 0.87244 1.000
#> 6 2.15 45 0 1 0.9392 0.0340 0.87244 1.000
#> 7 2.30 44 1 0 0.9178 0.0394 0.84061 0.995
#> 8 2.37 43 0 1 0.9178 0.0394 0.84061 0.995
#> 9 2.55 42 1 0 0.8960 0.0441 0.80952 0.982
#> 10 2.61 41 0 1 0.8960 0.0441 0.80952 0.982
#> 11 2.85 40 1 0 0.8736 0.0484 0.77879 0.968
#> 12 2.99 39 1 0 0.8512 0.0520 0.74916 0.953
#> 13 3.20 38 0 1 0.8512 0.0520 0.74916 0.953
#> 14 3.37 37 0 1 0.8512 0.0520 0.74916 0.953
#> 15 3.52 36 1 0 0.8275 0.0557 0.71833 0.937
#> 16 3.64 35 0 1 0.8275 0.0557 0.71833 0.937
#> 17 4.21 34 1 0 0.8032 0.0592 0.68725 0.919
#> 18 4.44 33 0 1 0.8032 0.0592 0.68725 0.919
#> 19 4.55 32 1 0 0.7781 0.0624 0.65578 0.900
#> 20 4.56 31 1 0 0.7530 0.0652 0.62512 0.881
#> 21 4.61 30 1 0 0.7279 0.0677 0.59515 0.861
#> 22 4.99 29 1 0 0.7028 0.0699 0.56582 0.840
#> 23 5.08 28 1 0 0.6777 0.0718 0.53705 0.818
#> 24 5.45 27 1 0 0.6526 0.0734 0.50881 0.796
#> 25 5.75 26 1 0 0.6275 0.0747 0.48107 0.774
#> 26 5.75 25 1 0 0.6024 0.0758 0.45379 0.751
#> 27 6.06 24 0 1 0.6024 0.0758 0.45379 0.751
#> 28 6.29 23 0 1 0.6024 0.0758 0.45379 0.751
#> 29 6.36 22 0 1 0.6024 0.0758 0.45379 0.751
#> 30 6.42 21 1 0 0.5737 0.0774 0.42192 0.725
#> 31 6.51 20 1 0 0.5450 0.0787 0.39076 0.699
#> 32 6.79 19 0 1 0.5450 0.0787 0.39076 0.699
#> 33 7.01 18 1 0 0.5147 0.0799 0.35805 0.671
#> 34 7.71 17 0 1 0.5147 0.0799 0.35805 0.671
#> 35 8.12 16 0 1 0.5147 0.0799 0.35805 0.671
#> 36 8.44 15 0 1 0.5147 0.0799 0.35805 0.671
#> 37 8.56 14 0 1 0.5147 0.0799 0.35805 0.671
#> 38 9.07 13 0 1 0.5147 0.0799 0.35805 0.671
#> 39 9.30 12 0 1 0.5147 0.0799 0.35805 0.671
#> 40 9.75 11 1 0 0.4679 0.0853 0.30080 0.635
#> 41 9.94 10 1 0 0.4212 0.0887 0.24737 0.595
#> 42 10.11 9 1 0 0.3744 0.0903 0.19733 0.551
#> 43 10.25 8 1 0 0.3276 0.0903 0.15049 0.505
#> 44 10.88 7 1 0 0.2808 0.0887 0.10686 0.455
#> 45 11.20 6 0 1 0.2808 0.0887 0.10686 0.455
#> 46 11.49 5 1 0 0.2246 0.0870 0.05418 0.395
#> 47 12.05 4 1 0 0.1685 0.0814 0.00901 0.328
#> 48 12.14 3 1 0 0.1123 0.0710 0.00000 0.252
#> 49 13.98 2 1 0 0.0562 0.0533 0.00000 0.161
#> 50 15.86 1 0 1 0.0562 0.0533 0.00000 0.161
## show survival probabilities as percentage and
## count number of events within intervals of a
## given time-grid:
summary(fit0,times=c(1,5,10,12),percent=TRUE,intervals=TRUE)
#> time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 0 1 50 0 0 100 0 100 100
#> 2 1 5 50 13 9 70 7 57 84
#> 3 5 10 28 9 10 42 9 25 59
#> 4 10 12 9 4 1 22 9 5 40
## the result of summary has a print function
## which passes ... to print and print.listof
sx <- summary(fit0,times=c(1,5,10,12),percent=TRUE,intervals=TRUE)
print(sx,digits=3)
#> time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 0 1 50 0 0 100.0 0.00 100.00 100.0
#> 2 1 5 50 13 9 70.3 6.99 56.58 84.0
#> 3 5 10 28 9 10 42.1 8.87 24.74 59.5
#> 4 10 12 9 4 1 22.5 8.70 5.42 39.5
## show absolute risks, i.e., cumulative incidences (1-survival)
summary(fit0,times=c(1,5,10,12),surv=FALSE,percent=TRUE,intervals=TRUE)
#> time0 time1 n.risk n.event n.lost cuminc se.cuminc lower upper
#> 1 0 1 50 0 0 0 0 100 100
#> 2 1 5 50 13 9 30 7 16 43
#> 3 5 10 28 9 10 58 9 41 75
#> 4 10 12 9 4 1 78 9 60 95
# Stratified Kaplan-Meier
# -----------------------------------------------------------------------------------------
fit1 <- prodlim(Hist(time,event)~sex,data=d)
print(summary(fit1,times=c(1,5,10),intervals=TRUE,percent=TRUE),digits=3)
#> sex time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 female 0 1 24 0 0 100.0 0.00 100.00 100.0
#> 2 female 1 5 24 8 3 63.5 10.37 43.18 83.8
#> 3 female 5 10 13 6 4 28.3 11.59 5.56 51.0
#> 4 male 0 1 26 0 0 100.0 0.00 100.00 100.0
#> 5 male 1 5 26 5 6 77.7 8.92 60.23 95.2
#> 6 male 5 10 15 3 6 56.7 12.68 31.83 81.5
# old behaviour
print(summary(fit1,times=c(1,5,10),intervals=TRUE,percent=TRUE,format="list"),digits=3)
#> sex=female :
#> time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 0 1 24 0 0 100.0 0.0 100.00 100.0
#> 2 1 5 24 8 3 63.5 10.4 43.18 83.8
#> 3 5 10 13 6 4 28.3 11.6 5.56 51.0
#>
#> sex=male :
#> time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 0 1 26 0 0 100.0 0.00 100.0 100.0
#> 2 1 5 26 5 6 77.7 8.92 60.2 95.2
#> 3 5 10 15 3 6 56.7 12.68 31.8 81.5
#>
summary(fit1,times=c(1,5,10),intervals=TRUE,percent=TRUE)
#> sex time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 female 0 1 24 0 0 100 0 100 100
#> 2 female 1 5 24 8 3 64 10 43 84
#> 3 female 5 10 13 6 4 28 12 6 51
#> 4 male 0 1 26 0 0 100 0 100 100
#> 5 male 1 5 26 5 6 78 9 60 95
#> 6 male 5 10 15 3 6 57 13 32 82
fit2 <- prodlim(Hist(time,event)~Z,data=d)
print(summary(fit2,times=c(1,5,10),intervals=TRUE,percent=TRUE),digits=3)
#> Z time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 B 0 1 17 0 0 100.0 0.0 100.00 100.0
#> 2 B 1 5 17 6 2 61.1 12.5 36.54 85.6
#> 3 B 5 10 9 1 3 54.3 12.8 29.13 79.5
#> 4 A 0 1 16 0 0 100.0 0.0 100.00 100.0
#> 5 A 1 5 16 3 3 76.9 11.7 54.02 99.8
#> 6 A 5 10 10 4 3 39.6 15.6 9.01 70.1
#> 7 C 0 1 17 0 0 100.0 0.0 100.00 100.0
#> 8 C 1 5 17 4 4 73.9 11.3 51.89 96.0
#> 9 C 5 10 9 4 4 23.5 18.1 0.00 58.9
## Continuous strata (Beran estimator)
# -----------------------------------------------------------------------------------------
fit3 <- prodlim(Hist(time,event)~age,data=d)
print(summary(fit3,
times=c(1,5,10),
newdata=data.frame(age=c(20,50,70)),
intervals=TRUE,
percent=TRUE),digits=3)
#> age time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 20 0 1 13 0 0 100.0 0.00 100.0 100.0
#> 2 20 1 5 13 1 2 90.9 8.67 73.9 100.0
#> 3 20 5 10 10 2 4 70.1 14.71 41.3 99.0
#> 4 50 0 1 25 0 0 100.0 0.00 100.0 100.0
#> 5 50 1 5 25 3 6 86.4 7.39 71.9 100.0
#> 6 50 5 10 16 5 7 45.8 15.11 16.2 75.4
#> 7 70 0 1 20 0 0 100.0 0.00 100.0 100.0
#> 8 70 1 5 20 9 2 52.2 11.67 29.4 75.1
#> 9 70 5 10 9 3 2 33.9 11.48 11.4 56.4
## stratified Beran estimator
# -----------------------------------------------------------------------------------------
fit4 <- prodlim(Hist(time,event)~age+sex,data=d)
print(summary(fit4,
times=c(1,5,10),
newdata=data.frame(age=c(20,50,70),sex=c("female","male","male")),
intervals=TRUE,
percent=TRUE),digits=3)
#> age sex time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 20 female 0 1 8 0 0 100.0 0.00 100.0 100.0
#> 2 20 female 1 5 8 1 1 85.7 13.23 59.8 100.0
#> 3 20 female 5 10 6 2 4 NA NA NA NA
#> 4 50 male 0 1 15 0 0 100.0 0.00 100.0 100.0
#> 5 50 male 1 5 15 1 5 93.3 6.44 80.7 100.0
#> 6 50 male 5 10 9 2 3 65.3 17.61 30.8 99.8
#> 7 70 male 0 1 10 0 0 100.0 0.00 100.0 100.0
#> 8 70 male 1 5 10 4 2 51.4 17.69 16.8 86.1
#> 9 70 male 5 10 4 1 2 34.3 18.30 0.0 70.2
print(summary(fit4,
times=c(1,5,10),
newdata=data.frame(age=c(20,50,70),sex=c("female","male","male")),
intervals=TRUE,
percent=TRUE),digits=3)
#> age sex time0 time1 n.risk n.event n.lost surv se.surv lower upper
#> 1 20 female 0 1 8 0 0 100.0 0.00 100.0 100.0
#> 2 20 female 1 5 8 1 1 85.7 13.23 59.8 100.0
#> 3 20 female 5 10 6 2 4 NA NA NA NA
#> 4 50 male 0 1 15 0 0 100.0 0.00 100.0 100.0
#> 5 50 male 1 5 15 1 5 93.3 6.44 80.7 100.0
#> 6 50 male 5 10 9 2 3 65.3 17.61 30.8 99.8
#> 7 70 male 0 1 10 0 0 100.0 0.00 100.0 100.0
#> 8 70 male 1 5 10 4 2 51.4 17.69 16.8 86.1
#> 9 70 male 5 10 4 1 2 34.3 18.30 0.0 70.2
## assess results from summary
x <- summary(fit4,times=10,newdata=expand.grid(age=c(60,40,50),sex=c("male","female")))
cbind(names(x$table),do.call("rbind",lapply(x$table,round,2)))
#> NULL
x <- summary(fit4,times=10,newdata=expand.grid(age=c(60,40,50),sex=c("male","female")))
## Competing risks: Aalen-Johansen
# -----------------------------------------------------------------------------------------
d <- SimCompRisk(30)
crfit <- prodlim(Hist(time,event)~X1,data=d)
summary(crfit,times=c(1,2,5))
#> X1 time cause n.risk n.event n.lost cuminc se.cuminc lower upper
#> 1 0 1 1 14 0 0 0.1176 0.0781 0.00000 0.271
#> 2 0 2 1 14 0 0 0.1176 0.0781 0.00000 0.271
#> 3 0 5 1 7 0 0 0.2353 0.1029 0.03365 0.437
#> 4 1 1 1 13 0 0 0.0000 0.0000 0.00000 0.000
#> 5 1 2 1 12 0 0 0.0000 0.0000 0.00000 0.000
#> 6 1 5 1 6 0 0 0.2593 0.1295 0.00548 0.513
#> 7 0 1 2 14 0 0 0.0588 0.0571 0.00000 0.171
#> 8 0 2 2 14 0 0 0.0588 0.0571 0.00000 0.171
#> 9 0 5 2 7 0 0 0.3072 0.1148 0.08215 0.532
#> 10 1 1 2 13 0 0 0.0000 0.0000 0.00000 0.000
#> 11 1 2 2 12 0 0 0.0000 0.0000 0.00000 0.000
#> 12 1 5 2 6 0 0 0.0926 0.0881 0.00000 0.265
summary(crfit,times=c(1,2,5),cause=1,intervals=TRUE)
#> X1 time0 time1 cause n.risk n.event n.lost cuminc se.cuminc lower upper
#> 1 0 0 1 1 17 2 0 0.118 0.0781 0.00000 0.271
#> 2 0 1 2 1 14 0 0 0.118 0.0781 0.00000 0.271
#> 3 0 2 5 1 14 2 1 0.235 0.1029 0.03365 0.437
#> 4 1 0 1 1 13 0 0 0.000 0.0000 0.00000 0.000
#> 5 1 1 2 1 13 0 1 0.000 0.0000 0.00000 0.000
#> 6 1 2 5 1 12 3 2 0.259 0.1295 0.00548 0.513
summary(crfit,times=c(1,2,5),cause=1)
#> X1 time cause n.risk n.event n.lost cuminc se.cuminc lower upper
#> 1 0 1 1 14 0 0 0.118 0.0781 0.00000 0.271
#> 2 0 2 1 14 0 0 0.118 0.0781 0.00000 0.271
#> 3 0 5 1 7 0 0 0.235 0.1029 0.03365 0.437
#> 4 1 1 1 13 0 0 0.000 0.0000 0.00000 0.000
#> 5 1 2 1 12 0 0 0.000 0.0000 0.00000 0.000
#> 6 1 5 1 6 0 0 0.259 0.1295 0.00548 0.513
summary(crfit,times=c(1,2,5),cause=1:2)
#> X1 time cause n.risk n.event n.lost cuminc se.cuminc lower upper
#> 1 0 1 1 14 0 0 0.1176 0.0781 0.00000 0.271
#> 2 0 2 1 14 0 0 0.1176 0.0781 0.00000 0.271
#> 3 0 5 1 7 0 0 0.2353 0.1029 0.03365 0.437
#> 4 1 1 1 13 0 0 0.0000 0.0000 0.00000 0.000
#> 5 1 2 1 12 0 0 0.0000 0.0000 0.00000 0.000
#> 6 1 5 1 6 0 0 0.2593 0.1295 0.00548 0.513
#> 7 0 1 2 14 0 0 0.0588 0.0571 0.00000 0.171
#> 8 0 2 2 14 0 0 0.0588 0.0571 0.00000 0.171
#> 9 0 5 2 7 0 0 0.3072 0.1148 0.08215 0.532
#> 10 1 1 2 13 0 0 0.0000 0.0000 0.00000 0.000
#> 11 1 2 2 12 0 0 0.0000 0.0000 0.00000 0.000
#> 12 1 5 2 6 0 0 0.0926 0.0881 0.00000 0.265
# extract the actual tables from the summary
sumfit <- summary(crfit,times=c(1,2,5),print=FALSE)
sumfit$table[[1]] # cause 1
#> NULL
sumfit$table[[2]] # cause 2
#> NULL
# '