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

Date: Sun 08 Apr 2007 - 07:47:07 GMT

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?

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

R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun Apr 08 17:50:56 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sun 08 Apr 2007 - 08:31:03 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.