revdep/checks.noindex/riskRegression/new/riskRegression.Rcheck/riskRegression-Ex.R

In gbm-developers/gbm: Generalized Boosted Regression Models

pkgname <- "riskRegression"
source(file.path(R.home("share"), "R", "examples-header.R"))
options(warn = 1)
library('riskRegression')

base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
cleanEx()
nameEx("CSC")
### * CSC

flush(stderr()); flush(stdout())

### Name: CSC
### Title: Cause-specific Cox proportional hazard regression
### Aliases: CSC
### Keywords: survival

### ** Examples


library(prodlim)
library(survival)
data(Melanoma)
## fit two cause-specific Cox models
## different formula for the two causes
fit1 <- CSC(list(Hist(time,status)~sex+age,Hist(time,status)~invasion+epicel+log(thick)),
            data=Melanoma)
print(fit1)
## Not run: 
##D library(Publish)
##D publish(fit1)
## End(Not run)

## model hazard of all cause mortality instead of hazard of type 2 
fit1a <- CSC(list(Hist(time,status)~sex+age,Hist(time,status)~invasion+epicel+log(thick)),
             data=Melanoma,
             surv.type="surv")

## the predicted probabilities are similar 
plot(predictRisk(fit1,times=500,cause=1,newdata=Melanoma),
     predictRisk(fit1a,times=500,cause=1,newdata=Melanoma))

## special case where cause 2 has no covariates
fit1b <- CSC(list(Hist(time,status)~sex+age,Hist(time,status)~1),
             data=Melanoma)
print(fit1b)
predict(fit1b,cause=1,times=100,newdata=Melanoma)


## same formula for both causes
fit2 <- CSC(Hist(time,status)~invasion+epicel+age,
            data=Melanoma)
print(fit2)

## combine a cause-specific Cox regression model for cause 2 
## and a Cox regression model for the event-free survival:
## different formula for cause 2 and event-free survival
fit3 <- CSC(list(Hist(time,status)~sex+invasion+epicel+age,
                 Hist(time,status)~invasion+epicel+age),
            surv.type="surv",
            data=Melanoma)
print(fit3)

## same formula for both causes
fit4 <- CSC(Hist(time,status)~invasion+epicel+age,
            data=Melanoma,
            surv.type="surv")
print(fit4)

## strata
fit5 <- CSC(Hist(time,status)~invasion+epicel+age+strata(sex),
            data=Melanoma,
            surv.type="surv")
print(fit5)

## sanity checks

cox1 <- coxph(Surv(time,status==1)~invasion+epicel+age+strata(sex),data=Melanoma)
cox2 <- coxph(Surv(time,status!=0)~invasion+epicel+age+strata(sex),data=Melanoma)
all.equal(coef(cox1),coef(fit5$models[[1]]))
all.equal(coef(cox2),coef(fit5$models[[2]]))

## predictions
##
## surv.type = "hazard": predictions for both causes can be extracted
## from the same fit
fit2 <- CSC(Hist(time,status)~invasion+epicel+age, data=Melanoma)
predict(fit2,cause=1,newdata=Melanoma[c(17,99,108),],times=c(100,1000,10000))
predictRisk(fit2,cause=1,newdata=Melanoma[c(17,99,108),],times=c(100,1000,10000))
predictRisk(fit2,cause=2,newdata=Melanoma[c(17,99,108),],times=c(100,1000,10000))
predict(fit2,cause=1,newdata=Melanoma[c(17,99,108),],times=c(100,1000,10000))
predict(fit2,cause=2,newdata=Melanoma[c(17,99,108),],times=c(100,1000,10000))

## surv.type = "surv" we need to change the cause of interest 
library(survival)
fit5.2 <- CSC(Hist(time,status)~invasion+epicel+age+strata(sex),
            data=Melanoma,
            surv.type="surv",cause=2)
## now this does not work
try(predictRisk(fit5.2,cause=1,newdata=Melanoma,times=4))

## but this does
predictRisk(fit5.2,cause=2,newdata=Melanoma,times=100)
predict(fit5.2,cause=2,newdata=Melanoma,times=100)
predict(fit5.2,cause=2,newdata=Melanoma[4,],times=100)




cleanEx()
nameEx("Ctree")
### * Ctree

flush(stderr()); flush(stdout())

### Name: Ctree
### Title: S3-Wrapper for ctree.
### Aliases: Ctree

### ** Examples

library(prodlim)
library(party)
library(survival)
set.seed(50)
d <- SimSurv(50)
nd <- data.frame(X1=c(0,1,0),X2=c(-1,0,1))
f <- Ctree(Surv(time,status)~X1+X2,data=d)
predictRisk(f,newdata=nd,times=c(3,8))




cleanEx()
nameEx("FGR")
### * FGR

flush(stderr()); flush(stdout())

### Name: FGR
### Title: Formula wrapper for crr from cmprsk
### Aliases: FGR
### Keywords: survival

### ** Examples


library(prodlim)
library(survival)
library(cmprsk)
library(lava)
d <- prodlim::SimCompRisk(100)
f1 <- FGR(Hist(time,cause)~X1+X2,data=d)
print(f1)

## crr allows that some covariates are multiplied by
## a function of time (see argument tf of crr)
## by FGR uses the identity matrix
f2 <- FGR(Hist(time,cause)~cov2(X1)+X2,data=d)
print(f2)

## same thing, but more explicit:
f3 <- FGR(Hist(time,cause)~cov2(X1)+cov1(X2),data=d)
print(f3)

## both variables can enter cov2:
f4 <- FGR(Hist(time,cause)~cov2(X1)+cov2(X2),data=d)
print(f4)

## change the function of time
qFun <- function(x){x^2}
noFun <- function(x){x}
sqFun <- function(x){x^0.5}

## multiply X1 by time^2 and X2 by time:
f5 <- FGR(Hist(time,cause)~cov2(X1,tf=qFun)+cov2(X2),data=d)
print(f5)
print(f5$crrFit)
## same results as crr
with(d,crr(ftime=time,
           fstatus=cause,
           cov2=d[,c("X1","X2")],
           tf=function(time){cbind(qFun(time),time)}))

## still same result, but more explicit
f5a <- FGR(Hist(time,cause)~cov2(X1,tf=qFun)+cov2(X2,tf=noFun),data=d)
f5a$crrFit

## multiply X1 by time^2 and X2 by sqrt(time)
f5b <- FGR(Hist(time,cause)~cov2(X1,tf=qFun)+cov2(X2,tf=sqFun),data=d,cause=1)

## additional arguments for crr
f6<- FGR(Hist(time,cause)~X1+X2,data=d, cause=1,gtol=1e-5)
f6
f6a<- FGR(Hist(time,cause)~X1+X2,data=d, cause=1,gtol=0.1)
f6a



cleanEx()
nameEx("IPA")
### * IPA

flush(stderr()); flush(stdout())

### Name: IPA
### Title: Explained variation for settings with binary, survival and
###   competing risk outcome
### Aliases: IPA rsquared rsquared.default rsquared.glm rsquared.coxph
###   rsquared.CauseSpecificCox IPA.default IPA.glm IPA.coxph
###   IPA.CauseSpecificCox

### ** Examples

library(prodlim)
library(data.table)
# binary outcome
library(lava)
set.seed(18)
learndat <- sampleData(48,outcome="binary")
lr1 = glm(Y~X1+X2+X7+X9,data=learndat,family=binomial)
IPA(lr1)

## validation data
valdat=sampleData(94,outcome="binary")
IPA(lr1,newdata=valdat)

## predicted risks externally given
p1=predictRisk(lr1,newdata=valdat)
IPA(p1,formula=Y~1,valdat)

# survival
library(survival)
data(pbc)
pbc=na.omit(pbc)
pbctest=(1:NROW(pbc)) %in% sample(1:NROW(pbc),size=.632*NROW(pbc))
pbclearn=pbc[pbctest,]
cox1= coxph(Surv(time,status!=0)~age+sex+log(bili)+log(albumin)+log(protime),
      data=pbclearn,x=TRUE)

## same data
IPA(cox1,formula=Surv(time,status!=0)~1,times=1000)

## validation data
pbcval=pbc[!pbctest,]
IPA(cox1,formula=Surv(time,status!=0)~1,newdata=pbcval,times=1000)

## predicted risks externally given
p2=predictRisk(cox1,newdata=pbcval,times=1000)
IPA(cox1,formula=Surv(time,status!=0)~1,newdata=pbcval,times=1000)
 
# competing risks
data(Melanoma)
Melanomatest=(1:NROW(Melanoma)) %in% sample(1:NROW(Melanoma),size=.632*NROW(Melanoma))
Melanomalearn=Melanoma[Melanomatest,]
fit1 <- CSC(list(Hist(time,status)~sex,
                 Hist(time,status)~invasion+epicel+age),
                 data=Melanoma)
IPA(fit1,times=1000,cause=2)

## validation data
Melanomaval=Melanoma[!Melanomatest,]
IPA(fit1,formula=Hist(time,status)~1,newdata=Melanomaval,times=1000)

## predicted risks externally given
p3= predictRisk(fit1,cause=1,newdata=Melanomaval,times=1000)
IPA(p3,formula=Hist(time,status)~1,cause=1,newdata=Melanomaval,times=1000)
 



cleanEx()
nameEx("Melanoma")
### * Melanoma

flush(stderr()); flush(stdout())

### Name: Melanoma
### Title: Malignant melanoma data
### Aliases: Melanoma
### Keywords: datasets

### ** Examples


data(Melanoma)



cleanEx()
nameEx("Paquid")
### * Paquid

flush(stderr()); flush(stdout())

### Name: Paquid
### Title: Paquid sample
### Aliases: Paquid
### Keywords: datasets

### ** Examples

data(Paquid)



cleanEx()
nameEx("Score.list")
### * Score.list

flush(stderr()); flush(stdout())

### Name: Score.list
### Title: Score risk predictions
### Aliases: Score.list Score

### ** Examples

# binary outcome
library(lava)
set.seed(18)
learndat <- sampleData(48,outcome="binary")
testdat <- sampleData(40,outcome="binary")

## score logistic regression models
lr1 = glm(Y~X1+X2+X7+X9,data=learndat,family=binomial)
lr2 = glm(Y~X3+X5,data=learndat,family=binomial)
Score(list("LR(X1+X2+X7+X9)"=lr1,"LR(X3+X5)"=lr2),formula=Y~1,data=testdat)

## ROC curve and calibration plot
xb=Score(list("LR(X1+X2+X7+X9)"=lr1,"LR(X3+X5+X6)"=lr2),formula=Y~1,
         data=testdat,plots=c("calibration","ROC"))
## Not run: 
##D plotROC(xb)
##D plotCalibration(xb)
## End(Not run)

## compute AUC for a list of continuous markers
markers = as.list(testdat[,.(X6,X7,X8,X9,X10)])
Score(markers,formula=Y~1,data=testdat,metrics=c("auc"))

# cross-validation
## Not run: 
##D     learndat=sampleData(400,outcome="binary")
##D     lr1a = glm(Y~X6,data=learndat,family=binomial)
##D     lr2a = glm(Y~X7+X8+X9,data=learndat,family=binomial)
##D     ## bootstrap cross-validation
##D     x1=Score(list("LR1"=lr1a,"LR2"=lr2a),formula=Y~1,data=learndat,split.method="bootcv",B=100)
##D     x1
##D     ## leave-one-out and leave-pair-out bootstrap
##D     x2=Score(list("LR1"=lr1a,"LR2"=lr2a),formula=Y~1,data=learndat,
##D              split.method="loob",
##D              B=100,plots="calibration")
##D     x2
## End(Not run)
# survival outcome

# Score Cox regression models
## Not run: 
##D library(survival)
##D library(rms)
##D library(prodlim)
##D set.seed(18)
##D trainSurv <- sampleData(100,outcome="survival")
##D testSurv <- sampleData(40,outcome="survival")
##D cox1 = coxph(Surv(time,event)~X1+X2+X7+X9,data=trainSurv, y=TRUE, x = TRUE)
##D cox2 = coxph(Surv(time,event)~X3+X5+X6,data=trainSurv, y=TRUE, x = TRUE)
##D xs=Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##D          formula=Surv(time,event)~1,data=testSurv,conf.int=FALSE,times=c(5,8))
##D xs
## End(Not run)

# Integrated Brier score
## Not run: 
##D xs=Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##D          formula=Surv(time,event)~1,data=testSurv,conf.int=FALSE,
##D          summary="ibs",
##D          times=sort(unique(testSurv$time)))
## End(Not run)

# time-dependent AUC for list of markers
## Not run: 
##D survmarkers = as.list(testSurv[,.(X6,X7,X8,X9,X10)])
##D Score(survmarkers,
##D       formula=Surv(time,event)~1,metrics="auc",data=testSurv,
##D       conf.int=TRUE,times=c(5,8))
##D 
##D # compare models on test data
##D Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##D       formula=Surv(time,event)~1,data=testSurv,conf.int=TRUE,times=c(5,8))
## End(Not run)
# crossvalidation models in traindata
## Not run: 
##D     library(survival)
##D     set.seed(18)
##D     trainSurv <- sampleData(400,outcome="survival")
##D     cox1 = coxph(Surv(time,event)~X1+X2+X7+X9,data=trainSurv, y=TRUE, x = TRUE)
##D     cox2 = coxph(Surv(time,event)~X3+X5+X6,data=trainSurv, y=TRUE, x = TRUE)
##D     x1 = Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##D                formula=Surv(time,event)~1,data=trainSurv,conf.int=TRUE,times=c(5,8),
##D                split.method="loob",B=100,plots="calibration")
##D 
##D     x2= Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##D               formula=Surv(time,event)~1,data=trainSurv,conf.int=TRUE,times=c(5,8),
##D               split.method="bootcv",B=100)
## End(Not run)

# restrict number of comparisons
## Not run: 
##D     Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##D           formula=Surv(time,event)~1,data=trainSurv,contrasts=TRUE,
##D           null.model=FALSE,conf.int=TRUE,times=c(5,8),split.method="bootcv",B=3)
##D 
##D     # competing risks outcome
##D     set.seed(18)
##D     trainCR <- sampleData(40,outcome="competing.risks")
##D     testCR <- sampleData(40,outcome="competing.risks")
##D     library(riskRegression)
##D     library(cmprsk)
##D     # Cause-specific Cox regression
##D     csc1 = CSC(Hist(time,event)~X1+X2+X7+X9,data=trainCR)
##D     csc2 = CSC(Hist(time,event)~X3+X5+X6,data=trainCR)
##D     # Fine-Gray regression
##D     fgr1 = FGR(Hist(time,event)~X1+X2+X7+X9,data=trainCR,cause=1)
##D     fgr2 = FGR(Hist(time,event)~X3+X5+X6,data=trainCR,cause=1)
##D     Score(list("CSC(X1+X2+X7+X9)"=csc1,"CSC(X3+X5+X6)"=csc2,
##D                "FGR(X1+X2+X7+X9)"=fgr1,"FGR(X3+X5+X6)"=fgr2),
##D           formula=Hist(time,event)~1,data=testCR,se.fit=1L,times=c(5,8))
## End(Not run)



## Not run: 
##D     # reproduce some results of Table IV of Blanche et al. Stat Med 2013
##D     data(Paquid)
##D     ResPaquid <- Score(list("DSST"=-Paquid$DSST,"MMSE"=-Paquid$MMSE),
##D                        formula=Hist(time,status)~1,
##D                        data=Paquid,
##D                        null.model = FALSE,
##D                        conf.int=TRUE,
##D                        metrics=c("auc"),
##D                        times=c(3,5,10),
##D                        plots="ROC")
##D     ResPaquid
##D     plotROC(ResPaquid,time=5)
## End(Not run)



cleanEx()
nameEx("SuperPredictor")
### * SuperPredictor

flush(stderr()); flush(stdout())

### Name: SuperPredictor
### Title: Formula interface for SuperLearner::SuperLearner
### Aliases: SuperPredictor

### ** Examples

## Not run: 
##D library(SuperLearner)
##D library(data.table)
##D d = sampleData(338, outcome="binary")
##D spfit = SuperPredictor(Y~X1+X2+X3+X4+X5+X6+X7+X8+X9+X10,data=d)
##D predictRisk(spfit)
##D x <- Score(list(spfit),data=d,formula=Y~1)
## End(Not run)



cleanEx()
nameEx("SurvResponseVar")
### * SurvResponseVar

flush(stderr()); flush(stdout())

### Name: SurvResponseVar
### Title: Extract the time and event variable from a Cox model
### Aliases: SurvResponseVar

### ** Examples

## Not run: 
##D SurvResponseVar(Surv(time,event)~X1+X2)
##D SurvResponseVar(Hist(time,event==0)~X1+X2)
##D SurvResponseVar(Surv(start,time, status,type="counting") ~ X3+X5)
##D SurvResponseVar(Surv(start,event=status, time2=time,type="counting") ~ X3+X5)
##D 
##D SurvResponseVar(survival::Surv(start,event=status, time2=time,type="counting") ~ X3+X5)
##D SurvResponseVar(status ~ X3+X5)
##D SurvResponseVar(I(status == 1) ~ X3+X5)
##D SurvResponseVar(list(Hist(time, event) ~ X1+X6,Hist(time, event) ~ X6))
## End(Not run)



cleanEx()
nameEx("ate")
### * ate

flush(stderr()); flush(stdout())

### Name: ate
### Title: Compute the Average Treatment Effects Via
### Aliases: ate

### ** Examples

library(survival)
library(rms)
library(prodlim)
set.seed(10)

#### Survival settings  ####
#### ATE with Cox model ####

## generate data
n <- 100
dtS <- sampleData(n, outcome="survival")
dtS$time <- round(dtS$time,1)
dtS$X1 <- factor(rbinom(n, prob = c(0.3,0.4) , size = 2), labels = paste0("T",0:2))

## estimate the Cox model
fit <- cph(formula = Surv(time,event)~ X1+X2,data=dtS,y=TRUE,x=TRUE)

## compute the ATE at times 5, 6, 7, and 8 using X1 as the treatment variable
## Not run: 
##D ## only point estimate (argument se = FALSE)
##D ateFit1a <- ate(fit, data = dtS, treatment = "X1", times = 5:8,
##D                se = FALSE)
##D 
##D ## standard error / confidence intervals computed using the influence function
##D ## (argument se = TRUE and B = 0)
##D ateFit1b <- ate(fit, data = dtS, treatment = "X1", times = 5:8,
##D                se = TRUE, B = 0)
##D 
##D ## same as before with in addition the confidence bands for the ATE
##D ## (argument band = TRUE)
##D ateFit1c <- ate(fit, data = dtS, treatment = "X1", times = 5:8,
##D                se = TRUE, band = TRUE, B = 0)
##D 
##D ## standard error / confidence intervals computed using 100 boostrap samples
##D ## (argument se = TRUE and B = 100) 
##D ateFit1d <- ate(fit, data = dtS, treatment = "X1",
##D                 times = 5:8, se = TRUE, B = 100)
##D ## NOTE: for real applications 100 bootstrap samples is not enougth 
##D 
##D ## same but using 2 cpus for generating and analyzing the boostrap samples
##D ## (parallel computation, argument mc.cores = 2) 
##D ateFit1e <- ate(fit, data = dtS, treatment = "X1",
##D                 times = 5:8, se = TRUE, B = 100, mc.cores = 2)
## End(Not run)

#### Survival settings without censoring ####
#### ATE with glm                        ####

## generate data
n <- 100
dtB <- sampleData(n, outcome="binary")
dtB[, X2 := as.numeric(X2)]

## estimate a logistic regression model
fit <- glm(formula = Y ~ X1+X2, data=dtB, family = "binomial")

## compute the ATE using X1 as the treatment variable
## only point estimate (argument se = FALSE)
ateFit1a <- ate(fit, data = dtB, treatment = "X1", se = FALSE)

## Not run: 
##D ## standard error / confidence intervals computed using the influence function
##D ateFit1b <- ate(fit, data = dtB, treatment = "X1",
##D                times = 5, ## just for having a nice output not used in computations
##D                se = TRUE, B = 0)
##D 
##D ## standard error / confidence intervals computed using 100 boostrap samples
##D ateFit1d <- ate(fit, data = dtB, treatment = "X1",
##D                 times = 5, se = TRUE, B = 100)
##D 
##D ## using the lava package
##D ateLava <- estimate(fit, function(p, data){
##D a <- p["(Intercept)"] ; b <- p["X11"] ; c <- p["X2"] ;
##D R.X11 <- expit(a + b + c * data[["X2"]])
##D R.X10 <- expit(a + c * data[["X2"]])
##D list(risk0=R.X10,risk1=R.X11,riskdiff=R.X11-R.X10)},
##D average=TRUE)
##D ateLava
##D 
##D ateFit1b$meanRisk
## End(Not run)

#### Competing risks settings               ####
#### ATE with cause specific Cox regression ####

## Not run: 
##D ## generate data
##D n <- 500
##D set.seed(10)
##D dt <- sampleData(n, outcome="competing.risks")
##D dt$time <- round(dt$time,1)
##D dt$X1 <- factor(rbinom(n, prob = c(0.2,0.3) , size = 2), labels = paste0("T",0:2))
##D 
##D ## estimate cause specific Cox model
##D fitCR <-  CSC(Hist(time,event)~ X1+X8,data=dt,cause=1)
##D 
##D ## compute the ATE at times 10, 15, 20 using X1 as the treatment variable
##D ateFit2a <- ate(fitCR, data = dt, treatment = "X1", times = c(10,15,20),
##D                 cause = 1, se = FALSE)
##D 
##D ## standard error / confidence intervals computed using the influence function
##D ## (argument se = TRUE and B = 0)
##D ateFit2b <- ate(fitCR, data = dt, treatment = "X1", times = c(10,15,20),
##D                 cause = 1, se = TRUE, B = 0)
##D 
##D ## same as before with in addition the confidence bands for the ATE
##D ## (argument band = TRUE)
##D ateFit2c <- ate(fitCR, data = dt, treatment = "X1", times = c(10,15,20), 
##D                cause = 1, se = TRUE, band = TRUE, B = 0)
##D 
##D ## standard error / confidence intervals computed using 100 boostrap samples
##D ## (argument se = TRUE and B = 100) 
##D ateFit2d <- ate(fitCR, data = dt, treatment = "X1", times = c(10,15,20), 
##D                 cause = 1, se = TRUE, B = 100)
##D ## NOTE: for real applications 100 bootstrap samples is not enougth 
##D 
##D ## same but using 2 cpus for generating and analyzing the boostrap samples
##D ## (parallel computation, argument mc.cores = 2) 
##D ateFit2e <- ate(fitCR, data = dt, treatment = "X1", times = c(10,15,20), 
##D                 cause = 1, se = TRUE, B = 100, mc.cores = 2)
## End(Not run)

#### time-dependent covariates ###
## Not run: 
##D library(survival)
##D fit <- coxph(Surv(time, status) ~ celltype+karno + age + trt, veteran)
##D vet2 <- survSplit(Surv(time, status) ~., veteran,
##D                        cut=c(60, 120), episode ="timegroup")
##D fitTD <- coxph(Surv(tstart, time, status) ~ celltype+karno + age + trt,
##D                data= vet2,x=1)
##D set.seed(16)
##D resVet <- ate(fitTD,formula=Hist(entry=tstart,time=time,event=status)~1,
##D           data = vet2, treatment = "celltype", contrasts = NULL,
##D         times=5,verbose=1,
##D         landmark = c(0,30,60,90), cause = 1, B = 10, se = 1,
##D         band = FALSE, mc.cores=1)
##D resVet
## End(Not run)

## Not run: 
##D set.seed(137)
##D d=sampleDataTD(127)
##D library(survival)
##D d[,status:=1*(event==1)]
##D d[,X3:=as.factor(X3)]
##D ## ignore competing risks
##D cox1TD <- coxph(Surv(start,time, status,type="counting") ~ X3+X5+X6+X8,
##D                 data=d, x = TRUE)
##D resTD1 <- ate(cox1TD,formula=Hist(entry=start,time=time,event=status)~1,
##D         data = d, treatment = "X3", contrasts = NULL,
##D         times=.5,verbose=1,
##D         landmark = c(0,0.5,1), B = 20, se = 1,
##D         band = FALSE, mc.cores=1)
##D resTD1
##D ## account for competing risks
##D cscTD <- CSC(Hist(time=time, event=event,entry=start) ~ X3+X5+X6+X8, data=d)
##D set.seed(16)
##D resTD <- ate(cscTD,formula=Hist(entry=start,time=time,event=event)~1,
##D         data = d, treatment = "X3", contrasts = NULL,
##D         times=.5,verbose=1,
##D         landmark = c(0,0.5,1), cause = 1, B = 20, se = 1,
##D         band = FALSE, mc.cores=1)
##D resTD
## End(Not run)



cleanEx()
nameEx("autoplot.Score")
### * autoplot.Score

flush(stderr()); flush(stdout())

### Name: autoplot.Score
### Title: ggplot AUC curve
### Aliases: autoplot.Score

### ** Examples

library(survival)
library(ggplot2)
d=sampleData(100,outcome="survival")
nd=sampleData(100,outcome="survival")
f1=coxph(Surv(time,event)~X1+X6+X8,data=d,x=TRUE,y=TRUE)
f2=coxph(Surv(time,event)~X2+X5+X9,data=d,x=TRUE,y=TRUE)
xx=Score(list(f1,f2), formula=Surv(time,event)~1,
data=nd, metrics="auc", null.model=FALSE, times=seq(3:10))
g <- autoplot(xx)
print(g)
aucgraph <- plotAUC(xx)
plotAUC(xx,conf.int=TRUE)
plotAUC(xx,which="contrasts")
plotAUC(xx,which="contrasts",conf.int=TRUE)





cleanEx()
nameEx("autoplot.ate")
### * autoplot.ate

flush(stderr()); flush(stdout())

### Name: autoplot.ate
### Title: Plot Average Risks
### Aliases: autoplot.ate

### ** Examples

## Not run: 
##D library(survival)
##D library(rms)
##D library(ggplot2)
##D #### simulate data ####
##D n <- 1e2
##D set.seed(10)
##D dtS <- sampleData(n,outcome="survival")
##D 
##D 
##D #### Cox model ####
##D fit <- cph(formula = Surv(time,event)~ X1+X2,data=dtS,y=TRUE,x=TRUE)
##D 
##D #### Average treatment effect ####
##D seqTimes <- sort(unique(fit$y[,1]))
##D seqTimes5 <- seqTimes[seqTimes>5 & seqTimes<10]
##D ateFit <- ate(fit, data = dtS, treatment = "X1", contrasts = NULL,
##D               times = seqTimes, B = 0, band = TRUE, nsim.band = 500, y = TRUE,
##D               mc.cores=1)
##D 
##D #### display #### 
##D ggplot2::autoplot(ateFit)
##D 
##D outGG <- autoplot(ateFit, band = TRUE, ci = TRUE, alpha = 0.1)
##D dd <- as.data.frame(outGG$data[treatment == 0])
##D outGG$plot + facet_wrap(~treatment, labeller = label_both)
## End(Not run)



cleanEx()
nameEx("autoplot.predictCSC")
### * autoplot.predictCSC

flush(stderr()); flush(stdout())

### Name: autoplot.predictCSC
### Title: Plot Predictions From a Cause-specific Cox Proportional Hazard
###   Regression
### Aliases: autoplot.predictCSC

### ** Examples

library(survival)
library(rms)
library(ggplot2)
library(prodlim)
#### simulate data ####
set.seed(10)
d <- sampleData(1e2, outcome = "competing.risks")

#### CSC model ####
m.CSC <- CSC(Hist(time,event)~ X1 + X2 + X6, data = d)

pred.CSC <- predict(m.CSC, newdata = d[1:2,], time = 1:5, cause = 1)#'
autoplot(pred.CSC)


#### stratified CSC model ####
m.SCSC <- CSC(Hist(time,event)~ strata(X1) + strata(X2) + X6,
              data = d)
pred.SCSC <- predict(m.SCSC, time = 1:3, newdata = d[1:4,],
                     cause = 1, keep.newdata = TRUE, keep.strata = TRUE)
autoplot(pred.SCSC, group.by = "strata")



cleanEx()
nameEx("autoplot.predictCox")
### * autoplot.predictCox

flush(stderr()); flush(stdout())

### Name: autoplot.predictCox
### Title: Plot Predictions From a Cox Model
### Aliases: autoplot.predictCox

### ** Examples

library(survival)
library(ggplot2)

#### simulate data ####
set.seed(10)
d <- sampleData(1e2, outcome = "survival")

#### Cox model ####
m.cox <- coxph(Surv(time,event)~ X1 + X2 + X3,
                data = d, x = TRUE, y = TRUE)

## display baseline hazard
e.basehaz <- predictCox(m.cox)

autoplot(e.basehaz, type = "cumhazard")

## display predicted survival
pred.cox <- predictCox(m.cox, newdata = d[1:4,],
  times = 1:5, type = "survival", keep.newdata = TRUE)
autoplot(pred.cox)
autoplot(pred.cox, group.by = "covariates")
autoplot(pred.cox, group.by = "covariates", reduce.data = TRUE)

## predictions with confidence interval/bands
pred.cox <- predictCox(m.cox, newdata = d[1,,drop=FALSE],
  times = 1:5, type = "survival", band = TRUE, se = TRUE, keep.newdata = TRUE)
autoplot(pred.cox, ci = TRUE, band = TRUE)
autoplot(pred.cox, ci = TRUE, band = TRUE, alpha = 0.1)

#### Stratified Cox model ####
m.cox.strata <- coxph(Surv(time,event)~ strata(X1) + strata(X2) + X3 + X6,
                      data = d, x = TRUE, y = TRUE)

pred.cox.strata <- predictCox(m.cox.strata, newdata = d[1:5,,drop=FALSE],
                              time = 1:5, keep.newdata = TRUE)

## display
res <- autoplot(pred.cox.strata, type = "survival", group.by = "strata")

## customize display
res$plot + facet_wrap(~strata, labeller = label_both)
res$plot %+% res$data[strata == "0, 1"]



cleanEx()
nameEx("boot2pvalue")
### * boot2pvalue

flush(stderr()); flush(stdout())

### Name: boot2pvalue
### Title: Compute the p.value from the distribution under H1
### Aliases: boot2pvalue

### ** Examples

set.seed(10)

#### no effect ####
x <- rnorm(1e3) 
boot2pvalue(x, null = 0, estimate = mean(x), alternative = "two.sided")
## expected value of 1
boot2pvalue(x, null = 0, estimate = mean(x), alternative = "greater")
## expected value of 0.5
boot2pvalue(x, null = 0, estimate = mean(x), alternative = "less")
## expected value of 0.5

#### positive effect ####
x <- rnorm(1e3, mean = 1) 
boot2pvalue(x, null = 0, estimate = 1, alternative = "two.sided")
## expected value of 0.32 = 2*pnorm(q = 0, mean = -1) = 2*mean(x<=0)
boot2pvalue(x, null = 0, estimate = 1, alternative = "greater")  
## expected value of 0.16 = pnorm(q = 0, mean = 1) = mean(x<=0)
boot2pvalue(x, null = 0, estimate = 1, alternative = "less")
## expected value of 0.84 = 1-pnorm(q = 0, mean = 1) = mean(x>=0)

#### negative effect ####
x <- rnorm(1e3, mean = -1) 
boot2pvalue(x, null = 0, estimate = -1, alternative = "two.sided") 
## expected value of 0.32 = 2*(1-pnorm(q = 0, mean = -1)) = 2*mean(x>=0)
boot2pvalue(x, null = 0, estimate = -1, alternative = "greater")
## expected value of 0.84 = pnorm(q = 0, mean = -1) = mean(x<=0)
boot2pvalue(x, null = 0, estimate = -1, alternative = "less") # pnorm(q = 0, mean = -1)
## expected value of 0.16 = 1-pnorm(q = 0, mean = -1) = mean(x>=0)



cleanEx()
nameEx("boxplot.Score")
### * boxplot.Score

flush(stderr()); flush(stdout())

### Name: boxplot.Score
### Title: Boxplot risk quantiles
### Aliases: boxplot.Score

### ** Examples

# binary outcome
library(data.table)
library(prodlim)
db=sampleData(40,outcome="binary")
fitconv=glm(Y~X3+X5,data=db,family=binomial)
fitnew=glm(Y~X1+X3+X5+X6+X7,data=db,family=binomial)
x=Score(list(new=fitnew,conv=fitconv),
        formula=Y~1,contrasts=list(c(2,1)),
               data=db,plots="box",null.model=FALSE)
boxplot(x)

# survival outcome
library(survival)
ds=sampleData(40,outcome="survival")
fit=coxph(Surv(time,event)~X6+X9,data=ds,x=TRUE,y=TRUE)
## Not run: 
##D  
##D scoreobj=Score(list("Cox"=fit),
##D                 formula=Hist(time,event)~1, data=ds,
##D                 metrics=NULL, plots="box",
##D                 times=c(1,5),null.model=FALSE)
##D boxplot(scoreobj,timepoint=5)
##D boxplot(scoreobj,timepoint=1)
##D 
## End(Not run)

# competing risks outcome
library(survival)
data(Melanoma, package = "riskRegression")
fit = CSC(Hist(time,event,cens.code="censored")~invasion+age+sex,data=Melanoma)
scoreobj=Score(list("CSC"=fit),
               formula=Hist(time,event,cens.code="censored")~1,
               data=Melanoma,plots="box",times=5*365.25,null.model=FALSE)
par(mar=c(4,12,4,4))
boxplot(scoreobj,timepoint=5*365.25)

# more than 2 competing risks
m=lava::lvm(~X1+X2+X3)
lava::distribution(m, "eventtime1") <- lava::coxWeibull.lvm(scale = 1/100)
lava::distribution(m, "eventtime2") <- lava::coxWeibull.lvm(scale = 1/100)
lava::distribution(m, "eventtime3") <- lava::coxWeibull.lvm(scale = 1/100)
lava::distribution(m, "censtime") <- lava::coxWeibull.lvm(scale = 1/100)
lava::regression(m,eventtime2~X3)=1.3
m <- lava::eventTime(m,
time ~ min(eventtime1 = 1, eventtime2 = 2, eventtime3 = 3, censtime = 0), "event")
set.seed(101)
dcr=as.data.table(lava::sim(m,101))
fit = CSC(Hist(time,event)~X1+X2+X3,data=dcr)
scoreobj=Score(list("my model"=fit),
               formula=Hist(time,event)~1,
               data=dcr,plots="box",times=5,null.model=FALSE)
boxplot(scoreobj)





graphics::par(get("par.postscript", pos = 'CheckExEnv'))
cleanEx()
nameEx("colCenter_cpp")
### * colCenter_cpp

flush(stderr()); flush(stdout())

### Name: colCenter_cpp
### Title: Apply - by column
### Aliases: colCenter_cpp

### ** Examples

x <- matrix(1,6,5)
sweep(x, MARGIN = 1, FUN = "-", STATS = 1:6)
colCenter_cpp(x, 1:6 )



cleanEx()
nameEx("colCumProd")
### * colCumProd

flush(stderr()); flush(stdout())

### Name: colCumProd
### Title: Apply cumprod in each column
### Aliases: colCumProd

### ** Examples

x <- matrix(1:8,ncol=2)
colCumProd(x)



cleanEx()
nameEx("colCumSum")
### * colCumSum

flush(stderr()); flush(stdout())

### Name: colCumSum
### Title: Apply cumsum in each column
### Aliases: colCumSum

### ** Examples

x <- matrix(1:8,ncol=2)
colCumSum(x)



cleanEx()
nameEx("colMultiply_cpp")
### * colMultiply_cpp

flush(stderr()); flush(stdout())

### Name: colMultiply_cpp
### Title: Apply * by column
### Aliases: colMultiply_cpp

### ** Examples

x <- matrix(1,6,5)
sweep(x, MARGIN = 1, FUN = "*", STATS = 1:6)
colMultiply_cpp(x, 1:6 )



cleanEx()
nameEx("colScale_cpp")
### * colScale_cpp

flush(stderr()); flush(stdout())

### Name: colScale_cpp
### Title: Apply / by column
### Aliases: colScale_cpp

### ** Examples

x <- matrix(1,6,5)
sweep(x, MARGIN = 1, FUN = "/", STATS = 1:6)
colScale_cpp(x, 1:6 )



cleanEx()
nameEx("colSumsCrossprod")
### * colSumsCrossprod

flush(stderr()); flush(stdout())

### Name: colSumsCrossprod
### Title: Apply crossprod and colSums
### Aliases: colSumsCrossprod

### ** Examples

x <- matrix(1:8,ncol=2)
y <- matrix(1:16,ncol=8)
colSumsCrossprod(x,y,0)

x <- matrix(1:8,ncol=2)
y <- matrix(1:16,ncol=2)
colSumsCrossprod(x,y,1)



cleanEx()
nameEx("confint.ate")
### * confint.ate

flush(stderr()); flush(stdout())

### Name: confint.ate
### Title: Confidence Intervals and Confidence Bands for the Predicted
###   Absolute Risk (Cumulative Incidence Function)
### Aliases: confint.ate

### ** Examples

library(survival)
library(data.table)

## ## generate data ####
set.seed(10)
d <- sampleData(70,outcome="survival")
d[, X1 := paste0("T",rbinom(.N, size = 2, prob = c(0.51)))]
## table(d$X1)

#### stratified Cox model ####
fit <- coxph(Surv(time,event)~X1 + strata(X2) + X6,
             data=d, ties="breslow", x = TRUE, y = TRUE)

#### average treatment effect ####
fit.ate <- ate(fit, treatment = "X1", times = 1:3, data = d,
               se = TRUE, iid = TRUE, band = TRUE)
print(fit.ate, type = "meanRisk")
dt.ate <- as.data.table(fit.ate)

## manual calculation of se
dd <- copy(d)
dd$X1 <- rep(factor("T0", levels = paste0("T",0:2)), NROW(dd))
out <- predictCox(fit, newdata = dd, se = TRUE, times = 1:3, average.iid = TRUE)
term1 <- -out$survival.average.iid
term2 <- sweep(1-out$survival, MARGIN = 2, FUN = "-", STATS = colMeans(1-out$survival))
sqrt(colSums((term1 + term2/NROW(d))^2)) 
## fit.ate$meanRisk[treatment=="T0",meanRisk.se]

## note
out2 <- predictCox(fit, newdata = dd, se = TRUE, times = 1:3, iid = TRUE)
mean(out2$survival.iid[,1,1])
out$survival.average.iid[1,1]

## check confidence intervals (no transformation)
dt.ate[,.(lower = pmax(0,value + qnorm(0.025) * se),
          lower2 = lower,
          upper = value + qnorm(0.975) * se,
          upper2 = upper)]

## add confidence intervals computed on the log-log scale
## and backtransformed
outCI <- confint(fit.ate,
                 meanRisk.transform = "loglog", diffRisk.transform = "atanh",
                 ratioRisk.transform = "log")
print(outCI, type = "meanRisk")

dt.ate[type == "ate", newse := se/(value*log(value))]
dt.ate[type == "ate", .(lower = exp(-exp(log(-log(value)) - 1.96 * newse)),
                        upper = exp(-exp(log(-log(value)) + 1.96 * newse)))]



cleanEx()
nameEx("confint.predictCSC")
### * confint.predictCSC

flush(stderr()); flush(stdout())

### Name: confint.predictCSC
### Title: Confidence Intervals and Confidence Bands for the Predicted
###   Absolute Risk (Cumulative Incidence Function)
### Aliases: confint.predictCSC

### ** Examples

library(survival)
library(prodlim)
#### generate data ####
set.seed(10)
d <- sampleData(100) 

#### estimate a stratified CSC model ###
fit <- CSC(Hist(time,event)~ X1 + strata(X2) + X6, data=d)

#### compute individual specific risks
fit.pred <- predict(fit, newdata=d[1:3], times=c(3,8), cause = 1,
                    se = TRUE, iid = TRUE, band = TRUE)
fit.pred

## check confidence intervals
newse <- fit.pred$absRisk.se/(-fit.pred$absRisk*log(fit.pred$absRisk))
cbind(lower = as.double(exp(-exp(log(-log(fit.pred$absRisk)) + 1.96 * newse))),
      upper = as.double(exp(-exp(log(-log(fit.pred$absRisk)) - 1.96 * newse)))
)

#### compute confidence intervals without transformation
confint(fit.pred, absRisk.transform = "none")
cbind(lower = as.double(fit.pred$absRisk - 1.96 * fit.pred$absRisk.se),
      upper = as.double(fit.pred$absRisk + 1.96 * fit.pred$absRisk.se)
)





cleanEx()
nameEx("confint.predictCox")
### * confint.predictCox

flush(stderr()); flush(stdout())

### Name: confint.predictCox
### Title: Confidence Intervals and Confidence Bands for the predicted
###   Survival/Cumulative Hazard
### Aliases: confint.predictCox

### ** Examples

library(survival)

#### generate data ####
set.seed(10)
d <- sampleData(40,outcome="survival") 

#### estimate a stratified Cox model ####
fit <- coxph(Surv(time,event)~X1 + strata(X2) + X6,
             data=d, ties="breslow", x = TRUE, y = TRUE)

#### compute individual specific survival probabilities  
fit.pred <- predictCox(fit, newdata=d[1:3], times=c(3,8), type = "survival",
                       se = TRUE, iid = TRUE, band = TRUE)
fit.pred

## check standard error
sqrt(rowSums(fit.pred$survival.iid[1,,]^2)) ## se for individual 1

## check confidence interval
newse <- fit.pred$survival.se/(-fit.pred$survival*log(fit.pred$survival))
cbind(lower = as.double(exp(-exp(log(-log(fit.pred$survival)) + 1.96 * newse))),
      upper = as.double(exp(-exp(log(-log(fit.pred$survival)) - 1.96 * newse)))
)

#### compute confidence intervals without transformation
confint(fit.pred, survival.transform = "none")
cbind(lower = as.double(fit.pred$survival - 1.96 * fit.pred$survival.se),
      upper = as.double(fit.pred$survival + 1.96 * fit.pred$survival.se)
)




cleanEx()
nameEx("getSplitMethod")
### * getSplitMethod

flush(stderr()); flush(stdout())

### Name: getSplitMethod
### Title: Input for data splitting algorithms
### Aliases: getSplitMethod

### ** Examples

# 3-fold crossvalidation
getSplitMethod("cv3",B=4,N=37)

# bootstrap with replacement
getSplitMethod("loob",B=4,N=37)

# bootstrap without replacement
getSplitMethod("loob",B=4,N=37,M=20)




cleanEx()
nameEx("iidCox")
### * iidCox

flush(stderr()); flush(stdout())

### Name: iidCox
### Title: Extract iid decomposition from a Cox model
### Aliases: iidCox iidCox.coxph iidCox.cph iidCox.phreg
###   iidCox.CauseSpecificCox

### ** Examples

library(survival)
library(data.table)
library(prodlim)
set.seed(10)
d <- sampleData(100, outcome = "survival")[,.(eventtime,event,X1,X6)]
setkey(d, eventtime)

m.cox <- coxph(Surv(eventtime, event) ~ X1+X6, data = d, y = TRUE, x = TRUE)
system.time(IF.cox <- iidCox(m.cox))
system.time(IF.cox_approx <- iidCox(m.cox, store.iid = "approx"))


IF.cox.all <- iidCox(m.cox, tau.hazard = sort(unique(c(7,d$eventtime))))
IF.cox.beta <- iidCox(m.cox, baseline.iid = FALSE)




cleanEx()
nameEx("influenceTest")
### * influenceTest

flush(stderr()); flush(stdout())

### Name: influenceTest
### Title: Influence test [Experimental!!]
### Aliases: influenceTest influenceTest.list influenceTest.default

### ** Examples

library(lava)
library(survival)
library(prodlim)
library(data.table)
n <- 100

#### Under H1
set.seed(1)
newdata <- data.frame(X1=0:1)

## simulate non proportional hazard using lava
m <- lvm()
regression(m) <- y ~ 1
regression(m) <- s ~ exp(-2*X1)
distribution(m,~X1) <- binomial.lvm()
distribution(m,~cens) <- coxWeibull.lvm(scale=1)
distribution(m,~y) <- coxWeibull.lvm(scale=1,shape=~s)
eventTime(m) <- eventtime ~ min(y=1,cens=0)
d <- as.data.table(sim(m,n))
setkey(d, eventtime)

## fit cox models
m.cox <- coxph(Surv(eventtime, status) ~ X1, 
               data = d, y = TRUE, x = TRUE)

mStrata.cox <- coxph(Surv(eventtime, status) ~ strata(X1), 
                     data = d, y = TRUE, x = TRUE)

## compare models
# one time point
outIF <- influenceTest(list(m.cox, mStrata.cox), 
              type = "survival", newdata = newdata, times = 0.5)
confint(outIF)
                                 
# several timepoints
outIF <- influenceTest(list(m.cox, mStrata.cox), 
              type = "survival", newdata = newdata, times = c(0.5,1,1.5))
confint(outIF)

#### Under H0 (Cox) ####
set.seed(1)
## simulate proportional hazard using lava
m <- lvm()
regression(m) <- y ~ 1
distribution(m,~X1) <- binomial.lvm()
distribution(m,~cens) <- coxWeibull.lvm()
distribution(m,~y) <- coxWeibull.lvm()
eventTime(m) <- eventtime ~ min(y=1,cens=0)
d <- as.data.table(sim(m,n))
setkey(d, eventtime)

## fit cox models
Utime <- sort(unique(d$eventtime))
m.cox <- coxph(Surv(eventtime, status) ~ X1, 
               data = d, y = TRUE, x = TRUE)

mStrata.cox <- coxph(Surv(eventtime, status) ~ strata(X1), 
                     data = d, y = TRUE, x = TRUE)

p.cox <- predictCox(m.cox, newdata = newdata, time = Utime, type = "survival")
p.coxStrata <- predictCox(mStrata.cox, newdata = newdata, time = Utime, type = "survival")

## display
library(ggplot2)
autoplot(p.cox)
autoplot(p.coxStrata)
 
## compare models
outIF <- influenceTest(list(m.cox, mStrata.cox), 
                       type = "survival", newdata = newdata, times = Utime[1:6])
confint(outIF)

#### Under H0 (CSC) ####
set.seed(1)
ff <- ~ f(X1,2) + f(X2,-0.033)
ff <- update(ff, ~ .+ f(X3,0) + f(X4,0) + f(X5,0))
ff <- update(ff, ~ .+ f(X6,0) + f(X7,0) + f(X8,0) + f(X9,0))
d <- sampleData(n, outcome = "competing.risk", formula = ff)
d[,X1:=as.numeric(as.character(X1))]
d[,X2:=as.numeric(as.character(X2))]
d[,X3:=as.numeric(as.character(X3))]
d[,X4:=as.numeric(as.character(X4))]
d[,X5:=as.numeric(as.character(X5))]
setkey(d, time)

Utime <- sort(unique(d$time))

## fit cox models
m.CSC <- CSC(Hist(time, event) ~ X1 + X2, data = d)
mStrata.CSC <- CSC(Hist(time, event) ~ strata(X1) + X2 + X3, data = d)

## compare models
outIF <- influenceTest(list(m.CSC, mStrata.CSC), 
             cause = 1, newdata = unique(d[,.(X1,X2,X3)]), times = Utime[1:5])
confint(outIF)



cleanEx()
nameEx("ipcw")
### * ipcw

flush(stderr()); flush(stdout())

### Name: ipcw
### Title: Estimation of censoring probabilities
### Aliases: ipcw ipcw.none ipcw.marginal ipcw.nonpar ipcw.cox ipcw.aalen
### Keywords: survival

### ** Examples


library(prodlim)
library(rms)
dat=SimSurv(30)

dat <- dat[order(dat$time),]

# using the marginal Kaplan-Meier for the censoring times

WKM=ipcw(Hist(time,status)~X2,
  data=dat,
  method="marginal",
  times=sort(unique(dat$time)),
  subject.times=dat$time,keep=c("fit"))
plot(WKM$fit)
WKM$fit

# using the Cox model for the censoring times given X2
library(survival)
WCox=ipcw(Hist(time=time,event=status)~X2,
  data=dat,
  method="cox",
  times=sort(unique(dat$time)),
  subject.times=dat$time,keep=c("fit"))
WCox$fit

plot(WKM$fit)
lines(sort(unique(dat$time)),
      1-WCox$IPCW.times[1,],
      type="l",
      col=2,
      lty=3,
      lwd=3)
lines(sort(unique(dat$time)),
      1-WCox$IPCW.times[5,],
      type="l",
      col=3,
      lty=3,
      lwd=3)

# using the stratified Kaplan-Meier
# for the censoring times given X2

WKM2=ipcw(Hist(time,status)~X2,
  data=dat,
  method="nonpar",
  times=sort(unique(dat$time)),
  subject.times=dat$time,keep=c("fit"))
plot(WKM2$fit,add=FALSE)





cleanEx()
nameEx("penalizedS3")
### * penalizedS3

flush(stderr()); flush(stdout())

### Name: penalizedS3
### Title: S3-wrapper for S4 function penalized
### Aliases: penalizedS3

### ** Examples

library(prodlim)
## Not run: 
##D ## too slow
##D library(penalized)
##D set.seed(8)
##D d <- sampleData(200,outcome="binary")
##D newd <- sampleData(80,outcome="binary")
##D fitridge <- penalizedS3(Y~X1+X2+pen(7:8), data=d, type="ridge",
##D standardize=TRUE, model="logistic",trace=FALSE)
##D fitlasso <- penalizedS3(Y~X1+X2+pen(7:8), data=d, type="lasso",
##D standardize=TRUE, model="logistic",trace=FALSE)
##D # fitnet <- penalizedS3(Y~X1+X2+pen(7:8), data=d, type="elastic.net",
##D # standardize=TRUE, model="logistic",trace=FALSE)
##D predictRisk(fitridge,newdata=newd)
##D predictRisk(fitlasso,newdata=newd)
##D # predictRisk(fitnet,newdata=newd)
##D Score(list(fitridge),data=newd,formula=Y~1)
##D Score(list(fitridge),data=newd,formula=Y~1,split.method="bootcv",B=2)
## End(Not run)
## Not run: 
##D  data(nki70) ## S4 fit
##D pen <- penalized(Surv(time, event), penalized = nki70[,8:77],
##D                  unpenalized = ~ER+Age+Diam+N+Grade, data = nki70,
##D lambda1 = 1)
##D penS3 <- penalizedS3(Surv(time,event)~ER+Age+Diam+pen(8:77)+N+Grade,
##D                      data=nki70, lambda1=1)
##D ## or
##D penS3 <- penalizedS3(Surv(time,event)~ER+pen(TSPYL5,Contig63649_RC)+pen(10:77)+N+Grade,
##D                      data=nki70, lambda1=1)
##D ## also this works
##D penS3 <- penalizedS3(Surv(time,event)~ER+Age+pen(8:33)+Diam+pen(34:77)+N+Grade,
##D                     data=nki70, lambda1=1)
## End(Not run)



cleanEx()
nameEx("plot.riskRegression")
### * plot.riskRegression

flush(stderr()); flush(stdout())

### Name: plot.riskRegression
### Title: Plotting predicted risk
### Aliases: plot.riskRegression
### Keywords: survival

### ** Examples


library(survival)
library(prodlim)
data(Melanoma)
fit.arr <- ARR(Hist(time,status)~invasion+age+strata(sex),data=Melanoma,cause=1)
plot(fit.arr,xlim=c(500,3000))





cleanEx()
nameEx("plotAUC")
### * plotAUC

flush(stderr()); flush(stdout())

### Name: plotAUC
### Title: Plot of time-dependent AUC curves
### Aliases: plotAUC

### ** Examples

library(survival)
library(prodlim)
d=sampleData(100,outcome="survival")
nd=sampleData(100,outcome="survival")
f1=coxph(Surv(time,event)~X1+X6+X8,data=d,x=TRUE,y=TRUE)
f2=coxph(Surv(time,event)~X2+X5+X9,data=d,x=TRUE,y=TRUE)
xx=Score(list("X1+X6+X8"=f1,"X2+X5+X9"=f2), formula=Surv(time,event)~1,
data=nd, metrics="auc", null.model=FALSE, times=seq(3:10))
aucgraph <- plotAUC(xx)
plotAUC(xx,conf.int=TRUE)
## difference between 
plotAUC(xx,which="contrasts",conf.int=TRUE)





cleanEx()
nameEx("plotBrier")
### * plotBrier

flush(stderr()); flush(stdout())

### Name: plotBrier
### Title: Plot Brier curve
### Aliases: plotBrier

### ** Examples

# survival
library(survival)
library(prodlim)
ds1=sampleData(40,outcome="survival")
ds2=sampleData(40,outcome="survival")
f1 <- coxph(Surv(time,event)~X1+X3+X5+X7+X9,data=ds1,x=TRUE)
f2 <- coxph(Surv(time,event)~X2+X4+X6+X8+X10,data=ds1,x=TRUE)
xscore <- Score(list(f1,f2),formula=Hist(time,event)~1,data=ds2,times=0:12,metrics="brier")
plotBrier(xscore)



cleanEx()
nameEx("plotCalibration")
### * plotCalibration

flush(stderr()); flush(stdout())

### Name: plotCalibration
### Title: Plot Calibration curve
### Aliases: plotCalibration

### ** Examples

library(prodlim)
# binary 
db=sampleData(100,outcome="binary")
fb1=glm(Y~X1+X5+X7,data=db,family="binomial")
fb2=glm(Y~X1+X3+X6+X7,data=db,family="binomial")
xb=Score(list(model1=fb1,model2=fb2),Y~1,data=db,
          plots="cal")
plotCalibration(xb,brier.in.legend=TRUE)
plotCalibration(xb,bars=TRUE,model="model1")
plotCalibration(xb,models=1,bars=TRUE,names.cex=1.3)

# survival
library(survival)
library(prodlim)
dslearn=sampleData(56,outcome="survival")
dstest=sampleData(100,outcome="survival")
fs1=coxph(Surv(time,event)~X1+X5+X7,data=dslearn,x=1)
fs2=coxph(Surv(time,event)~strata(X1)+X3+X6+X7,data=dslearn,x=1)
xs=Score(list(Cox1=fs1,Cox2=fs2),Surv(time,event)~1,data=dstest,
          plots="cal",metrics=NULL)
plotCalibration(xs)
plotCalibration(xs,cens.method="local",pseudo=1)
plotCalibration(xs,method="quantile")


# competing risks

## Not run: 
##D data(Melanoma)
##D f1 <- CSC(Hist(time,status)~age+sex+epicel+ulcer,data=Melanoma)
##D f2 <- CSC(Hist(time,status)~age+sex+logthick+epicel+ulcer,data=Melanoma)
##D x <- Score(list(model1=f1,model2=f2),Hist(time,status)~1,data=Melanoma,
##D            cause= 2,times=5*365.25,plots="cal")
##D plotCalibration(x)
## End(Not run)




cleanEx()
nameEx("plotEffects")
### * plotEffects

flush(stderr()); flush(stdout())

### Name: plotEffects
### Title: Plotting time-varying effects from a risk regression model.
### Aliases: plotEffects
### Keywords: survival

### ** Examples


library(survival)
library(prodlim)
data(Melanoma)

fit.tarr <- ARR(Hist(time,status)~strata(sex),
                data=Melanoma,
                cause=1)
plotEffects(fit.tarr)

fit.tarr <- ARR(Hist(time,status)~strata(sex)+strata(invasion),
                data=Melanoma,
                cause=1,
                times=seq(800,3000,20))
plotEffects(fit.tarr,formula=~sex)
plotEffects(fit.tarr,formula=~invasion)
plotEffects(fit.tarr,
            formula=~invasion,
            level="invasionlevel.1")

## legend arguments are transcluded:
plotEffects(fit.tarr,
            formula=~invasion,
            legend.bty="b",
            legend.cex=1)

## and other smart arguments too:
plotEffects(fit.tarr,
	    formula=~invasion,
	    legend.bty="b",
axis2.las=2,
	    legend.cex=1)





cleanEx()
nameEx("plotPredictRisk")
### * plotPredictRisk

flush(stderr()); flush(stdout())

### Name: plotPredictRisk
### Title: Plotting predicted risks curves.
### Aliases: plotPredictRisk
### Keywords: survival

### ** Examples

library(survival)
# generate survival data
# no effect
set.seed(8)
d <- sampleData(80,outcome="survival",formula = ~f(X6, 0) + f(X7, 0))
d[,table(event)]
f <- coxph(Surv(time,event)~X6+X7,data=d,x=1)
plotPredictRisk(f)

# large effect
set.seed(8)
d <- sampleData(80,outcome="survival",formula = ~f(X6, 0.1) + f(X7, -0.1))
d[,table(event)]
f <- coxph(Surv(time,event)~X6+X7,data=d,x=1)
plotPredictRisk(f)

# generate competing risk data
# small effect
set.seed(8)
d <- sampleData(40,formula = ~f(X6, 0.01) + f(X7, -0.01))
d[,table(event)]
f <- CSC(Hist(time,event)~X5+X6,data=d)
plotPredictRisk(f)

# large effect
set.seed(8)
d <- sampleData(40,formula = ~f(X6, 0.1) + f(X7, -0.1))
d[,table(event)]
f <- CSC(Hist(time,event)~X5+X6,data=d)
plotPredictRisk(f)



cleanEx()
nameEx("plotROC")
### * plotROC

flush(stderr()); flush(stdout())

### Name: plotROC
### Title: Plot ROC curves
### Aliases: plotROC

### ** Examples

## binary
set.seed(18)
library(randomForest)
library(prodlim)
bdl <- sampleData(40,outcome="binary")
bdt <- sampleData(58,outcome="binary")
bdl[,y:=factor(Y)]
bdt[,y:=factor(Y)]
fb1 <- glm(y~X1+X2+X3+X4+X5+X6+X7+X8+X9+X10,data=bdl,family="binomial")
fb2 <- randomForest(y~X1+X2+X3+X4+X5+X6+X7+X8+X9+X10,data=bdl)
xb <- Score(list("glm"=fb1,"rf"=fb2),y~1,data=bdt,
            plots="roc",metrics=c("auc","brier"))
plotROC(xb,brier.in.legend=1L)

# with cross-validation
## Not run: 
##D xb3 <- Score(list("glm"=fb1,"rf"=fb2),y~1,data=bdl,
##D             plots="roc",B=3,split.method="bootcv",
##D             metrics=c("auc"))
## End(Not run)
## survival
set.seed(18)
library(survival)
sdl <- sampleData(40,outcome="survival")
sdt <- sampleData(58,outcome="survival")
fs1 <- coxph(Surv(time,event)~X3+X5+X6+X7+X8+X10,data=sdl,x=TRUE)
fs2 <- coxph(Surv(time,event)~X1+X2+X9,data=sdl,x=TRUE)
xs <- Score(list(model1=fs1,model2=fs2),Hist(time,event)~1,data=sdt,
            times=5,plots="roc",metrics="auc")
plotROC(xs)
## competing risks
data(Melanoma)
f1 <- CSC(Hist(time,status)~age+sex+epicel+ulcer,data=Melanoma)
f2 <- CSC(Hist(time,status)~age+sex+logthick+epicel+ulcer,data=Melanoma)
x <- Score(list(model1=f1,model2=f2),Hist(time,status)~1,data=Melanoma,
            cause=1,times=5*365.25,plots="roc",metrics="auc")
plotROC(x)



cleanEx()
nameEx("plotRisk")
### * plotRisk

flush(stderr()); flush(stdout())

### Name: plotRisk
### Title: plot predicted risks
### Aliases: plotRisk

### ** Examples

library(prodlim)
## uncensored
learndat = sampleData(40,outcome="binary")
testdat = sampleData(40,outcome="binary")
lr1 = glm(Y~X1+X2+X7+X9,data=learndat,family="binomial")
lr2 = glm(Y~X3+X5+X6,data=learndat,family="binomial")
xb=Score(list("LR(X1+X2+X7+X9)"=lr1,"LR(X3+X5+X6)"=lr2),formula=Y~1,
         data=testdat,summary="risks",null.model=0L)
plotRisk(xb)
## survival
library(survival)
learndat = sampleData(40,outcome="survival")
testdat = sampleData(40,outcome="survival")
cox1 = coxph(Surv(time,event)~X1+X2+X7+X9,data=learndat,x=TRUE)
cox2 = coxph(Surv(time,event)~X3+X5+X6,data=learndat,x=TRUE)
xs=Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),formula=Surv(time,event)~1,
         data=testdat,summary="risks",null.model=0L,times=c(3,5,6))
plotRisk(xs,times=5)
## competing risk
## Not run: 
##D library(prodlim)
##D library(survival)
##D set.seed(8)
##D learndat = sampleData(80,outcome="competing.risk")
##D testdat = sampleData(140,outcome="competing.risk")
##D m1 = FGR(Hist(time,event)~X2+X7+X9,data=learndat,cause=1)
##D m2 = CSC(Hist(time,event)~X2+X7+X9,data=learndat,cause=1)
##D xcr=Score(list("FGR"=m1,"CSC"=m2),formula=Hist(time,event)~1,
##D          data=testdat,summary="risks",null.model=0L,times=c(3,5))
##D plotRisk(xcr,times=1)
## End(Not run)



cleanEx()
nameEx("predict.CauseSpecificCox")
### * predict.CauseSpecificCox

flush(stderr()); flush(stdout())

### Name: predict.CauseSpecificCox
### Title: Predicting Absolute Risk from Cause-Specific Cox Models
### Aliases: predict.CauseSpecificCox predictBig.CauseSpecificCox

### ** Examples

library(survival)
library(prodlim)
#### generate data ####
set.seed(5)
d <- sampleData(80,outcome="comp") ## training dataset
nd <- sampleData(4,outcome="comp") ## validation dataset
d$time <- round(d$time,1) ## create tied events
ttt <- sort(sample(x = unique(d$time), size = 10))

## estimate a CSC model based on the coxph function
CSC.fit <- CSC(Hist(time,event)~ X3+X8, data=d, method = "breslow")

## compute the absolute risk of cause 1, in the validation dataset
## at time 1:10
CSC.risk <-  predict(CSC.fit, newdata=nd, times=1:10, cause=1)
CSC.risk

## compute absolute risks with CI for cause 2
## (without displaying the value of the covariates)
predict(CSC.fit,newdata=nd,times=1:10,cause=2,se=TRUE,
        keep.newdata = FALSE)

## other example
library(survival)
CSC.fit.s <- CSC(list(Hist(time,event)~ strata(X1)+X2+X9,
 Hist(time,event)~ X2+strata(X4)+X8+X7),data=d, method = "breslow")
predict(CSC.fit.s,cause=1,times=ttt,se=1L) ## note: absRisk>1 due to small number of observations

## using the cph function instead of coxph
CSC.cph <- CSC(Hist(time,event)~ X1+X2,data=d, method = "breslow", fitter = "cph")#' 
predict(CSC.cph, newdata = d, cause = 2, times = ttt)

## landmark analysis
T0 <- 1
predCSC_afterT0 <- predict(CSC.fit, newdata = d, cause = 2, times = ttt[ttt>T0], landmark = T0)
predCSC_afterT0



cleanEx()
nameEx("predict.FGR")
### * predict.FGR

flush(stderr()); flush(stdout())

### Name: predict.FGR
### Title: Predict subject specific risks (cumulative incidence) based on
###   Fine-Gray regression model
### Aliases: predict.FGR

### ** Examples

library(prodlim)
library(survival)
set.seed(10)
d <- sampleData(101, outcome = "competing.risk")
tFun<-function(t) {t}
fgr<-FGR(Hist(time, event)~X1+strata(X2)+X6+cov2(X7, tf=tFun),
         data=d, cause=1)
predictRisk(fgr,times=5,newdata=d[1:10])



cleanEx()
nameEx("predict.riskRegression")
### * predict.riskRegression

flush(stderr()); flush(stdout())

### Name: predict.riskRegression
### Title: Predict individual risk.
### Aliases: predict.riskRegression
### Keywords: survival

### ** Examples


data(Melanoma)
library(prodlim)
library(survival)

fit.tarr <- ARR(Hist(time,status)~age+invasion+strata(sex),data=Melanoma,cause=1)
predict(fit.tarr,newdata=data.frame(age=48,
                     invasion=factor("level.1",
                         levels=levels(Melanoma$invasion)),
                     sex=factor("Female",levels=levels(Melanoma$sex))))
predict(fit.tarr,newdata=data.frame(age=48,
                     invasion=factor("level.1",
                         levels=levels(Melanoma$invasion)),
                     sex=factor("Male",levels=levels(Melanoma$sex))))
predict(fit.tarr,newdata=data.frame(age=c(48,58,68),
                     invasion=factor("level.1",
                         levels=levels(Melanoma$invasion)),
                     sex=factor("Male",levels=levels(Melanoma$sex))))
predict(fit.tarr,newdata=Melanoma[1:4,])




cleanEx()
nameEx("predictCox")
### * predictCox

flush(stderr()); flush(stdout())

### Name: predictCox
### Title: Fast computation of survival probabilities, hazards and
###   cumulative hazards from Cox regression models
### Aliases: predictCox

### ** Examples


library(survival)
library(data.table)

#### generate data ####
set.seed(10)
d <- sampleData(40,outcome="survival") ## training dataset
nd <- sampleData(4,outcome="survival") ## validation dataset
d$time <- round(d$time,1) ## create tied events
# table(duplicated(d$time))

#### stratified Cox model ####
fit <- coxph(Surv(time,event)~X1 + strata(X2) + X6,
             data=d, ties="breslow", x = TRUE, y = TRUE)

## compute the baseline cumulative hazard
fit.haz <- predictCox(fit)
cbind(survival::basehaz(fit), fit.haz$cumhazard)

## compute individual specific cumulative hazard and survival probabilities 
fit.pred <- predictCox(fit, newdata=nd, times=c(3,8), se = TRUE, band = TRUE)
fit.pred

####  other examples ####
# one strata variable
fitS <- coxph(Surv(time,event)~strata(X1)+X2,
              data=d, ties="breslow", x = TRUE, y = TRUE)

predictCox(fitS)
predictCox(fitS, newdata=nd, times = 1)

# two strata variables
set.seed(1)
d$U=sample(letters[1:5],replace=TRUE,size=NROW(d))
d$V=sample(letters[4:10],replace=TRUE,size=NROW(d))
nd$U=sample(letters[1:5],replace=TRUE,size=NROW(nd))
nd$V=sample(letters[4:10],replace=TRUE,size=NROW(nd))
fit2S <- coxph(Surv(time,event)~X1+strata(U)+strata(V)+X2,
              data=d, ties="breslow", x = TRUE, y = TRUE)

cbind(survival::basehaz(fit2S),predictCox(fit2S,type="cumhazard")$cumhazard)
predictCox(fit2S)
predictCox(fitS, newdata=nd, times = 3)

# left truncation
test2 <- list(start=c(1,2,5,2,1,7,3,4,8,8), 
              stop=c(2,3,6,7,8,9,9,9,14,17), 
              event=c(1,1,1,1,1,1,1,0,0,0), 
              x=c(1,0,0,1,0,1,1,1,0,0)) 
m.cph <- coxph(Surv(start, stop, event) ~ 1, test2, x = TRUE)
as.data.table(predictCox(m.cph))

basehaz(m.cph)



cleanEx()
nameEx("predictCoxPL")
### * predictCoxPL

flush(stderr()); flush(stdout())

### Name: predictCoxPL
### Title: Computation of survival probabilities from Cox regression models
###   using the product limit estimator.
### Aliases: predictCoxPL

### ** Examples

library(survival)

#### generate data ####
set.seed(10)
d <- sampleData(40,outcome="survival")
nd <- sampleData(4,outcome="survival")
d$time <- round(d$time,1)

#### Cox model ####
fit <- coxph(Surv(time,event)~ X1 + X2 + X6,
             data=d, ties="breslow", x = TRUE, y = TRUE)

## exponential approximation
predictCox(fit, newdata = d, times = 1:5)

## product limit
predictCoxPL(fit, newdata = d, times = 1:5)

#### stratified Cox model ####
fitS <- coxph(Surv(time,event)~ X1 + strata(X2) + X6,
             data=d, ties="breslow", x = TRUE, y = TRUE)

## exponential approximation
predictCox(fitS, newdata = d, times = 1:5)

## product limit
predictCoxPL(fitS, newdata = d, times = 1:5)

#### fully stratified Cox model ####
fitS <- coxph(Surv(time,event)~ 1,
             data=d, ties="breslow", x = TRUE, y = TRUE)

## product limit
GS <- survfit(Surv(time,event)~1, data = d)
range(predictCoxPL(fitS)$survival - GS$surv)

fitS <- coxph(Surv(time,event)~ strata(X2),
             data=d, ties="breslow", x = TRUE, y = TRUE)

## product limit
GS <- survfit(Surv(time,event)~X2, data = d)
range(predictCoxPL(fitS)$survival - GS$surv)




cleanEx()
nameEx("predictRisk")
### * predictRisk

flush(stderr()); flush(stdout())

### Name: predictRisk
### Title: Extrating predicting risks from regression models
### Aliases: predictRisk predictRisk.CauseSpecificCox
###   predictRisk.riskRegression predictRisk.FGR predictRisk.prodlim
###   predictRisk.rfsrc predictRisk.aalen predictRisk.ARR
###   predictRisk.cox.aalen predictRisk.coxph predictRisk.cph
###   predictRisk.default predictRisk.matrix predictRisk.pecCtree
###   predictRisk.pecCforest predictRisk.psm predictRisk.selectCox
###   predictRisk.survfit predictRisk.randomForest predictRisk.lrm
###   predictRisk.glm predictRisk.rpart predictRisk.gbm
###   predictRisk.flexsurvreg predictRisk.double predictRisk.integer
###   predictRisk.factor predictRisk.numeric predictRisk.formula
###   predictRisk.BinaryTree predictRisk.coxphTD predictRisk.CSCTD
###   predictRisk.coxph.penal predictRisk.ranger predictRisk.penfitS3
###   predictRisk.SuperPredictor
### Keywords: survival

### ** Examples

## binary outcome
library(rms)
set.seed(7)
d <- sampleData(80,outcome="binary")
nd <- sampleData(80,outcome="binary")
fit <- lrm(Y~X1+X8,data=d)
predictRisk(fit,newdata=nd)
## Not run: 
##D library(SuperLearner)
##D set.seed(1)
##D sl = SuperLearner(Y = d$Y, X = d[,-1], family = binomial(),
##D       SL.library = c("SL.mean", "SL.glmnet", "SL.randomForest"))
## End(Not run)

## survival outcome
# generate survival data
library(prodlim)
set.seed(100)
d <- sampleData(100,outcome="survival")
d[,X1:=as.numeric(as.character(X1))]
d[,X2:=as.numeric(as.character(X2))]
# then fit a Cox model
library(rms)
cphmodel <- cph(Surv(time,event)~X1+X2,data=d,surv=TRUE,x=TRUE,y=TRUE)
# or via survival
library(survival)
coxphmodel <- coxph(Surv(time,event)~X1+X2,data=d,x=TRUE,y=TRUE)

# Extract predicted survival probabilities 
# at selected time-points:
ttt <- quantile(d$time)
# for selected predictor values:
ndat <- data.frame(X1=c(0.25,0.25,-0.05,0.05),X2=c(0,1,0,1))
# as follows
predictRisk(cphmodel,newdata=ndat,times=ttt)
predictRisk(coxphmodel,newdata=ndat,times=ttt)

# stratified cox model
sfit <- coxph(Surv(time,event)~strata(X1)+X2,data=d,x=TRUE,y=TRUE)
predictRisk(sfit,newdata=d[1:3,],times=c(1,3,5,10))

## simulate learning and validation data
learndat <- sampleData(100,outcome="survival")
valdat <- sampleData(100,outcome="survival")
## use the learning data to fit a Cox model
library(survival)
fitCox <- coxph(Surv(time,event)~X1+X2,data=learndat,x=TRUE,y=TRUE)
## suppose we want to predict the survival probabilities for all subjects
## in the validation data at the following time points:
## 0, 12, 24, 36, 48, 60
psurv <- predictRisk(fitCox,newdata=valdat,times=seq(0,60,12))
## This is a matrix with event probabilities (1-survival)
## one column for each of the 5 time points
## one row for each validation set individual

# Do the same for a randomSurvivalForest model
# library(randomForestSRC)
# rsfmodel <- rfsrc(Surv(time,event)~X1+X2,data=learndat)
# prsfsurv=predictRisk(rsfmodel,newdata=valdat,times=seq(0,60,12))
# plot(psurv,prsfsurv)

## Cox with ridge option
f1 <- coxph(Surv(time,event)~X1+X2,data=learndat,x=TRUE,y=TRUE)
f2 <- coxph(Surv(time,event)~ridge(X1)+ridge(X2),data=learndat,x=TRUE,y=TRUE)
## Not run: 
##D plot(predictRisk(f1,newdata=valdat,times=10),
##D      riskRegression:::predictRisk.coxph(f2,newdata=valdat,times=10),
##D      xlim=c(0,1),
##D      ylim=c(0,1),
##D      xlab="Unpenalized predicted survival chance at 10",
##D      ylab="Ridge predicted survival chance at 10")
## End(Not run)

## competing risks

library(survival)
library(riskRegression)
library(prodlim)
train <- prodlim::SimCompRisk(100)
test <- prodlim::SimCompRisk(10)
cox.fit  <- CSC(Hist(time,cause)~X1+X2,data=train)
predictRisk(cox.fit,newdata=test,times=seq(1:10),cause=1)

## with strata
cox.fit2  <- CSC(list(Hist(time,cause)~strata(X1)+X2,Hist(time,cause)~X1+X2),data=train)
predictRisk(cox.fit2,newdata=test,times=seq(1:10),cause=1)




cleanEx()
nameEx("riskLevelPlot")
### * riskLevelPlot

flush(stderr()); flush(stdout())

### Name: riskLevelPlot
### Title: Level plots for risk prediction models
### Aliases: riskLevelPlot

### ** Examples


# ---------- logistic regression --------------------
expit <- function(x){exp(x)/(1+exp(x))}
partyData <- function(N){
  Age <- runif(N,.5,15)
  Parasites <- rnorm(N,mean=3.5-0.03*Age)
  Fever <- factor(rbinom(N,1,expit(-3.5-.3*Age+.55*Parasites+0.15*Age*Parasites)))
  data.frame(Fever,Age,Parasites)
}
d <- partyData(100)
f <- glm(Fever~Age+Parasites,data=d,family="binomial")
riskLevelPlot(f,Fever~Age+Parasites,d)
library(randomForest)
rf <- randomForest(Fever~Age+Parasites,data=d)
riskLevelPlot(f,Fever~Age+Parasites,d)
riskLevelPlot(rf,Fever~Age+Parasites,d)

# ---------- survival analysis --------------------

# --simulate an artificial data frame
# with survival response and three predictors

library(survival)
library(prodlim)
set.seed(140515)
sdat <- sampleData(43,outcome="survival")
# -- fit a Cox regression model 
survForm = Surv(time,event) ~ X8 + X9
cox <- coxph(survForm, data = sdat,x=TRUE)

# --choose a time horizon for the predictions and plot the risks
timeHorizon <- floor(median(sdat$time))
riskLevelPlot(cox, survForm, data = sdat, horizon = timeHorizon)

# ---------- competing risks --------------------

# -- simulate an artificial data frame
# with competing cause response and three predictors
library(cmprsk)
library(riskRegression)
set.seed(140515)
crdat <- sampleData(49)

# -- fit a cause-specific Cox regression model
crForm <- Hist(time,event)~X8+X9
csCox  <- CSC(crForm, data=crdat)

# -- choose a time horizon and plot the risk for a given cause
timeHorizon <- floor(median(crdat$time))
riskLevelPlot(csCox, crForm, data = crdat, horizon = timeHorizon, cause = 1)




cleanEx()
nameEx("riskRegression")
### * riskRegression

flush(stderr()); flush(stdout())

### Name: riskRegression
### Title: Risk Regression Fits a regression model for the risk of an event
###   - allowing for competing risks.
### Aliases: riskRegression ARR LRR
### Keywords: survival

### ** Examples


library(prodlim)
data(Melanoma,package="riskRegression")
## tumor thickness on the log-scale
Melanoma$logthick <- log(Melanoma$thick)

# Single binary factor

## absolute risk regression
library(survival)
library(prodlim)
fit.arr <- ARR(Hist(time,status)~sex,data=Melanoma,cause=1)
print(fit.arr)
# show predicted cumulative incidences
plot(fit.arr,col=3:4,newdata=data.frame(sex=c("Female","Male")))

## compare with non-parametric Aalen-Johansen estimate
library(prodlim)
fit.aj <- prodlim(Hist(time,status)~sex,data=Melanoma)
plot(fit.aj,conf.int=FALSE)
plot(fit.arr,add=TRUE,col=3:4,newdata=data.frame(sex=c("Female","Male")))

## with time-dependent effect
fit.tarr <- ARR(Hist(time,status)~strata(sex),data=Melanoma,cause=1)
plot(fit.tarr,newdata=data.frame(sex=c("Female","Male")))

## logistic risk regression
fit.lrr <- LRR(Hist(time,status)~sex,data=Melanoma,cause=1)
summary(fit.lrr)


# Single continuous factor

## tumor thickness on the log-scale
Melanoma$logthick <- log(Melanoma$thick)

## absolute risk regression 
fit2.arr <- ARR(Hist(time,status)~logthick,data=Melanoma,cause=1)
print(fit2.arr)
# show predicted cumulative incidences
plot(fit2.arr,col=1:5,newdata=data.frame(logthick=quantile(Melanoma$logthick)))

## comparison with nearest neighbor non-parametric Aalen-Johansen estimate
library(prodlim)
fit2.aj <- prodlim(Hist(time,status)~logthick,data=Melanoma)
plot(fit2.aj,conf.int=FALSE,newdata=data.frame(logthick=quantile(Melanoma$logthick)))
plot(fit2.arr,add=TRUE,col=1:5,lty=3,newdata=data.frame(logthick=quantile(Melanoma$logthick)))

## logistic risk regression
fit2.lrr <- LRR(Hist(time,status)~logthick,data=Melanoma,cause=1)
summary(fit2.lrr)

## change model for censoring weights
library(rms)
fit2a.lrr <- LRR(Hist(time,status)~logthick,
                 data=Melanoma,
                 cause=1,
                 cens.model="cox",
                 cens.formula=~sex+epicel+ulcer+age+logthick)
summary(fit2a.lrr)

##  compare prediction performance
Score(list(ARR=fit2.arr,AJ=fit2.aj,LRR=fit2.lrr),formula=Hist(time,status)~1,data=Melanoma)


# multiple regression
library(riskRegression)
library(prodlim)
# absolute risk model
multi.arr <- ARR(Hist(time,status)~logthick+sex+age+ulcer,data=Melanoma,cause=1)

# stratified model allowing different baseline risk for the two gender
multi.arr <- ARR(Hist(time,status)~thick+strata(sex)+age+ulcer,data=Melanoma,cause=1)

# stratify by a continuous variable: strata(age)
multi.arr <- ARR(Hist(time,status)~tp(thick,power=0)+strata(age)+sex+ulcer,
                 data=Melanoma,
                 cause=1)

fit.arr2a <- ARR(Hist(time,status)~tp(thick,power=1),data=Melanoma,cause=1)
summary(fit.arr2a)
fit.arr2b <- ARR(Hist(time,status)~timevar(thick),data=Melanoma,cause=1)
summary(fit.arr2b)

## logistic risk model
fit.lrr <- LRR(Hist(time,status)~thick,data=Melanoma,cause=1)
summary(fit.lrr)





## nearest neighbor non-parametric Aalen-Johansen estimate
library(prodlim)
fit.aj <- prodlim(Hist(time,status)~thick,data=Melanoma)
plot(fit.aj,conf.int=FALSE)

# prediction performance
x <- Score(list(fit.arr2a,fit.arr2b,fit.lrr),
             data=Melanoma,
             formula=Hist(time,status)~1,
             cause=1,
             split.method="none")





cleanEx()
nameEx("riskRegression.options")
### * riskRegression.options

flush(stderr()); flush(stdout())

### Name: riskRegression.options
### Title: Global options for 'riskRegression'
### Aliases: riskRegression.options

### ** Examples

options <- riskRegression.options()

## add new method.predictRiskIID
riskRegression.options(method.predictRiskIID = c(options$method.predictRiskIID,"xx"))

riskRegression.options()



cleanEx()
nameEx("rowCenter_cpp")
### * rowCenter_cpp

flush(stderr()); flush(stdout())

### Name: rowCenter_cpp
### Title: Apply - by row
### Aliases: rowCenter_cpp

### ** Examples

x <- matrix(1,6,5)
sweep(x, MARGIN = 2, FUN = "-", STATS = 1:5)
rowCenter_cpp(x, 1:5 )

rowCenter_cpp(x, colMeans(x) )



cleanEx()
nameEx("rowCumProd")
### * rowCumProd

flush(stderr()); flush(stdout())

### Name: rowCumProd
### Title: Apply cumprod in each row
### Aliases: rowCumProd

### ** Examples

x <- matrix(1:8,ncol=2)
rowCumProd(x)



cleanEx()
nameEx("rowCumSum")
### * rowCumSum

flush(stderr()); flush(stdout())

### Name: rowCumSum
### Title: Apply cumsum in each row
### Aliases: rowCumSum

### ** Examples

x <- matrix(1:8,ncol=2)
rowCumSum(x)



cleanEx()
nameEx("rowMultiply_cpp")
### * rowMultiply_cpp

flush(stderr()); flush(stdout())

### Name: rowMultiply_cpp
### Title: Apply * by row
### Aliases: rowMultiply_cpp

### ** Examples

x <- matrix(1,6,5)
sweep(x, MARGIN = 2, FUN = "*", STATS = 1:5)
rowMultiply_cpp(x, 1:5 )

rowMultiply_cpp(x, 1/colMeans(x) )




cleanEx()
nameEx("rowPaste")
### * rowPaste

flush(stderr()); flush(stdout())

### Name: rowPaste
### Title: Collapse Rows of Characters.
### Aliases: rowPaste

### ** Examples

## Not run: 
##D M <- matrix(letters,nrow = 26, ncol = 2)
##D rowPaste(M)
## End(Not run)



cleanEx()
nameEx("rowScale_cpp")
### * rowScale_cpp

flush(stderr()); flush(stdout())

### Name: rowScale_cpp
### Title: Apply / by row
### Aliases: rowScale_cpp

### ** Examples

x <- matrix(1,6,5)
sweep(x, MARGIN = 2, FUN = "/", STATS = 1:5)
rowScale_cpp(x, 1:5 )

rowScale_cpp(x, colMeans(x) )



cleanEx()
nameEx("rowSumsCrossprod")
### * rowSumsCrossprod

flush(stderr()); flush(stdout())

### Name: rowSumsCrossprod
### Title: Apply crossprod and rowSums
### Aliases: rowSumsCrossprod

### ** Examples

x <- matrix(1:10,nrow=5)
y <- matrix(1:20,ncol=4)
rowSumsCrossprod(x,y,0)

x <- matrix(1:10,nrow=5)
y <- matrix(1:20,ncol=5)
rowSumsCrossprod(x,y,1)



cleanEx()
nameEx("sampleData")
### * sampleData

flush(stderr()); flush(stdout())

### Name: sampleData
### Title: Simulate data with binary or time-to-event outcome
### Aliases: sampleData sampleDataTD

### ** Examples

sampleData(10,outcome="binary")
sampleData(10,outcome="survival")
sampleData(10,outcome="competing.risks")



cleanEx()
nameEx("selectCox")
### * selectCox

flush(stderr()); flush(stdout())

### Name: selectCox
### Title: Backward variable selection in the Cox regression model
### Aliases: selectCox
### Keywords: survival

### ** Examples


library(pec)
library(prodlim)
data(GBSG2)
library(survival)
f <- selectCox(Surv(time,cens)~horTh+age+menostat+tsize+tgrade+pnodes+progrec+estrec ,
	       data=GBSG2)




cleanEx()
nameEx("simActiveSurveillance")
### * simActiveSurveillance

flush(stderr()); flush(stdout())

### Name: simActiveSurveillance
### Title: Simulate data of a hypothetical active surveillance prostate
###   cancer study
### Aliases: simActiveSurveillance

### ** Examples

set.seed(71)
simActiveSurveillance(3)



cleanEx()
nameEx("simMelanoma")
### * simMelanoma

flush(stderr()); flush(stdout())

### Name: simMelanoma
### Title: Simulate data alike the Melanoma data
### Aliases: simMelanoma

### ** Examples

set.seed(71)
simMelanoma(3)



cleanEx()
nameEx("sliceMultiply_cpp")
### * sliceMultiply_cpp

flush(stderr()); flush(stdout())

### Name: sliceMultiply_cpp
### Title: Apply * by slice
### Aliases: sliceMultiply_cpp sliceMultiplyPointer_cpp

### ** Examples

x <- array(1, dim = c(2,6,5))
M <- matrix(1:12,2,6)
sweep(x, MARGIN = 1:2, FUN = "*", STATS = M)
sliceMultiply_cpp(x, M) 




cleanEx()
nameEx("sliceScale_cpp")
### * sliceScale_cpp

flush(stderr()); flush(stdout())

### Name: sliceScale_cpp
### Title: Apply / by slice
### Aliases: sliceScale_cpp sliceScalePointer_cpp

### ** Examples

x <- array(1, dim = c(2,6,5))
M <- matrix(1:12,2,6)
sweep(x, MARGIN = 1:2, FUN = "/", STATS = M)
sliceScale_cpp(x, M) 




cleanEx()
nameEx("subjectWeights")
### * subjectWeights

flush(stderr()); flush(stdout())

### Name: subjectWeights
### Title: Estimation of censoring probabilities at subject specific times
### Aliases: subjectWeights subjectWeights.none subjectWeights.km
###   subjectWeights.marginal subjectWeights.nonpar subjectWeights.cox
###   subjectWeights.aalen
### Keywords: survival

### ** Examples


library(prodlim)
library(survival)
dat=SimSurv(300)

dat <- dat[order(dat$time,-dat$status),]

# using the marginal Kaplan-Meier for the censoring times

WKM=subjectWeights(Hist(time,status)~X2,data=dat,method="marginal")
plot(WKM$fit)
WKM$fit
WKM$weights

# using the Cox model for the censoring times given X2

WCox=subjectWeights(Surv(time,status)~X2,data=dat,method="cox")
WCox
plot(WCox$weights,WKM$weights)

# using the stratified Kaplan-Meier for the censoring times given X2

WKM2 <- subjectWeights(Surv(time,status)~X2,data=dat,method="nonpar")
plot(WKM2$fit,add=FALSE)





cleanEx()
nameEx("subsetIndex")
### * subsetIndex

flush(stderr()); flush(stdout())

### Name: subsetIndex
### Title: Extract Specific Elements From An Object
### Aliases: subsetIndex subsetIndex.default subsetIndex.matrix

### ** Examples

M <- matrix(rnorm(50),5,10)
subsetIndex(M, index = c(0,0,1), default = 0)
subsetIndex(M, index = c(0,2,3,NA), default = 0)
subsetIndex(M, index = c(0,NA,2,3,NA), default = 0)

C <- 1:10
subsetIndex(C, index = c(0,0,1,5,NA), default = 0)



### * <FOOTER>
###
cleanEx()
options(digits = 7L)
base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
grDevices::dev.off()
###
### Local variables: ***
### mode: outline-minor ***
### outline-regexp: "\\(> \\)?### [*]+" ***
### End: ***
quit('no')