setwd("C:/Users/tsuesse/Dropbox/MERGM/Marginal")
rm(list=ls())
source("fcts.Lazega.r")
source("marginal.ERGM.fct.r")
dim(Y.MCMC)
beta0<-beta.new
Psi.new
gwesp.alpha.new
gwesp.tau.new
load("SolvePsi_22th_n36_M30000.RData")    # MCMC sample from MERGM Suesse (2012)
dim(Y.MCMC)
beta0<-beta.new
Psi.new
gwesp.alpha.new
gwesp.tau.new
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)
R0<-  Diagonal(N,1/sqrt(hilf$v.Y))%*% hilf$cov.Y%*%Diagonal(N,1/sqrt(hilf$v.Y))
p<-dim(X)[2]
p
dim(R0)
hilf
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)
}
results<-cbind(beta.hat[,2,],se.beta.hat[,2,],se.rob.beta.hat[,2,])
ind<-c(1,7,13)
results<-results[,c(ind,ind+1,ind+2,ind+3,ind+4,ind+5)]
output(round(results,3),nsmall=3)
results<-cbind(beta.hat[,2,],se.beta.hat[,2,],se.rob.beta.hat[,2,])
ind<-c(1,7,13)
results<-results[,c(ind,ind+1,ind+2,ind+3,ind+4,ind+5)]
output(round(results,4),nsmall=3)
dim(beta)
dim(beta.hat)
beta.hat[,1,]
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)]
output(round(results,4),nsmall=3)
output(round(results,4),nsmall=5)
results
round(results,4)
help(output)
output
output(round(results,4),nsmall=5,digits=3)
results
rm(list=ls())
source("fcts.Lazega.r")
source("marginal.ERGM.fct.r")
scenario<-1
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)
sim <- dim(Y.MCMC)[1]
hilf<-mean.var.w(Y.MCMC,w=rep(1/sim,sim),binary=TRUE)
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]
#        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
SessionInfo()
sessionInfo()