prodlim  R Documentation 
Nonparametric estimation in event history analysis. Featuring fast algorithms and user friendly syntax adapted from the survival package. The product limit algorithm is used for right censored data; the selfconsistency algorithm for interval censored data.
prodlim(
formula,
data = parent.frame(),
subset,
na.action = NULL,
reverse = FALSE,
conf.int = 0.95,
bandwidth = NULL,
caseweights,
discrete.level = 3,
x = TRUE,
maxiter = 1000,
grid,
tol = 7,
method = c("npmle", "one.step", "impute.midpoint", "impute.right"),
exact = TRUE,
type
)
formula 
A formula whose left hand side is a 
data 
A data.frame in which all the variables of

subset 
Passed as argument 
na.action 
All lines in data with any missing values in the variables of formula are removed. 
reverse 
For right censored data, if reverse=TRUE then the censoring distribution is estimated. 
conf.int 
The level (between 0 and 1) for twosided pointwise confidence intervals. Defaults to 0.95. Remark: only plain Waldtype confidence limits are available. 
bandwidth 
Smoothing parameter for nearest neighborhoods
based on the values of a continuous covariate. See function

caseweights 
Weights applied to the contribution of each
subject to change the number of events and the number at
risk. This can be used for bootstrap and survey analysis. Should
be a vector of the same length and the same order as 
discrete.level 
Numeric covariates are treated as factors
when their number of unique values exceeds not

x 
logical value: if 
maxiter 
For interval censored data only. Maximal number of iterations to obtain the nonparametric maximum likelihood estimate. Defaults to 1000. 
grid 
For interval censored data only. When method=one.step grid for onestep product limit estimate. Defaults to sorted list of unique left and right endpoints of the observed intervals. 
tol 
For interval censored data only. Numeric value whose negative exponential is used as convergence criterion for finding the nonparametric maximum likelihood estimate. Defaults to 7 meaning exp(7). 
method 
For interval censored data only. If equal to

exact 
If TRUE the grid of time points used for estimation includes all the L and R endpoints of the observed intervals. 
type 
In two state models either 
The response of formula
(ie the left hand side of the ‘~’ operator)
specifies the model.
In twostate models – the classical survival case – the standard
KaplanMeier method is applied. For this the response can be specified as a
Surv
or as a Hist
object. The Hist
function allows you to change the code for censored observations, e.g.
Hist(time,status,cens.code="4")
.
Besides a slight gain of computing efficiency, there are some extensions that are not included in the current version of the survival package:
(0) The KaplanMeier estimator for the censoring times reverse=TRUE
is correctly estimated when there are ties between event and censoring
times.
(1) A conditional version of the kernel smoothed KaplanMeier estimator for at most one continuous predictors using nearest neighborhoods (Beran 1981, Stute 1984, Akritas 1994).
(2) For clustercorrelated data the right hand side of formula
may
identify a cluster
variable. In that case Greenwood's variance
formula is replaced by the formula of Ying and Wei (1994).
(3) Competing risk models can be specified via Hist
response
objects in formula
.
The AalenJohansen estimator is applied for estimating the absolute risk of the competing causes, i.e., the cumulative incidence functions.
Under construction:
(U0) Interval censored event times specified via Hist
are used
to find the nonparametric maximum likelihood estimate. Currently this works
only for twostate models and the results should match with those from the
package ‘Icens’.
(U1) Extensions to more complex multistates models
(U2) The nonparametric maximum likelihood estimate for interval censored observations of competing risks models.
Object of class "prodlim". See print.prodlim
, predict.prodlim
, predict,
summary.prodlim
, plot.prodlim
.
Thomas A. Gerds tag@biostat.ku.dk
Andersen, Borgan, Gill, Keiding (1993) Springer 'Statistical Models Based on Counting Processes'
Akritas (1994) The Annals of Statistics 22, 12991327 Nearest neighbor estimation of a bivariate distribution under random censoring.
R Beran (1981) http://anson.ucdavis.edu/~beran/paper.html 'Nonparametric regression with randomly censored survival data'
Stute (1984) The Annals of Statistics 12, 917–926 'Asymptotic Normality of Nearest Neighbor Regression Function Estimates'
Ying, Wei (1994) Journal of Multivariate Analysis 50, 1729 The KaplanMeier estimate for dependent failure time observations
predictSurv
, predictSurvIndividual
,
predictAbsrisk
, Hist
, neighborhood
,
Surv
, survfit
, strata
,
##twostate survival model
dat < SimSurv(30)
with(dat,plot(Hist(time,status)))
fit < prodlim(Hist(time,status)~1,data=dat)
print(fit)
plot(fit)
summary(fit)
quantile(fit)
## Subset
fit1a < prodlim(Hist(time,status)~1,data=dat,subset=dat$X1==1)
fit1b < prodlim(Hist(time,status)~1,data=dat,subset=dat$X1==1 & dat$X2>0)
## clustered data
library(survival)
cdat < cbind(SimSurv(30),patnr=sample(1:5,size=30,replace=TRUE))
fit < prodlim(Hist(time,status)~cluster(patnr),data=cdat)
print(fit)
plot(fit)
summary(fit)
##compare KaplanMeier to survival package
dat < SimSurv(30)
pfit < prodlim(Surv(time,status)~1,data=dat)
pfit < prodlim(Hist(time,status)~1,data=dat) ## same thing
sfit < survfit(Surv(time,status)~1,data=dat,conf.type="plain")
## same result for the survival distribution function
all(round(pfit$surv,12)==round(sfit$surv,12))
summary(pfit,digits=3)
summary(sfit,times=quantile(unique(dat$time)))
##estimating the censoring survival function
rdat < data.frame(time=c(1,2,3,3,3,4,5,5,6,7),status=c(1,0,0,1,0,1,0,1,1,0))
rpfit < prodlim(Hist(time,status)~1,data=rdat,reverse=TRUE)
rsfit < survfit(Surv(time,1status)~1,data=rdat,conf.type="plain")
## When there are ties between times at which events are observed
## times at which subjects are right censored, then the convention
## is that events come first. This is not obeyed by the above call to survfit,
## and hence only prodlim delivers the correct reverse KaplanMeier:
cbind("Wrong:"=rsfit$surv,"Correct:"=rpfit$surv)
##stratified KaplanMeier
pfit.X2 < prodlim(Surv(time,status)~X2,data=dat)
summary(pfit.X2)
summary(pfit.X2,intervals=TRUE)
plot(pfit.X2)
##continuous covariate: StoneBeran estimate
prodlim(Surv(time,status)~X1,data=dat)
##both discrete and continuous covariates
prodlim(Surv(time,status)~X2+X1,data=dat)
##interval censored data
dat < data.frame(L=1:10,R=c(2,3,12,8,9,10,7,12,12,12),status=c(1,1,0,1,1,1,1,0,0,0))
with(dat,Hist(time=list(L,R),event=status))
dat$event=1
npmle.fitml < prodlim(Hist(time=list(L,R),event)~1,data=dat)
##competing risks
CompRiskFrame < data.frame(time=1:100,event=rbinom(100,2,.5),X=rbinom(100,1,.5))
crFit < prodlim(Hist(time,event)~X,data=CompRiskFrame)
summary(crFit)
plot(crFit)
summary(crFit,cause=2)
plot(crFit,cause=2)
# Changing the cens.code:
dat < data.frame(time=1:10,status=c(1,2,1,2,5,5,1,1,2,2))
fit < prodlim(Hist(time,status)~1,data=dat)
print(fit$model.response)
fit < prodlim(Hist(time,status,cens.code="2")~1,data=dat)
print(fit$model.response)
plot(fit)
plot(fit,cause="5")
##delayed entry
## lefttruncated event times with competing risk endpoint
dat < data.frame(entry=c(7,3,11,12,11,2,1,7,15,17,3),time=10:20,status=c(1,0,2,2,0,0,1,2,0,2,0))
fitd < prodlim(Hist(time=time,event=status,entry=entry)~1,data=dat)
summary(fitd)
plot(fitd)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.