# This R file does the simulations for the results presented 
# Tables 3-8

# set wordking directly to the folder including all files
# setwd("~/R/ERGM/Marginal")
# cd ~/R/ERGM/Marginal
rm(list=ls())


source("fcts.Lazega.r")
source("marginal.ERGM.fct.r")



# source("simu_marg.r")
# Rscript simu_marg.r
# load("Beta_Theta_Psi.RData")

# set the scenario (1-9), see Table 2 of the manuscript

scenario<-1


# 1 n=36 with data set X and original beta - GWESP
# 2 n=36 with data set X  beta=(-7,0,1,1,1,2) - GWESP
# 3 n=36 with data set X  beta=(-7,0,1,1,1,2) - degree(2)
# 4 n=36 with data set X  beta=(-7,0,1,1,1,2) - triangle
# 5 n=100 with new X  beta=(-6,0.02,0,1,0) - GWESP
# 6 n=100 with new X  beta=(-6,0.02,0,1,0) - triangle
# 7 n=100 with new X  beta=(-6,0.02,0,1,0) - degree(2)
# 8 n=36 with data set X  beta=(-7,0,1,1,1,2) - marg: sigma=1
# 9 n=100 with data set X   beta=(-6,0.02,0,1,0) - marg: sigma=1


# M is number of MCMC samples

switch(scenario,
       
{ # scenario = 1
  
  # n=36, GWESP  , original data set
# with 30,000 MCMC for ERGM in 2012 paper
# based on fitted MERGM to data set
# beta.old =  -7.372 0.039 0.764 0.883 0.78 1.75
load("SolvePsi_22th_n36_M30000.RData")    # OK
dim(Y.MCMC)
beta0<-beta.new
Psi.new
gwesp.alpha.new 
gwesp.tau.new   

n<-36
link<-"logit"
}
,
{  # scenario = 2
# n=36, GWESP  
# based on 
#gwesp.tau.old<-1
#gwesp.alpha.old<- 0.5  # rate parameter
  # beta.old =   -7    0    1    1    1    2
load("Psi_fixed_beta_GWESP_n36_M30000.RData")  # OK
dim(Y.MCMC)
Psi.new
# Psi.old don't use
gwesp.alpha.old   #0.5
gwesp.tau.old     # 1
n<-36
link<-"logit"
}
,
{  # scenario = 3
  #n=36 , degree(2)
  # beta.old =   -7    0    1    1    1    2
  load("Psi_fixed_beta_degree(2)_n36_M30000.RData")    # OK
  n<-36 
  link<-"logit"
}
,
{  # scenario = 4
  #n=36 , triangle
  # beta.old =   -6    0    0.5    1    1    1
  load("Psi_fixed_beta_triangle_n36_M30000.RData")
  n<-36
  link<-"logit"
}
,
{  # scenario = 5
  #n=100 , GWESP
  load("Psi_fixed_beta_GWESP_n100_M30000.RData")
  n<-100
  # beta.old =  -6 0.02    0    1    0
  link<-"logit"
}
,
{  # scenario = 6
  #n=100 , triangle
  load("Psi_fixed_beta_triangle_n100_M30000.RData")
  n<-100
  # beta.old =  -6 0.02    0    1    0
  link<-"logit"
}
,
{  # scenario = 7
  #n=100 , degree(2)
  load("Psi_fixed_beta_degree(2)_n100_M30000.RData")
  n<-100
  # beta.old =  -6 0.02    0    1    0
  link<-"logit"
}
,
{  # scenario = 8
  #n=36 , marg-model
  load("Probit_n36_M30000.RData")
  n<-36
   # beta.old =   -7    0    1    1    1    2
  link<-"probit"
}
,
{  # scenario = 8
  #n=100 , marg-model
  load("Probit_n100_M30000.RData")
  n<-100
  # beta.old =   -7    0    1    1    1    2
  link<-"probit"
}
)
beta0 <- beta.old




N<-n*(n-1)/2
A<-get.A(n)   # get Matrix A as discussed in Section 7.1 
sim <- dim(Y.MCMC)[1]


hilf<-mean.var.w(Y.MCMC,w=rep(1/sim,sim),binary=TRUE)

# extract "TRUE" correlation structure from MCMC data sets (MERGM)
R0<-  Diagonal(N,1/sqrt(hilf$v.Y))%*% hilf$cov.Y%*%Diagonal(N,1/sqrt(hilf$v.Y))


if(n==100){
  Y.MCMC<-Y.MCMC[seq(1,sim,by=3),]  # reduce from M=30,000 to M=10,000 for n=100, otherwise takes too long
  sim <- dim(Y.MCMC)[1]
  sim<-2e3
}#end 

# just make it smaller for testing purposes
# sim<-100

# R<- Diagonal(N,1/sqrt(hilf$v.Y))%*%COV%*%Diagonal(N,1/sqrt(hilf$v.Y))
# diag(R)<-0; mean(R)


p<-dim(X)[2]

# Fitting Methods as in paper
#        k   ==   0      1       2      3      4
method.names<-c("ind","sparse","a1a2","exch","TRUE")

se.rob.beta.hat<-se.beta.hat<-beta.hat<-array(NA,dim=c(sim,5,p),dimnames=list(NULL,method.names,1:p))

#max 2 parameters
alpha.hat<- array(NA,dim=c(sim,5,2),dimnames=list(NULL,method.names,1:2))


cat("Start Simulations")




for(i in 1:sim){

cat("Sim=",i,"\n")  
  
 for(k in 0:4){
   
   
   
   st<-system.time({
     
     hilf <- try(GEE(beta0=NULL,Y=Y.MCMC[i,],X,corr=k,A=A,R0=R0,link=link,iter.max=20))
   })
     
   cat("Time needed",st,"\n") 
   cat("Method=",method.names[k+1],"\n")
   
   
   
   if(!inherits(hilf, "try-error") && !is.null(hilf) && max(abs(hilf$beta.hat))<1e2 ){
     
     cat("Max abs beta",max(abs(hilf$beta.hat)),"\n") 
     
     beta.hat[i,k+1,] <-  as.numeric(hilf$beta.hat)
     se.beta.hat[i,k+1,] <- hilf$se.beta.hat
     se.rob.beta.hat[i,k+1,] <- hilf$se.robust.beta.hat
     
     if(!is.null(hilf$alpha1)){
     alpha.hat[i,k+1,1] <-   hilf$alpha1
     }
     if(!is.null(hilf$alpha2)){
       alpha.hat[i,k+1,1:2] <-   hilf$alpha2
     }
     if(!is.null(hilf$alpha)){
       alpha.hat[i,k+1,1] <-   hilf$alpha
     }
   
   # hilf$R[30:40,30:40]
     
   }else{
     
     cat("Not converged or some beta too large\n")
   
    }#end if(!inherits(hilf, "try-error")){
     
   }#end for k

  if(i%%10==0){
  cat("##########################################################################\n")  
  cat("Iteration",i,"\n")
  cat("##########################################################################\n")  
  Sys.sleep(3)
  }



  
  
}#end for sim

#######################################################################
# Collecting all information from the simulation and 
# summarise results

beta0.all <- array(beta0,dim=c(p,sim,5)) 
beta0.all<-aperm(beta0.all,c(2,3,1))

# mse or RMSE
RMSE<-sqrt(apply((beta.hat-beta0.all)^2,2:3,mean,na.rm=TRUE))



# bias
BIAS<-apply((beta.hat-beta0.all),2:3,mean,na.rm=TRUE)


# adjust robust variance estimators
sim<-dim(se.rob.beta.hat)[1]

se.rob.beta.hat.new<-se.rob.beta.hat
for(i in 1:sim){
  for(j in 1:5){
    if(is.na(se.rob.beta.hat.new[i,j,])>0){
       se.rob.beta.hat.new[i,j,]<-NA
    }
  }}


# coverage
z<-qnorm(0.975)
cov95.naive <-apply(beta.hat - z*se.beta.hat <= beta0.all & beta.hat + z*se.beta.hat >= beta0.all,2:3,mean,na.rm=TRUE)
cov95.robust<-apply(beta.hat - z*se.rob.beta.hat <= beta0.all & beta.hat + z*se.rob.beta.hat >= beta0.all,2:3,mean,na.rm=TRUE) 

length95.naive  <- apply(z*se.beta.hat,2:3,mean,na.rm=TRUE)
length95.robust <- apply(z*se.rob.beta.hat,2:3,mean,na.rm=TRUE)

z<-qnorm(0.95)
cov90.naive <-apply(beta.hat - z*se.beta.hat <= beta0.all & beta.hat + z*se.beta.hat >= beta0.all,2:3,mean,na.rm=TRUE)
cov90.robust<-apply(beta.hat - z*se.rob.beta.hat <= beta0.all & beta.hat + z*se.rob.beta.hat >= beta0.all,2:3,mean,na.rm=TRUE) 

length90.naive  <- apply(z*se.beta.hat,2:3,mean,na.rm=TRUE)
length90.robust <- apply(z*se.rob.beta.hat,2:3,mean,na.rm=TRUE)

z<-qnorm(0.995)
cov99.naive <-apply(beta.hat - z*se.beta.hat <= beta0.all & beta.hat + z*se.beta.hat >= beta0.all,2:3,mean,na.rm=TRUE)
cov99.robust<-apply(beta.hat - z*se.rob.beta.hat <= beta0.all & beta.hat + z*se.rob.beta.hat >= beta0.all,2:3,mean,na.rm=TRUE) 

length99.naive  <- apply(z*se.beta.hat,2:3,mean,na.rm=TRUE)
length99.robust <- apply(z*se.rob.beta.hat,2:3,mean,na.rm=TRUE)

# save all summary data, as CI coverage, CI length,  RMSE and BIAS

save(alpha.hat,beta0,beta0.all,R0,A,RMSE,BIAS,se.rob.beta.hat,se.beta.hat,
     beta.hat,X,
     cov95.robust,cov95.naive,length95.naive,length95.robust,
     cov90.robust,cov90.naive,length90.naive,length90.robust,
     cov99.robust,cov99.naive,length99.naive,length99.robust,
     file=paste("Results_n",n,"_sim",sim,"scenario",scenario,".RData",sep=""))

# break
# break
# stop

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

# some testing, don't run

#################################################################################################
if(0){
  
  # 1 n=36 with data set X and original beta - GWESP
  # 2 n=36 with data set X  beta=(-7,0,1,1,1,2) - GWESP
  # 3 n=36 with data set X  beta=(-7,0,1,1,1,2) - degree(2)
  # 4 n=36 with data set X  beta=(-7,0,1,1,1,2) - triangle
  # 5 n=100 with new X  beta=(-6,0.02,0,1,0) - GWESP
  # 6 n=100 with new X  beta=(-6,0.02,0,1,0) - triangle
  # 7 n=100 with new X  beta=(-6,0.02,0,1,0) - degree(2)
  
  source("marginal.ERGM.fct.r")
  
  n.index<-c(36,36,36,36,100,100,100,36,100)
  sim.index<-c(rep(3e4,4),rep(1e4,3),3e4,1e3)
  
  BIAS.all <- NULL
  RMSE.all <- NULL
  cov95.naive.all <- NULL
  cov95.robust.all <- NULL
  length95.naive.all <- NULL
  length95.robust.all <- NULL
  
  for(i in 1:9){
  
  file=paste("Results_n",n.index[i],"_sim",sim.index[i],"scenario",i,".RData",sep="")
  load(file)
  
  beta0<-c(beta0)
  p<-length(beta0)
  #beta0.all <- array(beta0,dim=c(p,sim,5)) 
  #beta0.all<-aperm(beta0.all,c(2,3,1))
  
  # adjust robust variance estimators
  sim<-dim(se.rob.beta.hat)[1]
  
  se.rob.beta.hat1<-se.rob.beta.hat
  se.beta.hat1<-se.beta.hat
  for(k in 1:sim){
    for(j in 1:5){
      if(max(is.na(se.rob.beta.hat1[k,j,]))){
        se.rob.beta.hat1[i,j,]<-NA
      }
      if(max(is.na(se.beta.hat1[k,j,]))){
        se.beta.hat1[i,j,]<-NA
      }
      
    }}
  z<-qnorm(0.975)
 
  # apply((beta.hat[,,]),c(2:3),mean,na.rm=T)
  
  if(i==9 | i==8){
    beta0<-beta0/sqrt(2);
    length95.naive<-length95.naive/sqrt(2);
    length95.robust<-length95.robust/sqrt(2);
    }
  
  beta0.all <- array(beta0,dim=c(p,sim,5)) 
  beta0.all<-aperm(beta0.all,c(2,3,1))
  
  # mse or RMSE
  RMSE<-sqrt(apply((beta.hat-beta0.all)^2,2:3,mean,na.rm=TRUE))
  
  
  
  # bias
  BIAS<-apply((beta.hat-beta0.all),2:3,mean,na.rm=TRUE)
  
  
  cov95.robust<-apply(beta.hat - z*se.rob.beta.hat1 <= beta0.all & beta.hat + z*se.rob.beta.hat1 >= beta0.all,2:3,mean,na.rm=TRUE) 
  cov95.naive<-apply(beta.hat - z*se.beta.hat <= beta0.all & beta.hat + z*se.beta.hat >= beta0.all,2:3,mean,na.rm=TRUE) 
  
 
 
 
 
 nconv.beta<-apply(is.na(beta.hat[,,1]),c(2),mean)
 conv.beta<-1-nconv.beta
 
 nconv.naive<-apply(is.na(se.beta.hat1),c(2),mean)
 nconv.robust<-apply(is.na(se.rob.beta.hat1),c(2),mean)
 conv.naive<-1-nconv.naive
 conv.robust<-1-nconv.robust
 #conv2<-apply(is.na(se.beta.hat),c(2),mean) 
  
 
 
  if(n.index[i]==100){
    BIAS<- cbind(BIAS,0)
    RMSE<- cbind(RMSE,0)
    cov95.naive<- cbind(cov95.naive,0)
    cov95.robust<- cbind(cov95.robust,0)
    length95.naive<- cbind(length95.naive,0)
    length95.robust<- cbind(length95.robust,0)
    beta0<-c(beta0,0)
    }
  
 
  BIAS<-cbind(round(100*nconv.beta,1),100*BIAS)
  RMSE<- cbind(round(100*nconv.beta,1),100*RMSE)
  cov95.naive<- cbind(round(100*(conv.beta-conv.naive)/conv.beta,1),100*cov95.naive)
  cov95.robust<- cbind(round(100*(conv.beta-conv.robust)/conv.beta,1),100*cov95.robust)
  length95.naive<- cbind(round(100*(conv.beta-conv.naive)/conv.beta,1),length95.naive)
  length95.robust<- cbind(round(100*(conv.beta-conv.robust)/conv.beta,1),length95.robust)
   
   
  method.names<-c("& ind","& sparse",paste(i," & a1a2",sep=""),"& exch","& TRUE")
  rownames(BIAS)<- rownames(RMSE)<- rownames(cov95.naive)<- rownames(cov95.robust)<- rownames(length95.naive)<- rownames(length95.robust)<-  method.names

  
  
  
  
  BIAS.all <- rbind(BIAS.all,BIAS)
  RMSE.all <- rbind(RMSE.all,RMSE)
  cov95.naive.all <- rbind(cov95.naive.all,cov95.naive)
  cov95.robust.all<- rbind(cov95.robust.all,cov95.robust)
  length95.naive.all<- rbind(length95.naive.all,length95.naive)
  length95.robust.all<- rbind( length95.robust.all,length95.robust)
  
  
  
  }#end  for(i in 1:7){
  
  
  BIAS.all<-round(BIAS.all,3)
  RMSE.all<-round(RMSE.all,3)
  cov95.naive.all<-round(cov95.naive.all,3)
  cov95.robust.all<-round(cov95.robust.all,3)
  length95.naive.all<-round(length95.naive.all,3)
  length95.robust.all<-round(length95.robust.all,3)
  
  
  source("marginal.ERGM.fct.r")
  output(BIAS.all,nsmall=3,special="\\hline",every=5)
  output(RMSE.all,nsmall=3,special="\\hline",every=5)
  output(cov95.naive.all,nsmall=1,special="\\hline",every=5)
  output(cov95.robust.all,nsmall=1,special="\\hline",every=5)
  output(length95.naive.all,nsmall=3,special="\\hline",every=5)
  output(length95.robust.all,nsmall=3,special="\\hline",every=5)
  
  
  
}

if(0){


n<-36
N<-n*(n-1)/2
ind<-indexfct(n)

if(0){
  st1<-system.time({
    # get correlation matrix for typical undirected network
    
    #N<-dim(ind)[2]
    edges<-matrix(0,0,2)
    for(i in 1:(N-1)){
      #cat("i=",i,"\n")
      #cat("i=",i," j=",j,"\n")
      for(j in (i+1):N){
        
        hilf<- intersect(ind[,i],ind[,j] )
        
        #print(hilf)
        #print(length(hilf))
        #print(is.null(hilf))
        
        if( length(hilf) ){
          edges<- rbind(edges,c(i,j))
        }
        
      }#end
    }#end 
    
    
    # for full model use Z1,Z2,Z3
    Z1 <-    sparseMatrix(i=c(edges[,1]),j=c(edges[,2]),x=rep(1,dim(edges)[1]),dims=c(N,N))
    Z1 <-   Z1+t(Z1)
    #Z1 <-    sparseMatrix(i=c(edges[,1],edges[,2]),j=c(edges[,2],edges[,1]),x=rep(1,2*dim(edges)[1]))
  })
}


st2<-system.time({
  A<-(outer(t(ind),(ind),"abs.diff"))  # calculates differences of pairs, if zero index must match
  #A[1,,,69]
  A<- 1*(apply(A,c(1,4),min)==0)
  #A[1,69]
  diag(A)<-0
  
  A<-as(A,"dgCMatrix")
})


R<-diag(N)+0.1*A
#image(R)
#image(solve(R))


sum(R==0)/2
sum(R==1)
sum(R==0.1)/2
N*(N-1)/2
N


#all.equal(Z1,A)



#play around with exchangeable structure

#

n<-10
B<- matrix(1,n,n);diag(B)<-0;
R<-Diagonal(n) + 0.1*B

solve(R)

#


Z2<-Diagonal(N)

Z3 <- matrix(1,N,N)
Z3<- Z3 -Z1 - Z2


# for sparse model use just Z1 and Z2 (ignore Z3 which will be usually small)


# exchangeable structure
# define 2 Matrices
# 
M1 <- Diagonal(N)
M2 <- matrix(1,N,N)
diag(M2)<-0


}#end tests



#################################################################
#   PLOTS
#################################################################

if(0){
# plots


ind.lt<-lower.tri(R11)
summary(Corr1[ind.lt])
quantile(Corr1[ind.lt],probs=c(0.8,0.9,0.95,0.98,0.99,0.999))




Corr1<-round(R11,2)

breaks.ini<-as.numeric(names(table(c(Corr1))))
#breaks<-c(breaks[breaks<0],-1e-4,0,1e-4,breaks[breaks>0])
#breaks<-c(1e-4,breaks)-1e-5
if(breaks.ini[1]==0){
  breaks<-c(0,breaks.ini[-1]-1e-5,1)
  breaks1<-c(0,breaks.ini[-1],1)
}else{
  breaks<- c(0,breaks.ini-1e-5,1)
  breaks1<- c(0,breaks.ini,1)
}

dyads<-N

l.breaks<-length(breaks)

cols<-grey.colors(l.breaks,start=0,end=1)
cols<-cols[-(l.breaks-1)]

#jpeg(filename = paste("Corr",k,"_dir.jpeg",sep=""),width = 3000, height =3000,  res = 300)
nf<-layout(matrix(c(1,2),1,2,byrow=T),widths =c(6,1) ,heights=c(4))
par(mar=c(3,3,3,2))
#par(mar=c(2,2,2,1))
image(x=1:dyads,y=1:dyads,Corr1,lwd=2,col=cols,breaks=breaks,axes=FALSE,xlab="dyads",ylab="dyads",main="Lawyers",cex.main=2)
axis(1,at=1:dyads,labels=names)
axis(2,at=1:dyads,labels=names)




par(mar=c(0,0,0,0))
plot.new()
image.plot(zlim=range(0,Corr1),legend.only=TRUE,col=cols,smallplot=c(.1,.4, .1,.9),
           cex.main=0.8,cex.axis=0.8,cex.lab=0.8,cex.sub=0.8,breaks=1:l.breaks,
           #axis.args=list(at=c(0.02,0.1,1), labels=c("0.02","0.1","1")) ) 
           axis.args=list(at=seq(1.5,l.breaks,by=1), labels=breaks1[-(l.breaks-1)] ) ) 
#dev.off()
}#end if 0/1