library(mvnfast)
set.seed(239)
expit<-function(x){return(1/(1+exp(-x)))}
n<-1000
Sig<-diag(c(0.1,0.5)) %*% (matrix(c(1,0.8,0.8,1),2,2)) %*% diag(c(0.1,0.5))
X<-rmvn(n,c(1,1),Sig)
etaZ<-1+X[,1]-2*X[,2]
eX<-expit(etaZ)
Z<-rbinom(n,1,eX)
etaY<-(2+3*X[,1]+X[,2]+X[,1]*X[,2])+Z
Y<-rnorm(n)+etaY

X1<-X[,1]
X2<-X[,2]
round(coef(summary(lm(Y~X1+X2+X1:X2+Z))),4)

round(coef(summary(lm(Y~X1+Z))),4)

round(coef(summary(lm(Y~X1+Z+eX))),4)

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

set.seed(239)

n<-1000
Sig<-diag(c(0.1,0.5)) %*% (matrix(c(1,0.8,0.8,1),2,2)) %*% diag(c(0.1,0.5))
X<-rmvn(n,c(1,1),Sig)
etaZ<-1+X[,1]-2*X[,2]
eX<-expit(etaZ)
Z<-rbinom(n,1,eX)
etaY<-(2+3*X[,1]+X[,2]+X[,1]*X[,2])+Z*(1+X[,1]+X[,2])
Y<-rnorm(n)+etaY

X1<-X[,1]
X2<-X[,2]
round(coef(summary(lm(Y~X1+X2+X1:X2+Z+Z:X1+Z:X2))),4)

round(coef(summary(lm(Y~X1+Z+Z:X1+Z:X2))),4)

round(coef(summary(lm(Y~X1+Z+Z:X1+Z:X2+eX))),4)


set.seed(239)
nreps<-5000
ests.mat<-array(0,c(nreps,4,3))

for(irep in 1:nreps){

	X<-rmvn(n,c(1,1),Sig)
	etaZ<-1+X[,1]-2*X[,2]
	eX<-expit(etaZ)
	Z<-rbinom(n,1,eX)
	etaY<-(2+3*X[,1]+X[,2]+X[,1]*X[,2])+Z*(1+X[,1]+X[,2])
	Y<-rnorm(n)+etaY

	X1<-X[,1]
	X2<-X[,2]

	ests.mat[irep,1,]<-coef(lm(Y~X1+X2+X1:X2+Z+Z:X1+Z:X2))[c(4,6,7)]
	ests.mat[irep,2,]<-coef(lm(Y~X1+Z+Z:X1+Z:X2))[c(3,4,5)]

	ests.mat[irep,3,]<-coef(lm(Y~X1+Z+Z:X1+Z:X2+eX))[c(3,5,6)]

	ests.mat[irep,4,]<-coef(lm(Y~X1+Z+Z:X1+Z:X2+eX+eX:X1+eX:X2))[c(3,5,6)]

}

pdf('PS1.pdf',paper='USr',width=11,height=9)
boxplot(ests.mat[,1:3,1],names=c('Correct','Incorrect','PS Regression'));abline(h=1,col='red',lty=2)
title(expression(psi[0]))
dev.off()
pdf('PS2.pdf',paper='USr',width=11,height=9)
boxplot(ests.mat[,1:3,2],names=c('Correct','Incorrect','PS Regression'));abline(h=1,col='red',lty=2)
title(expression(psi[1]))
dev.off()
pdf('PS3.pdf',paper='USr',width=11,height=9)
boxplot(ests.mat[,1:3,3],names=c('Correct','Incorrect','PS Regression'));abline(h=1,col='red',lty=2)
title(expression(psi[2]))
dev.off()


pdf('PS1x.pdf',paper='USr',width=11,height=9)
boxplot(ests.mat[,1:4,1],names=c('Correct','Incorrect','PS Regression','PS 2'));abline(h=1,col='red',lty=2)
title(expression(psi[0]))
dev.off()
pdf('PS2x.pdf',paper='USr',width=11,height=9)
boxplot(ests.mat[,1:4,2],names=c('Correct','Incorrect','PS Regression','PS 2'));abline(h=1,col='red',lty=2)
title(expression(psi[1]))
dev.off()
pdf('PS3x.pdf',paper='USr',width=11,height=9)
boxplot(ests.mat[,1:4,3],names=c('Correct','Incorrect','PS Regression','PS 2'));abline(h=1,col='red',lty=2)
title(expression(psi[2]))
dev.off()



vm<-apply(ests.mat,2:3,var)
vm[4,]/vm[1,]