# NOT RUN {
##### 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
# }
Run the code above in your browser using DataLab