##### 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. #####
# 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)
Run the code above in your browser using DataLab