#install.packages(c('ipw','gdata','car'))
library(ipw,gdata)
nhefs <- read.csv("nhefs_book.csv",header=TRUE)

## Prepare dataset
nhefs <- within(nhefs, {

    require(car)

    education <- car::recode(school,
                recodes = '0:8   = "1. 8th grade or less";
                           9:11  = "2. HS dropout";
                           12    = "3. HS";
                           13:15 = "4. College dropout";
                           16:30 = "5. College or more"')    
    education <- factor(education)

    university <- (school >= 16)

    age50 <- as.numeric(age > 50)

    ## Censoring indicator
    censored <- as.numeric(is.na(wt82_71))

    ## Reverse code to create non-censoring indicator
    uncensored <- as.numeric(censored == 0)
})

## Convert to factor
nhefs[,c("hbp","exercise","active","education")] <-
    lapply(nhefs[,c("hbp","exercise","active","education")], factor)

## Drop people without education
nhefs <- subset(nhefs, !is.na(education))

summary(nhefs)
###############################################################


## Load ipw package for IPW method
require(ipw)

boxplot(wt82_71~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Weight change (lbs) 1971-1982')
summary(lm(wt82_71~qsmk,data=nhefs))

attach(nhefs)
mean(wt82_71[qsmk==1],na.rm=T)-mean(wt82_71[qsmk==0],na.rm=T)
detach(nhefs)

boxplot(age~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Age')

quit.age.hist<-hist(nhefs$age[nhefs$qsmk==1])
no.quit.age.hist<-hist(nhefs$age[nhefs$qsmk==0])
plot(quit.age.hist, col=rgb(0,0,1,1/4),xlim=range(0,100),ylim=range(0,200))
plot(no.quit.age.hist, col=rgb(1,0,1,1/4),add=T)

summary(lm(age~qsmk,data=nhefs))

boxplot(income~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Income')


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

fq.age<-glm(qsmk~age+income+race,family=binomial,data=nhefs)
summary(fq.age)
ps.age<-rep(NA,nrow(nhefs))
ps.age[!is.na(nhefs$income)]<-fitted(fq.age)


quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,500),
	xlab='Propensity Score value',main='PS values')
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
legend(0,500,c('Quit=Y','Quit=N'),fill=c(rgb(0,0,1,1/4),rgb(1,0,1/2,1/4)))
box()

boxplot(ps.age~qsmk,data=nhefs)

psq<-quantile(ps.age,prob=c(0:5)/5)
ps.strata<-cut(ps.age,psq,include.lowest=TRUE,labels=FALSE)

f1<-lm(wt82_71~qsmk,subset=(ps.strata==1),data=nhefs);summary(f1)
f2<-lm(wt82_71~qsmk,subset=(ps.strata==2),data=nhefs);summary(f2)
f3<-lm(wt82_71~qsmk,subset=(ps.strata==3),data=nhefs);summary(f3)
f4<-lm(wt82_71~qsmk,subset=(ps.strata==4),data=nhefs);summary(f4)
f5<-lm(wt82_71~qsmk,subset=(ps.strata==5),data=nhefs);summary(f5)

(ps.adj.est<-(coef(f1)[2]+coef(f2)[2]+coef(f3)[2]+coef(f4)[2]+coef(f5)[2])/5)

f.all<-lm(wt82_71~qsmk*as.factor(ps.strata),subset=(!is.na(income)),data=nhefs);summary(f.all)