############################################################################################# ##Investigating hospital heterogeneity with a multi-state frailty model: ##application to nosocomial pneumonia disease in intensive care unit" ## by B. Liquet, J.F. Timsit and V. Rondeau ############################################################################################# # We propose two main functions : Simul.semi.markov and Analysis.multi.state # Simul.semi.markov: Create a simple simulated data set # Analysis.multi.state: perform a multi-state model based on the package "frailtypack" require(frailtypack) ############################################################################################# ## Create a simple simulated data set function: Simul.semi.markov ## With a Weibull distribution for each transition with two explanatory variables ## v1 and v2 follow Bernouilli distribution with probability 0.5 ############################################################################################# ############################################################################################# ### Example: one run of the simulation corresponding to table 4: ############################################################################################# G <- 30 # number of groups ni <- 150 # number of patients per group par01<-c(1.3,15);par02<-c(1.3,35);par03<-c(1.25,15);par12<-c(1.3,45);par13<-c(1.25,41) ### Weibull parameters for each transition function coef01<-c(0.8,1);coef02<-c(.6,1.2);coef03<-c(1.3,.3);coef12<-c(.7,1.1);coef13<-c(0.6,1.2) ### coefhk: true effects of the explanatory variable v1 and v2 p <- 0.5 # parameter of the Bernouilli for v1 and v2 gamma <- 0.15 # frailty variance paragamma <- c(1/gamma,gamma) #parameter for the gamma distribution for the frailties terms loi01<-list(par=par01,coef=coef01) loi02<-list(par=par02,coef=coef02) loi03<-list(par=par03,coef=coef03) loi12<-list(par=par12,coef=coef12) loi13<-list(par=par13,coef=coef13) ############################################################################################### ### ========================================== function Simul.semi.markov ==================### ### Simul.semi.markov: function which generate data from a semi-markov multi-state model Simul.semi.markov <-function(G=30,ni=150,paragamma=paragamma,p=0.5,loi01,loi02,loi03,loi12,loi13){ vecnj <- rep(ni,G) n <- sum(vecnj) # sample size of the generated data pass01<-rep(0,n) # pass01 indicates if subject transit from 0 to 1 at time "time0" pass02<-rep(0,n) # pass01 indicates if subject transit from 0 to 2 at time "time0" pass03<-rep(0,n) # pass01 indicates if subject transit from 0 to 3 at time "time0" pass12<-rep(0,n) # pass12 indicates if subject transit from 1 to 2 at time "time1" pass13<-rep(0,n) # pass13 indicates if subject transit from 1 to 3 at time "time1" v1<-rbinom(n,1,p) # generates the values of the variable v1 v2<-rbinom(n,1,p) # generates the values of the variable v2 obs<-c(1:n);time0<-rep(0,n);time1<-rep(0,n);centre <- rep(0,n) ## database: initialization at 0 database<-data.frame(obs,TIME0=time0,PASS01=pass01,v1=rep(0,n),v2=rep(0,n),PASS02=pass02,PASS03=pass03,TIME1=time1,PASS12=pass12,PASS13=pass13,IDCENTRE=centre) ##time0 represent the time of the first jump from 0 to k ##time1 represent the time spend in the state 1 before jumping from 1 to k centrej <- 0 i <- 0 for(nj in vecnj){ u01 <- rgamma(1,shape= paragamma[1],scale=paragamma[2]) u02 <- rgamma(1,shape= paragamma[1],scale=paragamma[2]) u03 <- rgamma(1,shape= paragamma[1],scale=paragamma[2]) u12 <- rgamma(1,shape= paragamma[1],scale=paragamma[2]) u13 <- rgamma(1,shape= paragamma[1],scale=paragamma[2]) centrej <- centrej+1 for (k in 1:nj){ i <- i+1 t01<- gene.weibull(loi01\$par,u01*exp(loi01\$coef[1]*v1[i]+loi01\$coef[2]*v2[i])) t02<- gene.weibull(loi02\$par,u02*exp(loi02\$coef[1]*v1[i]+loi02\$coef[2]*v2[i])) t03<- gene.weibull(loi03\$par,u03*exp(loi03\$coef[1]*v1[i]+loi03\$coef[2]*v2[i])) tk<-c(t01,t02,t03) ## 3 simulated times for the 3 transitions 01 02 03 if(which.min(tk)==1){ ## transition 01, t12<- gene.weibull(loi12\$par,u12*exp(loi12\$coef[1]*v1[i]+loi12\$coef[2]*v2[i])) t13<- gene.weibull(loi13\$par,u13*exp(loi13\$coef[1]*v1[i]+loi13\$coef[2]*v2[i])) rtt<-c(t12,t13) ## 2 simulated times for the transitions 12 13 if(which.min(rtt)==1){ ## then transition 12. bh<-c(i,tk[1],1,v1[i],v2[i],0,0,t12,1,0,centrej) }else{ ## then transition 13. bh<-c(i,tk[1],1,v1[i],v2[i],0,0,t13,0,1,centrej) } database[i,]<-bh } if(which.min(tk)==2){ ## transition 02. bh<-c(i,tk[2],0,v1[i],v2[i],1,0,0,0,0,centrej) database[i,]<-bh } if(which.min(tk)==3){ ## transition 03. bh<-c(i,tk[3],0,v1[i],v2[i],0,1,0,0,0,centrej) database[i,]<-bh } }} return(database) } ### ==== End function Simul.semi.markov =================================### ############################################################################################## ## gene.weibull generates survival time fom a proportional hazard weibull intensity ## X represents the predictor : exp(beta1*v1+beta2*v2) ## para: vector of the parameter of the weibull baseline intensity gene.weibull <- function(para,X){ b <- para[2] a <- para[1] return(b*((-(1/X)*log(1-runif(1)))**(1/a))) } ############################################################################################### ### ======================================= function analysis.multi-state ==================### ## Analysis.multi-state: function which perform a multi-state model ## non-homogeneous markov ("non-homogeneous") or semi-markov model ("semi") ## parametric (Weibull) estimation or penalized likelihood approach (penalized) Analysis.multi.state <- function(data=mydata,model="semi",estimator="penalized"){ n <- dim(data)[1] data1<-data[which(data\$PASS01==1),] ## data include observation which have transit from 0 to 1 n.VAP <- dim(data1)[1] if(estimator=="penalized"){ res01 <-frailtyPenal(Surv(TIME0,PASS01)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data,cross.validation=TRUE,n.knots=8,kappa1=1000,Frailty=TRUE,hazard="Splines") res02 <-frailtyPenal(Surv(TIME0,PASS02)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data,cross.validation=TRUE,n.knots=8,kappa1=1000,Frailty=TRUE,hazard="Splines") res03 <-frailtyPenal(Surv(TIME0,PASS03)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data,cross.validation=TRUE,n.knots=8,kappa1=1000,Frailty=TRUE,hazard="Splines") if(model=="semi"){ ### SEMI-MARKOV ### res12 <-frailtyPenal(Surv(TIME1,PASS12)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data1,cross.validation=TRUE,n.knots=8,kappa1=1000,Frailty=TRUE,hazard="Splines") res13 <-frailtyPenal(Surv(TIME1,PASS13)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data1,cross.validation=TRUE,n.knots=8,kappa1=1000,Frailty=TRUE,hazard="Splines") }else{ #### MARKOV NON-HOMOGENEOUS ### res12 <-frailtyPenal(Surv(TIME0,TIME0+TIME1,PASS12)~as.factor(v1)+v2+cluster(IDCENTRE),data=data1,cross.validation=TRUE,n.knots=8,kappa1=1000,Frailty=TRUE,hazard="Splines") res13 <-frailtyPenal(Surv(TIME0,TIME0+TIME1,PASS13)~as.factor(v1)+v2+cluster(IDCENTRE),data=data1,cross.validation=TRUE,n.knots=8,kappa1=1000,Frailty=TRUE,hazard="Splines") } LCV <- ((res01\$LCV+res01\$LCV+res01\$LCV)*n+(res12\$LCV+res13\$LCV)*n.VAP)/n }else{ res01 <-frailtyPenal(Surv(TIME0,PASS01)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data,Frailty=TRUE,hazard="Weibull") res02 <-frailtyPenal(Surv(TIME0,PASS02)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data,hazard="Weibull") res03 <-frailtyPenal(Surv(TIME0,PASS03)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data,hazard="Weibull") if(model=="semi"){ res12 <-frailtyPenal(Surv(TIME1,PASS12)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data1,Frailty=TRUE,hazard="Weibull") res13 <-frailtyPenal(Surv(TIME1,PASS13)~ as.factor(v1)+v2+cluster(IDCENTRE),data=data1,Frailty=TRUE,hazard="Weibull") }else{ res12 <-frailtyPenal(Surv(TIME0,TIME0+TIME1,PASS12)~as.factor(v1)+v2+cluster(IDCENTRE),data=data1,Frailty=TRUE,hazard="Weibull") res13 <-frailtyPenal(Surv(TIME0,TIME0+TIME1,PASS13)~as.factor(v1)+v2+cluster(IDCENTRE),data=data1,Frailty=TRUE,hazard="Weibull") } LCV <- ((res01\$AIC+res01\$AIC+res01\$AIC)*n+(res12\$AIC+res13\$AIC)*n.VAP)/n } res <- list(res01=res01,res02=res02,res03=res03,res12=res12,res13=res13,LCV.multi.state=LCV) return(res) } ### ==================End function Analysis.multi-state =======================================================### ################################################################################################################### ## Example on a simulated data set ## We perform the 4 model: semi-Markov (parametric and semi-parametric); ## Markov non-homogeneous (paremetric and semi-parametric ## We present also a table to choose the best model based on the LCV criterion ################################################################################################################### set.seed(133) mydata <- Simul.semi.markov(G=30,ni=150,paragamma=paragamma,p=0.5,loi01,loi02,loi03,loi12,loi13) model.semi.markov.penalized <- Analysis.multi.state(data=mydata,model="semi",estimator="penalized") model.non.homogeneous.penalized <- Analysis.multi.state(data=mydata,model="non-homogeneous",estimator="penalized") model.semi.markov.weibull <- Analysis.multi.state(data=mydata,model="semi",estimator="weibull") model.non.homogeneous.weibull <- Analysis.multi.state(data=mydata,model="non-homogeneous",estimator="weibull") ###table for the choice of the model tab <- matrix(c(model.non.homogeneous.weibull\$LCV,model.semi.markov.weibull\$LCV,model.non.homogeneous.weibull\$LCV-model.semi.markov.weibull\$LCV ,model.non.homogeneous.penalized\$LCV,model.semi.markov.penalized\$LCV ,model.non.homogeneous.penalized\$LCV-model.semi.markov.penalized\$LCV, model.non.homogeneous.weibull\$LCV-model.non.homogeneous.penalized\$LCV,model.semi.markov.weibull\$LCV-model.semi.markov.penalized\$LCV ,NA ),nrow=3,ncol=3,byrow=T) rownames(tab) <- c("parametric","semi-parametric","difference") tab <- as.data.frame(tab) colnames(tab) <- c("Non-homogeneous","Semi-Markov","difference") print(tab)