Title: | Estimation of Smooth Hazard Models for Interval-Censored Data |
---|---|
Description: | Estimation of two-state (survival) models and irreversible illness- death models with possibly interval-censored, left-truncated and right-censored data. Proportional intensities regression models can be specified to allow for covariates effects separately for each transition. We use either a parametric approach with Weibull baseline intensities or a semi-parametric approach with M-splines approximation of baseline intensities in order to obtain smooth estimates of the hazard functions. Parameter estimates are obtained by maximum likelihood in the parametric approach and by penalized maximum likelihood in the semi-parametric approach. |
Authors: | Celia Touraine [aut], Thomas Alexander Gerds [aut, cre], Pierre Joly [aut] (Author and maintainer of the Fortran code), Cecile Proust-Lima [aut] (Author of the Fortran code), Helene Jacqmin-Gadda [aut] (Author of the Fortran code), Amadou Diakite [aut] (Author of the Fortran code), W.D. Cody [aut] (Author of the Fortran code), A.H. Morris [aut] (Author of the Fortran code), B.W. Brown [aut] (Author of the Fortran code), Robin Genuer [ctb] |
Maintainer: | Thomas Alexander Gerds <[email protected]> |
License: | GPL (>= 2) |
Version: | 2024.04.10 |
Built: | 2024-11-07 06:36:55 UTC |
Source: | CRAN |
Fit an illness-death model using either a semi-parametric approach (penalized likelihood with an approximation of the transition intensity functions by linear combination of M-splines) or a parametric approach (specifying Weibull distributions on the transition intensities). Left-truncated, right-censored, and interval-censored data are allowed. State 0 corresponds to the initial state, state 1 to the transient one, state 2 to the absorbant one. The allowed transitions are: 0 –> 1, 0 –> 2 and 1 –> 2.
idm( formula01, formula02, formula12, data, maxiter = 200, eps = c(5, 5, 3), n.knots = c(7, 7, 7), knots = "equidistant", CV = FALSE, kappa = c(1000000, 500000, 20000), method = "Weib", conf.int = 0.95, print.iter = FALSE, subset = NULL, na.action = na.fail )
idm( formula01, formula02, formula12, data, maxiter = 200, eps = c(5, 5, 3), n.knots = c(7, 7, 7), knots = "equidistant", CV = FALSE, kappa = c(1000000, 500000, 20000), method = "Weib", conf.int = 0.95, print.iter = FALSE, subset = NULL, na.action = na.fail )
formula01 |
A formula specifying a regression model for the
|
formula02 |
A formula specifying a regression model for the
|
formula12 |
A formula specifying a regression model for the
|
data |
A data frame in which to interpret the variables of
|
maxiter |
Maximum number of iterations. The default is 200. |
eps |
A vector of 3 integers >0 used to define the power of
three convergence criteria: 1. for the regression parameters,
2. for the likelihood, 3. for the second derivatives. The default
is |
n.knots |
For |
knots |
Argument only active for the penalized likelihood approach
The algorithm needs at least 5 knots and allows no more than 20 knots. |
CV |
Binary variable equals to 1 when search (by approximated
cross validation) of the smoothing parameters |
kappa |
Argument only active for the penalized likelihood approach |
method |
type of estimation method: "Splines" for a penalized likelihood approach with approximation of the transition intensities by M-splines, "Weib" for a parametric approach with a Weibull distribution on the transition intensities. Default is "Weib". |
conf.int |
Level of confidence pointwise confidence intervals of the transition intensities, i.e.,
a value between 0 and 1, the default is |
print.iter |
boolean parameter. Equals to |
subset |
expression indicating the subset of the rows of data to be used in the fit. All observations are included by default. |
na.action |
how NAs are treated. The default is first, any na.action attribute of data, second a na.action setting of options, and third 'na.fail' if that is unset. The 'factory-fresh' default is na.omit. Another possible value is NULL. |
The estimated parameters are obtained using the robust Marquardt algorithm (Marquardt, 1963) which is a combination between a Newton-Raphson algorithm and a steepest descent algorithm.
call |
the call that produced the result. |
coef |
regression parameters. |
loglik |
vector containing the log-likelihood without and with covariate. |
cv |
vector containing the convergence criteria. |
niter |
number of iterations. |
converged |
integer equal to 1 when the model converged, 2, 3 or 4 otherwise. |
modelPar |
Weibull parameters. |
N |
number of subjects. |
events1 |
number of events 0 –> 1. |
events2 |
number of events 0 –> 2 or 0 –> 1 –> 2. |
NC |
vector containing the number of covariates on transitions 0 –> 1, 0 –> 2, 1 –> 2. |
responseTrans |
model response for the 0 –> 1
transition. |
responseAbs |
model
response for the 0 –> 2 transition. |
time |
times for which transition intensities have been evaluated for plotting. Vector in the Weibull approach. Matrix in the penalized likelihhod approach for which the colums corresponds to the transitions 0 –> 1, 1 –> 2, 0 –> 2. |
intensity01 |
matched values of the intensities for transition 0 –> 1. |
lowerIntensity01 |
lower confidence intervals for the values of the intensities for transition 0 –> 1. |
upperIntensity01 |
upper confidence intervals for the values of the intensities for transition 0 –> 1. |
intensity02 |
matched values of the intensities for transition 0 –> 2. |
lowerIntensity02 |
lower confidence intervals for the values of the intensities for transition 0 –> 2. |
upperIntensity02 |
upper confidence intervals for the values of the intensities for transition 0 –> 2. |
intensity12 |
matched values of the intensities for transition 1 –> 2. |
lowerIntensity12 |
lower confidence intervals for the values of the intensities for transition 1 –> 2. |
upperIntensity12 |
upper confidence intervals for the values of the intensities for transition 1 –> 2. |
RR |
vector of relative risks. |
V |
variance-covariance matrix derived from the Hessian of the log-likelihood if using method="Weib" or, from the Hessian of the penalized log-likelihood if using method="Splines". |
se |
standart errors of the regression parameters. |
Xnames01 |
names of covariates on 0 –> 1. |
Xnames02 |
names of covariates on 0 –> 2. |
Xnames12 |
names of covariates on 1 –> 2. |
knots01 |
knots to approximate by M-splines the intensity of the 0 –> 1 transition. |
knots02 |
knots to approximate by M-splines the intensity of the 0 –> 2 transition. |
knots12 |
knots to approximate by M-splines the intensity of the 1 –> 2 transition. |
nknots01 |
number of knots on transition 0 –> 1. |
nknots02 |
number of knots on transition 0 –> 2. |
nknots12 |
number of knots on transition 1 –> 2. |
theta01 |
square root of splines coefficients for transition 0 –> 1. |
theta02 |
square root of splines coefficients for transition 0 –> 2. |
theta12 |
square root of splines coefficients for transition 1 –> 2. |
CV |
a binary variable equals to 1 when search of the smoothing parameters kappa by approximated cross-validation, 1 otherwise. The default is 0. |
kappa |
vector containing the smoothing parameters for transition 0 –> 1, 0 –> 2, 1 –> 2 used to estimate the model by the penalized likelihood approach. |
CVcrit |
cross validation criteria. |
DoF |
degrees of freedom of the model. |
na.action |
observations deleted if missing values. |
R: Celia Touraine <[email protected]> Fortran: Pierre Joly <[email protected]>
D. Marquardt (1963). An algorithm for least-squares estimation of nonlinear parameters. SIAM Journal of Applied Mathematics, 431-441.
print.idm
summary.idm
predict.idm
library(lava) library(prodlim) set.seed(17) d <- simulateIDM(100) # right censored data fitRC <- idm(formula01=Hist(time=observed.illtime,event=seen.ill)~X1+X2, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2, formula12=Hist(time=observed.lifetime,event=seen.exit)~X1+X2,data=d, conf.int=FALSE) fitRC set.seed(17) d <- simulateIDM(300) fitRC.splines <- idm(formula01=Hist(time=observed.illtime,event=seen.ill)~X1+X2, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2, formula12=Hist(time=observed.lifetime,event=seen.exit)~1,data=d, conf.int=FALSE,method="splines") # interval censored data fitIC <- idm(formula01=Hist(time=list(L,R),event=seen.ill)~X1+X2, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2, formula12=Hist(time=observed.lifetime,event=seen.exit)~X1+X2,data=d, conf.int=FALSE) fitIC data(Paq1000) # Illness-death model with certif on the 3 transitions # Weibull parametrization and likelihood maximization fit.weib <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, data=Paq1000) # Illness-death model with certif on transitions 01 and 02 # Splines parametrization and penalized likelihood maximization fit.splines <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, formula12=~1, method="Splines", data=Paq1000) fit.weib summary(fit.splines)
library(lava) library(prodlim) set.seed(17) d <- simulateIDM(100) # right censored data fitRC <- idm(formula01=Hist(time=observed.illtime,event=seen.ill)~X1+X2, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2, formula12=Hist(time=observed.lifetime,event=seen.exit)~X1+X2,data=d, conf.int=FALSE) fitRC set.seed(17) d <- simulateIDM(300) fitRC.splines <- idm(formula01=Hist(time=observed.illtime,event=seen.ill)~X1+X2, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2, formula12=Hist(time=observed.lifetime,event=seen.exit)~1,data=d, conf.int=FALSE,method="splines") # interval censored data fitIC <- idm(formula01=Hist(time=list(L,R),event=seen.ill)~X1+X2, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2, formula12=Hist(time=observed.lifetime,event=seen.exit)~X1+X2,data=d, conf.int=FALSE) fitIC data(Paq1000) # Illness-death model with certif on the 3 transitions # Weibull parametrization and likelihood maximization fit.weib <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, data=Paq1000) # Illness-death model with certif on transitions 01 and 02 # Splines parametrization and penalized likelihood maximization fit.splines <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, formula12=~1, method="Splines", data=Paq1000) fit.weib summary(fit.splines)
Function to generate an illness-death model for simulation.
idmModel( scale.illtime = 1/100, shape.illtime = 1, scale.lifetime = 1/100, shape.lifetime = 1, scale.waittime = 1/100, shape.waittime = 1, scale.censtime = 1/100, shape.censtime = 1, n.inspections = 5, schedule = 10, punctuality = 5 )
idmModel( scale.illtime = 1/100, shape.illtime = 1, scale.lifetime = 1/100, shape.lifetime = 1, scale.waittime = 1/100, shape.waittime = 1, scale.censtime = 1/100, shape.censtime = 1, n.inspections = 5, schedule = 10, punctuality = 5 )
scale.illtime |
Weilbull scale for latent illness time |
shape.illtime |
Weilbull shape for latent illness time |
scale.lifetime |
Weilbull scale for latent life time |
shape.lifetime |
Weilbull shape for latent life time |
scale.waittime |
Weilbull scale for latent life time |
shape.waittime |
Weilbull shape for latent life time |
scale.censtime |
Weilbull scale for censoring time |
shape.censtime |
Weilbull shape for censoring time |
n.inspections |
Number of inspection times |
schedule |
Mean of the waiting time between adjacent inspections. |
punctuality |
Standard deviation of waiting time between inspections. |
Based on the functionality of the lava PACKAGE the function generates a latent variable model (latent illtime, waittime and lifetime) and censoring mechanism (censtime, inspection1,inspection2,...,inspectionK).
The function sim.idmModel
then simulates
right censored lifetimes and interval censored illness times.
A latent variable model object lvm
Thomas Alexander Gerds
library(lava) library(prodlim) # generate illness-death model based on exponentially # distributed times m <- idmModel(scale.illtime=1/70, shape.illtime=1.8, scale.lifetime=1/50, shape.lifetime=0.7, scale.waittime=1/30, shape.waittime=0.7) round(sim(m,6),1) # Estimate the parameters of the Weibull models # based on the uncensored exact event times # and the uncensored illstatus. set.seed(18) d <- sim(m,100,latent=FALSE) d$uncensored.status <- 1 f <- idm(formula01=Hist(time=illtime,event=illstatus)~1, formula02=Hist(time=lifetime,event=uncensored.status)~1, data=d, conf.int=FALSE) print(f) # Change the rate of the 0->2 and 0->1 transitions # also the rate of the 1->2 transition # and also lower the censoring rate m <- idmModel(scale.lifetime=1/2000, scale.waittime=1/30, scale.illtime=1/1000, scale.censtime=1/1000) set.seed(18) d <- sim(m,50,latent=TRUE) d$uncensored.status <- 1 f <- idm(formula01=Hist(time=observed.illtime,event=illstatus)~1, formula02=Hist(time=observed.lifetime,event=uncensored.status)~1, data=d, conf.int=FALSE) print(f) # Estimate based on the right censored observations fc <- idm(formula01=Hist(time=illtime,event=seen.ill)~1, formula02=Hist(time=observed.lifetime,event=seen.exit)~1, data=d, conf.int=FALSE) print(fc) # Estimate based on interval censored and right censored observations fi <- idm(formula01=Hist(time=list(L,R),event=seen.ill)~1, formula02=Hist(time=observed.lifetime,event=seen.exit)~1, data=d, conf.int=FALSE) print(fi) # Estimation of covariate effects: # X1, X2, X3 m <- idmModel(shape.waittime=2, scale.lifetime=1/2000, scale.waittime=1/300, scale.illtime=1/10000, scale.censtime=1/10000) distribution(m,"X1") <- binomial.lvm(p=0.3) distribution(m,"X2") <- normal.lvm(mean=120,sd=20) distribution(m,"X3") <- normal.lvm(mean=50,sd=20) regression(m,to="latent.illtime",from="X1") <- 1.7 regression(m,to="latent.illtime",from="X2") <- 0.07 regression(m,to="latent.illtime",from="X3") <- -0.1 regression(m,to="latent.waittime",from="X1") <- 1.8 regression(m,to="latent.lifetime",from="X1") <- 0.7 set.seed(28) d <- sim(m,100,latent=TRUE) head(d) table(ill=d$seen.ill,death=d$seen.exit) # Estimation based on uncensored data d$uncensored.status <- 1 # uncensored data F1 <- idm(formula01=Hist(time=illtime,event=illstatus)~X1+X2+X3, formula02=Hist(time=lifetime,event=uncensored.status)~X1+X2+X3, data=d,conf.int=FALSE) print(F1) # Estimation based on right censored data F2 <- idm(formula01=Hist(time=illtime,event=seen.ill)~X1+X2+X3, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2+X3, data=d,conf.int=FALSE) print(F2) # Estimation based on interval censored and right censored data F3 <- idm(formula01=Hist(time=list(L,R),event=seen.ill)~X1+X2+X3, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2+X3, data=d,conf.int=FALSE) print(F3) cbind(uncensored=F1$coef,right.censored=F2$coef,interval.censored=F3$coef)
library(lava) library(prodlim) # generate illness-death model based on exponentially # distributed times m <- idmModel(scale.illtime=1/70, shape.illtime=1.8, scale.lifetime=1/50, shape.lifetime=0.7, scale.waittime=1/30, shape.waittime=0.7) round(sim(m,6),1) # Estimate the parameters of the Weibull models # based on the uncensored exact event times # and the uncensored illstatus. set.seed(18) d <- sim(m,100,latent=FALSE) d$uncensored.status <- 1 f <- idm(formula01=Hist(time=illtime,event=illstatus)~1, formula02=Hist(time=lifetime,event=uncensored.status)~1, data=d, conf.int=FALSE) print(f) # Change the rate of the 0->2 and 0->1 transitions # also the rate of the 1->2 transition # and also lower the censoring rate m <- idmModel(scale.lifetime=1/2000, scale.waittime=1/30, scale.illtime=1/1000, scale.censtime=1/1000) set.seed(18) d <- sim(m,50,latent=TRUE) d$uncensored.status <- 1 f <- idm(formula01=Hist(time=observed.illtime,event=illstatus)~1, formula02=Hist(time=observed.lifetime,event=uncensored.status)~1, data=d, conf.int=FALSE) print(f) # Estimate based on the right censored observations fc <- idm(formula01=Hist(time=illtime,event=seen.ill)~1, formula02=Hist(time=observed.lifetime,event=seen.exit)~1, data=d, conf.int=FALSE) print(fc) # Estimate based on interval censored and right censored observations fi <- idm(formula01=Hist(time=list(L,R),event=seen.ill)~1, formula02=Hist(time=observed.lifetime,event=seen.exit)~1, data=d, conf.int=FALSE) print(fi) # Estimation of covariate effects: # X1, X2, X3 m <- idmModel(shape.waittime=2, scale.lifetime=1/2000, scale.waittime=1/300, scale.illtime=1/10000, scale.censtime=1/10000) distribution(m,"X1") <- binomial.lvm(p=0.3) distribution(m,"X2") <- normal.lvm(mean=120,sd=20) distribution(m,"X3") <- normal.lvm(mean=50,sd=20) regression(m,to="latent.illtime",from="X1") <- 1.7 regression(m,to="latent.illtime",from="X2") <- 0.07 regression(m,to="latent.illtime",from="X3") <- -0.1 regression(m,to="latent.waittime",from="X1") <- 1.8 regression(m,to="latent.lifetime",from="X1") <- 0.7 set.seed(28) d <- sim(m,100,latent=TRUE) head(d) table(ill=d$seen.ill,death=d$seen.exit) # Estimation based on uncensored data d$uncensored.status <- 1 # uncensored data F1 <- idm(formula01=Hist(time=illtime,event=illstatus)~X1+X2+X3, formula02=Hist(time=lifetime,event=uncensored.status)~X1+X2+X3, data=d,conf.int=FALSE) print(F1) # Estimation based on right censored data F2 <- idm(formula01=Hist(time=illtime,event=seen.ill)~X1+X2+X3, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2+X3, data=d,conf.int=FALSE) print(F2) # Estimation based on interval censored and right censored data F3 <- idm(formula01=Hist(time=list(L,R),event=seen.ill)~X1+X2+X3, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2+X3, data=d,conf.int=FALSE) print(F3) cbind(uncensored=F1$coef,right.censored=F2$coef,interval.censored=F3$coef)
M-spline estimate of the transition intensity function and the cumulative transition intensity function for survival and illness-death models
intensity(times, knots, number.knots, theta, linear.predictor = 0)
intensity(times, knots, number.knots, theta, linear.predictor = 0)
times |
Time points at which to estimate the intensity function |
knots |
Knots for the M-spline |
number.knots |
Number of knots for the M-splines (and I-splines see details) |
theta |
The coefficients for the linear combination of M-splines (and I-splines see details) |
linear.predictor |
Linear predictor beta*Z. When it is non-zero,
transition and cumulative transition are multiplied by |
The estimate of the transition intensity function is a linear
combination of M-splines and the estimate of the cumulative transition
intensity function is a linear combination of I-splines (the integral of a
M-spline is called I-spline). The coefficients theta
are the same for
the M-splines and I-splines.
Important: the theta parameters returned by idm
and shr
are in fact
the square root of the splines coefficients. See examples.
This function is a R-translation of a corresponding Fortran function called susp
. susp
is
used internally by idm
and shr
.
times |
The time points at which the following estimates are evaluated. |
intensity |
The transition intensity function evaluated at |
cumulative.intensity |
The cumulative transition intensity function evaluated at |
survival |
The "survival" function, i.e., exp(-cumulative.intensity) |
R: Celia Touraine <[email protected]> and Thomas Alexander Gerds <[email protected]> Fortran: Pierre Joly <[email protected]>
data(testdata) fit.su <- shr(Hist(time=list(l, r), id) ~ cov, data = testdata,method = "Splines",CV = TRUE) intensity(times = fit.su$time, knots = fit.su$knots, number.knots = fit.su$nknots, theta = fit.su$theta^2) data(Paq1000) fit.idm <- idm(formula02 = Hist(time = t, event = death, entry = e) ~ certif, formula01 = Hist(time = list(l,r), event = dementia) ~ certif, formula12 = ~ certif, method = "Splines", data = Paq1000) # Probability of survival in state 0 at age 80 for a subject with no cep given # that he is in state 0 at 70 su0 <- (intensity(times = 80, knots = fit.idm$knots01, number.knots = fit.idm$nknots01, theta = fit.idm$theta01^2)$survival *intensity(times = 80, knots = fit.idm$knots02, number.knots = fit.idm$nknots02, theta = fit.idm$theta02^2)$survival)/ (intensity(times = 70, knots = fit.idm$knots01, number.knots = fit.idm$nknots01, theta = fit.idm$theta01^2)$survival *intensity(times = 70, knots = fit.idm$knots02, number.knots = fit.idm$nknots02, theta = fit.idm$theta02^2)$survival) # Same result as: predict(fit.idm, s = 70, t = 80, conf.int = FALSE) # see first element
data(testdata) fit.su <- shr(Hist(time=list(l, r), id) ~ cov, data = testdata,method = "Splines",CV = TRUE) intensity(times = fit.su$time, knots = fit.su$knots, number.knots = fit.su$nknots, theta = fit.su$theta^2) data(Paq1000) fit.idm <- idm(formula02 = Hist(time = t, event = death, entry = e) ~ certif, formula01 = Hist(time = list(l,r), event = dementia) ~ certif, formula12 = ~ certif, method = "Splines", data = Paq1000) # Probability of survival in state 0 at age 80 for a subject with no cep given # that he is in state 0 at 70 su0 <- (intensity(times = 80, knots = fit.idm$knots01, number.knots = fit.idm$nknots01, theta = fit.idm$theta01^2)$survival *intensity(times = 80, knots = fit.idm$knots02, number.knots = fit.idm$nknots02, theta = fit.idm$theta02^2)$survival)/ (intensity(times = 70, knots = fit.idm$knots01, number.knots = fit.idm$nknots01, theta = fit.idm$theta01^2)$survival *intensity(times = 70, knots = fit.idm$knots02, number.knots = fit.idm$nknots02, theta = fit.idm$theta02^2)$survival) # Same result as: predict(fit.idm, s = 70, t = 80, conf.int = FALSE) # see first element
Paquid data set composed of 1000 subjects selected randomly from the Paquid data set of 3675 subjects.
A data frame with 1000 rows and the following 8 columns.
dementia status, 0=non-demented, 1=demented
death status, 0=alive, 1=dead
age at entry in the study
for demented subjects: age at the visit before the diagnostic visit; for non-demented subjects: age at the last visit (censoring age)
for demented subjects: age at the diagnostic visit; for non-demented subjects: age at the last visit (censoring age)
for dead subjects: age at death; for alive subject: age at the latest news
primary school
certificate:0=without certificate
, 1=with certificate
gender: 0=female
, 1=male
data(Paq1000)
data(Paq1000)
Plot estimated baseline transition intensities from an object of class
idm
optionally with confidence limits.
## S3 method for class 'idm' plot( x, conf.int = FALSE, citype = "shadow", add = FALSE, axes = TRUE, col, lwd, lty, xlim, ylim, xlab, ylab, legend = TRUE, transition = c("01", "02", "12"), ... )
## S3 method for class 'idm' plot( x, conf.int = FALSE, citype = "shadow", add = FALSE, axes = TRUE, col, lwd, lty, xlim, ylim, xlab, ylab, legend = TRUE, transition = c("01", "02", "12"), ... )
x |
a |
conf.int |
If TRUE show confidence limits |
citype |
Type of confidence limits, can be "shadow" or "bars" |
add |
If TRUE add to existing plot |
axes |
If TRUE axes are drawn |
col |
Color of the lines |
lwd |
Width of the lines |
lty |
Type of the lines |
xlim |
Limits for x-axis |
ylim |
Limits for y-axis |
xlab |
Label for x-axis |
ylab |
Label for y-axis |
legend |
If TRUE a legend is drawn, which can be further controlled via |
transition |
Choose one of the transition intensities: |
... |
Passed to |
Print a plot of the baseline transition intensities of an illness-death model estimated using a Weibull approach.
library(lava) library(prodlim) m <- idmModel(scale.lifetime=1/10,scale.illtime=1/8) distribution(m,"X") <- binomial.lvm() regression(m,latent.lifetime~X) <- 0.7 set.seed(30) d <- sim(m,100) fit.weib <- idm(formula02=Hist(observed.lifetime,event=seen.exit)~1, formula01=Hist(time=list(L,R),event=seen.ill)~1,data=d,conf.int=FALSE) plot(fit.weib) ## FIXME: the limits for the 01 transition are a bit wide!? ## with bootstrap confidence limits fit.weib <- idm(formula02=Hist(observed.lifetime,event=seen.exit)~1, formula01=Hist(time=list(L,R),event=seen.ill)~1,data=d,conf.int=TRUE) plot(fit.weib)
library(lava) library(prodlim) m <- idmModel(scale.lifetime=1/10,scale.illtime=1/8) distribution(m,"X") <- binomial.lvm() regression(m,latent.lifetime~X) <- 0.7 set.seed(30) d <- sim(m,100) fit.weib <- idm(formula02=Hist(observed.lifetime,event=seen.exit)~1, formula01=Hist(time=list(L,R),event=seen.ill)~1,data=d,conf.int=FALSE) plot(fit.weib) ## FIXME: the limits for the 01 transition are a bit wide!? ## with bootstrap confidence limits fit.weib <- idm(formula02=Hist(observed.lifetime,event=seen.exit)~1, formula01=Hist(time=list(L,R),event=seen.ill)~1,data=d,conf.int=TRUE) plot(fit.weib)
Plot estimated baseline survival function from an object of class
shr
. Pointwise confidence limits are available.
## S3 method for class 'shr' plot( x, type = "shr", add = FALSE, newdata = NULL, cause = NULL, col, lty, lwd, ylim, xlim, xlab = "Time", ylab, legend = TRUE, confint = TRUE, timeOrigin = 0, axes = TRUE, percent = TRUE, ... )
## S3 method for class 'shr' plot( x, type = "shr", add = FALSE, newdata = NULL, cause = NULL, col, lty, lwd, ylim, xlim, xlab = "Time", ylab, legend = TRUE, confint = TRUE, timeOrigin = 0, axes = TRUE, percent = TRUE, ... )
x |
a |
type |
type of function to plot. The default is "shr". |
add |
boolean. |
newdata |
newdata. |
cause |
cause. |
col |
col. |
lty |
lty. |
lwd |
lwd. |
ylim |
ylim. |
xlim |
xlim. |
xlab |
xlab. |
ylab |
ylab. |
legend |
legend. |
confint |
confint. |
timeOrigin |
timeOrigin. |
axes |
axes. |
percent |
percent. |
... |
other graphical parameters. |
Print a plot of a suvival model.
R: Celia Touraine <[email protected]> Fortran: Pierre Joly <[email protected]>
# Weibull survival model library(prodlim) data(testdata) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata) # pointwise confidence limits plot(fit.su) # no pointwise confidence limits plot(fit.su,confint=FALSE)
# Weibull survival model library(prodlim) data(testdata) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata) # pointwise confidence limits plot(fit.su) # no pointwise confidence limits plot(fit.su,confint=FALSE)
Predict transition probabilities and cumulative probabilities from an object
of class idmSplines
with confidence intervals are calculated.
## S3 method for class 'idm' predict( object, s, t, newdata, nsim = 200, seed = 21, conf.int = 0.95, lifeExpect = FALSE, maxtime, ... )
## S3 method for class 'idm' predict( object, s, t, newdata, nsim = 200, seed = 21, conf.int = 0.95, lifeExpect = FALSE, maxtime, ... )
object |
an |
s |
time point at which prediction is made. |
t |
time horizon for prediction. |
newdata |
A data frame with covariate values for prediction. |
nsim |
number of simulations for the confidence intervals calculations. The default is 200. |
seed |
Seed passed to |
conf.int |
Level of confidence, i.e., a value between 0 and 1,
the default is |
lifeExpect |
Logical. If |
maxtime |
The upper limit of integration for calculations of life expectancies from Weibull parametrizations. |
... |
other parameters. |
a list containing the following predictions with pointwise confidence intervals:
p00 |
the transition probability
|
p01 |
the transition probability
|
p11 |
the transition probability
|
p12 |
the transition probability
|
p02_0 |
the probability of direct transition from state 0 to state 2. |
p02_1 |
the probability of transition from state 0 to state 2 via state 1. |
p02 |
transition probability |
F01 |
the lifetime
risk of disease. |
F0. |
the probability of exit from state
0. |
R: Celia Touraine <[email protected]> and Thomas Alexander Gerds <[email protected]> Fortran: Pierre Joly <[email protected]>
set.seed(100) d=simulateIDM(n = 200) fit <- idm(formula01=Hist(time=list(L,R),event=seen.ill)~X1+X2+X3, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2+X3, data=d,conf.int=FALSE) predict(fit,s=0,t=80,conf.int=FALSE,lifeExpect=FALSE) predict(fit,s=0,t=80,nsim=4,conf.int=TRUE,lifeExpect=FALSE) predict(fit,s=0,t=80,nsim=4,conf.int=FALSE,lifeExpect=TRUE) data(Paq1000) library(prodlim) fit.paq <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif,data=Paq1000) predict(fit.paq,s=70,t=80,newdata=data.frame(certif=1)) predict(fit.paq,s=70,lifeExpect=TRUE,newdata=data.frame(certif=1)) fit.splines <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, formula12=~1, method="Splines", data=Paq1000) predict(fit.splines,s=70,t=80,newdata=data.frame(certif=1)) predict(fit.splines,s=70,t=80,lifeExpect=TRUE,newdata=data.frame(certif=1),nsim=20)
set.seed(100) d=simulateIDM(n = 200) fit <- idm(formula01=Hist(time=list(L,R),event=seen.ill)~X1+X2+X3, formula02=Hist(time=observed.lifetime,event=seen.exit)~X1+X2+X3, data=d,conf.int=FALSE) predict(fit,s=0,t=80,conf.int=FALSE,lifeExpect=FALSE) predict(fit,s=0,t=80,nsim=4,conf.int=TRUE,lifeExpect=FALSE) predict(fit,s=0,t=80,nsim=4,conf.int=FALSE,lifeExpect=TRUE) data(Paq1000) library(prodlim) fit.paq <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif,data=Paq1000) predict(fit.paq,s=70,t=80,newdata=data.frame(certif=1)) predict(fit.paq,s=70,lifeExpect=TRUE,newdata=data.frame(certif=1)) fit.splines <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, formula12=~1, method="Splines", data=Paq1000) predict(fit.splines,s=70,t=80,newdata=data.frame(certif=1)) predict(fit.splines,s=70,t=80,lifeExpect=TRUE,newdata=data.frame(certif=1),nsim=20)
idm
objectsPrint a summary of a fitted illness-death model
## S3 method for class 'idm' print(x, conf.int = 0.95, digits = 4, pvalDigits = 4, eps = 0.0001, ...)
## S3 method for class 'idm' print(x, conf.int = 0.95, digits = 4, pvalDigits = 4, eps = 0.0001, ...)
x |
Class |
conf.int |
The level of confidence for the hazard ratios. The default is |
digits |
Number of digits to print. |
pvalDigits |
Number of digits to print for p-values. |
eps |
Passed to |
... |
Not used. |
No return value.
Celia Touraine <[email protected]>, Thomas A. Gerds <[email protected]>
data(Paq1000) library(prodlim) fit.splines <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, formula12=~1, method="Splines", data=Paq1000) print(fit.splines)
data(Paq1000) library(prodlim) fit.splines <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, formula12=~1, method="Splines", data=Paq1000) print(fit.splines)
shrSplines
objectsPrint a summary of a fitted illness-death model using the penalized likelihood approach.
## S3 method for class 'shr' print(x, conf.int = 0.95, digits = 4, pvalDigits = 4, eps = 0.0001, ...)
## S3 method for class 'shr' print(x, conf.int = 0.95, digits = 4, pvalDigits = 4, eps = 0.0001, ...)
x |
a |
conf.int |
The level of confidence for the hazard ratios. The default is |
digits |
number of digits to print. |
pvalDigits |
number of digits to print for p-values. |
eps |
convergence criterion used for p-values. |
... |
other unusued arguments. |
No return value.
R: Celia Touraine <[email protected]> Fortran: Pierre Joly <[email protected]>
# a penalized survival model library(prodlim) data(testdata) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata,method="Splines") print(fit.su)
# a penalized survival model library(prodlim) data(testdata) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata,method="Splines") print(fit.su)
Fit a survival model using either a semi-parametric approach (penalized likelihood with an approximation of the hazard function by linear combination of M-splines) or a parametric approach (specifying a Weibull distribution on the hazard function). Left-truncated, right-censored, and interval-censored data are allowed.
shr( formula, data, eps = c(5, 5, 3), n.knots = 7, knots = "equidistant", CV = FALSE, kappa = 10000, conf.int = 0.95, maxiter = 200, method = "Weib", print.iter = FALSE, na.action = na.omit )
shr( formula, data, eps = c(5, 5, 3), n.knots = 7, knots = "equidistant", CV = FALSE, kappa = 10000, conf.int = 0.95, maxiter = 200, method = "Weib", print.iter = FALSE, na.action = na.omit )
formula |
a formula object with the response on the left hand side
and the terms on the right hand side. The
response must be a survival object or Hist object as returned by
the |
data |
a data frame in which to interpret the variables named
in the |
eps |
a vector of length 3 for the convergence criteria
(criterion for parameters, criterion for likelihood, criterion for
second derivatives). The default is |
n.knots |
Argument only active for the penalized likelihood approach |
knots |
Argument only active for the penalized likelihood approach
The algorithm reuqires at least 5 knots and allows no more than 20 knots. |
CV |
binary variable equals to 1 when search (by approximated cross validation) of the smoothing parameter kappa and 0 otherwise. Argument for the penalized likelihood approach. The default is 0. |
kappa |
Argument only active for the penalized likelihood approach |
conf.int |
Level of confidence pointwise confidence intervals of the survival and hazard functions, i.e.,
a value between 0 and 1, the default is |
maxiter |
maximum number of iterations. The default is 200. |
method |
type of estimation method: "Splines" for a penalized
likelihood approach with approximation of the hazard function by
M-splines, "Weib" for a parametric approach with a Weibull
distribution on the hazard function. Default is |
print.iter |
boolean parameter. Equals to |
na.action |
how NAs are treated. The default is first, any na.action attribute of data, second a na.action setting of options, and third 'na.fail' if that is unset. The 'factory-fresh' default is na.omit. Another possible value is NULL. |
The estimated parameters are obtained using the robust Marquardt algorithm (Marquardt, 1963) which is a combination between a Newton-Raphson algorithm and a steepest descent algorithm.
regression parameters.
vector containing the log-likelihood without and with covariate.
Weibull parameters.
number of subjects.
number of covariates.
number of events.
model response: Hist
or Surv
object.
integer equal to 1 when the model converged, 2, 3 or 4 otherwise.
times for which survival and hazard functions have been evaluated for plotting.
matched values of the hazard function.
lower confidence limits for hazard function.
upper confidence limits for hazard function.
matched values of the survival function.
lower confidence limits for survival function.
upper confidence limits for survival function.
vector of relative risks.
variance-covariance matrix.
standard errors.
knots of the M-splines estimate of the hazard function.
number of knots.
a binary variable equals to 1 when search of the smoothing parameter kappa by approximated cross-validation, 1 otherwise. The default is 0.
number of iterations.
vector containing the convergence criteria.
observations deleted if missing values.
R: Celia Touraine [email protected] Fortran: Pierre Joly [email protected]
D. Marquardt (1963). An algorithm for least-squares estimation of nonlinear parameters. SIAM Journal of Applied Mathematics, 431-441.
shr
, print.shr
,
summary.shr
, print.shr
,
# Weibull survival model library(prodlim) data(testdata) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata) fit.su summary(fit.su) shr.spline <- shr(Hist(time=list(l,r),id)~cov,data=testdata,method="splines",n.knots=6) shr.spline shr.spline.q <- shr(Hist(time=list(l,r),id)~cov,data=testdata, method="splines",n.knots=6,knots="quantiles") plot(shr.spline.q) ## manual placement of knots shr.spline.man <- shr(Hist(time=list(l,r),id)~cov,data=testdata,method="splines",knots=seq(0,7,1))
# Weibull survival model library(prodlim) data(testdata) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata) fit.su summary(fit.su) shr.spline <- shr(Hist(time=list(l,r),id)~cov,data=testdata,method="splines",n.knots=6) shr.spline shr.spline.q <- shr(Hist(time=list(l,r),id)~cov,data=testdata, method="splines",n.knots=6,knots="quantiles") plot(shr.spline.q) ## manual placement of knots shr.spline.man <- shr(Hist(time=list(l,r),id)~cov,data=testdata,method="splines",knots=seq(0,7,1))
Function to simulate illness-death model data
## S3 method for class 'idmModel' sim( x, n, illness.known.at.death = TRUE, compliance = 1, latent = FALSE, keep.inspectiontimes = FALSE, ... )
## S3 method for class 'idmModel' sim( x, n, illness.known.at.death = TRUE, compliance = 1, latent = FALSE, keep.inspectiontimes = FALSE, ... )
x |
An |
n |
Number of observations |
illness.known.at.death |
Affects the value of variable seen.ill |
compliance |
Probability of missing an inspection time. |
latent |
if TRUE keep the latent event times |
keep.inspectiontimes |
if |
... |
Extra arguments given to |
Based on the functionality of the lava PACKAGE
A data set with interval censored observations from an illness-death model
Thomas Alexander Gerds
example(idmModel) help(idmModel)
example(idmModel) help(idmModel)
Function to simulate interval censored survival data
## S3 method for class 'survIC' sim(x, n, compliance = 1, latent = TRUE, keep.inspectiontimes = FALSE, ...)
## S3 method for class 'survIC' sim(x, n, compliance = 1, latent = TRUE, keep.inspectiontimes = FALSE, ...)
x |
An |
n |
Number of observations |
compliance |
Probability of missing an inspection time. |
latent |
if TRUE keep the latent event times |
keep.inspectiontimes |
if |
... |
Extra arguments given to |
Based on the functionality of the lava PACKAGE
A data set with interval censored observations
Thomas Alexander Gerds
library(lava) example(survIC) help(survIC) ol <- survIC() dat.ol <- sim(ol,10)
library(lava) example(survIC) help(survIC) ol <- survIC() dat.ol <- sim(ol,10)
Simulate data from an illness-death model with interval censored event times and covariates
simulateIDM(n = 100)
simulateIDM(n = 100)
n |
number of observations |
Simulate data from an illness-death model with interval censored event times and covariates for the purpose of illustrating the help pages of the SmoothHazard package. See the body of the function for details, i.e., evaluate simulateIDM
Object with class data.frame
which contains the simulated data.
idmModel sim.idmModel
# simulateIDM simulateIDM(100)
# simulateIDM simulateIDM(100)
Summarize the event history data of an illness-death regression model and show regression coefficients for transition intensities
## S3 method for class 'idm' summary(object, conf.int = 0.95, digits = 4, pvalDigits = 4, eps = 0.0001, ...)
## S3 method for class 'idm' summary(object, conf.int = 0.95, digits = 4, pvalDigits = 4, eps = 0.0001, ...)
object |
a |
conf.int |
The level of confidence for the hazard ratios. The default is |
digits |
number of digits to print. |
pvalDigits |
number of digits to print for p-values. |
eps |
convergence criterion used for p-values. |
... |
other unusued arguments. |
No return value.
R: Celia Touraine <[email protected]> Fortran: Pierre Joly <[email protected]>
library(prodlim) data(Paq1000) fit.splines <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, formula12=~1, method="Splines", data=Paq1000) summary(fit.splines)
library(prodlim) data(Paq1000) fit.splines <- idm(formula02=Hist(time=t,event=death,entry=e)~certif, formula01=Hist(time=list(l,r),event=dementia)~certif, formula12=~1, method="Splines", data=Paq1000) summary(fit.splines)
Print a short summary of a fitted illness-death model using the penalized likelihood approach.
## S3 method for class 'shr' summary(object, conf.int = 0.95, digits = 4, pvalDigits = 4, eps = 0.0001, ...)
## S3 method for class 'shr' summary(object, conf.int = 0.95, digits = 4, pvalDigits = 4, eps = 0.0001, ...)
object |
a |
conf.int |
The level of confidence for the hazard ratios. The default is |
digits |
number of digits to print. |
pvalDigits |
number of digits to print for p-values. |
eps |
convergence criterion used for p-values. |
... |
other unusued arguments. |
No return value.
Celia Touraine <[email protected]>
# a penalized survival model data(testdata) library(prodlim) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata,method="Splines") summary(fit.su) # Weibull survival model data(testdata) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata) summary(fit.su)
# a penalized survival model data(testdata) library(prodlim) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata,method="Splines") summary(fit.su) # Weibull survival model data(testdata) fit.su <- shr(Hist(time=list(l,r),id)~cov,data=testdata) summary(fit.su)
Function to generate a latent variable model for interval censored survival times.
survIC( scale.time = 1/100, shape.time = 1, n.inspections = 5, schedule = 10, punctuality = 5 )
survIC( scale.time = 1/100, shape.time = 1, n.inspections = 5, schedule = 10, punctuality = 5 )
scale.time |
Weilbull scale for latent time |
shape.time |
Weilbull shape for latent time |
n.inspections |
Number of inspection times |
schedule |
Mean of the waiting time between adjacent inspections. |
punctuality |
Standard deviation of waiting time between inspections. |
Based on the functionality of the lava PACKAGE the function generates a latent variable model with a latent time and a censoring mechanism (censtime, inspection1,inspection2,...,inspectionK).
The function sim.survIC
then simulates
interval censored times.
A latent variable model object lvm
Thomas Alexander Gerds
library(lava) library(prodlim) # generate survival model based on exponentially # distributed times m <- survIC(scale.time=1/50, shape.time=0.7) round(sim(m,6),1) # Estimate the parameters of the Weibull models # based on the uncensored exact event times # and the uncensored illstatus. set.seed(18) d <- sim(m,100,latent=FALSE) d$uncensored.status <- 1 f <- shr(Hist(time=list(L,R),event=uncensored.status)~1, data=d, conf.int=FALSE) print(f)
library(lava) library(prodlim) # generate survival model based on exponentially # distributed times m <- survIC(scale.time=1/50, shape.time=0.7) round(sim(m,6),1) # Estimate the parameters of the Weibull models # based on the uncensored exact event times # and the uncensored illstatus. set.seed(18) d <- sim(m,100,latent=FALSE) d$uncensored.status <- 1 f <- shr(Hist(time=list(L,R),event=uncensored.status)~1, data=d, conf.int=FALSE) print(f)
A simulated data frame for survival models composed of right-censored and interval-censored data.
A data frame with 936 observations on the following 4 variables.
for diseased subjects: left endpoint of censoring interval; for non-diseased subjects: right censoring time
for diseased subjects: right endpoint of censoring interval; for non-diseased subjects: right censoring time for the disease event
disease status
covariate
data(testdata) head(testdata)
data(testdata) head(testdata)