############################################
###R-code for data examples in chapter 19###
############################################

#---------------------------------------------------------------------------------------------
# Read cohort and nested case-control data for table 19.1, 19.3, 19.4 and Figure 19.2 and 19.3
#---------------------------------------------------------------------------------------------

cohort = read.table("tb.cohort.txt", header = T)
ncc = read.table("tb.ncc.txt", header = T)

#-------------
# Libraries
#-------------

library(survival)
library(timereg)     #timereg used for Lin-Ying additive model
library(survey)      #survey used for calibration
library(dplyr)       #dplyr used for data management
library(multipleNCC) #multipleNCC used for weighted analysis of NCC data 



#----------------------------------------------------------------------------------
#Table 19.1
#----------------------------------------------------------------------------------

#Recode logdose
cohort$logdose = log2(cohort$dose + 1)
ncc$logdose = log2(ncc$dose + 1)
#-------------------------------------
#Create data set for weighted analysis
#-------------------------------------
#Each subject only included once in weighted analyses
#Remove duplicated subjects, but retain the case record if the subject is sampled 
#as a control, and later become a case

ncc.red = ncc %>% group_by(id) %>% arrange(desc(case)) %>% 
  mutate(count = row_number()) %>% filter(count==1)

##Create samplestat indicator 
ncc.red$samplestat = 1 + ncc.red$case
ncc.merge = select(ncc.red, id, samplestat)

#Merge with cohort
cohort = left_join(cohort, ncc.merge)
cohort = mutate(cohort, samplestat = ifelse(is.na(samplestat), 0, samplestat))


#---------
#Analyses
#---------
cohortfit = coxph(Surv(age.in, age.out, bc) ~ logdose + I(age.rx/10) + cluster(id),
                  data = cohort) 
summary(cohortfit)

NCCtrad = coxph(Surv(age.in, age.out, case) ~ logdose + I(age.rx/10) + strata(setno), 
                data = ncc)
summary(NCCtrad) 

fitkm = wpl(Surv(age.in, age.out, bc) ~ logdose + I(age.rx/10), data = cohort, 
            samplestat = cohort$samplestat, m = 2) 
summary(fitkm)

fitglm = wpl(Surv(age.in, age.out, bc) ~ logdose + I(age.rx/10), data = cohort,
             samplestat = cohort$samplestat, weight.method = "glm")
summary(fitglm)

fitgam = wpl(Surv(age.in, age.out, bc) ~ logdose + I(age.rx/10), data = cohort,
             samplestat = cohort$samplestat, weight.method = "gam")
summary(fitgam)

fitchen = wpl(Surv(age.in, age.out, bc) ~ logdose + I(age.rx/10), data = cohort,
              samplestat = cohort$samplestat, weight.method = "Chen")
summary(fitchen)


#----------------------------------------------------------------------------------
#Table 19.2
#----------------------------------------------------------------------------------

#Attaching data set from multipleNCC pakcage
data(CVD_Accidents)


#Creating two synthetic NCC-datasets. It could alternatively have been done with the 
#ccwc-function in the Epi-package.

CVD_Accidents$id = 1:nrow(CVD_Accidents)
CVD_Accidents$stratum = 0
CVD_Accidents$case = 0
CVD_Accidents$sbpKat = as.numeric(CVD_Accidents$sbp > 140)


#------------
#CVD-endpoint
#------------
CVD_case = filter(CVD_Accidents, dead2 == 1)
CVD_NCC =  CVD_case
ant.CVD = sum(CVD_Accidents$dead2)
CVD_NCC$stratum = 1:ant.CVD

j = 1
empty = 0

set.seed(354)

#sampling controls for CVD-cases
for(i in 1:ant.CVD)  {
  possible = which(CVD_Accidents$agestart < CVD_case$agestop[i] & 
                     CVD_Accidents$agestop > CVD_case$agestop[i])
  
  if(length(possible) > 0) {
    control = sample(possible, 1)
    
    CVD_NCC = rbind(CVD_NCC, CVD_Accidents[control, ])
    CVD_NCC$case[i] = 1
    CVD_NCC$stratum[(j + ant.CVD)] = i
    
    j = j + 1 
    
  }
  if(length(possible) == 0)  {
    empty = empty + 1
  }
}

#-----------------
#Accident-endpoint
#-----------------
Acc_case = filter(CVD_Accidents, dead4 == 1)
Acc_NCC = Acc_case
ant.Acc = sum(CVD_Accidents$dead4)
Acc_NCC$stratum = 1:ant.Acc

j = 1
empty = 1
#sampling controls for ALAV-cases
for(i in 1:ant.Acc)  {
  possible = which(CVD_Accidents$agestart < Acc_case$agestop[i] & 
                     CVD_Accidents$agestop > Acc_case$agestop[i])
  
  if(length(possible) > 0) {
    control = sample(possible, 1)
    
    Acc_NCC = rbind(Acc_NCC, CVD_Accidents[control, ])
    Acc_NCC$case[i] = 1
    Acc_NCC$stratum[(j + ant.Acc)] = i
    
    j = j + 1 
    
  }
  if(length(possible) == 0)  {
    empty = empty + 1
  }
}

#Combine the two NCC studies 
NCCdataW = rbind(CVD_NCC, Acc_NCC)

#Each subject only included once in weighted analyses
#Remove duplicated subjects, but retain the case record if the subject is sampled 
#as a control, and later become a case

NCCdataW = NCCdataW %>% group_by(id) %>% arrange(desc(case)) %>% 
  mutate(count = row_number()) %>% filter(count == 1)


#Creating a sampling indicator which will be merged to the cohort data
NCCdataW$samplestat2 = 1
NCCdataW$samplestat2[NCCdataW$dead2 == 1] = 2
NCCdataW$samplestat2[NCCdataW$dead4 == 1] = 3


nccw.merge = select(NCCdataW, id, samplestat2)
CVD_Accidents = left_join(CVD_Accidents, nccw.merge)
CVD_Accidents = mutate(CVD_Accidents, samplestat2 = ifelse(is.na(samplestat2), 0, samplestat2))

##cohort analysis CVD-endpoint
cohortCVD = coxph(Surv(agestart, agestop, dead2) ~ sbpKat, data = CVD_Accidents)
summary(cohortCVD)

##cohort analysis ALAV-endpoint
cohortALAV = coxph(Surv(agestart, agestop, dead4) ~ sbpKat, data = CVD_Accidents)
summary(cohortALAV)


##Traditional analysis CVD-endpoint
TradCVD = coxph(Surv(agestart, agestop, case) ~ sbpKat + strata(stratum), data = CVD_NCC)
summary(TradCVD)

##Traditional analysis ALAV-endpoint
TradALAV = coxph(Surv(agestart, agestop, case) ~ sbpKat + strata(stratum), data = Acc_NCC)
summary(TradALAV)

##WPL analysis with KM-weights both endpoints
wplKM = wpl(Surv(agestart, agestop, dead24) ~ sbpKat, data = CVD_Accidents,
            samplestat = CVD_Accidents$samplestat2, variance = "Modelbased")
summary(wplKM)


#----------------------------------------------------------------------------------
#Table 19.3
#----------------------------------------------------------------------------------

#Recode logdose
cohort$logdose = log2(cohort$dose + 1)
ncc$logdose = log2(ncc$dose + 1)

#-------------------------------------
#Create data set for weighted analysis
#-------------------------------------
#Each subject only included once in weighted analyses
#Remove duplicated subjects, but retain the case record if the subject is sampled 
#as a control, and later become a case

ncc.red = ncc %>% group_by(id) %>% arrange(desc(case)) %>% 
  mutate(count = row_number()) %>% filter(count==1)

##Create samplestat indicator 
ncc.red$samplestat = 1 + ncc.red$case
ncc.merge = select(ncc.red, id, samplestat)

#Merge with the cohort
cohort = left_join(cohort, ncc.merge)
cohort = mutate(cohort , samplestat=ifelse(is.na(samplestat), 0, samplestat))

#Function for printing a summary of coefficients, standard errors etc. from 
#Lin-Ying model

sumcoef = function(lyobj){
  coeff = lyobj$gamma
  naivese = sqrt(diag(lyobj$var.gamma))
  robse = sqrt(diag(lyobj$robvar.gamma))
  Z = coeff/robse
  pval = 2*pnorm(-abs(Z))
  res = cbind(round(coeff, 6), round(naivese, 6), round(robse, 6),round(Z, 2), round(pval, 3))
  colnames(res) = c("estimate", "naive se", "robust se", "Z", "p-value")
  res
}


#Cohort analysis
linying.coh = aalen(Surv(age.in, age.out, bc) ~ const(logdose) + 
                      const(age.rx/10), data = cohort)
summary(linying.coh)
sumcoef(linying.coh)

##Estimate sampling probabilities
cohort$pkm = KMprob(cohort$age.out, cohort$samplestat, 2, cohort$age.in)

##Add samplinging probabilities to NCCdata
cohpkm = select(cohort, id, pkm)
ncc.red = inner_join(ncc.red, cohpkm)

linying.ncc = aalen(Surv(age.in, age.out, bc) ~ const(logdose) + const(age.rx/10),
                    data = ncc.red, weights = 1/ncc.red$pkm, robust = T)
summary(linying.ncc)
sumcoef(linying.ncc)

#----------------------------------------------------------------------------------
#Figure 19.2
#----------------------------------------------------------------------------------

#Recode logdose
cohort$logdose = log2(cohort$dose + 1)
ncc$logdose = log2(ncc$dose + 1)

#-------------------------------------
#Create data set for weighted analysis
#-------------------------------------
#Each subject only included once in weighted analyses
#Remove duplicated subjects, but retain the case record if the subject is sampled 
#as a control, and later become a case

ncc.red = ncc %>% group_by(id) %>% arrange(desc(case)) %>% 
  mutate(count = row_number()) %>% filter(count==1)

##Create samplestat indicator 
ncc.red$samplestat = 1 + ncc.red$case
ncc.merge = select(ncc.red, id, samplestat)

#Merge with the cohort
cohort = left_join(cohort, ncc.merge)
cohort = mutate(cohort , samplestat=ifelse(is.na(samplestat), 0, samplestat))

#Estimate sampling weights
cohort$pkm = KMprob(cohort$age.out, cohort$samplestat, 2, cohort$age.in)

##Run cohort model
coxfit = coxph(Surv(age.in, age.out, bc)~logdose + I(age.rx/10)+cluster(id), data = cohort)

##Run weighted NCC-model
coxfitw = coxph(Surv(age.in, age.out, bc)~logdose + I(age.rx/10) + cluster(id), 
                data = cohort, subset = samplestat > 0, weights = 1/pkm)

#plot

par(mfrow = c(1, 2))
plot(survfit(coxfit, newdata = data.frame(logdose = 0, age.rx = 0)), fun = "cumhaz", 
     mark.time = F, ylim=c(0, 0.22), xlim = c(25,80), xlab = "Age", ylab = "Cumulative hazard") 
text(41, 0.2, "Cohort")
plot(survfit(coxfitw, newdata = data.frame(logdose = 0, age.rx = 0)), fun = "cumhaz", 
     mark.time = F, ylim = c(0, 0.22), xlim = c(25, 80), xlab = "Age", ylab = "Cumulative hazard") 
text(39, 0.2, "NCC")


#----------------------------------------------------------------------------------
#Figure 19.3
#----------------------------------------------------------------------------------

#Recode logdose
cohort$logdose = log2(cohort$dose + 1)
ncc$logdose = log2(ncc$dose + 1)

#-------------------------------------
#Create data set for weighted analysis
#-------------------------------------
#Each subject only included once in weighted analyses
#Remove duplicated subjects, but retain the case record if the subject is sampled 
#as a control, and later become a case

ncc.red = ncc %>% group_by(id) %>% arrange(desc(case)) %>% 
  mutate(count = row_number()) %>% filter(count==1)

##Create samplestat indicator 
ncc.red$samplestat = 1 + ncc.red$case
ncc.merge = select(ncc.red, id, samplestat)

#Merge with the cohort
cohort = left_join(cohort, ncc.merge)
cohort = mutate(cohort , samplestat=ifelse(is.na(samplestat), 0, samplestat))

#Cohort analysis
linying.coh = aalen(Surv(age.in, age.out, bc) ~ const(logdose) + 
                      const(age.rx/10), data = cohort)

##Estimate sampling probabilities
cohort$pkm = KMprob(cohort$age.out, cohort$samplestat, 2, cohort$age.in)

##Add samplinging probabilities to NCCdata
cohpkm = select(cohort, id, pkm)
ncc.red = inner_join(ncc.red, cohpkm)

linying.ncc = aalen(Surv(age.in, age.out, bc) ~ const(logdose) + const(age.rx/10),
                    data = ncc.red, weights = 1/ncc.red$pkm, robust = T)

stepfunS = function (datax, datay, type = "left") {
  if (missing(datay)) 
    x <- datax$x
  else x <- datax
  if (missing(datay)) 
    y <- datax$y
  else y <- datay
  n <- length(x)
  type <- charmatch(type, c("left", "right"))
  if (is.na(type)) 
    stop("The type must be 'left' or 'right' continuous")
  if (any(diff(x) < 0)) 
    stop("The x vector must be sorted")
  if (type == 2) {
    x <- rev(x)
    y <- rev(y)
  }
  if (n > 2) {
    dupy <- c(T, diff(y[-n]) != 0, T)
    dupx <- c(T, diff(x[-n]) != 0, T)
    x <- x[dupx & dupy]
    y <- y[dupx & dupy]
    n <- length(x)
  }
  xrep <- rep(x[2:n], rep(2, n - 1))
  yrep <- rep(y[1:(n - 1)], rep(2, n - 1))
  if (type == 1) 
    list(x = c(x[1], xrep), y = c(yrep, y[n]))
  else list(x = c(rev(xrep), x[1]), y = c(y[n], rev(yrep)))
}

#Plot
par(mfrow = c(1, 2))
plot(stepfunS(linying.coh$cum[, 1], linying.coh$cum[, 2]), xlim = c(25, 80), 
     ylim = c(-0.019, 0.18), xlab="Age", ylab = "Cumulative intercept",type = "l")
lines(stepfunS(linying.coh$cum[, 1], linying.coh$cum[, 2] + 
                 1.96 * linying.coh$robvar.cum[, 2] ^ 0.5), lty=2)
lines(stepfunS(linying.coh$cum[, 1], linying.coh$cum[, 2] - 
                 1.96 * linying.coh$robvar.cum[, 2] ^ 0.5), lty=2)
text(36, 0.18, "Cohort")

plot(stepfunS(linying.ncc$cum[, 1], linying.ncc$cum[, 2]), xlim=c(25,80), 
     ylim = c(-0.019, 0.18), xlab = "Age", ylab = "Cumulative intercept", type="l")
lines(stepfunS(linying.ncc$cum[, 1], linying.ncc$cum[, 2] + 
                 1.96 * linying.ncc$robvar.cum[, 2] ^ 0.5), lty=2)
lines(stepfunS(linying.ncc$cum[, 1], linying.ncc$cum[, 2] -
                 1.96 * linying.ncc$robvar.cum[, 2] ^ 0.5), lty=2)
text(35, 0.18, "NCC")

#----------------------------------------------------------------------------------
#Table 19.4
#----------------------------------------------------------------------------------

#Recode logdose
cohort$logdose = log2(cohort$dose + 1)
ncc$logdose = log2(ncc$dose + 1)

#-------------------------------------
#Create data set for weighted analysis
#-------------------------------------
#Each subject only included once in weighted analyses
#Remove duplicated subjects, but retain the case record if the subject is sampled 
#as a control, and later become a case

ncc.red = ncc %>% group_by(id) %>% arrange(desc(case)) %>% 
  mutate(count = row_number()) %>% filter(count==1)

##Create samplestat indicator 
ncc.red$samplestat = 1 + ncc.red$case
ncc.merge = select(ncc.red, id, samplestat)

#Merge with the cohort
cohort = left_join(cohort, ncc.merge)
cohort = mutate(cohort , samplestat=ifelse(is.na(samplestat), 0, samplestat))


##Estimate sampling probabilities
cohort$pkm = KMprob(cohort$age.out, cohort$samplestat, 2, cohort$age.in)

##Add samplinging probabilities to NCCdata
cohpkm = select(cohort, id, pkm)
ncc.red = inner_join(ncc.red, cohpkm)

tsi_cohort = coxph(Surv(age.out - age.in, bc) ~ logdose + I(age.rx/10), robust = T, 
                   data = cohort)
summary(tsi_cohort)

tsi_km = coxph(Surv(age.out - age.in, bc) ~ logdose + I(age.rx/10), robust = T,
               data = ncc.red, weights = 1/pkm)
summary(tsi_km)



#----------------------------------------------------------------------------------
#Table 19.5
#----------------------------------------------------------------------------------
#Recode logdose
cohort$logdose = log2(cohort$dose + 1)
ncc$logdose = log2(ncc$dose + 1)

#-------------------------------------
#Create data set for weighted analysis
#-------------------------------------
#Each subject only included once in weighted analyses
#Remove duplicated subjects, but retain the case record if the subject is sampled 
#as a control, and later become a case

ncc.red = ncc %>% group_by(id) %>% arrange(desc(case)) %>% 
  mutate(count = row_number()) %>% filter(count==1)

#Create samplestat indicator 
ncc.red$samplestat = 1 + ncc.red$case
ncc.merge = select(ncc.red, id, samplestat)

#Merge with the cohort
cohort = left_join(cohort, ncc.merge)
cohort = mutate(cohort , samplestat=ifelse(is.na(samplestat), 0, samplestat))



#Creating strata based on inclusion and censoring time
breaks.in = as.numeric(quantile(cohort$age.in[cohort$samplestat == 1],
                                probs = c(0, 0.30, 0.55, 1)))
breaks.in[1] = min(cohort$age.in) - 1
breaks.in[4] = max(cohort$age.in) + 1
breaks.out = as.numeric(quantile(cohort$age.out[cohort$samplestat == 1],
                                 probs = c(0, 0.40, 0.63, 1)) )
breaks.out[1] = min(cohort$age.out) - 1
breaks.out[4] = max(cohort$age.out) + 1
cohort$left = as.numeric(cut(cohort$age.in,breaks = breaks.in, include.lowest = T,
                             right = FALSE))
cohort$right = as.numeric(cut(cohort$age.out, breaks = breaks.out, include.lowest = T))

cohort$strat = as.numeric(interaction(cohort$left, cohort$right))

#Collapsing the two smallest strata  
cohort$strat[cohort$strat == 7] = 8

#The cases are in a separate stratum
cohort$strat[cohort$bc == 1] = 10

cohort$in.sample = cohort$samplestat > 0

#Creating a twophase design
dtwophs = twophase(id=list(~id, ~id), subset =~ in.sample, data = cohort,
                   strata = list(~strat, ~strat), method = "full")

#Get the dfbetas using nof as a surrogate of logdose and age.rx 
cox.imp = coxph(Surv(age.in, age.out, bc) ~ I(log2(1 + nof/100)) + 
                  I(age.rx/10), data = cohort) 
imp.dfb = resid(cox.imp, type="dfbeta")
colnames(imp.dfb) = paste("imp.dfb", 1:ncol(imp.dfb), sep="")

#Get the dfbetas using only age.rx
cox.imp.age = coxph(Surv(age.in, age.out, bc) ~ I(age.rx/10), data = cohort) 
imp.dfb.age = resid(cox.imp.age, type="dfbeta")

cohort.imp = cbind(cohort, imp.dfb, imp.dfb.age)

#Add dfbeta to the design
dtwophs.imp = twophase(id = list(~id, ~id), subset =~ in.sample, data = cohort.imp,
                       strata = list(~strat, ~strat), method="full")

#Calibrate
dcalibr = calibrate(dtwophs.imp, phase = 2, formula =~ imp.dfb1 + imp.dfb2,
                    calfun = "raking")
dcalibr.age = calibrate(dtwophs.imp,phase = 2, formula =~ imp.dfb.age, calfun = "raking")

#Using svycoxph for analysis
calibr = svycoxph(Surv(age.in, age.out, bc) ~ logdose + I(age.rx/10), design = dcalibr)
summary(calibr)

calibr.age = svycoxph(Surv(age.in, age.out, bc) ~ logdose + I(age.rx/10), 
                      design = dcalibr.age)
summary(calibr.age)