# Simulate data from a proportional hazards population model
require(survival)
n <- 1000
set.seed(731)
age <- 50 + 12*rnorm(n)
label(age) <- "Age"
cens <- 15*runif(n)
h <- .02*exp(.04*(age-50))
dt <- -log(runif(n))/h
label(dt) <- 'Follow-up Time'
e <- ifelse(dt <= cens,1,0)
dt <- pmin(dt, cens)
units(dt) <- "Year"
S <- Surv(dt,e)
f <- psm(S ~ lsp(age,c(40,70)))
survest(f, data.frame(age=seq(20,80,by=5)), times=2)
#Get predicted survival curve for 40 year old
survest(f, data.frame(age=40))
#Get hazard function for 40 year old
survest(f, data.frame(age=40), what="hazard")$surv #still called surv
Run the code above in your browser using DataLab