##### The following example is a TERGM adaptation of the #####
##### dyad-level example provided in figure 5(c) on page #####
##### 424 of Desmarais and Cranmer (2012) in the PSJ. At #####
##### each time step, it compares dyadic probabilities #####
##### (no tie, unidirectional tie, and reciprocal tie #####
##### probability) between a fitted model and a model #####
##### where the reciprocity effect is fixed at 0 based #####
##### on 20 randomly selected dyads per time step. The #####
##### results are visualized using a grouped bar plot. #####
## Not run:
# # create toy dataset and fit a model
# networks <- list()
# for (i in 1:3) { # create 3 random networks with 10 actors
# mat <- matrix(rbinom(100, 1, 0.25), nrow = 10, ncol = 10)
# diag(mat) <- 0 # loops are excluded
# nw <- network(mat) # create network object
# networks[[i]] <- nw # add network to the list
# }
# fit <- btergm(networks ~ edges + istar(2) + mutual, R = 200)
#
# # extract coefficients and create null hypothesis vector
# null <- coef(fit) # estimated coefs
# null[3] <- 0 # set mutual term = 0
#
# # sample 20 dyads per time step and compute probability ratios
# probabilities <- matrix(nrow = 9, ncol = length(networks))
# # nrow = 9 because three probabilities + upper and lower CIs
# colnames(probabilities) <- paste("t =", 1:length(networks))
# for (t in 1:length(networks)) {
# d <- dim(as.matrix(networks[[t]])) # how many row and column nodes?
# size <- d[1] * d[2] # size of the matrix
# nw <- matrix(1:size, nrow = d[1], ncol = d[2])
# nw <- nw[lower.tri(nw)] # sample only from lower triangle b/c
# samp <- sample(nw, 20) # dyadic probabilities are symmetric
# prob.est.00 <- numeric(0)
# prob.est.01 <- numeric(0)
# prob.est.11 <- numeric(0)
# prob.null.00 <- numeric(0)
# prob.null.01 <- numeric(0)
# prob.null.11 <- numeric(0)
# for (k in 1:20) {
# i <- arrayInd(samp[k], d)[1, 1] # recover 'i's and 'j's from sample
# j <- arrayInd(samp[k], d)[1, 2]
# # run interpretation function with estimated coefs and mutual = 0:
# int.est <- interpret(fit, type = "dyad", i = i, j = j, t = t)
# int.null <- interpret(fit, coefficients = null, type = "dyad",
# i = i, j = j, t = t)
# prob.est.00 <- c(prob.est.00, int.est[[1]][1, 1])
# prob.est.11 <- c(prob.est.11, int.est[[1]][2, 2])
# mean.est.01 <- (int.est[[1]][1, 2] + int.est[[1]][2, 1]) / 2
# prob.est.01 <- c(prob.est.01, mean.est.01)
# prob.null.00 <- c(prob.null.00, int.null[[1]][1, 1])
# prob.null.11 <- c(prob.null.11, int.null[[1]][2, 2])
# mean.null.01 <- (int.null[[1]][1, 2] + int.null[[1]][2, 1]) / 2
# prob.null.01 <- c(prob.null.01, mean.null.01)
# }
# prob.ratio.00 <- prob.est.00 / prob.null.00 # ratio of est. and null hyp
# prob.ratio.01 <- prob.est.01 / prob.null.01
# prob.ratio.11 <- prob.est.11 / prob.null.11
# probabilities[1, t] <- mean(prob.ratio.00) # mean estimated 00 tie prob
# probabilities[2, t] <- mean(prob.ratio.01) # mean estimated 01 tie prob
# probabilities[3, t] <- mean(prob.ratio.11) # mean estimated 11 tie prob
# ci.00 <- t.test(prob.ratio.00, conf.level = 0.99)$conf.int
# ci.01 <- t.test(prob.ratio.01, conf.level = 0.99)$conf.int
# ci.11 <- t.test(prob.ratio.11, conf.level = 0.99)$conf.int
# probabilities[4, t] <- ci.00[1] # lower 00 conf. interval
# probabilities[5, t] <- ci.01[1] # lower 01 conf. interval
# probabilities[6, t] <- ci.11[1] # lower 11 conf. interval
# probabilities[7, t] <- ci.00[2] # upper 00 conf. interval
# probabilities[8, t] <- ci.01[2] # upper 01 conf. interval
# probabilities[9, t] <- ci.11[2] # upper 11 conf. interval
# }
#
# # create barplots from probability ratios and CIs
# require("gplots")
# bp <- barplot2(probabilities[1:3, ], beside = TRUE, plot.ci = TRUE,
# ci.l = probabilities[4:6, ], ci.u = probabilities[7:9, ],
# col = c("tan", "tan2", "tan3"), ci.col = "grey40",
# xlab = "Dyadic tie values", ylab = "Estimated Prob./Null Prob.")
# mtext(1, at = bp, text = c("(0,0)", "(0,1)", "(1,1)"), line = 0, cex = 0.5)
#
#
# ##### The following examples illustrate the behavior of #####
# ##### the interpret function with undirected and/or #####
# ##### bipartite graphs with or without structural zeros. #####
#
# library("statnet")
# library("btergm")
#
# # micro-level interpretation for undirected network with structural zeros
# set.seed(12345)
# mat <- matrix(rbinom(400, 1, 0.1), nrow = 20, ncol = 20)
# mat[1, 5] <- 1
# mat[10, 7] <- 1
# mat[15, 3] <- 1
# mat[18, 4] < 1
# nw <- network(mat, directed = FALSE, bipartite = FALSE)
# cv <- matrix(rnorm(400), nrow = 20, ncol = 20)
# offsetmat <- matrix(rbinom(400, 1, 0.1), nrow = 20, ncol = 20)
# offsetmat[1, 5] <- 1
# offsetmat[10, 7] <- 1
# offsetmat[15, 3] <- 1
# offsetmat[18, 4] < 1
# model <- ergm(nw ~ edges + kstar(2) + edgecov(cv) + offset(edgecov(offsetmat)),
# offset.coef = -Inf)
# summary(model)
#
# # tie-level interpretation (note that dyad interpretation would not make any
# # sense in an undirected network):
# interpret(model, type = "tie", i = 1, j = 2) # 0.28 (= normal dyad)
# interpret(model, type = "tie", i = 1, j = 5) # 0.00 (= structural zero)
#
# # node-level interpretation; note the many 0 probabilities due to the
# # structural zeros; also note the warning message that the probabilities may
# # be slightly imprecise because -Inf needs to be approximated by some large
# # negative number (-9e8):
# interpret(model, type = "node", i = 1, j = 3:5)
#
# # repeat the same exercise for a directed network
# nw <- network(mat, directed = TRUE, bipartite = FALSE)
# model <- ergm(nw ~ edges + istar(2) + edgecov(cv) + offset(edgecov(offsetmat)),
# offset.coef = -Inf)
# interpret(model, type = "tie", i = 1, j = 2) # 0.13 (= normal dyad)
# interpret(model, type = "tie", i = 1, j = 5) # 0.00 (= structural zero)
# interpret(model, type = "dyad", i = 1, j = 2) # results for normal dyad
# interpret(model, type = "dyad", i = 1, j = 5) # results for i->j struct. zero
# interpret(model, type = "node", i = 1, j = 3:5)
#
# # micro-level interpretation for bipartite graph with structural zeros
# set.seed(12345)
# mat <- matrix(rbinom(200, 1, 0.1), nrow = 20, ncol = 10)
# mat[1, 5] <- 1
# mat[10, 7] <- 1
# mat[15, 3] <- 1
# mat[18, 4] < 1
# nw <- network(mat, directed = FALSE, bipartite = TRUE)
# cv <- matrix(rnorm(200), nrow = 20, ncol = 10) # some covariate
# offsetmat <- matrix(rbinom(200, 1, 0.1), nrow = 20, ncol = 10)
# offsetmat[1, 5] <- 1
# offsetmat[10, 7] <- 1
# offsetmat[15, 3] <- 1
# offsetmat[18, 4] < 1
# model <- ergm(nw ~ edges + b1star(2) + edgecov(cv)
# + offset(edgecov(offsetmat)), offset.coef = -Inf)
# summary(model)
#
# # tie-level interpretation; note the index for the second mode starts with 21
# interpret(model, type = "tie", i = 1, j = 21)
#
# # dyad-level interpretation does not make sense because network is undirected;
# # node-level interpretation prints warning due to structural zeros, but
# # computes the correct probabilities (though slightly imprecise because -Inf
# # is approximated by some small number:
# interpret(model, type = "node", i = 1, j = 21:25)
#
# # compute all dyadic probabilities
# dyads <- edgeprob(model)
# dyads
# ## End(Not run)
Run the code above in your browser using DataLab