# baseline hazard: Weibull
# N = sample size
# lambda = scale parameter in h0()
# rho = shape parameter in h0()
# beta = fixed effect parameter
# rateC = rate parameter of the exponential distribution of C
N=100
simulWeib <- function(N, lambda, rho, beta, rateC)
{
x=cbind(hba1c=rnorm(N,2,.5)-2,age=rnorm(N,40,5)-40,duration=rnorm(N,10,2)-10)
# Weibull latent event times
v <- runif(n=N)
Tlat <- (- log(v) / (lambda * exp(x %*% beta)))^(1 / rho)
# censoring times
C <- rexp(n=N, rate=rateC)
# follow-up times and event indicators
time <- pmin(Tlat, C)
status <- as.numeric(Tlat <= C)
# data set
data.frame(id=1:N,
           time=time,
           status=status,
           x=x)
}
library(survival)
set.seed(7003)
dat=simulWeib(N=100, lambda=0.01, rho=1, beta=c(1.1,1.2,1.3), rateC=0.001)
fit <- coxph(Surv(time, status) ~ x.hba1c + x.age + x.duration , data=dat)
summary(fit)

set.seed(7003)
betaHat1 <- rep(NA,1000)
betaHat2 <- rep(NA,1000)
betaHat3 <- rep(NA,1000)
for(k in 1:1000)
{
dat <- simulWeib(N=100, lambda=0.01, rho=1, beta=c(1.1,1.2,1.3), rateC=0.001)
fit <- coxph(Surv(time, status) ~ x.hba1c + x.age + x.duration , data=dat)
betaHat[k] <- c(fit$coef[1])
betaHat2[k] <- c(fit$coef[2])
betaHat3[k] <- c(fit$coef[3])
}

mean(betaHat[1:1000])
mean(betaHat2[1:1000])
mean(betaHat3[1:1000])