Simulation of the Frailty of the Cox PH model
Mohammad Ehsanul Karim
I am trying to do simulation of survival data to enable it to run
under frailty option. Below is the function a that I am using. My
questions are:
1. How do I modify it to get bigger (hopefully significant) value of
Variance of random effect?
2. What changes do I have to make in the function to run it under
correlated frailty model? (may be in kinship package)
3. Is there any program to run frailty by Inverse gaussian or stable
family in S-plus or R?
Thank you for your time.
Thanks in advance.
Mohammad Ehsanul Karim
wildscop at yahoo dot com
Institute of Statistical Research and Training
University of Dhaka
# ***************************************
"sim.data"<- function(g,m){
set.seed(123)
group <- rep(1:m,rep(g,m))
Frailty <- rep(rgamma(m,100,1),rep(g,m))
covariate <- rbinom(g*m,1,.05)
stimes <- rweibull(g*m,1.1,1/(5*Frailty*exp((covariate)/.5)))
cens <- 5 + 5*runif(25)
times <- pmin(stimes, cens)
censored <- as.numeric(cens > times)
data.mat <- cbind(group,covariate,times,censored,Frailty)
data.mat <- data.mat[rev(order(times)),1:length(data.mat[1,])]
data.fr <- data.frame(data.mat)
return(data.fr)
}
# ***************************************
# Example of 50 group each with 100 members
sim.fr<-sim.data(50,100)
# library(survival) # to run in R
fit.c <- coxph(Surv(times,censored) ~ covariate,data= sim.fr)
# fit.c gives the Usual cox proportional hazards model
fit.gm.em <- coxph(Surv(times,censored) ~ covariate + frailty(group,
dist='gamma', method='em'), data= sim.fr)
# fit.gm.em gives the gamma frailty model by EM algorithm
fit.c # result of Cox PH model
fit.gm.em # result of gamma frailty model