#Run replicates

nrep<-1000
ests.mat<-matrix(0,nrow=nrep,ncol=8)

for(irep in 1:nrep){
	#Predictors
	X<-mvrnorm(n,mu=rep(0,p),Sigma=Sig)

	#Data simulation
	Xmat<-cbind(rep(1,n),X)
	eta.Z<-Xmat %*% (c(1,al.select)*al)
	p.Z<-expit(eta.Z)
	Z<-rbinom(n,1,p.Z)
	eta.Y<-Z*theta+Xmat %*% (c(1,be.select)*be)

	if(irep == 1){
		fname<-paste(filecode,'-TruePS.pdf',sep='')
		pdf(file=fname,paper='USr',width=11,height=9)
		hist(p.Z,breaks=seq(0,1,by=0.05),main='True PS values');box()
		dev.off()
	}
	
	Y<-eta.Y+rnorm(n)*sig.Y
	#boxplot(Y~Z)

	#Complete outcome regression
	f0<-lm(Y~Z+X)
	ests.mat[irep,1]<-coef(f0)[2]

	#Correct outcome regression
	f1<-lm(Y~Z+X[,outcome.list])
	ests.mat[irep,2]<-coef(f1)[2]

	#Unadjusted outcome regression
	f2<-lm(Y~Z)
	ests.mat[irep,3]<-coef(f2)[2]

	#Incorrect outcome regression
	f3<-lm(Y~Z+X[,1:5])
	ests.mat[irep,4]<-coef(f3)[2]

	#Confounders only outcome regression
	f4<-lm(Y~Z+X[,confound.list])
	ests.mat[irep,5]<-coef(f4)[2]

	#Stratification by the propensity score: all predictors

	p1<-glm(Z~X,family=binomial)
	ps.vals<-fitted(p1)
	
	nstrata<-8
	n.in.stratum<-n/nstrata
	nreq<-n.in.stratum*0.1
	ps.q<-c(0,quantile(ps.vals,prob=c(1:(nstrata-1))/nstrata),1)
	ps.gp<-as.numeric(cut(ps.vals,ps.q))
	Zbygp<-table(Z,ps.gp)
	comparable.list<-c(1:nstrata)[Zbygp[1,] > nreq & Zbygp[2,] > nreq]

	n.common<-length(comparable.list)*n.in.stratum

	effect<-0
	for(i in 1:length(comparable.list)){
		igp<-comparable.list[i]
		Ysub<-Y[ps.gp==igp]	
		Zsub<-Z[ps.gp==igp]
		effect<-effect+(mean(Ysub[Zsub==1])-mean(Ysub[Zsub==0]))
	}
	ests.mat[irep,6]<-effect/length(comparable.list)

	#Stratification by the propensity score: using predictors of treatment only

	p1<-glm(Z~X[,treatment.list],family=binomial)
	ps.vals<-fitted(p1)
	
	nstrata<-8
	n.in.stratum<-n/nstrata
	nreq<-n.in.stratum*0.1
	ps.q<-c(0,quantile(ps.vals,prob=c(1:(nstrata-1))/nstrata),1)
	ps.gp<-as.numeric(cut(ps.vals,ps.q))
	Zbygp<-table(Z,ps.gp)
	comparable.list<-c(1:nstrata)[Zbygp[1,] > nreq & Zbygp[2,] > nreq]

	n.common<-length(comparable.list)*n.in.stratum

	effect<-0
	for(i in 1:length(comparable.list)){
		igp<-comparable.list[i]
		Ysub<-Y[ps.gp==igp]	
		Zsub<-Z[ps.gp==igp]
		effect<-effect+(mean(Ysub[Zsub==1])-mean(Ysub[Zsub==0]))
	}
	ests.mat[irep,7]<-effect/length(comparable.list)

	#Stratification by the propensity score: using confounders only

	p1<-glm(Z~X[,confound.list],family=binomial)
	ps.vals<-fitted(p1)
	
	nstrata<-8
	n.in.stratum<-n/nstrata
	nreq<-n.in.stratum*0.1
	ps.q<-c(0,quantile(ps.vals,prob=c(1:(nstrata-1))/nstrata),1)
	ps.gp<-as.numeric(cut(ps.vals,ps.q))
	Zbygp<-table(Z,ps.gp)
	comparable.list<-c(1:nstrata)[Zbygp[1,] > nreq & Zbygp[2,] > nreq]

	n.common<-length(comparable.list)*n.in.stratum

	effect<-0
	for(i in 1:length(comparable.list)){
		igp<-comparable.list[i]
		Ysub<-Y[ps.gp==igp]	
		Zsub<-Z[ps.gp==igp]
		effect<-effect+(mean(Ysub[Zsub==1])-mean(Ysub[Zsub==0]))
	}
	ests.mat[irep,8]<-effect/length(comparable.list)

}

mse.vals<-n*apply((ests.mat-theta)^2,2,mean)
name.list<-c('Comp.','Cor.','Unadj.','Inc.','Conf.','PS-All','PS-Tmt','PS-Conf')
boxplot(ests.mat,names=name.list,ylim=range(-2,6),main='Estimates')
abline(h=theta,col='red')
text(1:ncol(ests.mat),rep(6,ncol(ests.mat)),format(mse.vals,digits=4))

fname<-paste(filecode,'-PS-Strat.pdf',sep='')
dev.print(device = pdf, file=fname,width=11)

boxplot(Y~Z,names=c('Untreated','Treated'),main='Outcomes')