[R] Simulation of the Frailty of the Cox PH model

[Mohammad Ehsanul Karim]

Dear R-list users,

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 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)
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