## multiplicative model for genotype ###
set.seed(5)
library(MASS)
start=Sys.time()
k=1

#### define variable  ####

Case.T=rep(0,k)
Case2.T=rep(0,k)
bgstat.h=rep(0,k)
bgstat.b=rep(0,k)
SS.TT=rep(0,k)
bgestat.h=rep(0,k)
bgestat.b=rep(0,k)
Trio.tg=rep(0,k)
Trio.ge=rep(0,k)
hbgeg=rep(0,k)
bbgeg=rep(0,k)
Tdf.2=rep(0,k)
trio.b=rep(0,k)
trio.v=rep(0,k)
beta.hat=rep(0,k)
beta.bar=rep(0,k)
varb.hat=rep(0,k)
varb.bar=rep(0,k)
gamma1=rep(0,k)
gamma.var=rep(0,k)
r.rhat=matrix(0,k,2)
r.rbar=matrix(0,k,2)

cff1=matrix(0,k,2)
cff2=matrix(0,k,2)
cff3=matrix(0,k,2)

varr.hat=rep(0,k)
varr.bar=rep(0,k)

################### generate genotype of one-stage and two-stage #######

n=300
n3=600
b.g=0
q=0.5
tk=0.1

   
full.D=matrix(0,n,4)
Dom.D=matrix(0,n,4)
full2.D=matrix(0,n3,4)
Dom2.D=matrix(0,n3,4)
full3.c=matrix(0,(n3+n),4)
Dom3.c=matrix(0,(n3+n),4)

m1=1

while (m1<=n){
### stratification ####
u1=runif(1,0,1)
if(u1>q){p.Af=0.3;p.Am=0.3;b0=log(0.01)}
if(u1<=q){p.Af=0.3+tk;p.Am=0.3+tk;b0=log(0.05)}

p.af=1-p.Af
{
gf=runif(1,0,1);
gf[gf>0 & gf<=p.Af^2]=2;
gf[p.Af^2<gf&gf<=(p.Af^2+2*p.af*p.Af)]=1;
gf[gf>(p.Af^2+2*p.af*p.Af)& gf<1]=0
}
p.am=1-p.Am
{
gm=runif(1,0,1);
gm[gm>0 & gm<=p.Am^2]=2;
gm[p.Am^2<gm & gm<=(p.Am^2+2*p.am*p.Am)]=1;
gm[gm>(p.Am^2+2*p.am*p.Am)& gm<1]=0
}
if (gf==0 & gm==0) g0=0
if ((gf==0 & gm==1)|(gf==1 & gm==0)) g0=rbinom(1,1,0.5)
if ((gf==0 & gm==2)|(gf==2 & gm==0)) g0=1
if (gf==1 & gm==1) {g0=runif(1,0,1);
g0[0.25>=g0 & g0>0]=0 ;
g0[g0>0.25& g0<=0.75]=1; 
g0[g0>0.75& g0<1]=2}
if (gf==2 & gm==2) g0=2
if ((gf==1 & gm==2)|(gf==2 & gm==1)) {g0=runif(1,0,1);
g0[0.5>=g0 & g0>0]=1;
g0[g0>0.5& g0<1]=2 }

p3<-exp(b0+2*b.g)
if (g0==0 & p3<1){D=rbinom(1,1,exp(b0))}
if (g0==1 & p3<1){D=rbinom(1,1,exp(b0+b.g))}
if (g0==2 & p3<1){D=rbinom(1,1,exp(b0+2*b.g))}

if (D==1 & p3<1){

full.D[m1,]=cbind(D,g0,gf,gm);m1=m1+1 

}
 }

m2=1

while (m2<=n3){
### stratification ####
u1=runif(1,0,1)
if(u1>q){p.Af=0.3;p.Am=0.3;b0=log(0.01)}
if(u1<=q){p.Af=0.3+tk;p.Am=0.3+tk;b0=log(0.05)}

p.af=1-p.Af
{
gf=runif(1,0,1);
gf[gf>0 & gf<=p.Af^2]=2;
gf[p.Af^2<gf&gf<=(p.Af^2+2*p.af*p.Af)]=1;
gf[gf>(p.Af^2+2*p.af*p.Af)& gf<1]=0
}
p.am=1-p.Am
{
gm=runif(1,0,1);
gm[gm>0 & gm<=p.Am^2]=2;
gm[p.Am^2<gm & gm<=(p.Am^2+2*p.am*p.Am)]=1;
gm[gm>(p.Am^2+2*p.am*p.Am)& gm<1]=0
}
if (gf==0 & gm==0) g0=0
if ((gf==0 & gm==1)|(gf==1 & gm==0)) g0=rbinom(1,1,0.5)
if ((gf==0 & gm==2)|(gf==2 & gm==0)) g0=1
if (gf==1 & gm==1) {g0=runif(1,0,1);
g0[0.25>=g0 & g0>0]=0 ;
g0[g0>0.25& g0<=0.75]=1; 
g0[g0>0.75& g0<1]=2}
if (gf==2 & gm==2) g0=2
if ((gf==1 & gm==2)|(gf==2 & gm==1)) {g0=runif(1,0,1);
g0[0.5>=g0 & g0>0]=1;
g0[g0>0.5& g0<1]=2 }

p3<-exp(b0+2*b.g)
if (g0==0 & p3<1){D=rbinom(1,1,exp(b0))}
if (g0==1 & p3<1){D=rbinom(1,1,exp(b0+b.g))}
if (g0==2 & p3<1){D=rbinom(1,1,exp(b0+2*b.g))}

if (D==1 & p3<1){

full2.D[m2,]=cbind(D,g0,gf,gm);m2=m2+1 

}
 }

m3=1

while (m3<=(n3+n)){
### stratification ####
u1=runif(1,0,1)
if(u1>(1-q)){p.Af=0.3;p.Am=0.3;b0=log(0.01)}
if(u1<=(1-q)){p.Af=0.3+tk;p.Am=0.3+tk;b0=log(0.05)}
p.af=1-p.Af
{
gf=runif(1,0,1);
gf[gf>0 & gf<=p.Af^2]=2;
gf[p.Af^2<gf&gf<=(p.Af^2+2*p.af*p.Af)]=1;
gf[gf>(p.Af^2+2*p.af*p.Af)& gf<1]=0
}
p.am=1-p.Am
{
gm=runif(1,0,1);
gm[gm>0 & gm<=p.Am^2]=2;
gm[p.Am^2<gm & gm<=(p.Am^2+2*p.am*p.Am)]=1;
gm[gm>(p.Am^2+2*p.am*p.Am)& gm<1]=0
}
if (gf==0 & gm==0) g0=0
if ((gf==0 & gm==1)|(gf==1 & gm==0)) g0=rbinom(1,1,0.5)
if ((gf==0 & gm==2)|(gf==2 & gm==0)) g0=1
if (gf==1 & gm==1) {g0=runif(1,0,1);
g0[0.25>=g0 & g0>0]=0 ;
g0[g0>0.25& g0<=0.75]=1; 
g0[g0>0.75& g0<1]=2}
if (gf==2 & gm==2) g0=2
if ((gf==1 & gm==2)|(gf==2 & gm==1)) {g0=runif(1,0,1);
g0[0.5>=g0 & g0>0]=1;
g0[g0>0.5& g0<1]=2 }

p3<-exp(b0+2*b.g)
if (g0==0 & p3<1){D=rbinom(1,1,exp(b0))}
if (g0==1 & p3<1){D=rbinom(1,1,exp(b0+b.g))}
if (g0==2 & p3<1){D=rbinom(1,1,exp(b0+2*b.g))}

if (D==0 & p3<1){


full3.c[m3,]=cbind(D,g0,gf,gm);m3=m3+1 

}
 }

sub.D=full.D

SS.g=rep(0,n)


################### (CPG) family data likelihood and score ########

LL1= function(p) {
L<-rep(0,n)

 for(i in 1:n){
 lh<-0
 RR1=exp(p[1])
 RR2=exp(2*p[1])

 if((sub.D[i,3]==1 & sub.D[i,4]==2 &sub.D[i,2]==2 )|(sub.D[i,3]==2 & sub.D[i,4]==1 &sub.D[i,2]==2)) {lh=log(RR2/(RR1+RR2))}
 if((sub.D[i,3]==1 & sub.D[i,4]==2 &sub.D[i,2]==1 )|(sub.D[i,3]==2 & sub.D[i,4]==1 &sub.D[i,2]==1)){lh=log(RR1/(RR1+RR2))}
 if(sub.D[i,3]==1 & sub.D[i,4]==1 &sub.D[i,2]==2 ){lh=log(RR2/(1+2*RR1+RR2))}
 if(sub.D[i,3]==1 & sub.D[i,4]==1 &sub.D[i,2]==1 ) {lh=log(2*RR1/(1+2*RR1+RR2))}
 if(sub.D[i,3]==1 & sub.D[i,4]==1 &sub.D[i,2]==0 ) {lh=log(1/(1+2*RR1+RR2))}
 if((sub.D[i,3]==1 & sub.D[i,4]==0 &sub.D[i,2]==1 )|(sub.D[i,3]==0 & sub.D[i,4]==1 &sub.D[i,2]==1)) {lh=log(RR1/(1+RR1))}
 if((sub.D[i,3]==1 & sub.D[i,4]==0 &sub.D[i,2]==0 )|(sub.D[i,3]==0 & sub.D[i,4]==1 &sub.D[i,2]==0)){lh=log(1/(1+RR1))}
 L[i]<-lh
 }
  
-(sum(L))
}


out.nf <- nlm(LL1,p=c(0), hessian = TRUE,steptol = 1e-4, iterlim = 1000)


b.nf=out.nf$estimate

ck1=as.numeric(is.nan(det(out.nf$hessian)))
ck2=det(out.nf$hessian)
if (ck1==0 & ck2!=0 & ck2!=-Inf & ck2!=Inf){

bnf.v=ginv(out.nf$hessian)

for(i in 1:n){

phi1=exp(b.nf[1])
phi2=exp(2*b.nf[1])

if (sub.D[i,3]==2 & sub.D[i,4]==2 ){S.g=0}
if ((sub.D[i,3]==2 & sub.D[i,4]==1 & sub.D[i,2]==2)|(sub.D[i,3]==1 & sub.D[i,4]==2 & sub.D[i,2]==2)) {S.g=1/(1+phi1)}
if ((sub.D[i,3]==2 & sub.D[i,4]==1 & sub.D[i,2]==1)|(sub.D[i,3]==1 & sub.D[i,4]==2 & sub.D[i,2]==1)) {S.g=-phi1/(1+phi1)}
if ((sub.D[i,3]==2 & sub.D[i,4]==0 )|(sub.D[i,3]==0 & sub.D[i,4]==2 )) {S.g=0}
if (sub.D[i,3]==1 & sub.D[i,4]==1 &sub.D[i,2]==2 ){S.g=2/(1+phi1) }
if (sub.D[i,3]==1 & sub.D[i,4]==1 &sub.D[i,2]==1 ) {S.g=-(phi1-1)/(1+phi1) }
if (sub.D[i,3]==1 & sub.D[i,4]==1 &sub.D[i,2]==0 ) {S.g=-2*phi1/(1+phi1) }
if ((sub.D[i,3]==1 & sub.D[i,4]==0 &sub.D[i,2]==1 )|(sub.D[i,3]==0 & sub.D[i,4]==1 &sub.D[i,2]==1)) {S.g=1/(1+phi1)}
if ((sub.D[i,3]==1 & sub.D[i,4]==0 &sub.D[i,2]==0 )|(sub.D[i,3]==0 & sub.D[i,4]==1 &sub.D[i,2]==0)){S.g=-phi1/(1+phi1) }
if ((sub.D[i,3]==0 & sub.D[i,4]==0 )) {S.g=0}
SS.g[i]=S.g

}

 
#################### logistic score and hessian matrix##############

subGr=full.D[,2]
subGr.c=full3.c[,2]

sub.stat=c(rep(1,n),rep(0,(n3+n)))

G.cc=c(subGr,subGr.c)
G.rb=c(full.D[,2],full2.D[,2],full3.c[,2])
status=c(rep(1,(n+n3)),rep(0,(n+n3)))
fix.rb=glm(status~G.rb,family=binomial)
rr.bar=summary(fix.rb)
r.bar=as.matrix(fix.rb$coefficients)
rco.v=as.matrix(rr.bar$coefficients[,2])^2
gamma1[i1]=r.bar[2]
gamma.var[i1]=rco.v[2]
pr.bar=exp(r.bar[1]+r.bar[2]*G.rb)/(1+exp(r.bar[1]+r.bar[2]*G.rb))

####
fix.r=glm(sub.stat~G.cc,family=binomial)
rr.hat=summary(fix.r)
d2.2=solve(rr.hat$cov.unscaled)
r.hat=as.matrix(fix.r$coefficients)
pr.hat=exp(r.hat[1]+r.hat[2]*G.cc)/(1+exp(r.hat[1]+r.hat[2]*G.cc))
z=cbind(rep(1,(n+n+n3)),G.cc)
d2=t(z)%*%(pr.hat*(1-pr.hat)*z)
c22=t(z)%*%((sub.stat-pr.hat)^2*z)

####### adjust logistic regression Score####
SS.hat=((sub.stat-pr.hat)*z)
I.hat=t(z)%*%(pr.hat*(1-pr.hat)*z)
SShat.adj=SS.hat[,2]-I.hat[1,2]*I.hat[1,1]^(-1)*SS.hat[,1]
I.22=I.hat[2,2]-I.hat[1,2]*I.hat[1,1]^(-1)*I.hat[1,2]

z2=cbind(rep(1,(n+n3)*2),G.rb)
SS.bar=((status-pr.bar)*z2)
I.bar=t(z2)%*%(pr.bar*(1-pr.bar)*z2)
SSbar.adj=SS.bar[,2]-I.bar[1,2]*I.bar[1,1]^(-1)*SS.bar[,1]
I.33=I.bar[2,2]-I.bar[1,2]*I.bar[1,1]^(-1)*I.bar[1,2]

########## Two-stage analysis ######

d1.iv=ginv(out.nf$hessian)[1,1]
c12=d1.iv*(SS.g)%*%(SShat.adj[1:n])*solve(I.22)
c13=d1.iv*(SS.g)%*%(SSbar.adj[1:n])*solve(I.33)
c11=d1.iv*t(SS.g)%*%(SS.g)*d1.iv
cov22=solve(I.22)*SShat.adj%*%SShat.adj*solve(I.22)
cov33=solve(I.33)*SSbar.adj%*%SSbar.adj*solve(I.33)
cov23=solve(I.22)*(c(SShat.adj[1:n],rep(0,n3),SShat.adj[(n+1):(2*n+n3)]))%*%SSbar.adj*solve(I.33)

bg1=b.nf[1]-(c12-c13)*ginv(cov22+cov33-2*cov23)*(r.hat-r.bar)[2]
Q=(c12-c13)*ginv(cov22+cov33-2*cov23)
varb1.bar=d1.iv+Q*cov22*Q+Q*cov33*Q-2*c12*Q+2*c13*Q-2*Q*cov23*Q
SS.TT[i1]=bg1^2/varb1.bar

c11.trio=t(cbind(SS.g))%*%cbind(SS.g)

Rob.vnf=d1.iv
beta.hat[i1]=b.nf[1]
beta.bar[i1]=bg1
varb.hat[i1]=Rob.vnf
varb.bar[i1]=varb1.bar

###  Wald statistic #########
Case.T[i1]=(r.bar[2])^2/rco.v[2]
bgstat.h[i1]=(b.nf[1])^2/Rob.vnf
bgstat.b[i1]=(bg1)^2/varb1.bar


}



Sys.time()-start

###  Wald Test #########
t.c=Case.T[Case.T>qchisq(0.95,1)]
s.t=SS.TT[SS.TT>=qchisq(0.95,1)]
ct1=bgstat.h[bgstat.h>qchisq(0.95,1)]
ct2=bgstat.b[bgstat.b>qchisq(0.95,1)]
tt.c=length(t.c)/k
score.t=length(s.t)/k
hbg.t=length(ct1)/k
bbg.t=length(ct2)/k

### output #############
report1=matrix(c("","beta_0","beta_G","true",b0,b.g,
"beta.hat","",mean(beta.hat),
"beta.bar","",mean(beta.bar),
"Case_control_gamma","",mean(gamma1),
"Var.beta.hat","",mean(varb.hat),
"Var.beta.bar","",mean(varb.bar),
"Var.Case_control_gamma","",mean(gamma.var),
#"Simu var.hat","",var(beta.hat),
#"Simu var.bar","",var(beta.bar),
#"Simu var.CC","",var(gamma1),
"test hat","beta.hat=0",hbg.t,
"test bar","beta.bar=0",bbg.t,
"Case-control ","Test Gamma=0",tt.c),3,11)
print(t(report1))