####R code for generating data that follow a frailty model
####April, 2001


B <- c(2, -1) # define the slope for the two covariats
zij <- matrix(c( 1:50/10, runif(50)), nrow=2,ncol=50, byrow=TRUE)
#### define the two covariates
wi <- rgamma(5,1,1)   # generates the frailty 
lambdaij <- rep(0, 50)
for(k in 1:5) {
  for(j in (1+(k-1)*10):(10+(k-1)*10)) {lambdaij[j]<-wi[k]*exp(zij[,j]%*%B) }}
#### this generates the exponential parameters that includes 
#### the frailty and covariate
yi<-rexp(50, rate=lambdaij)   
#### generate the response variable using the lambdaij 
gyi<-100*sqrt(yi)     # use a transformation to give it a baseline hazard
clusters<-sort(rep(c(1,2,3,4,5),10))  # the frailty covariate
cov1<-zij[1,]
cov2<-zij[2,]
simudata<-cbind(clusters, gyi, cov1, cov2)
#### Now fit the model 
library(survival5)
coxph(Surv(gyi, rep(1, 50))~
               cov1+cov2+frailty(clusters, dist="gamma),data=simudata)
     # or this (same model)
coxph(Surv(gyi, rep(1, 50))~cov1+cov2+frailty(clusters),data=simudata)



####We can add censoring to the gyi if needed.