###############################################################
## Repeated Frailty
###############################################################
setwd("X://temp//Rdata")
# library(Design)
options(scipen=10)
library(Hmisc)
library(MASS)
library(KMsurv)
library(survival)
na.action(c(".",NA))

################################################################
#################################################################
# READ IN DATA
# survBladderAG has serious problem
# some subjects actually FU time > 4th recurrence
# if we include stratum (enum) == 5,
# then as.factor(enum) will caluse enum==4 have infinite small risk
# sice we treat artifact of the data
# subjects who have the 4th recurrent event can not have the 5th event
# DATA SHOULD USE survBladderAGrev -> delete 5th stratum
##################################################################

survBladder<-read.csv("survBladderORI.csv",
  header = TRUE, sep = ",", dec=".")
survBladder

survBladderAG<-read.csv("survBladderAG.csv",
  header = TRUE, sep = ",", dec=".")
survBladderAG

survBladderAGrev<-read.csv("survBladderAGrev.csv",
  header = TRUE, sep = ",", dec=".")
survBladderAGrev


survBladderWLW<-read.csv("survBladderWLW.csv",
  header = TRUE, sep = ",", dec=".")
survBladderWLW



#################################################################



##################################################################
# ORDERED Recurrent Events FRAILTY
# AG

###########################################################
###########################################################
##########################################################
###########################################################
###########################################################
#
# AG frailty
#

###########
# AG : treat + size + number
AG.full.mar<-coxph(Surv(time1, time2, status)~ treat + size + number
           + cluster(id), data=survBladderAGrev)
summary(AG.full.mar)

# AG + GAMMA : treat + size + number
AG.full.ga<-coxph(Surv(time1, time2, status)~ treat + size + number
           + frailty(id), data=survBladderAGrev)
summary(AG.full.ga)


##########
## DELETE SIZE
# AG : treat + number
AG.tr.num.mar<-coxph(Surv(time1, time2, status)~ treat + number
           + cluster(id), data=survBladderAGrev)
summary(AG.tr.num.mar)

# AG + GAMMA : treat + number
AG.tr.num.ga<-coxph(Surv(time1, time2, status)~ treat + number
           + frailty(id), data=survBladderAGrev)
summary(AG.tr.num.ga)


###########
## DELETE NUMBER
# AG : treat + size
AG.tr.si.mar<-coxph(Surv(time1, time2, status)~ treat + size
           + cluster(id), data=survBladderAGrev)
summary(AG.tr.si.mar)

# AG + GAMMA : treat + size
AG.tr.si.ga<-coxph(Surv(time1, time2, status)~ treat + size
           + frailty(id), data=survBladderAGrev)
summary(AG.tr.si.ga)




###########
## DELETE SIZE + NUMBER
# AG : treat
AG.tr.mar<-coxph(Surv(time1, time2, status)~ treat
           + cluster(id), data=survBladderAGrev)
summary(AG.tr.mar)

# AG + GAMMA : treat
AG.tr.ga<-coxph(Surv(time1, time2, status)~ treat
           + frailty(id), data=survBladderAGrev)
summary(AG.tr.ga)


###########################################################
###########################################################
###########################################################
###########################################################
###########################################################
###########################################################
#
# WLW frailty
#
###########################################################
## DELETE SIZE
# WLW : treat + number
WLW.tr.num.mar<-coxph(Surv(futime, status)~ treat + number
           + cluster(id), data=survBladderWLW)
summary(WLW.tr.num.mar)

# WLW + GAMMA : treat + number
WLW.tr.num.ga<-coxph(Surv(futime, status)~ treat + number
           + frailty(id), data=survBladderWLW)
summary(WLW.tr.num.ga)


# WLW + GAUSSIAN : treat + number
WLW.tr.num.norm<-coxph(Surv(futime, status)~ treat + number
           + frailty(id, dist='gauss', method="reml"),
           data=survBladderWLW)
summary(WLW.tr.num.norm)



###########
## DELETE SIZE + NUMBER
# WLW : treat
WLW.tr.mar<-coxph(Surv(futime, status)~ treat
           + strata(enum) + cluster(id), data=survBladderWLW)
summary(WLW.tr.mar)

# WLW + GAMMA : treat
WLW.tr.ga<-coxph(Surv(futime, status)~ treat
           + strata(enum) + frailty(id), data=survBladderWLW)
summary(WLW.tr.ga)

# WLW + GAUSSIAN : treat
WLW.tr.norm<-coxph(Surv(futime, status)~ treat
           + strata(enum) + frailty(id, dist='gauss', method="reml"),
           data=survBladderWLW)
summary(WLW.tr.norm)





###########################################################
###########################################################
###########################################################
###########################################################
###########################################################
###########################################################
#
# PWP frailty
#
###########################################################
## DELETE SIZE
# PWP : treat + number
PWP.tr.num.mar<-coxph(Surv(time1, time2, status)~ treat + number
           + strata(enum) + cluster(id), data=survBladderAGrev)
summary(PWP.tr.num.mar)

# PWP + GAMMA : treat + number
PWP.tr.num.ga<-coxph(Surv(time1, time2, status)~ treat + number
           + strata(enum) + frailty(id), data=survBladderAGrev)
summary(PWP.tr.num.ga)


# PWP + GAUSSIAN : treat + number
PWP.tr.num.norm<-coxph(Surv(time1, time2, status)~ treat + number
           + strata(enum) + frailty(id, dist='gauss', method="reml"),
           data=survBladderAGrev)
summary(PWP.tr.num.norm)



###########
## DELETE SIZE + NUMBER
# PWP : treat
PWP.tr.mar<-coxph(Surv(time1, time2, status)~ treat
           + strata(enum) + cluster(id), data=survBladderAGrev)
summary(PWP.tr.mar)

# PWP + GAMMA : treat
PWP.tr.ga<-coxph(Surv(time1, time2, status)~ treat
           + strata(enum) + frailty(id), data=survBladderAGrev)
summary(PWP.tr.ga)

# PWP + GAUSSIAN : treat
PWP.tr.norm<-coxph(Surv(time1, time2, status)~ treat
           + strata(enum) + frailty(id, dist='gauss', method="reml"),
           data=survBladderAGrev)
summary(PWP.tr.norm)