###################################
# Create by Bhramar Mukerjee      #
# Modified By Jaeil Ahn					  #
# G*E Interaction Model						#
# Power calculation by 5 methods  #
# Last update: 11.30.2007					#
###################################

result<-NULL;

# input parameters
# niter : the number of replications
# nca, nco : the number of case and control
# setge : OR before log-transformation
# setint : interaction before log-transformation
# base setting of Prevalence P_g and P_e, default :0.3  : pg ,pe
#	Default OR10=1; OR01=1;  
# ex :test<-GEsimulation(1000, 500, 500, 1.2, 0.05, 1.2, 0.3, 0.3, 1, 1)
GEsimulation<-function(niter,nca, nco, setge, alpha,  setint, pg, pe, OR10, OR01)
{

	#Following is an rcode for generating G-E data with random parameters
	# the number of simulations : default 5000 simulated datasets 
	nrep=niter;
	valpha=qnorm(1-alpha/2, 0,1 );
	
	# vector of the three estimates(CO : case only, CC : case control)
	vec_co=rep(0,nrep); 
	vec_cc=rep(0,nrep); 
	vec_prop=rep(0, nrep);
	
	# vector of the variance of the three estimates
	varvec_co=rep(0,nrep); 
	varvec_cc=rep(0,nrep);
	
	vec_ts=rep(0,nrep);
	varvec_ts=rep(0,nrep);
	 
	vec_prop=rep(0, nrep)
	varprop=rep(0, nrep)

	vec_propplus=rep(0, nrep)
	varpropplus=rep(0, nrep)

	varvec_tau=rep(0,nrep); 
	
	orgerec=rep(0,nrep); 
	orinterrec=rep(0,nrep); 
	
	vec_proptrue=rep(0,nrep);
	postprobm1=rep(0,nrep); 
	postprobm2=rep(0,nrep); 
	bmaestimate=rep(0,nrep);
	
	vec_propnew=rep(0,nrep); 
	vec_propnewplus=rep(0,nrep); 

	vec_re=rep(0,nrep); 
	vec_sigma=rep(0,nrep);
	vec_re2=rep(0,nrep); 
	
	varpropnew=rep(0,nrep); 
	varpropnewplus=rep(0,nrep); 
	ucipropnew=rep(0,nrep);
	lcipropnew=rep(0,nrep);
	
	vec_propnil=rep(0,nrep);
	msevec_propnil=rep(0,nrep);
	
	counter=rep(0,nrep); 
	nonnormal1=rep(0,nrep);
	nonnormalvar=rep(0,nrep);
	nonnormal2=rep(0,nrep);
	
	# vector of the variances of case-only and case-control estimate
	tau=rep(0,nrep); 
	k1=rep(0,nrep);
	k2=rep(0, nrep);

	#result vector
	GEresult_POWER<-NULL
	GEresult<-NULL
	
	# setting the number of case and control 
	# These numbers are changable (100, 200, 500)
	ncontrol=nco; 
	ncase=nca;
	
	# initialize the loop 
	i=0;
	# initialize the power setting 
	power_cc=0;
	power_co=0;
	power_bmaestimate=0;
	power_prop=0;
	power_propplus=0;

	power_propnew=0;
	power_propnewplus=0;
	power_re=0;
	power_ts=0;
	
	
	# base setting of OR: among(1, 1.25, 1.5, 2) 
	# later we will take log transformation. 
	orge=setge;


	#Following is the equation to change control probabilities
	#in each run
	if (orge==1) orge=1.00001;  # buffer to make denominator safe to run.
	f1=1-2*orge-pe+orge*pe-pg+orge*pg;
	f2=sqrt((-1+2*orge+pe-orge*pe+pg-orge*pg)^2-4*(1-orge)*(-orge+orge*pe+orge*pg-orge*pe*pg));
	p01=(f1+f2)/(2*(1-orge)); 
	p02=1-pg-p01; 
	p03=1-pe-p01; 
	p04=1-p01-p02-p03;
	prcontrol=c(p01,p02,p03,p04);


	# randomly generate the interaction. fi
	# ORinter varies from 1.809 to 2.21: default 2, we will take log later 
	ORinter=setint;
	#ORinter=exp(rnorm(1,mean=log(2),sd=0.05));
	#orinterrec[i]=log(ORinter)
	
	# p : p11+p12+p13+p14 = 1
	p=p01+p02*OR10+p03*OR01+p04*OR10*OR01*ORinter;

	# p1x's are expressed in terms of p0x
	p11=p01/p; 
	p12=(p02*OR10)/p; 
	p13=(p03*OR01)/p; 
	p14=(p04*OR10*OR01*ORinter)/p;
	prcase=c(p11,p12,p13,p14);
			
	#the tau involved in the weight function
	for(i in 1:nrep)
	{
		freqcontrol=rmultinom(1,ncontrol,prob=prcontrol);
		freqcase=rmultinom(1,ncase,prob=prcase);
	
		# setting cells of 2*4 table count : r01 = n0*p01;
		r01=as.double(freqcontrol[1,1]);
		r02=as.double(freqcontrol[2,1]);
		r03=as.double(freqcontrol[3,1]); 
		r04=as.double(freqcontrol[4,1]);
		r11=as.double(freqcase[1,1]);
		r12=as.double(freqcase[2,1]);
		r13=as.double(freqcase[3,1]); 
		r14=as.double(freqcase[4,1]);
		
		if((r01>0) & (r02>0) & (r03>0) & (r04>0))
		{	
		} else 
		{
			freqcontrol=rmultinom(1,ncontrol,prob=prcontrol);
			r01=as.double(freqcontrol[1,1]);
			r02=as.double(freqcontrol[2,1]);
			r03=as.double(freqcontrol[3,1]); 
			r04=as.double(freqcontrol[4,1]);
		}
		if((r11>0) & (r12>0) & (r13>0) & (r14>0))
		{	
		} else 
		{
			freqcase=rmultinom(1,ncase,prob=prcase);
			r11=as.double(freqcase[1,1]);
			r12=as.double(freqcase[2,1]);
			r13=as.double(freqcase[3,1]); 
			r14=as.double(freqcase[4,1]);
		}
	
		# Case Only
		vec_co[i]=log((r11*r14)/(r12*r13)); 
		varvec_co[i]=(1/r11+1/r12+1/r13+1/r14);
	
		# Case-Control 
		vec_cc[i]=log((r02*r03*r11*r14)/(r01*r04*r12*r13)); 
		varvec_cc[i]=(1/r01+1/r02+1/r03+1/r04)+(1/r11+1/r12+1/r13+1/r14);
		
		tau[i]=log((r01*r04)/(r02*r03));
		varvec_tau[i]=(1/r01+1/r02+1/r03+1/r04);
	
		# Two-stages
		if(abs(tau[i]/sqrt(varvec_tau[i]))>1.96)
		{
			# use case control
			vec_ts[i]=vec_cc[i];
			varvec_ts[i]=varvec_cc[i];
		
		}else
		{
			# use case only
			vec_ts[i]=vec_co[i];
			varvec_ts[i]=varvec_co[i];
		}
	
	
		# OLD EB
		vec_prop[i]=((varvec_cc[i])/((tau[i])^2+varvec_cc[i]))*(vec_co[i])+ (((tau[i])^2)/((tau[i])^2+varvec_cc[i]))*(vec_cc[i]);
	
		# EB : 
		m3=max(0,(tau[i])^2-varvec_tau[i]);
		k1[i]=((tau[i])^2)*(3*varvec_cc[i]+(tau[i])^2)/((varvec_cc[i]+(tau[i])^2)^2);
		varprop[i]=varvec_co[i]+((k1[i])^2)*varvec_tau[i];
		
		# OLD EB PLUS
		#vec_prop[i]=((varvec_cc[i])/(m3+varvec_cc[i]))*(vec_co[i])+ ((m3)/(m3+varvec_cc[i]))*(vec_cc[i]);
	
		# EB  PLUS: 
		# k2[i]=((m3)*(3*varvec_cc[i]+m3)/((varvec_cc[i]+m3)^2)^2);
		# varpropplus[i]=varvec_co[i]+((k1[i])^2)*varvec_tau[i];
		
		# End of Not true bayesian
	
   	# Start of Empirical  Bayesian
		#Following is the code to employ bayesian model averaged estimate
		p01hat=r01/ncontrol; 
		p02hat=r02/ncontrol; 
		p03hat=r03/ncontrol; 
		p04hat=r04/ncontrol; 
		p11hat=r11/ncase; 
		p12hat=r12/ncase; 
		p13hat=r13/ncase; 
		p14hat=r14/ncase;
		
		l1=dmultinom(freqcontrol,prob=c(p01hat,p02hat,p03hat,p04hat))*dmultinom(freqcase,prob=c(p11hat,p12hat,p13hat,p14hat));
	
		p01hatcaseonly=(r01+r03)*(r01+r02)/(ncontrol)^2;
		p02hatcaseonly=(r01+r02)*(r02+r04)/(ncontrol)^2;
		p03hatcaseonly=(r01+r03)*(r03+r04)/(ncontrol)^2;
		p04hatcaseonly=(r02+r04)*(r03+r04)/(ncontrol)^2;
		l2=dmultinom(freqcontrol,prob=c(p01hatcaseonly,p02hatcaseonly,p03hatcaseonly,p04hatcaseonly))*dmultinom(freqcase,prob=c(p11hat,p12hat,p13hat,p14hat));
	
		pi1=(((tau[i])^2)/((tau[i])^2+varvec_cc[i]));
		pi2=1-pi1
		bicm1=-2*log(l1)+6*log(ncase+ncontrol);
		bicm2=-2*log(l2)+5*log(ncase+ncontrol);
		postprobm1[i]=exp(-bicm1/2)/(exp(-bicm1/2)+exp(-bicm2/2));
		postprobm2[i]=1-postprobm1[i];
		bmaestimate[i]=postprobm1[i]*vec_cc[i]+postprobm2[i]*vec_co[i];
	
		# mpostprobm1=mean(postprobm1[vec_cc<5]);
		# Following is the NEWEB estimate suggested by the referee by plugging in 
		# the EB estimate of theta
		# NEW EB
		# Empirical Bayes 
		vec_propnew[i]=((varvec_tau[i])/((tau[i])^2+varvec_tau[i]))*(vec_co[i])+ 
				(((tau[i])^2)/((tau[i])^2+varvec_tau[i]))*(vec_cc[i]);
	
		k1[i]=((tau[i])^2)*(3*varvec_tau[i]+(tau[i])^2)/((varvec_tau[i]+(tau[i])^2)^2);
		varpropnew[i]=varvec_co[i]+((k1[i])^2)*varvec_tau[i];
		
		# power calculation
		# if absolute standardized value is larger than 1.96, we count that. 
		if(abs(vec_cc[i]/sqrt(varvec_cc[i]))>valpha)
		{power_cc=power_cc+1};
	
		if(abs(vec_co[i]/sqrt(varvec_co[i]))>valpha)
		{power_co=power_co+1};
	
		if(abs(vec_ts[i]/sqrt(varvec_ts[i]))>valpha)
		{power_ts=power_ts+1};

		if(abs(vec_prop[i]/sqrt(varprop[i]))>valpha)
		{power_prop=power_prop+1};
		
		if(abs(vec_propnew[i]/sqrt(varpropnew[i]))>valpha)
		{power_propnew=power_propnew+1};
	}
		
	# Gain powers from each approach
	pow_cc=power_cc / nrep;
	pow_co=power_co / nrep;
	pow_ts=power_ts / nrep;
	pow_prop = power_prop /nrep;
	pow_propnew=power_propnew / nrep;

	GEresult_POWER<-cbind(GEresult_POWER,  niter, alpha, nca, nco,setge, setint ,pow_cc, pow_co, pow_prop, pow_propnew, pow_ts);
	return(GEresult_POWER);
}