#Run replicates

nrep<-1000
ests.mat.psr<-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)
	
	Y<-eta.Y+rnorm(n)*sig.Y
	#boxplot(Y~Z)

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

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

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

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

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


#########################################################################
	#Regression on the propensity score: all predictors

	p1<-glm(Z~X,family=binomial)
	ps.vals<-fitted(p1)

	psr.fit1<-lm(Y~Z+ps.vals)
	ests.mat.psr[irep,6]<-coef(psr.fit1)[2]

	#Regression on the propensity score: using predictors of treatment only

	p1<-glm(Z~X[,treatment.list],family=binomial)
	ps.vals<-fitted(p1)
	
	psr.fit2<-lm(Y~Z+ps.vals)
	ests.mat.psr[irep,7]<-coef(psr.fit2)[2]


	#Regression on the propensity score: using confounders only

	p1<-glm(Z~X[,confound.list],family=binomial)
	ps.vals<-fitted(p1)
	
	psr.fit3<-lm(Y~Z+ps.vals)
	ests.mat.psr[irep,8]<-coef(psr.fit3)[2]


}

mse.vals<-n*apply((ests.mat.psr-theta)^2,2,mean)
name.list<-c('Comp.','Cor.','Unadj.','Inc.','Conf.','PSR-All','PSR-Tmt','PSR-Conf')
boxplot(ests.mat.psr,names=name.list,ylim=range(-2,6),main='Estimates')
abline(h=theta,col='red')
text(1:ncol(ests.mat.psr),rep(6,ncol(ests.mat.psr)),format(mse.vals,digits=4))
fname<-paste(filecode,'-PS-Reg.pdf',sep='')
dev.print(device = pdf, file=fname,width=11)


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