## 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)
# GLMnet example
fit <- GLMnet(Y~X1+X8,data=d) ## Uses CV as default
predictRisk(fit,newdata=nd)
if (FALSE) {
library(SuperLearner)
set.seed(1)
sl = SuperLearner(Y = d$Y, X = d[,-1], family = binomial(),
SL.library = c("SL.mean", "SL.glmnet", "SL.randomForest"))
}
## 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
set.seed(10)
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
# Use GLMnet to predict survvival probabilities
fitGLMnet <- GLMnet(Surv(time,event)~X1+X8,data=learndat) ## Use CV as standard.
psurv <- predictRisk(fitGLMnet,newdata=valdat,times=seq(0,60,12))
# Use hal9001 as an example
if (FALSE) {
fitHAL <- Hal9001(Surv(time,event)~X1+X8,data=learndat) ## Use CV as standard.
psurv <- predictRisk(fitHAL,newdata=valdat,times=seq(0,60,12))
}
if (require("randomForestSRC",quietly=TRUE)){
# 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)
if (FALSE) {
plot(predictRisk(f1,newdata=valdat,times=10),
riskRegression:::predictRisk.coxph(f2,newdata=valdat,times=10),
xlim=c(0,1),
ylim=c(0,1),
xlab="Unpenalized predicted survival chance at 10",
ylab="Ridge predicted survival chance at 10")
}
# aalen model
library(timereg)
data(sTRACE)
out <- aalen(Surv(time, status==9) ~ sex + diabetes + chf + vf,
data=sTRACE, max.time=7, n.sim=0, resample.iid=1)
print(methods(predictRisk))
predictRisk(object=out, newdata=sTRACE[1:5,], times=c(1, 2, 3))
# cox.aalen model
library(timereg)
data(sTRACE)
out <- cox.aalen(Surv(time,status==9) ~ prop(age) + prop(sex) +
prop(diabetes) + chf + vf,
data=sTRACE, max.time=7, n.sim=0, resample.iid=1)
predictRisk(object=out, newdata=sTRACE[1:5,], times=c(1, 2, 3))
## 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)
library(timereg)
data(bmt)
add <- comp.risk(Event(time, cause) ~ platelet + age + tcell, data=bmt, cause=1)
ndata <- data.frame(platelet=c(1, 0, 0), age=c(0, 1, 0), tcell=c(0, 0, 1))
predictRisk(object=add, newdata=ndata, times=c(1, 2, 3))
Run the code above in your browser using DataLab