R/predict.idm.R
In tagteam/SmoothHazard: Estimation of Smooth Hazard Models for Interval-Censored Data with Applications to Survival and Illness-Death Models
Defines functions
predict.idm
Documented in
predict.idm
# 0 : health state
# 1 : illness state
# 2 : death state
#' Predictions for an illness-death model using either a penalized likelihood
#' approach or a Weibull parametrization.
#'
#' Predict transition probabilities and cumulative probabilities from an object
#' of class \code{idmSplines} with confidence intervals are calculated.
#'
#' @param object an \code{idm} class objects returned by a call to the
#' \code{\link{idm}} function
#' @param s time point at which prediction is made.
#' @param t time horizon for prediction.
#' @param newdata A data frame with covariate values for prediction.
#' @param nsim number of simulations for the confidence intervals
#' calculations. The default is 200.
#' @param seed Seed passed to \code{set.seed} for Monte Carlo
#' simulation of confidence intervals.
#' @param conf.int Level of confidence, i.e., a value between 0 and 1,
#' the default is \code{0.95}. The default is also used when
#' \code{conf.int=TRUE}. To avoid computation of confidence
#' intervals, set \code{conf.int} to FALSE or NULL.
#' @param lifeExpect Logical. If \code{TRUE} compute life
#' expectancies, i.e., \code{t=Inf}.
#' @param maxtime The upper limit of integration for calculations of life expectancies from Weibull parametrizations.
#' @param ... other parameters.
#' @return a list containing the following predictions with pointwise
#' confidence intervals: \item{p00}{the transition probability
#' \eqn{p_{00}}.} \item{p01}{the transition probability
#' \eqn{p_{01}}.} \item{p11}{the transition probability
#' \eqn{p_{11}}.} \item{p12}{the transition probability
#' \eqn{p_{12}}.} \item{p02_0}{the probability of direct
#' transition from state 0 to state 2.} \item{p02_1}{the
#' probability of transition from state 0 to state 2 via state 1.}
#' \item{p02}{transition probability \eqn{p_{02}}. Note that
#' \code{p02}=\code{p_02_0}+\code{p02_1}.} \item{F01}{the lifetime
#' risk of disease. \code{F01}=\code{p01}+\code{p02_1}.}
#' \item{F0.}{the probability of exit from state
#' 0. \code{F0.}=\code{p02_0}+\code{p01}+\code{p02_1}.}
#' @author R: Celia Touraine <Celia.Touraine@@isped.u-bordeaux2.fr>
#' and Thomas Alexander Gerds <tag@@biostat.ku.dk> Fortran: Pierre
#' Joly <Pierre.Joly@@isped.u-bordeaux2.fr>
#' @seealso \code{\link{idm}}
#' @keywords methods
#' @examples
#'
#' \dontrun{
#' set.seed(100)
#' d=simulateIDM(n = 100)
#' 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)
#'
#'
#' }
#'
#' @export
predict.idm <- function(object,s,t,newdata,nsim=200,seed=21,conf.int=.95,lifeExpect=FALSE,maxtime,...) {
## if (lifeExpect==TRUE) t <- Inf
if (lifeExpect==TRUE) {
t <- Inf
if (!missing(maxtime) && is.numeric(maxtime)) {
maxtime <- min(maxtime,object$maxtime)
} else {
maxtime <- object$maxtime
}
}
if (any(s>t)) {stop("You must respect the condition 's<t' to calculate p(s,t)")}
do.conf.int <- !is.null(conf.int) && !is.na(conf.int) && !conf.int==FALSE
if (is.logical(conf.int)) conf.int <- .95
if (do.conf.int == TRUE){
stopifnot(0<conf.int && conf.int<1)
if (nsim < 2) stop("Need at least two simulations to construct confidence limits.")
}
if (missing(t) && lifeExpect==FALSE) stop("Argument t is missing.")
if (lifeExpect==TRUE) t <- Inf
if (missing(s)) stop("Argument s is missing.")
## if (missing(t) || is.infinite(t)) lifeExpect <- TRUE
# if covariates: cov=c(cov1,cov2,cov3,...)
nvar01 <- object$NC[1]
nvar02 <- object$NC[2]
nvar12 <- object$NC[3]
if (!missing(newdata)){
if (NROW(newdata)>1) stop("Argument newdata has more than one row\n.Currently this function works only for one covariate constallation at a time.")
if (length(object$Xnames01)>0)
Z01 <- as.matrix(model.frame(formula=update.formula(formula(object$terms$Formula01),NULL~.),data=newdata))
else
Z01 <- 0
if (length(object$Xnames02)>0)
Z02 <- as.matrix(model.frame(formula=update.formula(formula(object$terms$Formula02),NULL~.),data=newdata))
else
Z02 <- 0
if (length(object$Xnames12)>0)
Z12 <- as.matrix(model.frame(formula=update.formula(formula(object$terms$Formula12),NULL~.),data=newdata))
else
Z12 <- 0
}else{
vars <- unique(c(object$Xnames01,object$Xnames02,object$Xnames12))
newdata <- data.frame(matrix(0,ncol=pmax(1,length(vars))))
names(newdata) <- vars
Z01 <- matrix(rep(0,length(object$Xnames01)),nrow=1)
Z02 <- matrix(rep(0,length(object$Xnames02)),nrow=1)
Z12 <- matrix(rep(0,length(object$Xnames12)),nrow=1)
}
if(nvar01 > 0){
beta01 <- object$coef[1:nvar01]
names(beta01) <- paste0("beta01.",names(beta01))
bZ01 <- sum(Z01 * beta01)
}else{
bZ01 <- 0
beta01 <- NULL
}
if (nvar02 != 0) {
beta02 <- object$coef[(nvar01+1):(nvar01+nvar02)]
names(beta02) <- paste0("beta02.",names(beta02))
bZ02 <- sum(Z02 * beta02)
}else{
beta02 <- NULL
bZ02 <- 0
}
if (nvar12 != 0) {
beta12 <- object$coef[(nvar01+nvar02+1):(nvar01+nvar02+nvar12)]
names(beta12) <- paste0("beta12.",names(beta12))
bZ12 <- sum(Z12 * beta12)
}else{
beta12 <- NULL
bZ12 <- 0
}
## Splines
if (object$method=="Splines"){
nknots01 <- object$nknots01
nknots02 <- object$nknots02
nknots12 <- object$nknots12
knots01 <- object$knots01
knots02 <- object$knots02
knots12 <- object$knots12
the01 <- object$theta01
names(the01) <- paste0("the01.",1:length(the01))
the02 <- object$theta02
names(the02) <- paste0("the02.",1:length(the02))
the12 <- object$theta12
names(the12) <- paste0("the12.",1:length(the12))
if (do.conf.int == TRUE){
### conf.int prediction by Monte-Carlo
Vmean <- c(the01,the02,the12,beta01,beta02,beta12) # vector of estimates
set.seed(seed)
X <- mvtnorm::rmvnorm(nsim,Vmean,object$V)
colnames(X) <- names(Vmean)
# 1 set of simulated parameters for each element of the list
Xtheta01=X[,names(the01)]^2
Xtheta02=X[,names(the02)]^2
Xtheta12=X[,names(the12)]^2
if (!is.null(beta01))
linPred01=X[,names(beta01),drop=FALSE] %*% t(Z01)
else
linPred01=matrix(0,nrow=nsim,ncol=1)
if (!is.null(beta02))
linPred02=X[,names(beta02),drop=FALSE] %*% t(Z02)
else
linPred02=matrix(0,nrow=nsim,ncol=1)
if (!is.null(beta12))
linPred12=X[,names(beta12),drop=FALSE] %*% t(Z12)
else
linPred12=matrix(0,nrow=nsim,ncol=1)
if (lifeExpect==TRUE){
simResults <- do.call("rbind",
lapply(1:nsim,function(i){
lifexpect0.idmPl(s,
knots01,
nknots01,
Xtheta01[i,],
knots12,
nknots12,
Xtheta12[i,],
knots02,
nknots02,
Xtheta02[i,],
linPred01[i,],
linPred12[i,],
linPred02[i,])
}))
}else{
simResults <- do.call("rbind",
lapply(1:nsim,function(i){
Predict0.idmPl(s,
t,
knots01,
nknots01,
Xtheta01[i,],
knots12,
nknots12,
Xtheta12[i,],
knots02,
nknots02,
Xtheta02[i,],
linPred01[i,],
linPred12[i,],
linPred02[i,])
}))
}
q.lower <- (1-conf.int)/2
q.upper <- 1-q.lower
ci <- apply(simResults,2,function(x)quantile(unlist(x),c(q.lower,q.upper)))
}
if (lifeExpect==TRUE){
transprob <- unlist(lifexpect0.idmPl(s,
knots01,
nknots01,
the01^2,
knots12,
nknots12,
the12^2,
knots02,
nknots02,
the02^2,
bZ01,
bZ12,
bZ02))
}else{
transprob <- unlist(Predict0.idmPl(s,t,knots01,nknots01,the01^2,knots12,nknots12,the12^2,knots02,nknots02,the02^2,bZ01,bZ12,bZ02))
}
}else {
a01 <- object$modelPar[1]
b01 <- object$modelPar[2]
a02 <- object$modelPar[3]
b02 <- object$modelPar[4]
a12 <- object$modelPar[5]
b12 <- object$modelPar[6]
if (do.conf.int==TRUE) {
## conf.int prediction by Monte-Carlo
## vector of parameter estimates
Vmean <- c(sqrt(a01),sqrt(b01),sqrt(a02),sqrt(b02),sqrt(a12),sqrt(b12),beta01,beta02,beta12)
names(Vmean) <- c("a01","b01","a02","b02","a12","b12",names(beta01),names(beta02),names(beta12))
set.seed(seed)
X <- mvtnorm::rmvnorm(nsim,Vmean,object$V)
colnames(X) <- names(Vmean)
# set of simulated parameters for each element of the list
Xa01=X[,"a01"]^2
Xb01=1/(X[,"b01"]^2)
Xa02=X[,"a02"]^2
Xb02=1/(X[,"b02"]^2)
Xa12=X[,"a12"]^2
Xb12=1/(X[,"b12"]^2)
if (!is.null(beta01))
linPred01=X[,names(beta01),drop=FALSE] %*% t(Z01)
else
linPred01=matrix(0,nrow=nsim,ncol=1)
if (!is.null(beta02))
linPred02=X[,names(beta02),drop=FALSE] %*% t(Z02)
else
linPred02=matrix(0,nrow=nsim,ncol=1)
if (!is.null(beta12))
linPred12=X[,names(beta12),drop=FALSE] %*% t(Z12)
else
linPred12=matrix(0,nrow=nsim,ncol=1)
if (lifeExpect==TRUE){
simResults <- do.call("rbind",
lapply(1:nsim,function(i){
lifexpect0.idmWeib(s,
a01=Xa01[[i]],
b01=Xb01[[i]],
a02=Xa02[[i]],
b02=Xb02[[i]],
a12=Xa12[[i]],
b12=Xb12[[i]],
bZ01=linPred01[[i]],
bZ02=linPred02[[i]],
bZ12=linPred12[[i]],max=maxtime)
}))
}else{
simResults <- do.call("rbind",
lapply(1:nsim,function(i){
Predict0.idmWeib(s,
t,
a01=Xa01[[i]],
b01=Xb01[[i]],
a02=Xa02[[i]],
b02=Xb02[[i]],
a12=Xa12[[i]],
b12=Xb12[[i]],
bZ01=linPred01[[i]],
bZ02=linPred02[[i]],
bZ12=linPred12[[i]])
}))
}
q.lower <- (1-conf.int)/2
q.upper <- 1-q.lower
ci <- apply(simResults,2,function(x)quantile(unlist(x),c(q.lower,q.upper)))
}
if (lifeExpect==TRUE){
transprob <- unlist(lifexpect0.idmWeib(s,
a01,
1/b01,
a02,
1/b02,
a12,
1/b12,
bZ01,
bZ02,
bZ12,max=maxtime))
}else{
transprob <- unlist(Predict0.idmWeib(s,
t,
a01,
1/b01,
a02,
1/b02,
a12,
1/b12,
bZ01,
bZ02,
bZ12))
}
}
if (do.conf.int==TRUE){
transprob <- data.frame(cbind(transprob,t(ci)))
names(transprob) <- c("Estimate",paste("Lower",round(100*conf.int),sep="."),paste("Upper",round(100*conf.int),sep="."))
transprob <- cbind("Parameter"=rownames(transprob),transprob)
rownames(transprob) <- NULL
}else{
transprob <- data.frame(cbind(transprob))
names(transprob) <- c("Estimate")
transprob <- cbind("Parameter"=rownames(transprob),transprob)
rownames(transprob) <- NULL
}
out <- list(transprob=transprob)
out <- c(out,list(newdata=newdata))
out <- c(out,list(s=s,t=t,conf.int=ifelse(do.conf.int,conf.int,FALSE)))
class(out) <- "predict.idm"
out
}
Predict0.idmPl <- function(s,t,knots01,nknots01,the01,knots12,nknots12,the12,knots02,nknots02,the02,bZ01=0,bZ12=0,bZ02=0) {
if (s>(min(knots01[nknots01+6],knots02[nknots02+6],knots12[nknots12+6]))) {stop("argument s is off")}
if (any(t>knots12[nknots12+6])) {stop("argument t is off")}
if (any(s<knots01[1])) {stop("argument s is off")}
p11 <- S.pl(s,t,knots12,nknots12,the12,bZ12)
p12 <- 1-p11
p00 <- S.pl(s,t,knots01,nknots01,the01,bZ01)*S.pl(s,t,knots02,nknots02,the02,bZ02)
p02_0 <- sapply(t,function(t) {integrate(f=function(x){S.pl(s,x,knots01,nknots01,the01,bZ01)*S.pl(s,x,knots02,nknots02,the02,bZ02)*intensity(times=x,knots=knots02,number.knots=nknots02,theta=the02,linear.predictor=bZ02)$intensity},lower=s,upper=t)$value})
p01 <- sapply(t,function(t) {integrate(f=function(x){S.pl(s,x,knots01,nknots01,the01,bZ01)*S.pl(s,x,knots02,nknots02,the02,bZ02)*intensity(times=x,knots=knots01,number.knots=nknots01,theta=the01,linear.predictor=bZ01)$intensity*S.pl(x,t,knots12,nknots12,the12,bZ12)},lower=s,upper=t)$value})
p02_1 <- 1-p00-p02_0-p01
p02 <- p02_0+p02_1
RM<- integrate(f = function(x) {S.pl(s, x, knots01, nknots01, the01, bZ01) * S.pl(s, x, knots02, nknots02, the02, bZ02) }, lower = s, upper = t)$value
F01<- integrate(f = function(x) {S.pl(s, x, knots01, nknots01, the01, bZ01) * S.pl(s, x, knots02, nknots02, the02, bZ02) * intensity(times = x, knots = knots01, number.knots = nknots01, theta = the01,linear.predictor = bZ01)$intensity }, lower = s, upper = t)$value
list(p00=p00,p01=p01,p11=p11,p12=p12,p02_0=p02_0,p02_1=p02_1,p02=p02,F01=F01,F0.=p02_0+p01+p02_1, RM=RM)
}
Predict0.idmWeib <- function(s,t,a01,b01,a02,b02,a12,b12,bZ01=0,bZ02=0,bZ12=0) {
p11 = S.weib(s,t,a12,b12,bZ12)
p12 = 1-p11
p00 = S.weib(s,t,a01,b01,bZ01)*S.weib(s,t,a02,b02,bZ02)
p02_0 = sapply(t,function(t) {integrate(f=function(x){S.weib(s,x,a01,b01,bZ01)*S.weib(s,x,a02,b02,bZ02)*iweibull(x,a02,b02,bZ02)},lower=s,upper=t)$value })
p01 = sapply(t,function(t) {integrate(f=function(x){S.weib(s,x,a01,b01,bZ01)*S.weib(s,x,a02,b02,bZ02)*iweibull(x,a01,b01,bZ01)*S.weib(x,t,a12,b12,bZ12)},lower=s,upper=t)$value})
p02_1 = 1-p00-p02_0-p01
p02 = p02_0+p02_1
RM= integrate(f = function(x) {S.weib(s, x, a01, b01, bZ01) * S.weib(s, x, a02, b02, bZ02)}, lower = s, upper = t)$value
F01= integrate(f = function(x) {S.weib(s, x, a01, b01, bZ01) * S.weib(s, x, a02, b02, bZ02)*iweibull(x,a01,b01,bZ01)}, lower = s, upper = t)$value
## return(list(p00=p00,p01=p01,p11=p11,p12=p12,p02_0=p02_0,p02_1=p02_1,p02=p02,F01=p01+p02_1,F0.=p02_0+p01+p02_1))
list(p00=p00,p01=p01,p11=p11,p12=p12,p02_0=p02_0,p02_1=p02_1,p02=p02,F01=F01,F0.=p02_0+p01+p02_1, RM=RM)
}
# a = shape parameter
# b = scale parameter
iweibull <- function(x,a,b,bZ=0) {
res = (a/b) * (x/b)**(a-1) * exp(bZ)
return(res)
}
# S(s,t) = S(t)/S(s)
# if S(s)=0, S(s,t)=0
S.weib <- function(s,t,a,b,bZ=0) {
res <- 0
St <- (1-pweibull(t,shape=a,scale=b))^(exp(bZ))
Ss <- (1-pweibull(s,shape=a,scale=b))^(exp(bZ))
if (length(s)==1){
if (Ss==0){res <- 0}
else{res <- St/Ss}
}else{
idx0 <- which(Ss==0)
idx <- which(Ss!=0)
res[idx0] <- 0
res[idx] <- St/Ss[idx]
}
return(res)
}
A <- function(s,t,zi,nknots,the,bZ=0) {
res=rep(0,length(t))
TF = (t>=zi[length(zi)])
ind = which(TF)
if (sum(TF)!=0) {res[ind]=intensity((zi[nknots+6]-10^-5),zi,nknots,the,bZ)$cumul.intensity-intensity(s,zi,nknots,the,bZ)$cumul.intensity}
TF = (t<zi[length(zi)])
ind = which(TF)
if (sum(TF)!=0) {res[ind]=intensity(t[ind],zi,nknots,the,bZ)$cumul.intensity-intensity(s,zi,nknots,the,bZ)$cumul.intensity}
return(res)
}
### Survival function with two time s, t
# S(s,t) = S(t)/S(s)
# = exp(-A(s,t))
S.pl <- function(s,t,zi,nknots,the,bZ=0) {
if (length(t)>=length(s)){
res=rep(0,length(t))
TF = (t>zi[length(zi)])
ind = which(TF)
if (sum(TF)!=0) {res[ind]=0}
TF = (t<=zi[length(zi)])
ind = which(TF)
if (sum(TF)!=0) {res[ind]=intensity(t[ind],zi,nknots,the,bZ)$survival/intensity(s,zi,nknots,the,bZ)$survival}
}else{
res=rep(0,length(s))
if (t>zi[length(zi)]) {res=0}
else {res=intensity(t,zi,nknots,the,bZ)$survival/intensity(s,zi,nknots,the,bZ)$survival}
}
return(res)
}
lifexpect0.idmWeib <- function(s,a01,b01,a02,b02,a12,b12,bZ01=0,bZ02=0,bZ12=0,max) {
## print("lifexpect0.idmWeib")
# max <- 100
ET12 = integrate(
f=function(x) {
S.weib(s,x,a12,b12,bZ12)
},s,max)
ET0dot = integrate(f=function(x) {
S.weib(s,x,a01,b01,bZ01)*S.weib(s,x,a02,b02,bZ02)
},s,max)
ET01 = integrate(f=function(x){
sapply(x,function(x){
integrate(f=function(y){
S.weib(s,y,a01,b01,bZ01)*S.weib(s,y,a02,b02,bZ02)*iweibull(y,a01,b01,bZ01)*S.weib(y,x,a12,b12,bZ12)},
lower=s,
upper=x)$value})},s,max)
LTR=integrate(f=function(x){
S.weib(s,x,a01,b01,bZ01)*S.weib(s,x,a02,b02,bZ02)*iweibull(x,a01,b01,bZ01)},s,max)
list(LE.00=ET0dot$value,
LE.0.=ET01$value+ET0dot$value,
LE.01=ET01$value,
LE.11=ET12$value,
LTR=LTR$value)
}
lifexpect0.idmPl <- function(s,knots01,nknots01,the01,knots12,nknots12,the12,knots02,nknots02,the02,bZ01=0,bZ12=0,bZ02=0) {
ET12 = integrate(f=function(x) {
S.pl(s,x,knots12,nknots12,the12,bZ12)},s,knots12[nknots12+6])
ET0dot = integrate(f=function(x) {
S.pl(s,x,knots01,nknots01,the01,bZ01)*S.pl(s,x,knots02,nknots02,the02,bZ02) },s,knots02[nknots02+6])
ET01 = integrate(f=function(x) {
sapply(x,function(x) {integrate(f=function(y){
(S.pl(s,y,knots01,nknots01,the01,bZ01)
*S.pl(s,y,knots02,nknots02,the02,bZ02)*
intensity(times=y,knots=knots01,number.knots=nknots01,theta=the01,linear.predictor=bZ01)$intensity
*S.pl(y,x,knots12,nknots12,the12,bZ12))},
lower=s,upper=x)$value})},s,knots01[nknots01+6])
LTR=integrate(f = function(x) {S.pl(s, x, knots01, nknots01, the01, bZ01) * S.pl(s, x, knots02, nknots02, the02, bZ02) * intensity(times = x, knots = knots01, number.knots = nknots01, theta = the01,linear.predictor = bZ01)$intensity }, lower = s, upper = knots01[nknots01+6])$value
list(LE.00=ET0dot$value,
LE.0.=ET01$value+ET0dot$value,
LE.01=ET01$value,
LE.11=ET12$value,
LTR=LTR)
}
tagteam/SmoothHazard documentation built on April 5, 2024, 6:32 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions predict.idm
Documented in predict.idm
# 0 : health state
# 1 : illness state
# 2 : death state
#' Predictions for an illness-death model using either a penalized likelihood
#' approach or a Weibull parametrization.
#'
#' Predict transition probabilities and cumulative probabilities from an object
#' of class \code{idmSplines} with confidence intervals are calculated.
#'
#' @param object an \code{idm} class objects returned by a call to the
#' \code{\link{idm}} function
#' @param s time point at which prediction is made.
#' @param t time horizon for prediction.
#' @param newdata A data frame with covariate values for prediction.
#' @param nsim number of simulations for the confidence intervals
#' calculations. The default is 200.
#' @param seed Seed passed to \code{set.seed} for Monte Carlo
#' simulation of confidence intervals.
#' @param conf.int Level of confidence, i.e., a value between 0 and 1,
#' the default is \code{0.95}. The default is also used when
#' \code{conf.int=TRUE}. To avoid computation of confidence
#' intervals, set \code{conf.int} to FALSE or NULL.
#' @param lifeExpect Logical. If \code{TRUE} compute life
#' expectancies, i.e., \code{t=Inf}.
#' @param maxtime The upper limit of integration for calculations of life expectancies from Weibull parametrizations.
#' @param ... other parameters.
#' @return a list containing the following predictions with pointwise
#' confidence intervals: \item{p00}{the transition probability
#' \eqn{p_{00}}.} \item{p01}{the transition probability
#' \eqn{p_{01}}.} \item{p11}{the transition probability
#' \eqn{p_{11}}.} \item{p12}{the transition probability
#' \eqn{p_{12}}.} \item{p02_0}{the probability of direct
#' transition from state 0 to state 2.} \item{p02_1}{the
#' probability of transition from state 0 to state 2 via state 1.}
#' \item{p02}{transition probability \eqn{p_{02}}. Note that
#' \code{p02}=\code{p_02_0}+\code{p02_1}.} \item{F01}{the lifetime
#' risk of disease. \code{F01}=\code{p01}+\code{p02_1}.}
#' \item{F0.}{the probability of exit from state
#' 0. \code{F0.}=\code{p02_0}+\code{p01}+\code{p02_1}.}
#' @author R: Celia Touraine <Celia.Touraine@@isped.u-bordeaux2.fr>
#' and Thomas Alexander Gerds <tag@@biostat.ku.dk> Fortran: Pierre
#' Joly <Pierre.Joly@@isped.u-bordeaux2.fr>
#' @seealso \code{\link{idm}}
#' @keywords methods
#' @examples
#'
#' \dontrun{
#' set.seed(100)
#' d=simulateIDM(n = 100)
#' 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)
#'
#'
#' }
#'
#' @export
predict.idm <- function(object,s,t,newdata,nsim=200,seed=21,conf.int=.95,lifeExpect=FALSE,maxtime,...) {
## if (lifeExpect==TRUE) t <- Inf
if (lifeExpect==TRUE) {
t <- Inf
if (!missing(maxtime) && is.numeric(maxtime)) {
maxtime <- min(maxtime,object$maxtime)
} else {
maxtime <- object$maxtime
}
}
if (any(s>t)) {stop("You must respect the condition 's<t' to calculate p(s,t)")}
do.conf.int <- !is.null(conf.int) && !is.na(conf.int) && !conf.int==FALSE
if (is.logical(conf.int)) conf.int <- .95
if (do.conf.int == TRUE){
stopifnot(0<conf.int && conf.int<1)
if (nsim < 2) stop("Need at least two simulations to construct confidence limits.")
}
if (missing(t) && lifeExpect==FALSE) stop("Argument t is missing.")
if (lifeExpect==TRUE) t <- Inf
if (missing(s)) stop("Argument s is missing.")
## if (missing(t) || is.infinite(t)) lifeExpect <- TRUE
# if covariates: cov=c(cov1,cov2,cov3,...)
nvar01 <- object$NC[1]
nvar02 <- object$NC[2]
nvar12 <- object$NC[3]
if (!missing(newdata)){
if (NROW(newdata)>1) stop("Argument newdata has more than one row\n.Currently this function works only for one covariate constallation at a time.")
if (length(object$Xnames01)>0)
Z01 <- as.matrix(model.frame(formula=update.formula(formula(object$terms$Formula01),NULL~.),data=newdata))
else
Z01 <- 0
if (length(object$Xnames02)>0)
Z02 <- as.matrix(model.frame(formula=update.formula(formula(object$terms$Formula02),NULL~.),data=newdata))
else
Z02 <- 0
if (length(object$Xnames12)>0)
Z12 <- as.matrix(model.frame(formula=update.formula(formula(object$terms$Formula12),NULL~.),data=newdata))
else
Z12 <- 0
}else{
vars <- unique(c(object$Xnames01,object$Xnames02,object$Xnames12))
newdata <- data.frame(matrix(0,ncol=pmax(1,length(vars))))
names(newdata) <- vars
Z01 <- matrix(rep(0,length(object$Xnames01)),nrow=1)
Z02 <- matrix(rep(0,length(object$Xnames02)),nrow=1)
Z12 <- matrix(rep(0,length(object$Xnames12)),nrow=1)
}
if(nvar01 > 0){
beta01 <- object$coef[1:nvar01]
names(beta01) <- paste0("beta01.",names(beta01))
bZ01 <- sum(Z01 * beta01)
}else{
bZ01 <- 0
beta01 <- NULL
}
if (nvar02 != 0) {
beta02 <- object$coef[(nvar01+1):(nvar01+nvar02)]
names(beta02) <- paste0("beta02.",names(beta02))
bZ02 <- sum(Z02 * beta02)
}else{
beta02 <- NULL
bZ02 <- 0
}
if (nvar12 != 0) {
beta12 <- object$coef[(nvar01+nvar02+1):(nvar01+nvar02+nvar12)]
names(beta12) <- paste0("beta12.",names(beta12))
bZ12 <- sum(Z12 * beta12)
}else{
beta12 <- NULL
bZ12 <- 0
}
## Splines
if (object$method=="Splines"){
nknots01 <- object$nknots01
nknots02 <- object$nknots02
nknots12 <- object$nknots12
knots01 <- object$knots01
knots02 <- object$knots02
knots12 <- object$knots12
the01 <- object$theta01
names(the01) <- paste0("the01.",1:length(the01))
the02 <- object$theta02
names(the02) <- paste0("the02.",1:length(the02))
the12 <- object$theta12
names(the12) <- paste0("the12.",1:length(the12))
if (do.conf.int == TRUE){
### conf.int prediction by Monte-Carlo
Vmean <- c(the01,the02,the12,beta01,beta02,beta12) # vector of estimates
set.seed(seed)
X <- mvtnorm::rmvnorm(nsim,Vmean,object$V)
colnames(X) <- names(Vmean)
# 1 set of simulated parameters for each element of the list
Xtheta01=X[,names(the01)]^2
Xtheta02=X[,names(the02)]^2
Xtheta12=X[,names(the12)]^2
if (!is.null(beta01))
linPred01=X[,names(beta01),drop=FALSE] %*% t(Z01)
else
linPred01=matrix(0,nrow=nsim,ncol=1)
if (!is.null(beta02))
linPred02=X[,names(beta02),drop=FALSE] %*% t(Z02)
else
linPred02=matrix(0,nrow=nsim,ncol=1)
if (!is.null(beta12))
linPred12=X[,names(beta12),drop=FALSE] %*% t(Z12)
else
linPred12=matrix(0,nrow=nsim,ncol=1)
if (lifeExpect==TRUE){
simResults <- do.call("rbind",
lapply(1:nsim,function(i){
lifexpect0.idmPl(s,
knots01,
nknots01,
Xtheta01[i,],
knots12,
nknots12,
Xtheta12[i,],
knots02,
nknots02,
Xtheta02[i,],
linPred01[i,],
linPred12[i,],
linPred02[i,])
}))
}else{
simResults <- do.call("rbind",
lapply(1:nsim,function(i){
Predict0.idmPl(s,
t,
knots01,
nknots01,
Xtheta01[i,],
knots12,
nknots12,
Xtheta12[i,],
knots02,
nknots02,
Xtheta02[i,],
linPred01[i,],
linPred12[i,],
linPred02[i,])
}))
}
q.lower <- (1-conf.int)/2
q.upper <- 1-q.lower
ci <- apply(simResults,2,function(x)quantile(unlist(x),c(q.lower,q.upper)))
}
if (lifeExpect==TRUE){
transprob <- unlist(lifexpect0.idmPl(s,
knots01,
nknots01,
the01^2,
knots12,
nknots12,
the12^2,
knots02,
nknots02,
the02^2,
bZ01,
bZ12,
bZ02))
}else{
transprob <- unlist(Predict0.idmPl(s,t,knots01,nknots01,the01^2,knots12,nknots12,the12^2,knots02,nknots02,the02^2,bZ01,bZ12,bZ02))
}
}else {
a01 <- object$modelPar[1]
b01 <- object$modelPar[2]
a02 <- object$modelPar[3]
b02 <- object$modelPar[4]
a12 <- object$modelPar[5]
b12 <- object$modelPar[6]
if (do.conf.int==TRUE) {
## conf.int prediction by Monte-Carlo
## vector of parameter estimates
Vmean <- c(sqrt(a01),sqrt(b01),sqrt(a02),sqrt(b02),sqrt(a12),sqrt(b12),beta01,beta02,beta12)
names(Vmean) <- c("a01","b01","a02","b02","a12","b12",names(beta01),names(beta02),names(beta12))
set.seed(seed)
X <- mvtnorm::rmvnorm(nsim,Vmean,object$V)
colnames(X) <- names(Vmean)
# set of simulated parameters for each element of the list
Xa01=X[,"a01"]^2
Xb01=1/(X[,"b01"]^2)
Xa02=X[,"a02"]^2
Xb02=1/(X[,"b02"]^2)
Xa12=X[,"a12"]^2
Xb12=1/(X[,"b12"]^2)
if (!is.null(beta01))
linPred01=X[,names(beta01),drop=FALSE] %*% t(Z01)
else
linPred01=matrix(0,nrow=nsim,ncol=1)
if (!is.null(beta02))
linPred02=X[,names(beta02),drop=FALSE] %*% t(Z02)
else
linPred02=matrix(0,nrow=nsim,ncol=1)
if (!is.null(beta12))
linPred12=X[,names(beta12),drop=FALSE] %*% t(Z12)
else
linPred12=matrix(0,nrow=nsim,ncol=1)
if (lifeExpect==TRUE){
simResults <- do.call("rbind",
lapply(1:nsim,function(i){
lifexpect0.idmWeib(s,
a01=Xa01[[i]],
b01=Xb01[[i]],
a02=Xa02[[i]],
b02=Xb02[[i]],
a12=Xa12[[i]],
b12=Xb12[[i]],
bZ01=linPred01[[i]],
bZ02=linPred02[[i]],
bZ12=linPred12[[i]],max=maxtime)
}))
}else{
simResults <- do.call("rbind",
lapply(1:nsim,function(i){
Predict0.idmWeib(s,
t,
a01=Xa01[[i]],
b01=Xb01[[i]],
a02=Xa02[[i]],
b02=Xb02[[i]],
a12=Xa12[[i]],
b12=Xb12[[i]],
bZ01=linPred01[[i]],
bZ02=linPred02[[i]],
bZ12=linPred12[[i]])
}))
}
q.lower <- (1-conf.int)/2
q.upper <- 1-q.lower
ci <- apply(simResults,2,function(x)quantile(unlist(x),c(q.lower,q.upper)))
}
if (lifeExpect==TRUE){
transprob <- unlist(lifexpect0.idmWeib(s,
a01,
1/b01,
a02,
1/b02,
a12,
1/b12,
bZ01,
bZ02,
bZ12,max=maxtime))
}else{
transprob <- unlist(Predict0.idmWeib(s,
t,
a01,
1/b01,
a02,
1/b02,
a12,
1/b12,
bZ01,
bZ02,
bZ12))
}
}
if (do.conf.int==TRUE){
transprob <- data.frame(cbind(transprob,t(ci)))
names(transprob) <- c("Estimate",paste("Lower",round(100*conf.int),sep="."),paste("Upper",round(100*conf.int),sep="."))
transprob <- cbind("Parameter"=rownames(transprob),transprob)
rownames(transprob) <- NULL
}else{
transprob <- data.frame(cbind(transprob))
names(transprob) <- c("Estimate")
transprob <- cbind("Parameter"=rownames(transprob),transprob)
rownames(transprob) <- NULL
}
out <- list(transprob=transprob)
out <- c(out,list(newdata=newdata))
out <- c(out,list(s=s,t=t,conf.int=ifelse(do.conf.int,conf.int,FALSE)))
class(out) <- "predict.idm"
out
}
Predict0.idmPl <- function(s,t,knots01,nknots01,the01,knots12,nknots12,the12,knots02,nknots02,the02,bZ01=0,bZ12=0,bZ02=0) {
if (s>(min(knots01[nknots01+6],knots02[nknots02+6],knots12[nknots12+6]))) {stop("argument s is off")}
if (any(t>knots12[nknots12+6])) {stop("argument t is off")}
if (any(s<knots01[1])) {stop("argument s is off")}
p11 <- S.pl(s,t,knots12,nknots12,the12,bZ12)
p12 <- 1-p11
p00 <- S.pl(s,t,knots01,nknots01,the01,bZ01)*S.pl(s,t,knots02,nknots02,the02,bZ02)
p02_0 <- sapply(t,function(t) {integrate(f=function(x){S.pl(s,x,knots01,nknots01,the01,bZ01)*S.pl(s,x,knots02,nknots02,the02,bZ02)*intensity(times=x,knots=knots02,number.knots=nknots02,theta=the02,linear.predictor=bZ02)$intensity},lower=s,upper=t)$value})
p01 <- sapply(t,function(t) {integrate(f=function(x){S.pl(s,x,knots01,nknots01,the01,bZ01)*S.pl(s,x,knots02,nknots02,the02,bZ02)*intensity(times=x,knots=knots01,number.knots=nknots01,theta=the01,linear.predictor=bZ01)$intensity*S.pl(x,t,knots12,nknots12,the12,bZ12)},lower=s,upper=t)$value})
p02_1 <- 1-p00-p02_0-p01
p02 <- p02_0+p02_1
RM<- integrate(f = function(x) {S.pl(s, x, knots01, nknots01, the01, bZ01) * S.pl(s, x, knots02, nknots02, the02, bZ02) }, lower = s, upper = t)$value
F01<- integrate(f = function(x) {S.pl(s, x, knots01, nknots01, the01, bZ01) * S.pl(s, x, knots02, nknots02, the02, bZ02) * intensity(times = x, knots = knots01, number.knots = nknots01, theta = the01,linear.predictor = bZ01)$intensity }, lower = s, upper = t)$value
list(p00=p00,p01=p01,p11=p11,p12=p12,p02_0=p02_0,p02_1=p02_1,p02=p02,F01=F01,F0.=p02_0+p01+p02_1, RM=RM)
}
Predict0.idmWeib <- function(s,t,a01,b01,a02,b02,a12,b12,bZ01=0,bZ02=0,bZ12=0) {
p11 = S.weib(s,t,a12,b12,bZ12)
p12 = 1-p11
p00 = S.weib(s,t,a01,b01,bZ01)*S.weib(s,t,a02,b02,bZ02)
p02_0 = sapply(t,function(t) {integrate(f=function(x){S.weib(s,x,a01,b01,bZ01)*S.weib(s,x,a02,b02,bZ02)*iweibull(x,a02,b02,bZ02)},lower=s,upper=t)$value })
p01 = sapply(t,function(t) {integrate(f=function(x){S.weib(s,x,a01,b01,bZ01)*S.weib(s,x,a02,b02,bZ02)*iweibull(x,a01,b01,bZ01)*S.weib(x,t,a12,b12,bZ12)},lower=s,upper=t)$value})
p02_1 = 1-p00-p02_0-p01
p02 = p02_0+p02_1
RM= integrate(f = function(x) {S.weib(s, x, a01, b01, bZ01) * S.weib(s, x, a02, b02, bZ02)}, lower = s, upper = t)$value
F01= integrate(f = function(x) {S.weib(s, x, a01, b01, bZ01) * S.weib(s, x, a02, b02, bZ02)*iweibull(x,a01,b01,bZ01)}, lower = s, upper = t)$value
## return(list(p00=p00,p01=p01,p11=p11,p12=p12,p02_0=p02_0,p02_1=p02_1,p02=p02,F01=p01+p02_1,F0.=p02_0+p01+p02_1))
list(p00=p00,p01=p01,p11=p11,p12=p12,p02_0=p02_0,p02_1=p02_1,p02=p02,F01=F01,F0.=p02_0+p01+p02_1, RM=RM)
}
# a = shape parameter
# b = scale parameter
iweibull <- function(x,a,b,bZ=0) {
res = (a/b) * (x/b)**(a-1) * exp(bZ)
return(res)
}
# S(s,t) = S(t)/S(s)
# if S(s)=0, S(s,t)=0
S.weib <- function(s,t,a,b,bZ=0) {
res <- 0
St <- (1-pweibull(t,shape=a,scale=b))^(exp(bZ))
Ss <- (1-pweibull(s,shape=a,scale=b))^(exp(bZ))
if (length(s)==1){
if (Ss==0){res <- 0}
else{res <- St/Ss}
}else{
idx0 <- which(Ss==0)
idx <- which(Ss!=0)
res[idx0] <- 0
res[idx] <- St/Ss[idx]
}
return(res)
}
A <- function(s,t,zi,nknots,the,bZ=0) {
res=rep(0,length(t))
TF = (t>=zi[length(zi)])
ind = which(TF)
if (sum(TF)!=0) {res[ind]=intensity((zi[nknots+6]-10^-5),zi,nknots,the,bZ)$cumul.intensity-intensity(s,zi,nknots,the,bZ)$cumul.intensity}
TF = (t<zi[length(zi)])
ind = which(TF)
if (sum(TF)!=0) {res[ind]=intensity(t[ind],zi,nknots,the,bZ)$cumul.intensity-intensity(s,zi,nknots,the,bZ)$cumul.intensity}
return(res)
}
### Survival function with two time s, t
# S(s,t) = S(t)/S(s)
# = exp(-A(s,t))
S.pl <- function(s,t,zi,nknots,the,bZ=0) {
if (length(t)>=length(s)){
res=rep(0,length(t))
TF = (t>zi[length(zi)])
ind = which(TF)
if (sum(TF)!=0) {res[ind]=0}
TF = (t<=zi[length(zi)])
ind = which(TF)
if (sum(TF)!=0) {res[ind]=intensity(t[ind],zi,nknots,the,bZ)$survival/intensity(s,zi,nknots,the,bZ)$survival}
}else{
res=rep(0,length(s))
if (t>zi[length(zi)]) {res=0}
else {res=intensity(t,zi,nknots,the,bZ)$survival/intensity(s,zi,nknots,the,bZ)$survival}
}
return(res)
}
lifexpect0.idmWeib <- function(s,a01,b01,a02,b02,a12,b12,bZ01=0,bZ02=0,bZ12=0,max) {
## print("lifexpect0.idmWeib")
# max <- 100
ET12 = integrate(
f=function(x) {
S.weib(s,x,a12,b12,bZ12)
},s,max)
ET0dot = integrate(f=function(x) {
S.weib(s,x,a01,b01,bZ01)*S.weib(s,x,a02,b02,bZ02)
},s,max)
ET01 = integrate(f=function(x){
sapply(x,function(x){
integrate(f=function(y){
S.weib(s,y,a01,b01,bZ01)*S.weib(s,y,a02,b02,bZ02)*iweibull(y,a01,b01,bZ01)*S.weib(y,x,a12,b12,bZ12)},
lower=s,
upper=x)$value})},s,max)
LTR=integrate(f=function(x){
S.weib(s,x,a01,b01,bZ01)*S.weib(s,x,a02,b02,bZ02)*iweibull(x,a01,b01,bZ01)},s,max)
list(LE.00=ET0dot$value,
LE.0.=ET01$value+ET0dot$value,
LE.01=ET01$value,
LE.11=ET12$value,
LTR=LTR$value)
}
lifexpect0.idmPl <- function(s,knots01,nknots01,the01,knots12,nknots12,the12,knots02,nknots02,the02,bZ01=0,bZ12=0,bZ02=0) {
ET12 = integrate(f=function(x) {
S.pl(s,x,knots12,nknots12,the12,bZ12)},s,knots12[nknots12+6])
ET0dot = integrate(f=function(x) {
S.pl(s,x,knots01,nknots01,the01,bZ01)*S.pl(s,x,knots02,nknots02,the02,bZ02) },s,knots02[nknots02+6])
ET01 = integrate(f=function(x) {
sapply(x,function(x) {integrate(f=function(y){
(S.pl(s,y,knots01,nknots01,the01,bZ01)
*S.pl(s,y,knots02,nknots02,the02,bZ02)*
intensity(times=y,knots=knots01,number.knots=nknots01,theta=the01,linear.predictor=bZ01)$intensity
*S.pl(y,x,knots12,nknots12,the12,bZ12))},
lower=s,upper=x)$value})},s,knots01[nknots01+6])
LTR=integrate(f = function(x) {S.pl(s, x, knots01, nknots01, the01, bZ01) * S.pl(s, x, knots02, nknots02, the02, bZ02) * intensity(times = x, knots = knots01, number.knots = nknots01, theta = the01,linear.predictor = bZ01)$intensity }, lower = s, upper = knots01[nknots01+6])$value
list(LE.00=ET0dot$value,
LE.0.=ET01$value+ET0dot$value,
LE.01=ET01$value,
LE.11=ET12$value,
LTR=LTR)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.