library("Bergm")

data("lazega", package = "Bergm")

# Lazega (1)
par(mfrow = c(1, 2))
CC <- hcl.colors(3, "Teal")
#CC <- c("#DC778D", "#7595DD", "#00AD73")
set.seed(22)
plot(lazega,
     vertex.col = CC[lazega %v% "Office"], 
     vertex.cex = 2)
legend("topright",
       pch    = 21,
       pt.bg  = CC,
       legend = c("Boston", "Hartford", "Providence"),
       title  = "OFFICE")

# Lazega (2)
CC <- hcl.colors(3, "Oranges")
set.seed(22)
plot(lazega, 
     vertex.col = CC[lazega %v% "Practice"], 
     vertex.cex = 2)
legend("topright",
       pch    = 21,
       pt.bg  = CC,
       legend = c("Litigation", "Corporate"),
       title  = "PRACTICE")

# --- DIXON --- #

data(faux.dixon.high)
dixon <- faux.dixon.high
dixon %v% 'race' <- as.numeric(factor(dixon %v% 'race'))

# Dixon (1)
par(mfrow = c(1, 3))
CC <- hcl.colors(2, "Geyser")
#CC <- c("#5FB9F2", "#FA92A8")
set.seed(21)
plot(dixon,
     vertex.col = CC[dixon %v% 'sex'],
     vertex.cex = 2)
legend("topright",
       col    = CC,
       legend = c("Male", "Female"),
       pch    = 21,
       pt.bg  = CC,
       title  = "SEX")

# Dixon (2)
CC <- hcl.colors(5, "Purp")
#CC <- c("#0E3F5C", "#387D84", "#7BBCAB", "#D1FBD4")
set.seed(21)
plot(dixon,
     vertex.col = CC[dixon %v% 'race'],
     vertex.cex = 2)
legend("topright",
       legend = c("Black", "Hispanic", "Other", "White"),
       pch    = 21,
       pt.bg  = CC,
       title  = "RACE")

# Dixon (3)
CC <- hcl.colors(6, "Blues 3")
#CC <- c("#A93154", "#B84484", "#BD60A9", "#BA7EC5", "#B39BD9", "#AEB6E5")
set.seed(21)
plot(dixon,
     vertex.col = CC,
     vertex.cex = 2)
legend("topright",
       legend = paste('Grade', 7:12),
       pch    = 21,
       pt.bg  = CC,
       title  = "GRADE")

# Parameter inference using the bergm() function
# Model specification
m1 <- lazega ~ edges +
  nodematch("Office") +
  nodematch("Practice") +
  gwesp(0.5, fixed = TRUE)

# Prior specification
M.prior <- c(-4, 0.5, 0.5, 1)
S.prior <- diag(4, 4)

# Posterior estimation
set.seed(1)
p.m1 <- bergm(m1,
              prior.mean  = M.prior,
              prior.sigma = S.prior,
              burn.in     = 500,
              main.iters  = 3000,
              aux.iters   = 2500,
              nchains     = 8,
              gamma       = 0.6)
p.m1$specs <- c("edges", "nodematch.Office", "nodematch.Practice",
                "gwesp.fixed.0.5")

summary(p.m1)

plot(p.m1, lag = 100)

set.seed(1)
bgof(p.m1,
     aux.iters = 5000,
     n.deg     = 15,
     n.dist    = 9,
     n.esp     = 8)

## Parameter inference with missing data using the bergmM() function

set.seed(1)
missV <- sample(1:36, 4)
lazega[missV, ] <- lazega[, missV] <- NA
p.m1.M <- bergmM(m1,
                 prior.mean  = M.prior,
                 prior.sigma = S.prior,
                 burn.in     = 200,
                 main.iters  = 3000,
                 aux.iters   = 3000,
                 nchains     = 8,
                 gamma       = 0.6,
                 nImp        = 10,
                 seed        = 1)
p.m1.M$specs <- c("edges", "nodematch.Office", "nodematch.Practice",
                  "gwesp.fixed.0.5")

summary(p.m1.M)

#p.m1.M$impNets

## Evidence estimation using the evidence() function
data(faux.dixon.high)
dixon <- faux.dixon.high

m1 <- dixon ~ edges + mutual + absdiff("grade") + 
  nodefactor("race") + nodefactor("grade") + nodefactor("sex") +
  nodematch("race", diff = TRUE, levels = c("B", "O", "W")) + 
  nodematch("grade", diff = TRUE) + 
  nodematch("sex",  diff = FALSE) + 
  idegree(0:1) + odegree(0:1) + gwesp(0.1, fixed = TRUE)

M.prior1 <- c(-5, rep(0, 26))
S.prior1 <- diag(5, 27)

cj1 <- evidence(
  evidence.method = "CJ",
  formula         = m1,
  prior.mean      = M.prior1,
  prior.sigma     = S.prior1,
  aux.iters       = 2500,
  n.aux.draws     = 50,
  aux.thin        = 50,
  ladder          = 200,
  V.proposal      = 0.5,
  burn.in         = 5000,
  main.iters      = 30000,
  num.samples     = 25000, 
  estimate        = "CD",
  seed            = 1)
cj1$specs <- c("edges", "mutual", "absdiff.grade",
               "nodefactor.race.H", "nodefactor.race.O", "nodefactor.race.W",
               "nodefactor.grade.8", "nodefactor.grade.9",
               "nodefactor.grade.10", "nodefactor.grade.11",
               "nodefactor.grade.12", "nodefactor.sex.2", "nodematch.race.B",
               "nodematch.race.O", "nodematch.race.W", "nodematch.grade.7",
               "nodematch.grade.8", "nodematch.grade.9", "nodematch.grade.10",
               "nodematch.grade.11", "nodematch.grade.12", "nodematch.sex",
               "idegree0", "idegree1", "odegree0", "odegree1",
               "gwesp.fixed.0.1")

summary(cj1)

plot(cj1, lag = 200)

set.seed(1)
bgof(cj1,  
     sample.size = 100, 
     aux.iters   = 5000,
     n.ideg      = 20, 
     n.odeg      = 20, 
     n.dist      = 10, 
     n.esp       = 7)

m2 <- dixon ~ edges + absdiff("grade") + 
  nodefactor("race") + nodefactor("grade") + nodefactor("sex") +
  nodematch("race", diff = TRUE, levels = c("B", "O", "W")) + 
  nodematch("grade", diff = TRUE) + 
  nodematch("sex",  diff = FALSE) + 
  idegree(0:1) + odegree(0:1)

M.prior2 <- c(-5, rep(0, 24))
S.prior2 <- diag(5, 25)
cj2 <- evidence(
  evidence.method = "CJ",
  formula         = m2,
  prior.mean      = M.prior2,
  prior.sigma     = S.prior2,
  aux.iters       = 2500,
  n.aux.draws     = 50,
  aux.thin        = 50,
  ladder          = 200,
  V.proposal      = 0.5,
  burn.in         = 5000,
  main.iters      = 30000,
  num.samples     = 25000,
  estimate        = "CD",
  seed            = 1)

cj1$log.evidence
cj2$log.evidence

BF_12 <- exp(cj1$log.evidence - cj2$log.evidence)
BF_12