## ---- cache=FALSE, echo=FALSE----------------------------------------------------
library("ergm")


## ---- faux.mesa.high-------------------------------------------------------------
data("faux.mesa.high", package = "ergm")
(faux.mesa.high %v% "Grade") [1:20] 


## ---- node-attr-spec, tidy=FALSE-------------------------------------------------
summary(faux.mesa.high ~ nodecov("Grade")
                       + nodecov(~ Grade)
                       + nodecov(function(nw) nw %v% "Grade"))


## ---- node-attr-spec2, tidy=FALSE------------------------------------------------
summary(faux.mesa.high ~ nodecov(~ abs(Grade - mean(Grade)) / network.size(.)))


## ---- node-attr-spec3, tidy=FALSE------------------------------------------------
summary(faux.mesa.high ~ nodecov(~ (Grade - mean(Grade)) / network.size(.)))


## ---- node-attr-mat, tidy = FALSE------------------------------------------------
coef(summary(ergm(faux.mesa.high ~ edges + nodecov(~ cbind(Grade, Grade2 = Grade^2)))))


## ---- node-attr-mat2-------------------------------------------------------------
coef(summary(ergm(faux.mesa.high ~ edges + nodecov(~ poly(Grade, 2)))))


## ---- node-attr-rand-------------------------------------------------------------
set.seed(123)
randomcov <- structure(rbinom(network.size(faux.mesa.high), 1, 0.5), name = "random")
summary(faux.mesa.high ~ nodefactor(I(randomcov)))


## ---- node-attr-cats-------------------------------------------------------------
table(faux.mesa.high %v% "Grade")


## ---- node-attr-cats2------------------------------------------------------------
summary(faux.mesa.high ~ nodefactor(~ Grade, levels = -SMALLEST(3)))
summary(faux.mesa.high ~ nodefactor(~ Grade, levels = I(7:9)))
summary(faux.mesa.high ~ nodefactor(~ Grade, levels = c("7","8","9")))
summary(faux.mesa.high ~ nodefactor("Grade", levels = function(a) a %in% 7:9))
summary(faux.mesa.high ~ nodefactor("Grade", levels = ~ . %in% 7:9))


## ---- COLLAPSE_SMALLEST----------------------------------------------------------
summary(faux.mesa.high ~ nodefactor("Race", levels = TRUE))


## ---- COLLAPSE_SMALLEST2---------------------------------------------------------
library("magrittr")
summary(faux.mesa.high ~ nodefactor((~ Race) %>% COLLAPSE_SMALLEST(3, "BWO"), levels = TRUE))


## ---- mm-------------------------------------------------------------------------
summary(faux.mesa.high ~ mm(~ Grade >= 10))
summary(faux.mesa.high ~ mm(~ Grade >= 10, levels2 = NULL))


## ---- mm2, tidy=FALSE------------------------------------------------------------
summary(faux.mesa.high ~ mm("Grade",
  levels2 = ~ sapply(.levels, function(pair) pair[[1]] %in% c(7, 8) && pair[[2]] %in% c(7, 8))))

## ---- mm2b, tidy=FALSE-----------------------------------------------------------
summary(faux.mesa.high ~ mm(~ Grade, levels = TRUE ~ c("7","8"), levels2 = NULL ))


## ---- mm3, tidy=FALSE------------------------------------------------------------
summary(faux.mesa.high ~ mm(Grade >= 10 ~ Race, levels = TRUE ~ c("Hisp", "NatAm", "White")))


## ----nodemix---------------------------------------------------------------------
m <- matrix(c("homophilous", "", "", "homophilous"), 2, 2)
summary(faux.mesa.high ~ nodemix("Sex", levels2 = m))


## ---- sampson, out.width="100%", fig.cap="The monks dataset, with edges indicating directed liking relationships at any of three time points and nodes numbered from 1 to 18 and with group membership as assigned by Sampson indicated by L for Loyalists, O for Outcasts, and T for Young Turks."----
set.seed(2345)
data("sampson", package = "ergm")
lab <- paste0(1:18, " ", substr(samplike %v% "group", 1, 1), ": ", samplike %v% "vertex.names")
plot(samplike, displaylabels=TRUE, label = lab)


## ---- filter, cache=FALSE, tidy=FALSE--------------------------------------------
summary(samplike ~ nodematch("group", diff = TRUE, levels = "Turks")
                 + F(~ nodematch("group"), ~ nodefactor("group", levels = "Turks"))
                 + F(~ edges, ~ nodefactor("group", levels = "Turks") == 2)
                 + F(~ edges, ~ !nodefactor(~ group != "Turks")))


## ----nodematchTwice, tidy=FALSE--------------------------------------------------
cbind(summary(faux.mesa.high ~ nodematch("Grade")),
      summary(faux.mesa.high ~ F(~ edges, ~ absdiff("Grade") < 1)))


## ---- cloisterville, tidy=FALSE--------------------------------------------------
summary(samplike ~ ttriple + F(~ ttriple, ~ nodefactor("cloisterville") == 0))


## ---- undir, cache=FALSE, tidy=FALSE---------------------------------------------
cbind(summary(samplike ~ Symmetrize(~ edges, "weak") + Symmetrize(~ edges, "strong")  
          + Symmetrize(~ edges, "upper") + Symmetrize(~ edges, "lower")))


## ---- subgraph-induced, tidy=FALSE-----------------------------------------------
coef(summary(ergm(samplike ~ edges + mutual + S(~ edges + mutual, ~ (group == "Turks")),
                  control = snctrl(seed = 123))))


## ---- subgraph-between-----------------------------------------------------------
summary(samplike ~ S(~ cycle(4), (group != "Turks") ~ (group == "Turks")))


## ---- subgraph-between-undir-----------------------------------------------------
summary(samplike ~ Symmetrize(~ S(~ cycle(4), (group != "Turks") ~ (group == "Turks")), "weak"))


## ---- MutualityWithCovariate, tidy=FALSE-----------------------------------------
data("faux.dixon.high", package = "ergm")
FDHfit <- ergm(faux.dixon.high ~ edges + mutual + absdiff("grade")
               + Symmetrize(~ absdiff("grade"), "strong"), control = snctrl(seed = 321))
coef(summary(FDHfit))


## ---- interactions, tidy=FALSE---------------------------------------------------
summary(faux.mesa.high ~ nodefactor("Grade") : nodefactor("Sex", levels = TRUE))


## ---- interactions2--------------------------------------------------------------
m <- ergm(faux.mesa.high ~ edges + nodefactor("Grade") * nodefactor("Sex"))
print(summary(m), digits = 3)


## ---- DirectedInteract-----------------------------------------------------------
data("florentine", package = "ergm")
summary(flomarriage ~ nodecov("wealth") : nodecov("wealth") + nodecov(~ wealth^2))


## ---- lincomb-summ---------------------------------------------------------------
summary(samplike ~ nodeifactor("group", levels = TRUE))


## ---- lincomb-sum, tidy=FALSE----------------------------------------------------
summary(samplike ~ Sum(cbind(1, 0, 1) ~ nodeifactor("group", levels = TRUE), "nf.L_T") +
                    Sum("sum" ~ nodeifactor("group", levels = -2), "nf.L_T"))


## ---- lincomb-sums1, tidy=FALSE--------------------------------------------------
f1 <- samplike ~ edges + nodeifactor(~ group != "Outcasts")
summary(f1)
f2 <- samplike ~ edges + Sum(cbind(1, 0, 1) ~ nodeifactor("group", levels = TRUE), "nf.L_T")
summary(f2)


## ---- lincomb-sums1b, tidy=FALSE-------------------------------------------------
f3 <- samplike ~ edges + Parametrize(~ nodeifactor("group", levels = TRUE), "nf.L_T",
                               function(x, n, ...) c(x, 0, x), gradient = "linear")
summary(f3)


## ---- lincomb-sums1c, tidy=FALSE-------------------------------------------------
f4 <- samplike ~ edges + Parametrize(~ nodeifactor("group", levels = -2), "nf.L_T", "rep")
summary(f4)


## ---- lincomb-sums2, tidy=FALSE--------------------------------------------------
cbind( c(
  coef(ergm(f1))[2], coef(ergm(f2))[2], coef(ergm(f3))[2], coef(ergm(f4))[2]
))  


## ---- ProdIllustration, tidy=FALSE-----------------------------------------------
summary(faux.dixon.high ~ edges + mutual + Sum(list(~ edges, ~ mutual), "EdgesAndMutual")
        + Prod(list(~ edges, ~ mutual), "EdgesAndMutual"))


## ---- ColemanDataset-------------------------------------------------------------
library("sna")
data("coleman", package = "sna")
cole <- matrix(0, 2 * 73, 2 * 73)
cole[1 : 73, 1 : 73] <- coleman[1, , ]; cole[73 + (1 : 73), 73 + (1 : 73) ] <- coleman[2, , ]
diag(cole[1 : 73, 73 + (1 : 73) ] ) <- diag(cole[73 + (1 : 73), 1 : 73] ) <- 1
ncole <- network(cole)
ncole %v% "Semester" <- rep(c("Fall", "Spring"), each = 73)
ncole


## ---- ncole1---------------------------------------------------------------------
table(ncole %v% "Semester")


## ---- ncole2---------------------------------------------------------------------
summary(ncole ~ edges + nodematch("Semester"))


## ---- ncole3, tidy=FALSE---------------------------------------------------------
logit <- function(p) log(p / (1 - p))
cbind(logit((652 - 506) / (21170 - 10512)),
      coef(ergm(ncole ~ edges, constraints = ~ Dyads(fix = ~ nodematch("Semester")))))


## ---- ncole4, tidy=FALSE---------------------------------------------------------
cbind(logit(506 / 10512),
      coef(ergm(ncole ~ edges, constraints = ~ Dyads(vary = ~ nodematch("Semester")))))


## ---- ncole5---------------------------------------------------------------------
cbind(logit(652 / 21170), coef(ergm(ncole ~ edges)))


## ---- nodematchGroup-------------------------------------------------------------
summary(ncole ~ nodemix("Semester", levels = TRUE, levels2 = TRUE))


## ---- blocksOperator-------------------------------------------------------------
coef(ergm(ncole ~ edges, constraints = ~ blocks("Semester", levels2 = c(1, 4))))


## ---- blocksOperator2------------------------------------------------------------
coef(ergm(ncole ~ edges, constraints = ~ blocks("Semester", levels2 = c(2, 3))))


## ---- ValuedExampleConstruction, tidy=FALSE--------------------------------------
data("samplk", package = "ergm")
samplk.tot.m <- as.matrix(samplk1) + as.matrix(samplk2) + as.matrix(samplk3)
samplk.tot <- as.network(samplk.tot.m, directed = TRUE, matrix.type = "a",
                         ignore.eval = FALSE, names.eval = "nominations")


## ---- ValuedExampleSummary, tidy=FALSE-------------------------------------------
summary(samplk.tot ~ B(~ edges, ~ atleast(1)) + B(~ edges, ~ atleast(2))
                   + B(~ edges, ~ atleast(3)), response = "nominations")


## ---- ValuedExampleModelFit, tidy=FALSE------------------------------------------
mod <- ergm(samplk.tot ~ B(~ edges, ~ atleast(1)) + B(~ edges, ~ atleast(2))
          + B(~ edges, ~ atleast(3)), response = "nominations",
          reference = ~ DiscUnif(0, 3), control = snctrl(seed = 123))
coef(mod)
true <- c(EdgeVal0 = 218, EdgeVal1 = 38, EdgeVal2 = 20, EdgeVal3 = 30)
est <- c(1, exp(cumsum(coef(mod))), use.names = FALSE)
rbind(True_Proportions = true / sum(true), Estimated_Proportions = est / sum(est))


## ---- monks-full-fit, tidy=FALSE-------------------------------------------------
print(samplike)
summary(full.fit <- ergm(samplike ~ edges + mutual + transitiveties + cyclicalties,
                         eval.loglik = TRUE), control = snctrl(seed = 321))


## ---- monks-miss-----------------------------------------------------------------
samplike1 <- samplike
samplike1[1, ] <- NA
print(samplike1)


## ---- monks-miss-fit, tidy=FALSE-------------------------------------------------
summary(m1.fit <- ergm(samplike1 ~ edges + mutual + transitiveties + cyclicalties,
                       eval.loglik = TRUE), control = snctrl(seed = 321))


## ---- monks-miss-obs-fit, tidy = FALSE-------------------------------------------
samplike2 <- samplike
samplike2[1,] <- 0
samplike2 %v% "responded" <- rep(c(FALSE,TRUE),c(1,17))
print(samplike2)
summary(m2.fit <- ergm(samplike2 ~ edges + mutual + transitiveties + cyclicalties,
                       obs.constraints = ~ egocentric(~ responded, "out"),
                       control = snctrl(seed = 123)))


## ---- monks-miss-obs-ergmlhs-fit, tidy=FALSE-------------------------------------
samplike2 %ergmlhs% "obs.constraints" <- ~ egocentric(~ responded, "out")
summary(m3.fit <- ergm(samplike2 ~ edges + mutual + transitiveties + cyclicalties),
        control = snctrl(seed = 231))


## ---- allstatsBigExample, tidy=FALSE---------------------------------------------
system.time({
  EmptyNW <- network.initialize(8, directed = FALSE)
  a <- ergm.allstats(EmptyNW ~ edges + triangle + isolates + degree(4), force = TRUE)
})


## ---- meanval, warning=TRUE, tidy=FALSE------------------------------------------
ts <- summary(samplike ~ edges + mutual + transitiveties + cyclicalties)
emptynw <- network.initialize(network.size(samplike), directed = TRUE)
ts.fit <- ergm(emptynw ~ edges + mutual + transitiveties + cyclicalties,
               target.stats = ts, control = snctrl(seed = 123))
rbind(coef(full.fit), coef(ts.fit))


## ---- predict1-------------------------------------------------------------------
data("g4", package = "ergm")
g4 %v% "First" <- c(TRUE, FALSE, FALSE, FALSE)
SimpleERGM <- ergm(g4 ~ edges + nodecov("First"))
as.matrix(g4)


## ---- predict2-------------------------------------------------------------------
predict(SimpleERGM, conditional = TRUE, output = "matrix")


## ---- predict3-------------------------------------------------------------------
set.seed(123)
predict(SimpleERGM, conditional = FALSE, output = "matrix", nsim = 1000)


## ---- forceTRUE------------------------------------------------------------------
coef(ergm(g4 ~ edges, control = snctrl(force.main = TRUE, seed = 321)))


## ----glob-opts, cache=FALSE------------------------------------------------------
options(ergm.eval.loglik = FALSE)

sessionInfo()