# This files fits the methods in the Paper to the Lazega Data set, see Table 9 of the manuscript


# set working directly to folder with all files
#setwd("~/R/ERGM/Marginal")
rm(list=ls())

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

# load MCMC sample from MERGM fitted to Lazega data set published in Suesse (2012)
load("SolvePsi_22th_n36_M30000.RData")   
dim(Y.MCMC)
beta0<-beta.new
Psi.new
gwesp.alpha.new 
gwesp.tau.new   

n<-36



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)

#correlation matrix of all 630 dyads
R0<-  Diagonal(N,1/sqrt(hilf$v.Y))%*% hilf$cov.Y%*%Diagonal(N,1/sqrt(hilf$v.Y))


p<-dim(X)[2] # number fixed effects

# 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(5,2,p),dimnames=list(method.names,c("logit","probit"),1:p))

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


cat("Start fitting")

link<-c("logit","probit")

for(l in 1:2){

  for(k in 0:4){
    
    
    
    st<-system.time({
      
      hilf <- try(GEE(beta0=NULL,Y=Y.obs,X,corr=k,A=A,R0=R0,link=link[l]))
    })
    
    cat("Time needed",st,"\n") 
    cat("Method=",method.names[k+1],"\n")
    
    if(!inherits(hilf, "try-error") && !is.null(hilf)){
      beta.hat[k+1,l,] <-  as.numeric(hilf$beta.hat)
      se.beta.hat[k+1,l,] <- hilf$se.beta.hat
      se.rob.beta.hat[k+1,l,] <- hilf$se.robust.beta.hat
      
      if(!is.null(hilf$alpha1)){
        alpha.hat[k+1,l,1] <-   hilf$alpha1
      }
      if(!is.null(hilf$alpha2)){
        alpha.hat[k+1,l,1:2] <-   hilf$alpha2
      }
      if(!is.null(hilf$alpha)){
        alpha.hat[k+1,l,1] <-   hilf$alpha
      }
    
      
    }#end if(!inherits(hilf, "try-error")){
    
  }#end for k
  
  }#end for(l in 1:2){
  
  
  if(i%%10==0){
    cat("##########################################################################\n")  
    cat("Iteration",i,"\n")
    cat("##########################################################################\n")  
    Sys.sleep(3)
  }
  

#set l<-1  for logit and l<-2 for probit

l<-1

results<-cbind(beta.hat[,l,],se.beta.hat[,l,],se.rob.beta.hat[,l,])
ind<-c(1,7,13)
results<-results[,c(ind,ind+1,ind+2,ind+3,ind+4,ind+5)]

results

#output(round(results,4),nsmall=5,digits=3)  



if(0){ # DO NOT RUN
  
  # this R code produces the correlation structures in Figure 3
  # for the Lazega data set
  
  library(graphics)
  #library(scatterplot3d)
  #library(lattice)
  library(fields)
  # PLOTS CORRELATION STRUCTURE ERGM
  
  full<-T
  
  
  
  if(full){
    title<-"Lazega - Full Correlation Matrix"
    names<-NULL
    use.Psi<-T
    
    Corr1<-as.matrix(R0)
    dyads<-630
    
    
  }else{
    title<-"Lazega - Partial Correlation Matrix"
    use.Psi<-F
    ind1<-indexfct(n)
    
    
    ind<-indexfct(6)
    
    
    dyads<-15  #as Corr1 is 15 times 15 matrix
    
    
    names<-rep(NA,dyads)
    for(i in 1:dyads){
      names[i]<-paste("(",ind[1,i],",",ind[2,i],")",sep="")
    }
    
    indices<-NULL
    for(i in 1:15){
    indices<-c(indices,(1:4950)[apply(t(ind1), 1, function(x) all(x == ind[,i]))]) 
    }
    
    
    Corr1<-as.matrix(R0[indices,indices])
    
    
    
  }#end f
  
  
  Corr1<-abs(round(Corr1,3))
  
  
  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)
  }
  
  
  
  
  
  l.breaks<-length(breaks)
  
  cols<-grey.colors(l.breaks,start=0,end=1)
  cols<-cols[-(l.breaks-1)]
  
  
  #(100-l.breaks+2):100
  jpeg(filename = paste("R0_Lazega_full",full,".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=1,col=cols,breaks=breaks,axes=FALSE,xlab="",ylab="",main=title,cex.main=2)
  
  if(!full){
  axis(1,at=1:dyads,labels=names)
  axis(2,at=1:dyads,labels=names)
  }else{
  mtext("Dyads", side=1, line=1, cex=1,col="black")
  mtext("Dyads", side=2, line=1, cex=1,col="black")
  }
  
  
  if(use.Psi){
    list.arg<-list( at=0.5+c(1,ceiling(l.breaks/2),l.breaks-10,l.breaks-1),
                    labels=breaks1[c(1,ceiling(l.breaks/2),l.breaks-2,l.breaks-1)] )
  }else{
    list.arg<-list( at=seq(1.5,l.breaks,by=1) , labels=breaks1[-(l.breaks-1)] )
  }
  
  
  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=, labels=breaks1[-(l.breaks-1)] ) 
             ,axis.args=list.arg 
  ) 
  
  
  dev.off()
  
}