if (FALSE) {
# remove NAs
lung <- lung[!is.na(lung$inst),]
# transform inst into factor variable
lung$inst <- as.factor(lung$inst)
# just for illustration, create factor with only three ph.ecog classes
lung$ph.ecog[is.na(lung$ph.ecog)] <- 2
lung$ph.ecog[lung$ph.ecog==3] <- 2
lung$ph.ecog <- as.factor(lung$ph.ecog)
fix.form <- as.formula("Surv(time, status) ~ 1 + age + ph.ecog + sex")
ridge.obj <- coxridge(fix=fix.form, data=lung, xi.ridge=10,
control=list(print.iter=TRUE, exact = 1))
coef(ridge.obj)
# now add random institutional effect
ridge.obj2 <- coxridge(fix=fix.form, rnd = list(inst=~1),
data=lung, xi.ridge=10,control=list(print.iter=TRUE, exact = 1))
coef(ridge.obj2)
# print frailty Std.Dev.
print(ridge.obj2$Q)
# print frailties
print(ridge.obj2$ranef)
# now fit a time-varying effect for age
fix.form <- as.formula("Surv(time, status) ~ 1 + ph.ecog + sex")
vary.coef <- as.formula("~ age")
ridge.obj3 <- ridgelasso(fix=fix.form,vary.coef=vary.coef,
data=lung, xi.ridge=10,control=list(print.iter=TRUE))
summary(ridge.obj3)
# show fit
plot(ridge.obj3)
# predict survival curve of new subject, institution 1 and up to time 300
pred.obj <- predict(ridge.obj2, newdata=data.frame(inst=1, time=NA, status=NA, age=26,
ph.ecog=2,sex=1), time.grid=seq(0,300,by=1))
# plot predicted hazard function
plot(pred.obj$time.grid,pred.obj$haz,type="l",xlab="time",ylab="hazard")
# plot predicted survival function
plot(pred.obj$time.grid,pred.obj$survival,type="l",xlab="time",ylab="survival")
## specify a larger new data set
new.data <- data.frame(inst=c(1,1,6), time=c(20,40,200),
status=c(NA,NA,NA), age=c(26,26,54), ph.ecog=c(0,0,2),sex=c(1,1,1))
## as here no frailties have been specified, id.var needs to be given!
pred.obj2 <- predict(lasso.obj3, newdata=new.data,id.var = "inst")
# plot predicted hazard functions (for the available time intervals)
plot(pred.obj2$time.grid[!is.na(pred.obj2$haz[,1])],
pred.obj2$haz[,1][!is.na(pred.obj2$haz[,1])],
type="l",xlab="time",ylab="hazard",xlim=c(0,200),
ylim=c(0,max(pred.obj2$haz,na.rm=T)))
lines(pred.obj2$time.grid[!is.na(pred.obj2$haz[,3])],
pred.obj2$haz[,3][!is.na(pred.obj2$haz[,3])],
col="red",lty=2,)
# plot predicted survival functions (for the available time intervals)
plot(pred.obj2$time.grid[!is.na(pred.obj2$survival[,1])],
pred.obj2$survival[,1][!is.na(pred.obj2$survival[,1])],
type="l",xlab="time",ylab="hazard",xlim=c(0,200),
ylim=c(0,max(pred.obj2$survival,na.rm=T)))
lines(pred.obj2$time.grid[!is.na(pred.obj2$survival[,3])],
pred.obj2$survival[,3][!is.na(pred.obj2$survival[,3])],
col="red",lty=2,)
# see also demo("coxlasso-lung")
}
Run the code above in your browser using DataLab