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

###############################################################
## Homogeneous Poisson
# Simple POisson

# Extra-Poisson

# Negative-Binomial


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


#############
# SIMPLE COX MODEL for srivBladderAG DATA
# with STRATA (enum) == 5 ERROR
# BIG coefficient
# large factor(enum)5 -17.26606 0.0000000317
#
coxph(Surv(time1, time2, status)~ treat + size + number + factor(enum),
          data=survBladderAG)

#######
# CORRECT COX DATA and MODEL
# survBladderAGrev DATA
#
coxph(Surv(time1, time2, status)~ treat + size + number + factor(enum),
          data=survBladderAGrev)
          

#######
# SIMPLE COX MODEL COMPARISON
#
coxph(Surv(time1, time2, status)~ treat + size + number, data=survBladderAG)
coxph(Surv(time1, time2, status)~ treat + size + number, data=survBladderAGrev)

#######
# CORRECT AG MODEL
# survBladderAGrev DATA
#
bladder.AG<-coxph(Surv(time1, time2, status)~ treat + size + number
           + cluster(id), data=survBladderAGrev)
summary(bladder.AG)

#######
# CORRECT WLW MODEL
# survBladderWLW DATA
#
bladder.WLW<-coxph(Surv(futime, status)~ treat + size + number
           + strata(enum) + cluster(id), data=survBladderWLW)
summary(bladder.WLW)


#######
# CORRECT PWP MODEL
# survBladderAGrev DATA
#
bladder.PWP<-coxph(Surv(time1, time2, status)~ treat + size + number
           + strata(enum) + cluster(id), data=survBladderAGrev)
summary(bladder.PWP)



#########################################################3
#######
# CORRECT WLW TREAT*STRATA MODEL
# survBladderWLW DATA
#
#
coxph(Surv(futime, status)~  treat
           + size + number
           + treat*strata(enum) + cluster(id), data=survBladderWLW)

# beta_treat + beta_treat*enum2 = beta treat2
coxph(Surv(futime, status)~  treat1 + treat2 + treat3 + treat4
           + size + number
           + strata(enum) + cluster(id), data=survBladderWLW)

###########################################3
#######
# CORRECT PWP TREAT*STRATA MODEL
# survBladderAGrev DATA
#
#
coxph(Surv(time1, time2, status)~  treat
           + size + number
           + treat*strata(enum) + cluster(id), data=survBladderAGrev)

# beta_treat + beta_treat*enum2 = beta treat2
coxph(Surv(time1, time2, status)~  treat1 + treat2 + treat3 + treat4
           + size + number
           + strata(enum) + cluster(id), data=survBladderAGrev)


###########
# PWP
# CF AG DATA
#
coxph(Surv(time1, time2, status)~  treat
           + size + number
           + treat*strata(enum) + cluster(id), data=survBladderAG)

# beta_treat + beta_treat*enum2 = beta treat2
coxph(Surv(time1, time2, status)~  treat1 + treat2 + treat3 + treat4
           + size + number
           + strata(enum) + cluster(id), data=survBladderAG)






##########################################################
# simple gain used in Book
##########################################################
# CONSIDER recurrent events
# Assume INDENEPDENT
# calculate efficiency gain
# SIZE is not important, so delete first
#
PWPfitStrata1<-coxph(Surv(time1, time2, status)~ treat + size + number
           + strata(enum) + cluster(id), data=survBladderAGrev,
           subset=(enum==1))
PWPfitStrata1

## refit ALL recurrent events
## different from Therneau book page 202
PWPfitStrataAL<-coxph(Surv(time1, time2, status)~ treat + size + number
           + strata(enum) + cluster(id), data=survBladderAGrev,
           init=PWPfitStrata1$coef, iter=0)
PWPfitStrataAL

# simple gain
112-47
(EFFgain<-sqrt(47/112)*(0.3158))  # approximate ~= uncorrelated s.e. 0.2010

## different from Therneau book page 202
# expect gain --> 112-47=65 subjects
47/(((0.2045748/0.3158)**2)) # ->  112.0000

# real gain
47/(((0.2445/0.3158)**2))  # ->  78.4
78.4-47
31.4/65