rm(list=ls())
file.remove(list.files(pattern='.pdf'))

list.of.packages <- c("NHANES", "tableone", "Matching", "MatchIt","survey",
                       "twang","SuperLearner","glmnet","polspline","randomForest","SIS","Hmisc")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages,verbose=FALSE)

library(NHANES,verbose=FALSE)
library(SuperLearner,verbose=FALSE)
library(nnet,verbose=FALSE)
library(nnls,verbose=FALSE)
library(glmnet,verbose=FALSE)
library(polspline,verbose=FALSE)
library(randomForest,verbose=FALSE)
library(SIS,verbose=FALSE)
library(twang,verbose=FALSE)
library(tableone,verbose=FALSE)
library(survey,verbose=FALSE)
library(Matching,verbose=FALSE)
library(MatchIt,verbose=FALSE)
library(Hmisc,verbose=FALSE)

######################################################################
set.seed(37)

summary(NHANES)


NHANES$SmokeNow <- as.numeric(NHANES$SmokeNow)-1

small.nhanes <- na.omit(NHANES[NHANES$SurveyYr=="2011_12"& NHANES$Age > 17,c(3,4,8:11,13,25,61)])
dim(small.nhanes)  ## 1377 x 9
summary(small.nhanes)

#table(small.nhanes$Race3,small.nhanes$Education,small.nhanes$SmokeNow)
#table(small.nhanes$Race3,small.nhanes$MaritalStatus,small.nhanes$SmokeNow)
#table(small.nhanes$Education,small.nhanes$MaritalStatus,small.nhanes$SmokeNow)

temp.mat <- model.matrix(~Gender*Age+Gender*Race3+Gender*Education+Gender*MaritalStatus+Gender*HHIncome+Gender*Poverty+
                 Age*Race3+Age*Education+Age*MaritalStatus+Age*HHIncome+Age*Poverty+
                 Race3*Education+Race3*MaritalStatus+Race3*HHIncome+Race3*Poverty+
                 Education*MaritalStatus+Education*HHIncome+Education*Poverty+
                 MaritalStatus*HHIncome+MaritalStatus*Poverty+HHIncome*Poverty,data=small.nhanes)
colnames(temp.mat[,(apply(temp.mat,2,sum)<6)])
interact.data <- model.matrix(~Gender*Age+Gender*Race3+Gender*Education+Gender*MaritalStatus+Gender*HHIncome+Gender*Poverty+
                 Age*Race3+Age*Education+Age*MaritalStatus+Age*HHIncome+Age*Poverty+
                 Race3*Poverty+Education*Poverty+MaritalStatus*Poverty+HHIncome*Poverty,data=small.nhanes)
dim(interact.data)
interact.data <- data.frame(interact.data)
interact.data$SmokeNow <- small.nhanes$SmokeNow

vars.interact <- colnames(interact.data)[30:107]

## Check for balance
vars <- c("Gender", "Age", "Race3", "Education", "MaritalStatus", "Poverty")

tabUnmatched <- CreateTableOne(vars = vars, strata = "SmokeNow", data = small.nhanes, test = FALSE)
print(tabUnmatched, smd = TRUE)

plot(ecdf(small.nhanes$Age[small.nhanes$SmokeNow==0]),ylab="ECDF(Age)",xlab="Age",main="")
lines(ecdf(small.nhanes$Age[small.nhanes$SmokeNow==1]),col="red")

temp0 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==0])
temp1 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==1])
pdf("Age-eCFD-NHANES.pdf")
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
lines(temp1$x,temp1$y,col="red",lwd=3)
legend(20,1,c('SmokeNow=0','SmokeNow=1'),col=c('black','red'),lty=1,lwd=2)
dev.off()


##########
## on interactions?

SMDinteract <- CreateTableOne(vars = vars.interact, strata = "SmokeNow", data=interact.data,test = FALSE)
summary(ExtractSmd(SMDinteract))
boxplot(ExtractSmd(SMDinteract))



## Construct PS with logistic regression

ps.mod <- glm(SmokeNow~Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,data=small.nhanes,family="binomial")
ps.lr <- predict(ps.mod,type="response")
mean(as.numeric(ps.lr > 0.5) == small.nhanes$SmokeNow)

small.nhanes$ps.lr <- ps.lr
summary(ps.lr)
pdf("Overlap-LR.pdf")
boxplot(ps.lr[small.nhanes$SmokeNow==0],ps.lr[small.nhanes$SmokeNow==1],
             ylab="PS",xlab="Treatment Group",names=c(0,1),ylim=range(0,1),col='gray')
quints <- c(0,quantile(ps.lr,seq(.2,1,.2)))
quints
abline(h=quints[2:5],col="red",lwd=2)
title('PS Quintiles')
dev.off()
small.nhanes$ps.lr <- ps.lr

rbind(table(cut(ps.lr[small.nhanes$SmokeNow==0],quints)),table(cut(ps.lr[small.nhanes$SmokeNow==1],quints)))

## Check within quintiles
ps.lr.quints <- cut(ps.lr,quints,labels=1:5)
sum(table(cut(ps.lr,quints)))
summary(ps.lr.quints)
small.nhanes$ps.lr.quints <- ps.lr.quints

pdf("Age-eCFD-NHANES-PSquints-v1.pdf")
par(mfrow=c(2,3),mar=c(4,4,2,1))
for(j in 1:5) {
  plot(ecdf(small.nhanes$Age[small.nhanes$SmokeNow==0 & ps.lr.quints==j]),ylab="ECDF(Age)",xlab="Age",main="")
  lines(ecdf(small.nhanes$Age[small.nhanes$SmokeNow==1 & ps.lr.quints==j]),col="red")
}
dev.off()


pdf("Age-eCFD-NHANES-PSquints.pdf")
par(mfrow=c(2,3),mar=c(4,4,2,1))
for(j in 1:5) {
  nonsmj <- small.nhanes$Age[small.nhanes$SmokeNow==0 & small.nhanes$ps.lr.quints==j]
  temp0 <- Ecdf(nonsmj,pl=FALSE)
  smj <- small.nhanes$Age[small.nhanes$SmokeNow==1 & small.nhanes$ps.lr.quints==j]
  temp1 <- Ecdf(smj,pl=FALSE)
  plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
  lines(temp1$x,temp1$y,col="red",lwd=3)
}
dev.off()



SMD.table <- ExtractSmd(tabUnmatched)
for(j in 1:5) {
  tabPSquints <- CreateTableOne(vars = vars, strata = "SmokeNow", data = small.nhanes[ps.lr.quints==j,], test = FALSE)
  SMD.table <- cbind(SMD.table,ExtractSmd(tabPSquints))
}
round(SMD.table,3)

Max.SMD <- max(SMD.table)
Mean.SMD <- mean(SMD.table)
Med.SMD <- median(SMD.table)


run.this <- TRUE
if(run.this) {
	## deciles?
	pdf("Overlap-LR-dec.pdf")
	boxplot(ps.lr[small.nhanes$SmokeNow==0],ps.lr[small.nhanes$SmokeNow==1],
             ylab="PS",xlab="Treatment Group",names=c(0,1),col='gray')
		dec <- c(0,quantile(ps.lr,seq(.1,1,.1)))
		dec
		abline(h=dec[2:10],col="red",lwd=2)
		title('PS Deciles')
	dev.off()

	rbind(table(cut(ps.lr[small.nhanes$SmokeNow==0],dec)),
	table(cut(ps.lr[small.nhanes$SmokeNow==1],dec)))

	ps.lr.dec <- cut(ps.lr,dec,labels=1:10)
	SMD.10.table <- ExtractSmd(tabUnmatched)
	for(j in 1:10) {
  		tabPSdec <- CreateTableOne(vars = vars, strata = "SmokeNow", data = small.nhanes[ps.lr.dec==j,], test = FALSE)
  		SMD.10.table <- cbind(SMD.10.table,ExtractSmd(tabPSdec))
	}
	round(SMD.10.table,3)
	summary(as.vector(SMD.10.table))
}

## Check for IPW
ps.lr.weight <- small.nhanes$SmokeNow/ps.lr + (1-small.nhanes$SmokeNow)/(1-ps.lr)
temp0 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==0],weights=ps.lr.weight[small.nhanes$SmokeNow==0])
temp1 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==1],weights=ps.lr.weight[small.nhanes$SmokeNow==1])
pdf("Age-eCFD-NHANES-IPW.pdf")
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
lines(temp1$x,temp1$y,col="red",lwd=3)
legend(20,1,c('SmokeNow=0','SmokeNow=1'),col=c('black','red'),lty=1,lwd=2)
dev.off()

nhanes.IPW.lr <- svydesign(ids=~0, data=small.nhanes, weights=ps.lr.weight)
tabIPW <- svyCreateTableOne(vars = vars, strata = "SmokeNow", data = nhanes.IPW.lr, test = FALSE)
print(tabIPW, smd = TRUE)
round(cbind(SMD.table,ExtractSmd(tabIPW)),3)

##########
## on interactions?
interact.data$ps.lr.weight <- ps.lr.weight
nhanes.IPW.lr.interact <- svydesign(ids=~0, data=interact.data, weights=ps.lr.weight)
SMDinteract.IPW.a <- svyCreateTableOne(vars = vars.interact[1:25],
       strata = "SmokeNow", data = nhanes.IPW.lr.interact, test = FALSE)
SMDinteract.IPW.b <- svyCreateTableOne(vars = vars.interact[26:50],
       strata = "SmokeNow", data = nhanes.IPW.lr.interact, test = FALSE)
SMDinteract.IPW.c <- svyCreateTableOne(vars = vars.interact[51:78],
       strata = "SmokeNow", data = nhanes.IPW.lr.interact, test = FALSE)
SMDinteract.IPW <- c(ExtractSmd(SMDinteract.IPW.a),ExtractSmd(SMDinteract.IPW.b),ExtractSmd(SMDinteract.IPW.c))
summary(SMDinteract.IPW)
boxplot(cbind(ExtractSmd(SMDinteract),SMDinteract.IPW))
interact.table <- cbind(ExtractSmd(SMDinteract),SMDinteract.IPW)



## Matching sample
ps.lr.match <- Match(Tr=small.nhanes$SmokeNow,X=small.nhanes$ps.lr,estimand="ATE",ties=FALSE)
matched.samp <- small.nhanes[c(ps.lr.match$index.control,ps.lr.match$index.treated),]
summary(lm(matched.samp$BPSysAve~matched.samp$SmokeNow))  ## Same result
table(table(c(ps.lr.match$index.control,ps.lr.match$index.treated)))  

temp0 <- Ecdf(matched.samp$Age[matched.samp$SmokeNow==0])
temp1 <- Ecdf(matched.samp$Age[matched.samp$SmokeNow==1])
pdf("Age-eCFD-NHANES-PSmatch.pdf")
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
lines(temp1$x,temp1$y,col="red",lwd=3)
legend(20,1,c('SmokeNow=0','SmokeNow=1'),col=c('black','red'),lty=1,lwd=2)
dev.off()

MatchBalance(SmokeNow~Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,
                   data=small.nhanes,match.out=ps.lr.match)

tabMatched <- CreateTableOne(vars = vars, strata = "SmokeNow", data = matched.samp, test = FALSE)
round(cbind(SMD.table,ExtractSmd(tabMatched),ExtractSmd(tabIPW)),3)

Max.SMD <- c(NA,Max.SMD,max(ExtractSmd(tabMatched)),NA,max(ExtractSmd(tabIPW)))
Mean.SMD <- c(NA,Mean.SMD,mean(ExtractSmd(tabMatched)),NA,mean(ExtractSmd(tabIPW)))
Med.SMD <- c(NA,Med.SMD,median(ExtractSmd(tabMatched)),NA,median(ExtractSmd(tabIPW)))

ps.SMDtable <- round(cbind(SMD.table,ExtractSmd(tabMatched),ExtractSmd(tabIPW)),3)



####################################################
## GBM with twang
gbm.fit <- ps(SmokeNow~Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty, estimand = "ATE",data=as.data.frame(small.nhanes))
ps.gbm <- gbm.fit$ps$ks.mean.ATE 
small.nhanes$ps.gbm <- ps.gbm
mean(as.numeric(ps.gbm > 0.5) == small.nhanes$SmokeNow)

summary(ps.gbm)
boxplot(ps.gbm[small.nhanes$SmokeNow==0],ps.gbm[small.nhanes$SmokeNow==1],col='gray')
plot(ps.lr,ps.gbm,pch=19)  ## correlate pretty well
points(ps.lr[small.nhanes$SmokeNow==1],ps.gbm[small.nhanes$SmokeNow==1],col="red")  ## correlate very badly!!

pdf("Overlap-GBM.pdf")
boxplot(ps.gbm[small.nhanes$SmokeNow==0],ps.gbm[small.nhanes$SmokeNow==1],
             ylab="PS (GBM)",xlab="Treatment Group",names=c(0,1),col='gray')
quints.gbm <- c(0,quantile(ps.gbm,seq(.2,1,.2)))
quints.gbm
abline(h=quints.gbm[2:5],col="red",lwd=2)
dev.off()
small.nhanes$ps.gbm <- ps.gbm

rbind(table(cut(ps.gbm[small.nhanes$SmokeNow==0],quints.gbm)),
table(cut(ps.gbm[small.nhanes$SmokeNow==1],quints.gbm)))


ps.gbm.weight <- small.nhanes$SmokeNow/ps.gbm + (1-small.nhanes$SmokeNow)/(1-ps.gbm)

hist(ps.lr.weight);hist(ps.gbm.weight) 

nhanes.IPW.gbm <- svydesign(ids=~0, data=small.nhanes, weights=ps.gbm.weight)
tabIPW.gbm <- svyCreateTableOne(vars = vars, strata = "SmokeNow", data = nhanes.IPW.gbm, test = FALSE)
#round(cbind(ExtractSmd(tabUnmatched),ExtractSmd(tabIPW),ExtractSmd(tabIPW.gbm),ExtractSmd(tabIPW.SL)),3)
round(cbind(ExtractSmd(tabUnmatched),ExtractSmd(tabIPW),ExtractSmd(tabIPW.gbm)),3)

##########
## on interactions?
interact.data$ps.gbm.weight <- ps.gbm.weight
nhanes.IPW.gbm.interact <- svydesign(ids=~0, data=interact.data, weights=ps.gbm.weight)
SMDinteract.IPW.gbm.a <- svyCreateTableOne(vars = vars.interact[1:25],
       strata = "SmokeNow", data = nhanes.IPW.gbm.interact, test = FALSE)
SMDinteract.IPW.gbm.b <- svyCreateTableOne(vars = vars.interact[26:50],
       strata = "SmokeNow", data = nhanes.IPW.gbm.interact, test = FALSE)
SMDinteract.IPW.gbm.c <- svyCreateTableOne(vars = vars.interact[51:78],
       strata = "SmokeNow", data = nhanes.IPW.gbm.interact, test = FALSE)
SMDinteract.IPW.gbm <- c(ExtractSmd(SMDinteract.IPW.gbm.a),ExtractSmd(SMDinteract.IPW.gbm.b),ExtractSmd(SMDinteract.IPW.gbm.c))
summary(SMDinteract.IPW.gbm)
interact.table <- cbind(interact.table,SMDinteract.IPW.gbm)
boxplot(interact.table)



ps.gbm.quints <- cut(ps.gbm,quints.gbm,labels=1:5)
SMD.gbm.table <- NULL
for(j in 1:5) {
  tabPSquints.gbm <- CreateTableOne(vars = vars, strata = "SmokeNow", data = small.nhanes[ps.gbm.quints==j,], test = FALSE)
  SMD.gbm.table <- cbind(SMD.gbm.table,ExtractSmd(tabPSquints.gbm))
}
round(SMD.gbm.table,3)

ps.gbm.match <- Match(Tr=small.nhanes$SmokeNow,X=small.nhanes$ps.gbm,estimand="ATE",ties=FALSE)
matched.gbm.samp <- small.nhanes[c(ps.gbm.match$index.control,ps.gbm.match$index.treated),]
MatchBalance(SmokeNow~Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,
                   data=small.nhanes,match.out=ps.gbm.match)
tabMatched.gbm <- CreateTableOne(vars = vars, strata = "SmokeNow", data = matched.gbm.samp, test = FALSE)


temp0 <- Ecdf(matched.gbm.samp$Age[matched.gbm.samp$SmokeNow==0])
temp1 <- Ecdf(matched.gbm.samp$Age[matched.gbm.samp$SmokeNow==1])
pdf("Age-eCFD-NHANES-GBM-PSmatch.pdf")
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
lines(temp1$x,temp1$y,col="red",lwd=3)
legend(20,1,c('SmokeNow=0','SmokeNow=1'),col=c('black','red'),lty=1,lwd=2)
dev.off()

temp0 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==0],weights=ps.gbm.weight[small.nhanes$SmokeNow==0])
temp1 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==1],weights=ps.gbm.weight[small.nhanes$SmokeNow==1])
pdf("Age-eCFD-NHANES-GBM-IPW.pdf")
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
lines(temp1$x,temp1$y,col="red",lwd=3)
legend(20,1,c('SmokeNow=0','SmokeNow=1'),col=c('black','red'),lty=1,lwd=2)
dev.off()

SMD.gbm.table <- cbind(SMD.gbm.table,ExtractSmd(tabMatched.gbm),ExtractSmd(tabIPW.gbm))
round(SMD.gbm.table,3)


## SuperLearner
X.mat <- data.frame(cbind(small.nhanes$Gender,small.nhanes$Age,small.nhanes$Race3,
    small.nhanes$Education,small.nhanes$MaritalStatus,small.nhanes$HHIncome,small.nhanes$Poverty))
my.library <- c("SL.knn","SL.randomForest","SL.glmnet","SL.mean")
SL.fit <- SuperLearner(Y = small.nhanes$SmokeNow, X = X.mat, 
                 SL.library = my.library,verbose = TRUE, method ="method.NNLS", 
                 family=binomial())
small.nhanes$ps.SL <- ps.SL <- SL.fit$SL.predict
mean(as.numeric(ps.SL > 0.5) == small.nhanes$SmokeNow)

summary(ps.SL)
boxplot(ps.SL[small.nhanes$SmokeNow==0],ps.SL[small.nhanes$SmokeNow==1],col='gray')
title('PS using SuperLearner')
plot(ps.lr,ps.SL,pch=19)  ## correlate very badly!!
points(ps.lr[small.nhanes$SmokeNow==1],ps.SL[small.nhanes$SmokeNow==1],col="red",pch=19)  ## correlate very badly!!

plot(ps.gbm,ps.SL,pch=19)  ## correlate badly
points(ps.gbm[small.nhanes$SmokeNow==1],ps.SL[small.nhanes$SmokeNow==1],col="red",pch=19)  ## correlate very badly!!

pdf("Overlap-SL.pdf")
boxplot(ps.SL[small.nhanes$SmokeNow==0],ps.SL[small.nhanes$SmokeNow==1],
             ylab="PS (SL)",xlab="Treatment Group",names=c(0,1),col='gray')
quints.SL <- c(0,quantile(ps.SL,seq(.2,1,.2)))
quints.SL
abline(h=quints.SL[2:5],col="red",lwd=2)
dev.off()

rbind(table(cut(ps.SL[small.nhanes$SmokeNow==0],quints.SL)),
table(cut(ps.SL[small.nhanes$SmokeNow==1],quints.SL)))

ps.SL[ps.SL==0] <- min(ps.SL[ps.SL!=0])
ps.SL.weight <- small.nhanes$SmokeNow/ps.SL + (1-small.nhanes$SmokeNow)/(1-ps.SL)
nhanes.IPW.SL <- svydesign(ids=~0, data=small.nhanes, weights=ps.SL.weight)
tabIPW.SL <- svyCreateTableOne(vars = vars, strata = "SmokeNow", data = nhanes.IPW.SL, test = FALSE)
ExtractSmd(tabIPW.SL)

ps.SL.quints <- cut(ps.SL,quints.SL,labels=1:5)
SMD.SL.table <- NULL
for(j in 3) {
  tabPSquints.SL <- CreateTableOne(vars = vars, strata = "SmokeNow", data = small.nhanes[ps.SL.quints==j,], test = FALSE)
  SMD.SL.table <- cbind(SMD.SL.table,ExtractSmd(tabPSquints.SL))
}
round(SMD.SL.table,3)

ps.SL.match <- Match(Tr=small.nhanes$SmokeNow,X=small.nhanes$ps.SL,estimand="ATE",ties=FALSE)
matched.SL.samp <- small.nhanes[c(ps.SL.match$index.control,ps.SL.match$index.treated),]
MatchBalance(SmokeNow~Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,
                   data=small.nhanes,match.out=ps.SL.match)
tabMatched.SL <- CreateTableOne(vars = vars, strata = "SmokeNow", data = matched.SL.samp, test = FALSE)

temp0 <- Ecdf(matched.SL.samp$Age[matched.SL.samp$SmokeNow==0])
temp1 <- Ecdf(matched.SL.samp$Age[matched.SL.samp$SmokeNow==1])
pdf("Age-eCFD-NHANES-SL-PSmatch.pdf")
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
lines(temp1$x,temp1$y,col="red",lwd=3)
legend(20,1,c('SmokeNow=0','SmokeNow=1'),col=c('black','red'),lty=1,lwd=2)
dev.off()

temp0 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==0],weights=ps.SL.weight[small.nhanes$SmokeNow==0])
temp1 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==1],weights=ps.SL.weight[small.nhanes$SmokeNow==1])
pdf("Age-eCFD-NHANES-SL-IPW.pdf")
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
lines(temp1$x,temp1$y,col="red",lwd=3)
dev.off()


SMD.SL.table <- cbind(SMD.SL.table,ExtractSmd(tabMatched.SL),ExtractSmd(tabIPW.SL))
round(SMD.SL.table,3)

round(cbind(ExtractSmd(tabUnmatched),ExtractSmd(tabIPW),ExtractSmd(tabIPW.gbm),ExtractSmd(tabIPW.SL)),3)


##########
## on interactions?
interact.data$ps.SL.weight <- ps.SL.weight
nhanes.IPW.SL.interact <- svydesign(ids=~0, data=interact.data, weights=ps.SL.weight)
SMDinteract.IPW.SL.a <- svyCreateTableOne(vars = vars.interact[1:25],
       strata = "SmokeNow", data = nhanes.IPW.SL.interact, test = FALSE)
SMDinteract.IPW.SL.b <- svyCreateTableOne(vars = vars.interact[26:50],
       strata = "SmokeNow", data = nhanes.IPW.SL.interact, test = FALSE)
SMDinteract.IPW.SL.c <- svyCreateTableOne(vars = vars.interact[51:78],
       strata = "SmokeNow", data = nhanes.IPW.SL.interact, test = FALSE)
SMDinteract.IPW.SL <- c(ExtractSmd(SMDinteract.IPW.SL.a),ExtractSmd(SMDinteract.IPW.SL.b),ExtractSmd(SMDinteract.IPW.SL.c))
summary(SMDinteract.IPW.SL)
interact.table <- cbind(interact.table,SMDinteract.IPW.SL)
pdf("InteractBalance-IPW.pdf")
boxplot(interact.table,xlab="Approach", names=c("Unadjusted","IPW: LR", "IPW: GBM", "IPW: SL"),col='gray')
title('Balance on Interactions')
dev.off()

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

## Effect estimates ##
# Naive: simple model
summary(lm(BPSysAve~SmokeNow,data=small.nhanes))  ## -3.68

#Regression
summary(lm(BPSysAve~SmokeNow+Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,data=small.nhanes))  ## -1.10

#PS regression: quintiles
summary(lm(BPSysAve~SmokeNow+ps.lr.quints,data=small.nhanes))  ## -1.41

#PS Regression
summary(lm(BPSysAve~SmokeNow+ps.lr,data=small.nhanes))  ## -1.11

#PS Regression with quadratic term
summary(lm(BPSysAve~SmokeNow+ps.lr+I(ps.lr^2),data=small.nhanes))  ## -1.11

#IPW
summary(lm(BPSysAve~SmokeNow,weights=ps.lr.weight,data=small.nhanes))  ## -1.99

#Matching
matched.anal <- Match(Y=small.nhanes$BPSysAve, Tr=small.nhanes$SmokeNow, X=X.mat, estimand = "ATE", ties=FALSE)
matched.anal <- Match(Y=small.nhanes$BPSysAve, Tr=small.nhanes$SmokeNow, X=ps.lr, estimand = "ATE", ties=FALSE)
summary(matched.anal) # Note that this value changes from match to match due to randomness in matching algorithm


nhanes.allsmoke <- small.nhanes
nhanes.allsmoke$SmokeNow <- 1 
nhanes.nosmoke <- small.nhanes
nhanes.nosmoke$SmokeNow <- 0

all.ests <- NULL

## Regression
coef(lm(BPSysAve~SmokeNow,data=small.nhanes))[2] 
coef(lm(BPSysAve~SmokeNow+Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,data=small.nhanes))[2]

## ATE via regression
mod1.lm <- lm(BPSysAve~SmokeNow+Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,data=small.nhanes)
APO.lm.1 <- mean(predict(mod1.lm,nhanes.allsmoke))
APO.lm.0 <- mean(predict(mod1.lm,nhanes.nosmoke))
APO.lm.1 - APO.lm.0
all.ests <- c(all.ests,APO.lm.1 - APO.lm.0)

mod1.lmX <- lm(BPSysAve~SmokeNow+Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty+SmokeNow:HHIncome+SmokeNow:Gender+SmokeNow:Age,data=small.nhanes)
APO.lmX.1 <- mean(predict(mod1.lmX,nhanes.allsmoke))
APO.lmX.0 <- mean(predict(mod1.lmX,nhanes.nosmoke))
APO.lmX.1 - APO.lmX.0

## ATE via PS stratification
ps.lr.quints <- cut(ps.lr,quints,labels=1:5)
sum(table(cut(ps.lr,quints)))
summary(ps.lr.quints)
small.nhanes$ps.lr.quints <- ps.lr.quints
p.strat <- table(ps.lr.quints)/length(ps.lr.quints)
p.strat 

ATE.strat <- rep(NA,5)
for(j in 1:5) {
   ATE.strat[j] <- mean(small.nhanes$BPSysAve[small.nhanes$SmokeNow == 1 & small.nhanes$ps.lr.quints==j]) - 
     				mean(small.nhanes$BPSysAve[small.nhanes$SmokeNow == 0 & small.nhanes$ps.lr.quints==j])
}
ATE.strat
sum(ATE.strat*p.strat)
all.ests <- c(all.ests,sum(ATE.strat*p.strat))


## ATE via matching
dim(matched.samp)
mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 1]) - 
   mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 0])
all.ests <- c(all.ests,mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 1]) - 
   mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 0]))

## ATE via PS regression
mod1.PSlm1 <- lm(BPSysAve~SmokeNow+ps.lr,data=small.nhanes)
APO.PSlm1.1 <- mean(predict(mod1.PSlm1,nhanes.allsmoke))
APO.PSlm1.0 <- mean(predict(mod1.PSlm1,nhanes.nosmoke))
APO.PSlm1.1 - APO.PSlm1.0

mod1.PSlm2 <- lm(BPSysAve~SmokeNow+ps.lr+I(ps.lr^2),data=small.nhanes)
APO.PSlm2.1 <- mean(predict(mod1.PSlm2,nhanes.allsmoke))
APO.PSlm2.0 <- mean(predict(mod1.PSlm2,nhanes.nosmoke))
APO.PSlm2.1 - APO.PSlm2.0

mod1.PSlm3 <- lm(BPSysAve~SmokeNow+bs(ps.lr,df=4),data=small.nhanes)
APO.PSlm3.1 <- mean(predict(mod1.PSlm3,nhanes.allsmoke))
APO.PSlm3.0 <- mean(predict(mod1.PSlm3,nhanes.nosmoke))
APO.PSlm3.1 - APO.PSlm3.0
all.ests <- c(all.ests,APO.PSlm3.1 - APO.PSlm3.0)


## ATE via IPW
small.nhanes$ps.lr.weight <- small.nhanes$SmokeNow/small.nhanes$ps.lr + (1-small.nhanes$SmokeNow)/(1-small.nhanes$ps.lr)
IPW.est<-mean(small.nhanes$SmokeNow*small.nhanes$BPSysAve*small.nhanes$ps.lr.weight) - 
    mean((1-small.nhanes$SmokeNow)*small.nhanes$BPSysAve*small.nhanes$ps.lr.weight) 
all.ests <- c(all.ests,IPW.est)


coef(lm(BPSysAve ~ SmokeNow, weights = ps.lr.weight,data=small.nhanes))

mean(small.nhanes$SmokeNow*small.nhanes$BPSysAve/small.nhanes$ps.lr) -   
 mean((1-small.nhanes$SmokeNow)*small.nhanes$BPSysAve/(1-small.nhanes$ps.lr)) 


# Final table:
round(rbind(Max.SMD, Mean.SMD, Med.SMD, all.ests), 3)

# ATT
matched.anal.ATT <- Match(Y=small.nhanes$BPSysAve, Tr=small.nhanes$SmokeNow, X=ps.lr, estimand = "ATT", ties=FALSE)
summary(matched.anal.ATT) 

matched.samp.ATT <- small.nhanes[c(matched.anal.ATT$index.control,matched.anal.ATT$index.treated),]
dim(matched.samp.ATT)
mean(matched.samp.ATT$BPSysAve[matched.samp.ATT$SmokeNow == 1]) - 
   mean(matched.samp.ATT$BPSysAve[matched.samp.ATT$SmokeNow == 0])

temp0 <- Ecdf(matched.samp.ATT$Age[matched.samp.ATT$SmokeNow==0])
temp1 <- Ecdf(matched.samp.ATT$Age[matched.samp.ATT$SmokeNow==1])
pdf("Age-eCFD-NHANES-ATT-match.pdf")
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
lines(temp1$x,temp1$y,col="red",lwd=3)
dev.off()

ATT.match <- CreateTableOne(vars = vars, strata = "SmokeNow", data = matched.samp.ATT, test = FALSE)
SMD.ATT <- ExtractSmd(ATT.match)

## ATT via IPW
small.nhanes$ATT.lr.weight <- small.nhanes$SmokeNow + (1-small.nhanes$SmokeNow)*ps.lr/(1-ps.lr)
temp0 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==0],weights=ATT.lr.weight[small.nhanes$SmokeNow==0])
temp1 <- Ecdf(small.nhanes$Age[small.nhanes$SmokeNow==1],weights=ATT.lr.weight[small.nhanes$SmokeNow==1])
pdf("Age-eCFD-NHANES-ATT-IPW.pdf")
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",main="",type="l",lwd=3)
lines(temp1$x,temp1$y,col="red",lwd=3)
dev.off()

nhanes.ATT.IPW <- svydesign(ids=~0, data=small.nhanes, weights=ATT.lr.weight)
ATT.IPW <- svyCreateTableOne(vars = vars, strata = "SmokeNow", data = nhanes.ATT.IPW, test = FALSE)
print(ATT.IPW, smd = TRUE)
round(cbind(SMD.ATT,ExtractSmd(ATT.IPW)),3)

mean(small.nhanes$SmokeNow*small.nhanes$BPSysAve*small.nhanes$ATT.lr.weight) - 
   mean((1-small.nhanes$SmokeNow)*small.nhanes$BPSysAve*small.nhanes$ATT.lr.weight) 


#### with bootstrap

## ATE via IPW

small.nhanes$ps.lr.weight <- small.nhanes$SmokeNow/small.nhanes$ps.lr + (1-small.nhanes$SmokeNow)/(1-small.nhanes$ps.lr)
MainEst.lr <- mean(small.nhanes$SmokeNow*small.nhanes$BPSysAve*small.nhanes$ps.lr.weight) - 
   mean((1-small.nhanes$SmokeNow)*small.nhanes$BPSysAve*small.nhanes$ps.lr.weight) 

small.nhanes$ps.gbm.weight <- small.nhanes$SmokeNow/small.nhanes$ps.gbm + (1-small.nhanes$SmokeNow)/(1-small.nhanes$ps.gbm)
MainEst.gbm <- mean(small.nhanes$SmokeNow*small.nhanes$BPSysAve*small.nhanes$ps.gbm.weight) - 
   mean((1-small.nhanes$SmokeNow)*small.nhanes$BPSysAve*small.nhanes$ps.gbm.weight) 

small.nhanes$ps.SL[small.nhanes$ps.SL < 0.01] <- 0.01
small.nhanes$ps.SL.weight <- small.nhanes$SmokeNow/small.nhanes$ps.SL + (1-small.nhanes$SmokeNow)/(1-small.nhanes$ps.SL)
MainEst.SL <- mean(small.nhanes$SmokeNow*small.nhanes$BPSysAve*small.nhanes$ps.SL.weight) - 
   mean((1-small.nhanes$SmokeNow)*small.nhanes$BPSysAve*small.nhanes$ps.SL.weight) 


B <- 1000
n <- dim(small.nhanes)[1]

boot.lr <- boot.gbm <- boot.SL <- NULL
for(b in 1:B) {
  who <- sample(1:n,replace=TRUE)
  re.samp <- small.nhanes[who,]

  re.samp$ps.lr <- predict(glm(SmokeNow~Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,data=re.samp,family="binomial"),type="response")
  re.samp$ps.lr.weight <- re.samp$SmokeNow/re.samp$ps.lr + (1-re.samp$SmokeNow)/(1-re.samp$ps.lr)
  boot.lr <- c(boot.lr,mean(re.samp$SmokeNow*re.samp$BPSysAve*re.samp$ps.lr.weight) - 
     mean((1-re.samp$SmokeNow)*re.samp$BPSysAve*re.samp$ps.lr.weight)) 

  re.samp$gbm.fit <- ps(SmokeNow~Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,data=as.data.frame(re.samp),verbose=FALSE)$ps$ks.mean.ATE
  re.samp$ps.gbm.weight <- re.samp$SmokeNow/re.samp$ps.gbm + (1-re.samp$SmokeNow)/(1-re.samp$ps.gbm)
  boot.gbm <- c(boot.gbm,mean(re.samp$SmokeNow*re.samp$BPSysAve*re.samp$ps.gbm.weight) - 
     mean((1-re.samp$SmokeNow)*re.samp$BPSysAve*re.samp$ps.gbm.weight)) 

  X.mat.resamp <- X.mat[who,] 
  re.samp$ps.SL <- SuperLearner(Y = re.samp$SmokeNow, X = X.mat.resamp, 
                 SL.library = my.library,verbose = FALSE, method ="method.NNLS", 
                 family=binomial())$SL.predict
  re.samp$ps.SL[re.samp$ps.SL < 0.01] <- 0.01
  re.samp$ps.SL.weight <- re.samp$SmokeNow/re.samp$ps.SL + (1-re.samp$SmokeNow)/(1-re.samp$ps.SL)
  boot.SL <- c(boot.SL,mean(re.samp$SmokeNow*re.samp$BPSysAve*re.samp$ps.SL.weight) - 
     mean((1-re.samp$SmokeNow)*re.samp$BPSysAve*re.samp$ps.SL.weight))

if(b %in% seq(25,B,25)) print(paste("b=",b,sep=""))
}

par(mfrow=c(3,1))
hist(boot.lr)
hist(boot.gbm)
hist(boot.SL)