#-----------------------------------
# National Wilms' Tumour Group data
#----------------------------------

set.seed(2017-12-3)
nwts <- read.table("nwts-share.txt", header=TRUE)
summary(nwts)


#--------------
# New variables
#-------------

# A linear spline in age
nwts$age1 <- with(nwts, pmin(age, 1))
nwts$age2 <- with(nwts, pmax(age, 1))

#-----------------------
# The full-cohort model
#-----------------------

# Histology:age interaction and stage:tumour-diameter interaction
fullmodel <- glm(relaps~histol*(age1+age2)+ I(stage>2)*tumdiam, family=binomial, data=nwts)
summary(fullmodel)


#---------------------
# Case-control sample
#--------------------

# Here we take all the cases and a random sample of controls
nwts$id <- 1:nrow(nwts)
cases <- subset(nwts, relaps==1)
noncases <- subset(nwts, relaps==0)
controlsample <- sample(noncases$id, nrow(cases))
ccsample<- rbind(cases, noncases[noncases$id %in% controlsample,])
ccsample$weight<-with(ccsample, ifelse(relaps==1, 1, nrow(noncases)/nrow(cases)))

# We can compare the maximum likelihood estimator and the survey estimator:
library(survey)
ccmle <- glm(relaps~offset(log(weight))+histol*(age1+age2)+ I(stage>2)*tumdiam, family=binomial, data=ccsample)
summary(ccmle)
survey_cc <- svydesign(id=~1, weights=~weight, strata=~relaps, data=ccsample)
ccest <- svyglm(relaps~histol*(age1+age2)+ I(stage>2)*tumdiam, family=quasibinomial, design=survey_cc)
summary(ccest)

# The two seem fairly comparable: we know that asymptotically the maximum likelihood estimator must be better, 
# but the difference is small enough to not show up in a single comparison
round(cbind(coef(ccmle), coef(ccest))-coef(fullmodel),3)

# The simple survey estimator does not use the full cohort; we can declare a `twophase` object that does. 
# We do not need to specify weights because the software can work out what they are.
nwts_twophase <- twophase(id=list(~1,~1), strata=list(NULL, ~relaps), subset=~I((relaps==1)| id %in% controlsample), data=nwts)
twophaseest <- svyglm(relaps~histol*(age1+age2)+ I(stage>2)*tumdiam, family=quasibinomial, design=nwts_twophase)
summary(twophaseest)

# We still aren't using the whole cohort for anything, so the two analyses are almost identical


#-----------------------
# Using the whole cohort
#-----------------------

# We'll try to use the whole cohort now. First, just use `instit` instead of `histol`
# First, fit the model to the full data
phase1model <- glm(relaps~instit*(age1+age2)+ I(stage>2)*tumdiam, family=binomial, data=nwts)

# Extract the influence functions and create a new design object
inffun<-model.matrix(phase1model)*resid(phase1model, type="response")
colnames(inffun)<-paste0("if",1:ncol(inffun))
aug_twophase <- twophase(id=list(~1,~1), strata=list(NULL, ~relaps), subset=~I((relaps==1)| id %in% controlsample), data=cbind(nwts,inffun), method="simple")

# Calibrate, and fit the model of interest (ie, with `histol`) to the calibrated sample
calformula <- make.formula(colnames(inffun))
cal_twophase <- calibrate(aug_twophase, calformula, phase=2)
svyest_instit<-svyglm(relaps~histol*(age1+age2)+ I(stage>2)*tumdiam, family=quasibinomial, design=cal_twophase)
summary(svyest_instit)

# Comparing the uncalibrated and calibrated estimates, the coefficients have nearly all moved closer to the true full cohort value.
round(cbind(coef(twophaseest), coef(svyest_instit))-coef(fullmodel),3)


#---------------------------
# Calibration by imputation
#--------------------------

# It is always valid to just use a surrogate such as `instit` in calibration, but it is probably not optimal. 
# The attenuation bias in using a mismeasured predictor translates into a loss of precision in the calibrated estimate.  
# We can try to construct a regression imputation of histology instead:
impmodel<-svyglm(histol~instit*(relaps+I(stage>3))+I(age>10)+factor(study),family=quasibinomial,design=nwts_twophase)
nwts$imphistol <-as.vector(predict(impmodel,newdata=nwts,type="response",se.fit=FALSE))
with(nwts, by(imphistol, histol, summary))

# We now proceed as before. In particular, note that it is important to use `imphistol` for all observations in the phase-1 model, 
# even those where `histol` is available -- it was not available at phase 1.
phase1model_imp <- glm(relaps~imphistol*(age1+age2)+ I(stage>2)*tumdiam, family=binomial, data=nwts)

# Extract the influence functions and create a new design object
inffun_imp<-model.matrix(phase1model_imp)*resid(phase1model_imp, type="response")
colnames(inffun_imp)<-paste0("if",1:ncol(inffun_imp))
aug_twophase_imp <- twophase(id=list(~1,~1), strata=list(NULL, ~relaps), subset=~I((relaps==1)| id %in% controlsample), data=cbind(nwts,inffun_imp), method="simple")

# Calibrate, and fit the model of interest
calformula <- make.formula(colnames(inffun_imp))
cal_twophase_imp <- calibrate(aug_twophase_imp, calformula, phase=2)
svyest_imp<-svyglm(relaps~histol*(age1+age2)+ I(stage>2)*tumdiam, family=quasibinomial, design=cal_twophase_imp)
summary(svyest_imp)

# There has been a slight additional improvement; slight, because `instit` is overwhelming the best predictor of histology.
round(cbind(coef(twophaseest), coef(svyest_instit), coef(svyest_imp))-coef(fullmodel),3)

#----------------------
# Multiple imputation
#----------------------

# We can take multiple predictions from the fitted imputation model, produce influence functions for each of them, 
#then average the influence functions over the multiple imputations.  We need to use the full influence function, 
#not just the score function as before, because the matrix it's multiplied by can vary between imputations.
inffun_mi<-matrix(0,nrow=nrow(nwts),ncol=ncol(inffun_imp))
M<-50
for(i in 1:M){
nwts$mihistol <-rbinom(nrow(nwts),1, as.vector(predict(impmodel,newdata=nwts,type="response",se.fit=FALSE)))
phase1model_imp <- glm(relaps~mihistol*(age1+age2)+ I(stage>2)*tumdiam, family=binomial, data=nwts)
inffun_mi<-inffun_mi+(model.matrix(phase1model_imp)*resid(phase1model_imp, type="response"))%*%vcov(phase1model)
}
inffun_mi<-inffun_mi/M
colnames(inffun_mi)<-paste0("if",1:ncol(inffun_mi))

# Now, do the calibration and fit the model of interest
aug_twophase_mi <- twophase(id=list(~1,~1), strata=list(NULL, ~relaps), subset=~I((relaps==1)| id %in% controlsample), data=cbind(nwts,inffun_mi), method="simple")
calformula <- make.formula(colnames(inffun_mi))
cal_twophase_mi <- calibrate(aug_twophase_mi, calformula, phase=2)
svyest_mi<-svyglm(relaps~histol*(age1+age2)+ I(stage>2)*tumdiam, family=quasibinomial, design=cal_twophase_mi)
summary(svyest_mi)
round(cbind(coef(twophaseest), coef(svyest_instit),coef(svyest_imp),coef(svyest_mi))-coef(fullmodel),3)