R/interpretation.R
In leifeld/btergm: Temporal Exponential Random Graph Models by Bootstrapped Pseudolikelihood
Defines functions
marginalplot
interpret.btergm
interpret.ergm
dointerpret
Documented in
marginalplot
#' Micro-Level Interpretation of (T)ERGMs
#'
#' Micro-level interpretation of (T)ERGMs.
#'
#' The \code{interpret} function facilitates interpretation of ERGMs and TERGMs
#' at the micro level, as described in Desmarais and Cranmer (2012). There are
#' methods for \code{ergm} objects, \code{btergm} objects, and \code{mtergm}
#' objects. The function can be used to interpret these models at the tie or
#' edge level, dyad level, and block level. For example, what is the probability
#' that two specific nodes \code{i} (the sender) and \code{j} (the receiver) are
#' connected given the rest of the network and given the model? Or what is the
#' probability that any two nodes are tied at \code{t = 2} if they were tied (or
#' disconnected) at \code{t = 1} (i.e., what is the amount of tie stability)?
#' These tie- or edge-level questions can be answered if the \code{type = "tie"}
#' argument is used.
#'
#' Another example: What is the probability that node \code{i} has a tie to node
#' \code{j} but not vice-versa? Or that \code{i} and \code{j} maintain a
#' reciprocal tie? Or that they are disconnected? How much more or less likely
#' are \code{i} and \code{j} reciprocally connected if the \code{mutual} term in
#' the model is fixed at \code{0} (compared to the model that includes the
#' estimated parameter for reciprocity)? See example below. These dyad-level
#' questions can be answered if the \code{type = "dyad"} argument is used.
#'
#' Or what is the probability that a specific node \code{i} is connected to
#' nodes \code{j1} and \code{j2} but not to \code{j5} and \code{j7}? And how
#' likely is any node \code{i} to be connected to exactly four \code{j} nodes?
#' These node-level questions (focusing on the ties of node \code{i} or node
#' \code{j}) can be answered by using the \code{type = "node"} argument.
#'
#' The typical procedure is to manually enumerate all dyads or
#' sender-receiver-time combinations with certain properties and repeat the same
#' thing with some alternative properties for contrasting the two groups. Then
#' apply the \code{interpret} function to the two groups of dyads and compute a
#' measure of central tendency (e.g., mean or median) and possibly some
#' uncertainy measure (i.e., confidence intervals) from the distribution of
#' dyadic probabilities in each group. For example, if there is a gender
#' attribute, one can sample male-male or female-female dyads, compute the
#' distributions of edge probabilities for the two sets of dyads, and create
#' boxplots or barplots with confidence intervals for the two types of dyads in
#' order to contrast edge probabilities for male versus female same-sex dyads.
#'
#' See also the \code{\link{edgeprob}} function for automatic computation of all
#' dyadic edge probabilities.
#'
#' @param object An \code{ergm}, \code{btergm}, or \code{mtergm} object.
#' @param formula The formula to be used for computing probabilities. By
#' default, the formula embedded in the model object is retrieved and used.
#' @param coefficients The estimates on which probabilities should be based. By
#' default, the coefficients from the model object are retrieved and used.
#' Custom coefficients can be handed over, for example, in order to compare
#' versions of the model where the reciprocity term is fixed at \code{0}
#' versus versions of the model where the reciprocity term is left as in the
#' empirical result. This is one of the examples described in Desmarais and
#' Cranmer (2012).
#' @param target The response network on which probabilities are based.
#' Depending on whether the function is applied to an \code{ergm} or
#' \code{btergm}/\code{mtergm} object, this can be either a single network or
#' a list of networks. By default, the (list of) network(s) provided as the
#' left-hand side of the (T)ERGM formula is used.
#' @param type If \code{type = "tie"} is used, probabilities at the edge level
#' are computed. For example, what is the probability of a specific node
#' \code{i} to be connected to a specific node \code{j} given the rest of the
#' network and given the model? If \code{type = "dyad"} is used, probabilities
#' at the dyad level are computed. For example, what is the probability that
#' node \code{i} is connected to node \code{j} but not vice-versa, or what is
#' the probability that nodes \code{i} and \code{j} and mutually connected in
#' a directed network? If \code{type = "node"} is used, probabilities at the
#' node level are computed. For example, what is the probability that node
#' \code{i} is connected to a set of three other \code{j} nodes given the rest
#' of the network and the model?
#' @param i A single (sender) node \code{i} or a set of (sender) nodes \code{i}.
#' If \code{type = "node"} is used, this can be more than one node and should
#' be provided as a vector. The \code{i} argument can be either provided as
#' the index of the node in the sociomatrix (e.g., the fourth node would be
#' \code{i = 4}) or the row name of the node in the sociomatrix (e.g.,
#' \code{i = "Peter"}). If more than one node is provided and
#' \code{type = "node"}, there can be only one (receiver) node \code{j}. The
#' \code{i} and \code{j} arguments are used to specify for which nodes
#' probabilities should be computed. For example, what is the probability that
#' \code{i = 4} is connected to \code{j = 7}?
#' @param j A single (receiver) node \code{j} or a set of (receiver) nodes
#' \code{j}. If \code{type = "node"} is used, this can be more than one node
#' and should be provided as a vector. The \code{j} argument can be either
#' provided as the index of the node in the sociomatrix (e.g., the fourth node
#' would be \code{j = 4}) or the column name of the node in the sociomatrix
#' (e.g., \code{j = "Mary"}). If more than one node is provided and
#' \code{type = "node"}, there can be only one (sender) node \code{i}. The
#' \code{i} and \code{j} arguments are used to specify for which nodes
#' probabilities should be computed. For example, what is the probability that
#' \code{i = 4} is connected to \code{j = 7}?
#' @param t A vector of (numerical) time steps for which the probabilities
#' should be computed. This only applies to \code{btergm} amd \code{mtergm}
#' objects because \code{ergm} objects are by definition based on a single
#' time step. By default, all available time steps are used. It is, for
#' example, possible to compute probabilities only for a single time step by
#' specifying, e.g., \code{t = 5} in order to compute probabilities for the
#' fifth response network.
#' @param ... Further arguments to be passed on to subroutines.
#'
#' @references
#' Desmarais, Bruce A. and Skyler J. Cranmer (2012): Micro-Level Interpretation
#' of Exponential Random Graph Models with Application to Estuary Networks.
#' \emph{Policy Studies Journal} 40(3): 402--434.
#' \doi{10.1111/j.1541-0072.2012.00459.x}.
#'
#' Leifeld, Philip, Skyler J. Cranmer and Bruce A. Desmarais (2017): Temporal
#' Exponential Random Graph Models with btergm: Estimation and Bootstrap
#' Confidence Intervals. \emph{Journal of Statistical Software} 83(6): 1--36.
#' \doi{10.18637/jss.v083.i06}.
#'
#' Czarna, Anna Z., Philip Leifeld, Magdalena Smieja, Michael Dufner and Peter
#' Salovey (2016): Do Narcissism and Emotional Intelligence Win Us Friends?
#' Modeling Dynamics of Peer Popularity Using Inferential Network Analysis.
#' \emph{Personality and Social Psychology Bulletin} 42(11): 1588--1599.
#' \doi{10.1177/0146167216666265}.
#'
#' @examples
#' ##### 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. #####
#'
#' \dontrun{
#' # 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
#' }
#'
#' @docType methods
#' @aliases interpret-methods
#' @family interpretation
#' @export
setGeneric("interpret", function(object, ...) standardGeneric("interpret"),
package = "btergm")
#' Internal helper function for doing the actual interpretation computations
#'
#' Internal helper function for doing the actual interpretation computations.
#'
#' Takes the output of \code{\link{tergmprepare}} and a few parameters to
#' produce micro-level interpretation results. This is an internal worker
#' function. Users should call the high-level functions \code{interpret} and
#' \code{edgeprob} instead.
#'
#' @param l List of networks, with possibly just one network in this list, and a
#' formula (\code{l$form}, with temporal indices, as prepared by
#' \code{\link{tergmprepare}}).
#' @param type Type string. Can be \code{"tie"}, \code{"dyad"}, or
#' \code{"node"}.
#'
#' @importFrom utils combn
#' @importFrom network is.network is.directed
#'
#' @noRd
dointerpret <- function(l, coefficients, type, i, j, t) {
for (cv in 1:length(l$covnames)) {
assign(l$covnames[cv], l[[l$covnames[cv]]])
}
assign("offsmat", l$offsmat)
form <- as.formula(l$form)
# prepare i and j
if (!is.list(i)) {
i <- rep(list(i), length(l$networks))
num.actors <- numeric()
for (k in t) { # note that t can be a vector of time steps
num.actors[k] <- nrow(as.matrix(l$networks[[k]]))
}
if (length(table(num.actors)) > 1) {
warning(paste("'i' does not vary across time steps, but the number of",
"actors does. 'i' can be provided as a list or as a name."))
}
}
if (!is.list(j)) {
j <- rep(list(j), length(l$networks))
num.actors <- numeric()
for (k in t) {
num.actors[k] <- nrow(as.matrix(l$networks[[k]]))
}
if (length(table(num.actors)) > 1) {
warning(paste("'j' does not vary across time steps, but the number of",
"actors does. 'j' can be provided as a list or as a name."))
}
}
for (ll in 1:length(j)) {
if (length(j[[ll]]) > 1 && (type == "tie" || type == "dyad")) {
stop(paste("For computing dyadic or tie probabilities, only a single 'j'",
"node can be specified per time step."))
}
}
node_i <- i
node_j <- j
if (type == "tie") {
results <- numeric()
for (i in t) {
networks[[i]][node_i[[i]], node_j[[i]]] <- 0
stat0 <- summary(form, response = NULL)
networks[[i]][node_i[[i]], node_j[[i]]] <- 1
stat1 <- summary(form, response = NULL)
chgstat <- stat1 - stat0
if (length(chgstat) != length(coefficients)) {
stop(paste("Number of coefficients and statistics differ.",
"Did you fit a curved model? Curved models with non-fixed",
"parameters are currently not supported."))
}
for (k in 1:length(coefficients)) {
if (coefficients[k] == -Inf && chgstat[k] == 0) {
coefficients[k] <- 0 # no struct zero dyad; replace -Inf to avoid NaN
}
}
inf <- FALSE # next couple of lines: prevent errors with structural ones
for (k in 1:length(coefficients)) {
if (coefficients[k] == Inf && chgstat[k] != 0) {
inf <- TRUE
break
}
}
if (inf == TRUE) {
result <- 1
} else { # do the actual probability computations
coefficients[coefficients == Inf] <- 0
lp <- t(chgstat) %*% cbind(coefficients)
result <- c(1 / (1 + exp(-lp)))
}
names(result) <- "i->j = 1"
results[i] <- result
}
results <- results[!is.na(results)]
names(results) <- paste("t =", t)
} else if (type == "dyad") {
if (-Inf %in% coefficients) {
coefficients[coefficients == -Inf] <- -9e8
warning(paste("There are -Inf coefficients (possibly due to",
"offset terms). To yield interpretable results, -Inf is",
"approximated by a large negative number (-9e8). Note that",
"this may be imprecise."))
}
if (Inf %in% coefficients) {
coefficients[coefficients == Inf] <- 9e8
warning(paste("There are +Inf coefficients (possibly due to",
"offset terms). To yield interpretable results, +Inf is",
"approximated by a large positive number (9e8). Note that",
"this may be imprecise."))
}
results <- list()
for (i in t) {
# print error if undirected
if ((is.network(l$networks[[i]]) && !is.directed(l$networks[[i]])) ||
(is.matrix(l$networks[[i]]) && is.mat.directed(l$networks[[i]]) == FALSE)) {
stop(paste0("Network at t=", i, " is undirected. Dyadic ",
"probabilities do not make sense in undirected networks. Try ",
"type = \"tie\" instead!"))
}
eta_mat <- matrix(NA, 2, 2)
for (xi in 0:1) {
for (xj in 0:1) {
networks[[i]][node_i[[i]], node_j[[i]]] <- xi
networks[[i]][node_j[[i]], node_i[[i]]] <- xj
stat <- summary(form, response = NULL)
if (length(stat) != length(coefficients)) {
stop(paste("Number of coefficients and statistics differ.",
"Did you fit a curved model? Curved models with non-fixed",
"parameters are currently not supported."))
}
eta_mat[xi + 1, xj + 1] <- t(coefficients) %*% cbind(stat)
}
}
prob_mat <- matrix(NA, 2, 2)
for (xi in 0:1) {
for (xj in 0:1) {
etas <- c(eta_mat) - eta_mat[xi + 1, xj + 1]
prob_mat[xi + 1, xj + 1] <- 1 / (sum(exp(etas)))
}
}
rownames(prob_mat) <- c("i->j = 0", "i->j = 1")
colnames(prob_mat) <- c("j->i = 0", "j->i = 1")
results[[i]] <- prob_mat
}
results <- results[!sapply(results, is.null)]
names(results) <- paste("t =", t)
} else if (type == "node") {
if (-Inf %in% coefficients) {
coefficients[coefficients == -Inf] <- -9e8
warning(paste("There are -Inf coefficients (possibly due to",
"offset terms). To yield interpretable results, -Inf is",
"approximated by a large negative number (-9e8). Note that",
"this may be imprecise."))
}
if (Inf %in% coefficients) {
coefficients[coefficients == Inf] <- 9e8
warning(paste("There are +Inf coefficients (possibly due to",
"offset terms). To yield interpretable results, +Inf is",
"approximated by a large positive number (9e8). Note that",
"this may be imprecise."))
}
results <- list()
for (i in t) {
m <- length(node_i[[i]])
n <- length(node_j[[i]])
if (m == 1 && n > 1) {
labels <- c("Sender", "Receiver")
} else if (m > 1 && n == 1) {
labels <- c("Receiver", "Sender")
j.old <- node_j[[i]]
node_j[[i]] <- node_i[[i]]
node_i[[i]] <- j.old
m <- length(node_i[[i]])
n <- length(node_j[[i]])
} else {
stop(paste("Either 'i' or 'j' must contain more than one node (per",
"time step)."))
}
vecs <- rbind(rep(0, n), rep(1, n))
base <- rep(0, n)
for (l in 1:(n - 1)) {
places <- t(combn(1:n, l))
for (r in 1:nrow(places)) {
veci <- base
veci[places[r, ]] <- 1
vecs <- rbind(vecs, veci)
}
}
eta <- numeric(nrow(vecs))
for (l in 1:nrow(vecs)) {
ik <- node_i[[i]]
jk <- node_j[[i]]
networks[[i]][ik, jk] <- vecs[l, ]
stat <- summary(form, response = NULL)
if (length(stat) != length(coefficients)) {
stop(paste("Number of coefficients and statistics differ.",
"Did you fit a curved model? Curved models with non-fixed",
"parameters are currently not supported."))
}
eta[l] <- t(coefficients) %*% cbind(stat)
}
prob <- numeric(nrow(vecs))
for (l in 1:nrow(vecs)) {
prob[l] <- 1 / sum(exp(eta - eta[l]))
}
colnames(vecs) <- paste(labels[2], node_j[[i]])
rownames(vecs) <- rep(paste(labels[1], node_i[[i]]), nrow(vecs))
result <- cbind(prob, vecs)
colnames(result)[1] <- "probability"
results[[i]] <- result
}
results <- results[!sapply(results, is.null)]
names(results) <- paste("t =", t)
} else {
stop("'type' argument undefined or not recognized.")
}
return(results)
}
#' @noRd
interpret.ergm <- function(object, formula = getformula(object),
coefficients = coef(object), target = NULL, type = "tie", i, j) {
if (!is.null(target) && is.numeric(target)) {
stop(paste("'target' argument (short: 't') cannot be numeric."))
}
l <- tergmprepare(formula = formula, offset = FALSE, blockdiag = FALSE,
verbose = FALSE)
# extract response network and adjust formula
if (!is.null(target)) {
l$networks <- list(target)
}
# warn about indices with bipartite networks
if (is.network(l$networks[[1]]) && is.bipartite(l$networks[[1]]) ||
is.matrix(l$networks[[1]]) && !is.mat.onemode(l$networks[[1]])) {
if (is.numeric(j) && j < nrow(as.matrix(l$networks[[1]])) + 1) {
stop(paste0("In this bipartite network, the indices of the 'j' ",
"index start with ", nrow(as.matrix(l$networks[[1]])) + 1,
". Smaller indices denote within-block entries and do not make ",
"any sense."))
}
}
dointerpret(l, coefficients = coefficients, type = type, i = i, j = j,
t = 1)[[1]]
}
#' @noRd
interpret.btergm <- function(object, formula = getformula(object),
coefficients = coef(object), target = NULL, type = "tie", i, j,
t = 1:object@time.steps) {
l <- tergmprepare(formula = formula, offset = FALSE, blockdiag = FALSE,
verbose = FALSE)
# extract response networks and adjust formula
if (!is.null(target)) {
l$networks <- target
}
# warn about indices with bipartite networks
for (h in t) {
if (is.network(l$networks[[h]]) && is.bipartite(l$networks[[h]]) ||
is.matrix(l$networks[[h]]) && !is.mat.onemode(l$networks[[h]])) {
if (is.numeric(j) && j < nrow(as.matrix(l$networks[[h]])) + 1) {
stop(paste0("In this bipartite network, the indices of the 'j' ",
"index start with ", nrow(as.matrix(l$networks[[h]])) + 1,
" at t=", h, ". Smaller indices denote within-block entries and ",
"do not make any sense."))
}
}
}
dointerpret(l, coefficients = coefficients, type = type, i = i, j = j,
t = t)
}
#' @describeIn interpret Interpret method for \code{ergm} objects
setMethod("interpret", signature = className("ergm", "ergm"),
definition = interpret.ergm)
#' @describeIn interpret Interpret method for \code{btergm} objects
setMethod("interpret", signature = className("btergm", "btergm"),
definition = interpret.btergm)
#' @describeIn interpret Interpret method for \code{mtergm} objects
setMethod("interpret", signature = className("mtergm", "btergm"),
definition = interpret.btergm)
#' Create all predicted tie probabilities using MPLE
#'
#' Create all predicted tie probabilities using MPLE.
#'
#' For a given (T)ERGM, return a data frame with all predicted edge
#' probabilities along with the design matrix of the MPLE logit model, based
#' on the estimated coefficients and the design matrix, for all time points,
#' along with \code{i}, \code{j}, and \code{t} variables indicating where the
#' respective dyad is located.
#'
#' \code{edgeprob} is a convenience function that creates a data frame with all
#' dyads in the ERGM or TERGM along with their edge probabilities and their
#' predictor values (i.e., change statistics). This is useful for creating
#' marginal effects plots or contrasting multiple groups of dyads. This function
#' works faster than the \code{\link{interpret}} function.
#'
#' @param verbose Print details?
#' @inheritParams interpret
#'
#' @return The first variable in the resulting data frame contains the edge
#' value (i.e., the dependent variable, which is usually binary). The next
#' variables contain all the predictors from the ERGM or TERGM (i.e., the change
#' statistics). The next five variables contain the indices of the sender (i),
#' the receiver (j), the time step (t), the vertex id of i (i.name), and the
#' vertex id of j (j.name). These five variables serve to identify the dyad. The
#' last variable contains the computed edge probabilities.
#'
#' @aliases edgeprob
#' @family interpretation
#'
#' @importFrom network network is.network is.bipartite
#' @export
edgeprob <- function (object, verbose = FALSE) {
if ("ergm" %in% class(object)) {
tergm <- FALSE
} else if ("btergm" %in% class(object) || "mtergm" %in% class(object)) {
tergm <- TRUE
} else {
stop(paste("The edgeprob function is only applicable to ergm, btergm, and",
"mtergm objects."))
}
if (any(grepl("^gw.+#\\d{1,}$", names(coef(object))))) {
stop("MPLE-based (T)ERGMs with variable GW* decay are currently not supported.")
}
l <- tergmprepare(formula = getformula(object), offset = FALSE,
blockdiag = FALSE, verbose = FALSE)
for (cv in 1:length(l$covnames)) {
assign(l$covnames[cv], l[[l$covnames[cv]]])
}
assign("offsmat", l$offsmat)
form <- stats::as.formula(l$form)
covnames <- l$covnames[-1]
coefs <- coef(object)
if (verbose == TRUE) {
message("Creating data frame with predictors...")
}
Y <- NULL
dyads <- NULL
for (i in 1:length(l$networks)) {
mat <- as.matrix(l$networks[[i]])
imat <- matrix(rep(1:nrow(mat), ncol(mat)), nrow = nrow(mat))
if ((is.network(l$networks[[i]]) && is.bipartite(l$networks[[i]])) ||
(is.matrix(l$networks[[i]]) && is.mat.onemode(l$networks[[i]]) == FALSE)) {
mn <- nrow(mat) + 1
mx <- nrow(mat) + ncol(mat)
jmat <- matrix(rep(mn:mx, nrow(mat)), nrow = nrow(mat),
byrow = TRUE)
} else {
jmat <- matrix(rep(1:ncol(mat), nrow(mat)), nrow = nrow(mat),
byrow = TRUE)
}
f <- stats::as.formula(paste(l$form, " + edgecov(imat) + edgecov(jmat)"))
mpli <- ergm::ergmMPLE(f)
Y <- c(Y, mpli$response)
if (any(grepl("^gw.+#\\d{1,}$", colnames(mpli$predictor)))) {
stop("MPLE-based (T)ERGMs with variable GW* decay are currently not supported.")
}
mpli$predictor <- cbind(mpli$predictor, i)
# TODO: the previous line might say something like: "Note: Term nodeofactor("edu") skipped because it contributes no statistics." when a covariate is full of zeros.
# retrieve row and column names
mat.i <- as.matrix(l$networks[[i]])
i_ind <- mpli$predictor[, which(colnames(mpli$predictor) == "edgecov.imat")]
j_ind <- mpli$predictor[, which(colnames(mpli$predictor) == "edgecov.jmat")]
i.names <- rownames(mat.i)[i_ind]
if (isTRUE(l$bipartite)) {
j.names <- colnames(mat.i)[j_ind - nrow(mat.i)]
} else {
j.names <- colnames(mat.i)[j_ind]
}
if (is.null(dyads)) {
dyads <- cbind(as.data.frame(mpli$predictor), i.names, j.names)
} else { # before rbind can be use, we need to check if some covariate levels were absent either before or now
newColumns <- which(sapply(colnames(mpli$predictor), function(x) x %in% colnames(dyads)) == FALSE)
if (length(newColumns) > 0) {
for (j in 1:length(newColumns)) {
dyads <- cbind(dyads[, 1:(newColumns[j] - 1)], 0, dyads[, newColumns[j]:ncol(dyads)])
colnames(dyads)[newColumns[j]] <- names(newColumns)[j]
}
}
notPresentAnymore <- which(sapply(colnames(dyads)[1:(ncol(dyads) - 2)], function(x) x %in% colnames(mpli$predictor)) == FALSE)
if (length(notPresentAnymore) > 0) {
for (j in 1:length(notPresentAnymore)) {
mpli$predictor <- cbind(mpli$predictor[, 1:(notPresentAnymore[j] - 1)], 0, mpli$predictor[, notPresentAnymore[j]:ncol(mpli$predictor)])
colnames(mpli$predictor)[notPresentAnymore[j]] <- names(notPresentAnymore)[j]
}
}
dyads <- rbind(dyads, cbind(as.data.frame(mpli$predictor), i.names, j.names))
}
}
term.names <- colnames(dyads)[1:(ncol(dyads) - 5)]
term.names <- c(term.names, "i", "j", "t", "i.name", "j.name")
colnames(dyads) <- term.names
dyads <- cbind(Y, dyads)
colnames(dyads)[1] <- "tie"
class(dyads$i) <- "integer"
class(dyads$j) <- "integer"
class(dyads$t) <- "integer"
cf <- coef(object)
cf.length <- length(cf)
cf <- cf[!cf %in% c(Inf, -Inf)]
if (length(cf) != cf.length) {
warning(paste("There are structural zeros or ones. For these dyads, the",
"predicted probabilities are not valid and must be manually replaced",
"by 0 or 1, respectively."))
}
cbcoef <- cbind(cf)
chgstat <- dyads[, 2:(ncol(dyads) - 5)]
# handle decay term in curved ERGMs
if ("mtergm" %in% class(object) && ergm::is.curved(object@ergm)) {
curved.term <- vector(length = length(object@ergm$etamap$curved))
for (i in 1:length(object@ergm$etamap$curved)) {
curved.term[i] <- object@ergm$etamap$curved[[i]]$from[2]
}
cbcoef <- cbcoef[-c(curved.term)]
} else if ("ergm" %in% class(object) && ergm::is.curved(object)) {
curved.term <- vector(length = length(object$etamap$curved))
for (i in 1:length(object$etamap$curved)) {
curved.term[i] <- object$etamap$curved[[i]]$from[2]
}
cbcoef <- cbcoef[-c(curved.term)]
}
lp <- apply(chgstat, 1, function(x) t(x) %*% cbcoef)
result <- c(1 / (1 + exp(-lp)))
dyads <- cbind(dyads, result)
colnames(dyads)[ncol(dyads)] <- "probability"
dyads <- dyads[order(dyads$t, dyads$i, dyads$j), ]
rownames(dyads) <- NULL
return(dyads)
}
#' Plot marginal effects for two-way interactions in (T)ERGMs
#'
#' Plot marginal effects for two-way interactions in (T)ERGMs.
#'
#' The \code{marginalplot} function creates marginal effects plots for ERGMs
#' with interaction effects. The user has to supply the \code{ergm} object and
#' the coefficient names of the first main variable, the second main variable,
#' and the interaction term as stored in the coefficients vector inside the
#' \code{ergm} object. It is possible to draw continuous curves or discrete
#' error bars depending on the nature of the data (using the \code{point}
#' argument). The distribution of the second (conditioning) variable can be
#' plotted at the bottom of the viewport using the \code{rug} argument.
#'
#' The resulting marginal effects plot is a \code{ggplot2} plot. This means it
#' can be extended by plotting additional elements and using themes.
#'
#' @param model An \code{ergm} object as generated by the \pkg{ergm} package.
#' Note that marginal effects plots cannot be created for \code{btergm}
#' objects because the variance-covariance matrix is not valid. However, it
#' should be possible to apply the \code{marginalplot} function to
#' MCMC-MLE-estimated TERGMs because the \code{ergm} object is stored in the
#' \code{ergm} slot of an \code{mtergm} object. To do this, supply the
#' \code{ergm} object instead of the \code{mtergm} object (e.g.,
#' \code{marginalplot(mtergmobject@ergm)}).
#' @param var1 Name of the first main variable. This is the focal variable.
#' @param var2 Name of the second main variable. This is the conditioning
#' variable.
#' @param inter Name of the interaction effect.
#' @param ci Significance level.
#' @param rug Display the distribution of the conditioning variable at the
#' bottom of the plot?
#' @param point Display error bars for the levels of the conditioning variable
#' (instead of a continuous curve)?
#' @param structzeromat An optional matrix object which indicates dyads that
#' should be deleted prior to the calculation of the confidence interval for
#' the marginal effect curve. This is useful when such a matrix was used to
#' indicate structural zeros during estimation. In this event, the dyads
#' characterized by structural zeros are not allowed to be tied, therefore
#' they should be removed from the set of dyads used for the calculation of
#' marginal effects. The matrix should contain ones for structural zeros and
#' zeros for entries that should be used.
#' @param zeroline Draw a horizontal line to indicate zero for the first main
#' variable?
#' @param color Color of the curve, confidence interval, and distribution.
#' @param xlab Axis label for the second (conditioning) variable.
#' @param ylab Axis label for the first (focal) variable.
#'
#' @examples
#' \dontrun{
#' # data preparation
#' data("florentine")
#' n <- network.size(flobusiness)
#' wealth <- get.vertex.attribute(flobusiness, "wealth")
#' priorates <- get.vertex.attribute(flobusiness, "priorates")
#' wealth.icov <- matrix(rep(wealth, n), ncol = n, byrow = TRUE)
#' priorates.icov <- matrix(rep(priorates, n), ncol = n, byrow = TRUE)
#' interac <- wealth.icov * priorates.icov
#'
#' # estimate model with interaction effect
#' model <- ergm(flobusiness ~ edges + esp(1) + edgecov(wealth.icov)
#' + edgecov(priorates.icov) + edgecov(interac))
#'
#' # plot the interaction (note the additional optional ggplot2 elements)
#' marginalplot(model, var1 = "edgecov.wealth.icov",
#' var2 = "edgecov.priorates.icov", inter = "edgecov.interac",
#' color = "darkred", rug = TRUE, point = FALSE,
#' xlab = "Priorates", ylab = "Wealth") +
#' ggplot2::theme_bw() +
#' ggplot2::ggtitle("Interaction effect")
#' }
#'
#' @family interpretation
#' @export
marginalplot <- function(model, var1, var2, inter, ci = 0.95,
rug = FALSE, point = FALSE, structzeromat = NULL,
zeroline = TRUE, color = "black", xlab = NULL, ylab = NULL) {
# check arguments
if (!var1 %in% names(coef(model))) {
stop("'var1' not found.")
}
if (!var2 %in% names(coef(model))) {
stop("'var2' not found.")
}
# change statistics (needed for second variable)
ep <- edgeprob(model)
# remove '[[i]]' from edge covariate names in edgeprob
for (i in 1:ncol(ep)) {
if (grepl("((edge)|(dyad))cov", colnames(ep)[i])) {
colnames(ep)[i] <- substr(colnames(ep)[i], 1, nchar(colnames(ep)[i]) - 5)
}
}
# delete structural zeros from edgeprob output
if (!is.null(structzeromat)) {
include <- logical(nrow(ep))
for (r in 1:length(include)) {
if (structzeromat[ep$i[r], ep$j[r] - nrow(structzeromat)] == 1) {
include[r] <- FALSE
} else {
include[r] <- TRUE
}
}
ep <- ep[include == TRUE, ]
}
# unique values of the second variable
v2 <- sort(unlist(unique(ep[var2])))
names(v2) <- NULL
# coefficients
co <- coef(model)
beta1 <- co[match(var1, names(co))]
beta3 <- co[match(inter, names(co))]
# marginal effects
delta1 = beta1 + beta3 * v2
# variance-covariance matrix
if ("mtergm" %in% class(model)) {
model <- model@ergm
}
vc <- stats::vcov(model)
# variances
variances1 = vc[var1, var1] + (v2^2) * vc[inter, inter] +
(2 * v2 * vc[var1, inter])
# standard errors
se1 = sqrt(variances1)
# confidence intervals
z_score = qnorm(1 - ((1 - ci) / 2))
upper = delta1 + z_score * se1
lower = delta1 - z_score * se1
# aggregate in data frame
dta <- data.frame(v2 = v2, delta1 = delta1, upper = upper, lower = lower)
# axis labels
if (is.null(xlab)) {
xlab <- var2
}
if (is.null(ylab)) {
ylab <- var1
}
# plot error bars or line
if (point == TRUE) {
if (length(v2) == 2 && v2[1] == 0 && v2[2] == 1 && rug == FALSE) {
dta$v2 <- as.factor(dta$v2)
}
gp <- ggplot2::ggplot(data = dta, ggplot2::aes(x = v2, y = delta1)) +
ggplot2::geom_errorbar(data = dta,
ggplot2::aes(x = v2, ymin = lower, ymax = upper),
color = color) +
ggplot2::geom_point(data = dta, ggplot2::aes(x = v2, y = delta1)) +
ggplot2::ylab(ylab) +
ggplot2::xlab(xlab)
} else {
gp <- ggplot2::ggplot(data = dta, ggplot2::aes(x = v2, y = delta1)) +
ggplot2::geom_line(color = color) +
ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "lower", ymax = "upper"),
alpha = 0.15,
fill = color) +
ggplot2::ylab(ylab) +
ggplot2::xlab(xlab)
}
# add distribution to the plot
if (rug == TRUE) {
gp <- gp + ggplot2::geom_rug(data = ep[var2],
ggplot2::aes_string(x = var2),
inherit.aes = FALSE,
sides = "b",
col = color,
alpha = 0.1)
}
# add horizontal line for 0
if (zeroline == TRUE) {
gp <- gp + ggplot2::geom_hline(yintercept = 0, linetype = "dashed")
}
return(gp)
}
leifeld/btergm documentation built on April 3, 2024, 4:49 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions marginalplot interpret.btergm interpret.ergm dointerpret
Documented in marginalplot
#' Micro-Level Interpretation of (T)ERGMs
#'
#' Micro-level interpretation of (T)ERGMs.
#'
#' The \code{interpret} function facilitates interpretation of ERGMs and TERGMs
#' at the micro level, as described in Desmarais and Cranmer (2012). There are
#' methods for \code{ergm} objects, \code{btergm} objects, and \code{mtergm}
#' objects. The function can be used to interpret these models at the tie or
#' edge level, dyad level, and block level. For example, what is the probability
#' that two specific nodes \code{i} (the sender) and \code{j} (the receiver) are
#' connected given the rest of the network and given the model? Or what is the
#' probability that any two nodes are tied at \code{t = 2} if they were tied (or
#' disconnected) at \code{t = 1} (i.e., what is the amount of tie stability)?
#' These tie- or edge-level questions can be answered if the \code{type = "tie"}
#' argument is used.
#'
#' Another example: What is the probability that node \code{i} has a tie to node
#' \code{j} but not vice-versa? Or that \code{i} and \code{j} maintain a
#' reciprocal tie? Or that they are disconnected? How much more or less likely
#' are \code{i} and \code{j} reciprocally connected if the \code{mutual} term in
#' the model is fixed at \code{0} (compared to the model that includes the
#' estimated parameter for reciprocity)? See example below. These dyad-level
#' questions can be answered if the \code{type = "dyad"} argument is used.
#'
#' Or what is the probability that a specific node \code{i} is connected to
#' nodes \code{j1} and \code{j2} but not to \code{j5} and \code{j7}? And how
#' likely is any node \code{i} to be connected to exactly four \code{j} nodes?
#' These node-level questions (focusing on the ties of node \code{i} or node
#' \code{j}) can be answered by using the \code{type = "node"} argument.
#'
#' The typical procedure is to manually enumerate all dyads or
#' sender-receiver-time combinations with certain properties and repeat the same
#' thing with some alternative properties for contrasting the two groups. Then
#' apply the \code{interpret} function to the two groups of dyads and compute a
#' measure of central tendency (e.g., mean or median) and possibly some
#' uncertainy measure (i.e., confidence intervals) from the distribution of
#' dyadic probabilities in each group. For example, if there is a gender
#' attribute, one can sample male-male or female-female dyads, compute the
#' distributions of edge probabilities for the two sets of dyads, and create
#' boxplots or barplots with confidence intervals for the two types of dyads in
#' order to contrast edge probabilities for male versus female same-sex dyads.
#'
#' See also the \code{\link{edgeprob}} function for automatic computation of all
#' dyadic edge probabilities.
#'
#' @param object An \code{ergm}, \code{btergm}, or \code{mtergm} object.
#' @param formula The formula to be used for computing probabilities. By
#' default, the formula embedded in the model object is retrieved and used.
#' @param coefficients The estimates on which probabilities should be based. By
#' default, the coefficients from the model object are retrieved and used.
#' Custom coefficients can be handed over, for example, in order to compare
#' versions of the model where the reciprocity term is fixed at \code{0}
#' versus versions of the model where the reciprocity term is left as in the
#' empirical result. This is one of the examples described in Desmarais and
#' Cranmer (2012).
#' @param target The response network on which probabilities are based.
#' Depending on whether the function is applied to an \code{ergm} or
#' \code{btergm}/\code{mtergm} object, this can be either a single network or
#' a list of networks. By default, the (list of) network(s) provided as the
#' left-hand side of the (T)ERGM formula is used.
#' @param type If \code{type = "tie"} is used, probabilities at the edge level
#' are computed. For example, what is the probability of a specific node
#' \code{i} to be connected to a specific node \code{j} given the rest of the
#' network and given the model? If \code{type = "dyad"} is used, probabilities
#' at the dyad level are computed. For example, what is the probability that
#' node \code{i} is connected to node \code{j} but not vice-versa, or what is
#' the probability that nodes \code{i} and \code{j} and mutually connected in
#' a directed network? If \code{type = "node"} is used, probabilities at the
#' node level are computed. For example, what is the probability that node
#' \code{i} is connected to a set of three other \code{j} nodes given the rest
#' of the network and the model?
#' @param i A single (sender) node \code{i} or a set of (sender) nodes \code{i}.
#' If \code{type = "node"} is used, this can be more than one node and should
#' be provided as a vector. The \code{i} argument can be either provided as
#' the index of the node in the sociomatrix (e.g., the fourth node would be
#' \code{i = 4}) or the row name of the node in the sociomatrix (e.g.,
#' \code{i = "Peter"}). If more than one node is provided and
#' \code{type = "node"}, there can be only one (receiver) node \code{j}. The
#' \code{i} and \code{j} arguments are used to specify for which nodes
#' probabilities should be computed. For example, what is the probability that
#' \code{i = 4} is connected to \code{j = 7}?
#' @param j A single (receiver) node \code{j} or a set of (receiver) nodes
#' \code{j}. If \code{type = "node"} is used, this can be more than one node
#' and should be provided as a vector. The \code{j} argument can be either
#' provided as the index of the node in the sociomatrix (e.g., the fourth node
#' would be \code{j = 4}) or the column name of the node in the sociomatrix
#' (e.g., \code{j = "Mary"}). If more than one node is provided and
#' \code{type = "node"}, there can be only one (sender) node \code{i}. The
#' \code{i} and \code{j} arguments are used to specify for which nodes
#' probabilities should be computed. For example, what is the probability that
#' \code{i = 4} is connected to \code{j = 7}?
#' @param t A vector of (numerical) time steps for which the probabilities
#' should be computed. This only applies to \code{btergm} amd \code{mtergm}
#' objects because \code{ergm} objects are by definition based on a single
#' time step. By default, all available time steps are used. It is, for
#' example, possible to compute probabilities only for a single time step by
#' specifying, e.g., \code{t = 5} in order to compute probabilities for the
#' fifth response network.
#' @param ... Further arguments to be passed on to subroutines.
#'
#' @references
#' Desmarais, Bruce A. and Skyler J. Cranmer (2012): Micro-Level Interpretation
#' of Exponential Random Graph Models with Application to Estuary Networks.
#' \emph{Policy Studies Journal} 40(3): 402--434.
#' \doi{10.1111/j.1541-0072.2012.00459.x}.
#'
#' Leifeld, Philip, Skyler J. Cranmer and Bruce A. Desmarais (2017): Temporal
#' Exponential Random Graph Models with btergm: Estimation and Bootstrap
#' Confidence Intervals. \emph{Journal of Statistical Software} 83(6): 1--36.
#' \doi{10.18637/jss.v083.i06}.
#'
#' Czarna, Anna Z., Philip Leifeld, Magdalena Smieja, Michael Dufner and Peter
#' Salovey (2016): Do Narcissism and Emotional Intelligence Win Us Friends?
#' Modeling Dynamics of Peer Popularity Using Inferential Network Analysis.
#' \emph{Personality and Social Psychology Bulletin} 42(11): 1588--1599.
#' \doi{10.1177/0146167216666265}.
#'
#' @examples
#' ##### 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. #####
#'
#' \dontrun{
#' # 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
#' }
#'
#' @docType methods
#' @aliases interpret-methods
#' @family interpretation
#' @export
setGeneric("interpret", function(object, ...) standardGeneric("interpret"),
package = "btergm")
#' Internal helper function for doing the actual interpretation computations
#'
#' Internal helper function for doing the actual interpretation computations.
#'
#' Takes the output of \code{\link{tergmprepare}} and a few parameters to
#' produce micro-level interpretation results. This is an internal worker
#' function. Users should call the high-level functions \code{interpret} and
#' \code{edgeprob} instead.
#'
#' @param l List of networks, with possibly just one network in this list, and a
#' formula (\code{l$form}, with temporal indices, as prepared by
#' \code{\link{tergmprepare}}).
#' @param type Type string. Can be \code{"tie"}, \code{"dyad"}, or
#' \code{"node"}.
#'
#' @importFrom utils combn
#' @importFrom network is.network is.directed
#'
#' @noRd
dointerpret <- function(l, coefficients, type, i, j, t) {
for (cv in 1:length(l$covnames)) {
assign(l$covnames[cv], l[[l$covnames[cv]]])
}
assign("offsmat", l$offsmat)
form <- as.formula(l$form)
# prepare i and j
if (!is.list(i)) {
i <- rep(list(i), length(l$networks))
num.actors <- numeric()
for (k in t) { # note that t can be a vector of time steps
num.actors[k] <- nrow(as.matrix(l$networks[[k]]))
}
if (length(table(num.actors)) > 1) {
warning(paste("'i' does not vary across time steps, but the number of",
"actors does. 'i' can be provided as a list or as a name."))
}
}
if (!is.list(j)) {
j <- rep(list(j), length(l$networks))
num.actors <- numeric()
for (k in t) {
num.actors[k] <- nrow(as.matrix(l$networks[[k]]))
}
if (length(table(num.actors)) > 1) {
warning(paste("'j' does not vary across time steps, but the number of",
"actors does. 'j' can be provided as a list or as a name."))
}
}
for (ll in 1:length(j)) {
if (length(j[[ll]]) > 1 && (type == "tie" || type == "dyad")) {
stop(paste("For computing dyadic or tie probabilities, only a single 'j'",
"node can be specified per time step."))
}
}
node_i <- i
node_j <- j
if (type == "tie") {
results <- numeric()
for (i in t) {
networks[[i]][node_i[[i]], node_j[[i]]] <- 0
stat0 <- summary(form, response = NULL)
networks[[i]][node_i[[i]], node_j[[i]]] <- 1
stat1 <- summary(form, response = NULL)
chgstat <- stat1 - stat0
if (length(chgstat) != length(coefficients)) {
stop(paste("Number of coefficients and statistics differ.",
"Did you fit a curved model? Curved models with non-fixed",
"parameters are currently not supported."))
}
for (k in 1:length(coefficients)) {
if (coefficients[k] == -Inf && chgstat[k] == 0) {
coefficients[k] <- 0 # no struct zero dyad; replace -Inf to avoid NaN
}
}
inf <- FALSE # next couple of lines: prevent errors with structural ones
for (k in 1:length(coefficients)) {
if (coefficients[k] == Inf && chgstat[k] != 0) {
inf <- TRUE
break
}
}
if (inf == TRUE) {
result <- 1
} else { # do the actual probability computations
coefficients[coefficients == Inf] <- 0
lp <- t(chgstat) %*% cbind(coefficients)
result <- c(1 / (1 + exp(-lp)))
}
names(result) <- "i->j = 1"
results[i] <- result
}
results <- results[!is.na(results)]
names(results) <- paste("t =", t)
} else if (type == "dyad") {
if (-Inf %in% coefficients) {
coefficients[coefficients == -Inf] <- -9e8
warning(paste("There are -Inf coefficients (possibly due to",
"offset terms). To yield interpretable results, -Inf is",
"approximated by a large negative number (-9e8). Note that",
"this may be imprecise."))
}
if (Inf %in% coefficients) {
coefficients[coefficients == Inf] <- 9e8
warning(paste("There are +Inf coefficients (possibly due to",
"offset terms). To yield interpretable results, +Inf is",
"approximated by a large positive number (9e8). Note that",
"this may be imprecise."))
}
results <- list()
for (i in t) {
# print error if undirected
if ((is.network(l$networks[[i]]) && !is.directed(l$networks[[i]])) ||
(is.matrix(l$networks[[i]]) && is.mat.directed(l$networks[[i]]) == FALSE)) {
stop(paste0("Network at t=", i, " is undirected. Dyadic ",
"probabilities do not make sense in undirected networks. Try ",
"type = \"tie\" instead!"))
}
eta_mat <- matrix(NA, 2, 2)
for (xi in 0:1) {
for (xj in 0:1) {
networks[[i]][node_i[[i]], node_j[[i]]] <- xi
networks[[i]][node_j[[i]], node_i[[i]]] <- xj
stat <- summary(form, response = NULL)
if (length(stat) != length(coefficients)) {
stop(paste("Number of coefficients and statistics differ.",
"Did you fit a curved model? Curved models with non-fixed",
"parameters are currently not supported."))
}
eta_mat[xi + 1, xj + 1] <- t(coefficients) %*% cbind(stat)
}
}
prob_mat <- matrix(NA, 2, 2)
for (xi in 0:1) {
for (xj in 0:1) {
etas <- c(eta_mat) - eta_mat[xi + 1, xj + 1]
prob_mat[xi + 1, xj + 1] <- 1 / (sum(exp(etas)))
}
}
rownames(prob_mat) <- c("i->j = 0", "i->j = 1")
colnames(prob_mat) <- c("j->i = 0", "j->i = 1")
results[[i]] <- prob_mat
}
results <- results[!sapply(results, is.null)]
names(results) <- paste("t =", t)
} else if (type == "node") {
if (-Inf %in% coefficients) {
coefficients[coefficients == -Inf] <- -9e8
warning(paste("There are -Inf coefficients (possibly due to",
"offset terms). To yield interpretable results, -Inf is",
"approximated by a large negative number (-9e8). Note that",
"this may be imprecise."))
}
if (Inf %in% coefficients) {
coefficients[coefficients == Inf] <- 9e8
warning(paste("There are +Inf coefficients (possibly due to",
"offset terms). To yield interpretable results, +Inf is",
"approximated by a large positive number (9e8). Note that",
"this may be imprecise."))
}
results <- list()
for (i in t) {
m <- length(node_i[[i]])
n <- length(node_j[[i]])
if (m == 1 && n > 1) {
labels <- c("Sender", "Receiver")
} else if (m > 1 && n == 1) {
labels <- c("Receiver", "Sender")
j.old <- node_j[[i]]
node_j[[i]] <- node_i[[i]]
node_i[[i]] <- j.old
m <- length(node_i[[i]])
n <- length(node_j[[i]])
} else {
stop(paste("Either 'i' or 'j' must contain more than one node (per",
"time step)."))
}
vecs <- rbind(rep(0, n), rep(1, n))
base <- rep(0, n)
for (l in 1:(n - 1)) {
places <- t(combn(1:n, l))
for (r in 1:nrow(places)) {
veci <- base
veci[places[r, ]] <- 1
vecs <- rbind(vecs, veci)
}
}
eta <- numeric(nrow(vecs))
for (l in 1:nrow(vecs)) {
ik <- node_i[[i]]
jk <- node_j[[i]]
networks[[i]][ik, jk] <- vecs[l, ]
stat <- summary(form, response = NULL)
if (length(stat) != length(coefficients)) {
stop(paste("Number of coefficients and statistics differ.",
"Did you fit a curved model? Curved models with non-fixed",
"parameters are currently not supported."))
}
eta[l] <- t(coefficients) %*% cbind(stat)
}
prob <- numeric(nrow(vecs))
for (l in 1:nrow(vecs)) {
prob[l] <- 1 / sum(exp(eta - eta[l]))
}
colnames(vecs) <- paste(labels[2], node_j[[i]])
rownames(vecs) <- rep(paste(labels[1], node_i[[i]]), nrow(vecs))
result <- cbind(prob, vecs)
colnames(result)[1] <- "probability"
results[[i]] <- result
}
results <- results[!sapply(results, is.null)]
names(results) <- paste("t =", t)
} else {
stop("'type' argument undefined or not recognized.")
}
return(results)
}
#' @noRd
interpret.ergm <- function(object, formula = getformula(object),
coefficients = coef(object), target = NULL, type = "tie", i, j) {
if (!is.null(target) && is.numeric(target)) {
stop(paste("'target' argument (short: 't') cannot be numeric."))
}
l <- tergmprepare(formula = formula, offset = FALSE, blockdiag = FALSE,
verbose = FALSE)
# extract response network and adjust formula
if (!is.null(target)) {
l$networks <- list(target)
}
# warn about indices with bipartite networks
if (is.network(l$networks[[1]]) && is.bipartite(l$networks[[1]]) ||
is.matrix(l$networks[[1]]) && !is.mat.onemode(l$networks[[1]])) {
if (is.numeric(j) && j < nrow(as.matrix(l$networks[[1]])) + 1) {
stop(paste0("In this bipartite network, the indices of the 'j' ",
"index start with ", nrow(as.matrix(l$networks[[1]])) + 1,
". Smaller indices denote within-block entries and do not make ",
"any sense."))
}
}
dointerpret(l, coefficients = coefficients, type = type, i = i, j = j,
t = 1)[[1]]
}
#' @noRd
interpret.btergm <- function(object, formula = getformula(object),
coefficients = coef(object), target = NULL, type = "tie", i, j,
t = 1:object@time.steps) {
l <- tergmprepare(formula = formula, offset = FALSE, blockdiag = FALSE,
verbose = FALSE)
# extract response networks and adjust formula
if (!is.null(target)) {
l$networks <- target
}
# warn about indices with bipartite networks
for (h in t) {
if (is.network(l$networks[[h]]) && is.bipartite(l$networks[[h]]) ||
is.matrix(l$networks[[h]]) && !is.mat.onemode(l$networks[[h]])) {
if (is.numeric(j) && j < nrow(as.matrix(l$networks[[h]])) + 1) {
stop(paste0("In this bipartite network, the indices of the 'j' ",
"index start with ", nrow(as.matrix(l$networks[[h]])) + 1,
" at t=", h, ". Smaller indices denote within-block entries and ",
"do not make any sense."))
}
}
}
dointerpret(l, coefficients = coefficients, type = type, i = i, j = j,
t = t)
}
#' @describeIn interpret Interpret method for \code{ergm} objects
setMethod("interpret", signature = className("ergm", "ergm"),
definition = interpret.ergm)
#' @describeIn interpret Interpret method for \code{btergm} objects
setMethod("interpret", signature = className("btergm", "btergm"),
definition = interpret.btergm)
#' @describeIn interpret Interpret method for \code{mtergm} objects
setMethod("interpret", signature = className("mtergm", "btergm"),
definition = interpret.btergm)
#' Create all predicted tie probabilities using MPLE
#'
#' Create all predicted tie probabilities using MPLE.
#'
#' For a given (T)ERGM, return a data frame with all predicted edge
#' probabilities along with the design matrix of the MPLE logit model, based
#' on the estimated coefficients and the design matrix, for all time points,
#' along with \code{i}, \code{j}, and \code{t} variables indicating where the
#' respective dyad is located.
#'
#' \code{edgeprob} is a convenience function that creates a data frame with all
#' dyads in the ERGM or TERGM along with their edge probabilities and their
#' predictor values (i.e., change statistics). This is useful for creating
#' marginal effects plots or contrasting multiple groups of dyads. This function
#' works faster than the \code{\link{interpret}} function.
#'
#' @param verbose Print details?
#' @inheritParams interpret
#'
#' @return The first variable in the resulting data frame contains the edge
#' value (i.e., the dependent variable, which is usually binary). The next
#' variables contain all the predictors from the ERGM or TERGM (i.e., the change
#' statistics). The next five variables contain the indices of the sender (i),
#' the receiver (j), the time step (t), the vertex id of i (i.name), and the
#' vertex id of j (j.name). These five variables serve to identify the dyad. The
#' last variable contains the computed edge probabilities.
#'
#' @aliases edgeprob
#' @family interpretation
#'
#' @importFrom network network is.network is.bipartite
#' @export
edgeprob <- function (object, verbose = FALSE) {
if ("ergm" %in% class(object)) {
tergm <- FALSE
} else if ("btergm" %in% class(object) || "mtergm" %in% class(object)) {
tergm <- TRUE
} else {
stop(paste("The edgeprob function is only applicable to ergm, btergm, and",
"mtergm objects."))
}
if (any(grepl("^gw.+#\\d{1,}$", names(coef(object))))) {
stop("MPLE-based (T)ERGMs with variable GW* decay are currently not supported.")
}
l <- tergmprepare(formula = getformula(object), offset = FALSE,
blockdiag = FALSE, verbose = FALSE)
for (cv in 1:length(l$covnames)) {
assign(l$covnames[cv], l[[l$covnames[cv]]])
}
assign("offsmat", l$offsmat)
form <- stats::as.formula(l$form)
covnames <- l$covnames[-1]
coefs <- coef(object)
if (verbose == TRUE) {
message("Creating data frame with predictors...")
}
Y <- NULL
dyads <- NULL
for (i in 1:length(l$networks)) {
mat <- as.matrix(l$networks[[i]])
imat <- matrix(rep(1:nrow(mat), ncol(mat)), nrow = nrow(mat))
if ((is.network(l$networks[[i]]) && is.bipartite(l$networks[[i]])) ||
(is.matrix(l$networks[[i]]) && is.mat.onemode(l$networks[[i]]) == FALSE)) {
mn <- nrow(mat) + 1
mx <- nrow(mat) + ncol(mat)
jmat <- matrix(rep(mn:mx, nrow(mat)), nrow = nrow(mat),
byrow = TRUE)
} else {
jmat <- matrix(rep(1:ncol(mat), nrow(mat)), nrow = nrow(mat),
byrow = TRUE)
}
f <- stats::as.formula(paste(l$form, " + edgecov(imat) + edgecov(jmat)"))
mpli <- ergm::ergmMPLE(f)
Y <- c(Y, mpli$response)
if (any(grepl("^gw.+#\\d{1,}$", colnames(mpli$predictor)))) {
stop("MPLE-based (T)ERGMs with variable GW* decay are currently not supported.")
}
mpli$predictor <- cbind(mpli$predictor, i)
# TODO: the previous line might say something like: "Note: Term nodeofactor("edu") skipped because it contributes no statistics." when a covariate is full of zeros.
# retrieve row and column names
mat.i <- as.matrix(l$networks[[i]])
i_ind <- mpli$predictor[, which(colnames(mpli$predictor) == "edgecov.imat")]
j_ind <- mpli$predictor[, which(colnames(mpli$predictor) == "edgecov.jmat")]
i.names <- rownames(mat.i)[i_ind]
if (isTRUE(l$bipartite)) {
j.names <- colnames(mat.i)[j_ind - nrow(mat.i)]
} else {
j.names <- colnames(mat.i)[j_ind]
}
if (is.null(dyads)) {
dyads <- cbind(as.data.frame(mpli$predictor), i.names, j.names)
} else { # before rbind can be use, we need to check if some covariate levels were absent either before or now
newColumns <- which(sapply(colnames(mpli$predictor), function(x) x %in% colnames(dyads)) == FALSE)
if (length(newColumns) > 0) {
for (j in 1:length(newColumns)) {
dyads <- cbind(dyads[, 1:(newColumns[j] - 1)], 0, dyads[, newColumns[j]:ncol(dyads)])
colnames(dyads)[newColumns[j]] <- names(newColumns)[j]
}
}
notPresentAnymore <- which(sapply(colnames(dyads)[1:(ncol(dyads) - 2)], function(x) x %in% colnames(mpli$predictor)) == FALSE)
if (length(notPresentAnymore) > 0) {
for (j in 1:length(notPresentAnymore)) {
mpli$predictor <- cbind(mpli$predictor[, 1:(notPresentAnymore[j] - 1)], 0, mpli$predictor[, notPresentAnymore[j]:ncol(mpli$predictor)])
colnames(mpli$predictor)[notPresentAnymore[j]] <- names(notPresentAnymore)[j]
}
}
dyads <- rbind(dyads, cbind(as.data.frame(mpli$predictor), i.names, j.names))
}
}
term.names <- colnames(dyads)[1:(ncol(dyads) - 5)]
term.names <- c(term.names, "i", "j", "t", "i.name", "j.name")
colnames(dyads) <- term.names
dyads <- cbind(Y, dyads)
colnames(dyads)[1] <- "tie"
class(dyads$i) <- "integer"
class(dyads$j) <- "integer"
class(dyads$t) <- "integer"
cf <- coef(object)
cf.length <- length(cf)
cf <- cf[!cf %in% c(Inf, -Inf)]
if (length(cf) != cf.length) {
warning(paste("There are structural zeros or ones. For these dyads, the",
"predicted probabilities are not valid and must be manually replaced",
"by 0 or 1, respectively."))
}
cbcoef <- cbind(cf)
chgstat <- dyads[, 2:(ncol(dyads) - 5)]
# handle decay term in curved ERGMs
if ("mtergm" %in% class(object) && ergm::is.curved(object@ergm)) {
curved.term <- vector(length = length(object@ergm$etamap$curved))
for (i in 1:length(object@ergm$etamap$curved)) {
curved.term[i] <- object@ergm$etamap$curved[[i]]$from[2]
}
cbcoef <- cbcoef[-c(curved.term)]
} else if ("ergm" %in% class(object) && ergm::is.curved(object)) {
curved.term <- vector(length = length(object$etamap$curved))
for (i in 1:length(object$etamap$curved)) {
curved.term[i] <- object$etamap$curved[[i]]$from[2]
}
cbcoef <- cbcoef[-c(curved.term)]
}
lp <- apply(chgstat, 1, function(x) t(x) %*% cbcoef)
result <- c(1 / (1 + exp(-lp)))
dyads <- cbind(dyads, result)
colnames(dyads)[ncol(dyads)] <- "probability"
dyads <- dyads[order(dyads$t, dyads$i, dyads$j), ]
rownames(dyads) <- NULL
return(dyads)
}
#' Plot marginal effects for two-way interactions in (T)ERGMs
#'
#' Plot marginal effects for two-way interactions in (T)ERGMs.
#'
#' The \code{marginalplot} function creates marginal effects plots for ERGMs
#' with interaction effects. The user has to supply the \code{ergm} object and
#' the coefficient names of the first main variable, the second main variable,
#' and the interaction term as stored in the coefficients vector inside the
#' \code{ergm} object. It is possible to draw continuous curves or discrete
#' error bars depending on the nature of the data (using the \code{point}
#' argument). The distribution of the second (conditioning) variable can be
#' plotted at the bottom of the viewport using the \code{rug} argument.
#'
#' The resulting marginal effects plot is a \code{ggplot2} plot. This means it
#' can be extended by plotting additional elements and using themes.
#'
#' @param model An \code{ergm} object as generated by the \pkg{ergm} package.
#' Note that marginal effects plots cannot be created for \code{btergm}
#' objects because the variance-covariance matrix is not valid. However, it
#' should be possible to apply the \code{marginalplot} function to
#' MCMC-MLE-estimated TERGMs because the \code{ergm} object is stored in the
#' \code{ergm} slot of an \code{mtergm} object. To do this, supply the
#' \code{ergm} object instead of the \code{mtergm} object (e.g.,
#' \code{marginalplot(mtergmobject@ergm)}).
#' @param var1 Name of the first main variable. This is the focal variable.
#' @param var2 Name of the second main variable. This is the conditioning
#' variable.
#' @param inter Name of the interaction effect.
#' @param ci Significance level.
#' @param rug Display the distribution of the conditioning variable at the
#' bottom of the plot?
#' @param point Display error bars for the levels of the conditioning variable
#' (instead of a continuous curve)?
#' @param structzeromat An optional matrix object which indicates dyads that
#' should be deleted prior to the calculation of the confidence interval for
#' the marginal effect curve. This is useful when such a matrix was used to
#' indicate structural zeros during estimation. In this event, the dyads
#' characterized by structural zeros are not allowed to be tied, therefore
#' they should be removed from the set of dyads used for the calculation of
#' marginal effects. The matrix should contain ones for structural zeros and
#' zeros for entries that should be used.
#' @param zeroline Draw a horizontal line to indicate zero for the first main
#' variable?
#' @param color Color of the curve, confidence interval, and distribution.
#' @param xlab Axis label for the second (conditioning) variable.
#' @param ylab Axis label for the first (focal) variable.
#'
#' @examples
#' \dontrun{
#' # data preparation
#' data("florentine")
#' n <- network.size(flobusiness)
#' wealth <- get.vertex.attribute(flobusiness, "wealth")
#' priorates <- get.vertex.attribute(flobusiness, "priorates")
#' wealth.icov <- matrix(rep(wealth, n), ncol = n, byrow = TRUE)
#' priorates.icov <- matrix(rep(priorates, n), ncol = n, byrow = TRUE)
#' interac <- wealth.icov * priorates.icov
#'
#' # estimate model with interaction effect
#' model <- ergm(flobusiness ~ edges + esp(1) + edgecov(wealth.icov)
#' + edgecov(priorates.icov) + edgecov(interac))
#'
#' # plot the interaction (note the additional optional ggplot2 elements)
#' marginalplot(model, var1 = "edgecov.wealth.icov",
#' var2 = "edgecov.priorates.icov", inter = "edgecov.interac",
#' color = "darkred", rug = TRUE, point = FALSE,
#' xlab = "Priorates", ylab = "Wealth") +
#' ggplot2::theme_bw() +
#' ggplot2::ggtitle("Interaction effect")
#' }
#'
#' @family interpretation
#' @export
marginalplot <- function(model, var1, var2, inter, ci = 0.95,
rug = FALSE, point = FALSE, structzeromat = NULL,
zeroline = TRUE, color = "black", xlab = NULL, ylab = NULL) {
# check arguments
if (!var1 %in% names(coef(model))) {
stop("'var1' not found.")
}
if (!var2 %in% names(coef(model))) {
stop("'var2' not found.")
}
# change statistics (needed for second variable)
ep <- edgeprob(model)
# remove '[[i]]' from edge covariate names in edgeprob
for (i in 1:ncol(ep)) {
if (grepl("((edge)|(dyad))cov", colnames(ep)[i])) {
colnames(ep)[i] <- substr(colnames(ep)[i], 1, nchar(colnames(ep)[i]) - 5)
}
}
# delete structural zeros from edgeprob output
if (!is.null(structzeromat)) {
include <- logical(nrow(ep))
for (r in 1:length(include)) {
if (structzeromat[ep$i[r], ep$j[r] - nrow(structzeromat)] == 1) {
include[r] <- FALSE
} else {
include[r] <- TRUE
}
}
ep <- ep[include == TRUE, ]
}
# unique values of the second variable
v2 <- sort(unlist(unique(ep[var2])))
names(v2) <- NULL
# coefficients
co <- coef(model)
beta1 <- co[match(var1, names(co))]
beta3 <- co[match(inter, names(co))]
# marginal effects
delta1 = beta1 + beta3 * v2
# variance-covariance matrix
if ("mtergm" %in% class(model)) {
model <- model@ergm
}
vc <- stats::vcov(model)
# variances
variances1 = vc[var1, var1] + (v2^2) * vc[inter, inter] +
(2 * v2 * vc[var1, inter])
# standard errors
se1 = sqrt(variances1)
# confidence intervals
z_score = qnorm(1 - ((1 - ci) / 2))
upper = delta1 + z_score * se1
lower = delta1 - z_score * se1
# aggregate in data frame
dta <- data.frame(v2 = v2, delta1 = delta1, upper = upper, lower = lower)
# axis labels
if (is.null(xlab)) {
xlab <- var2
}
if (is.null(ylab)) {
ylab <- var1
}
# plot error bars or line
if (point == TRUE) {
if (length(v2) == 2 && v2[1] == 0 && v2[2] == 1 && rug == FALSE) {
dta$v2 <- as.factor(dta$v2)
}
gp <- ggplot2::ggplot(data = dta, ggplot2::aes(x = v2, y = delta1)) +
ggplot2::geom_errorbar(data = dta,
ggplot2::aes(x = v2, ymin = lower, ymax = upper),
color = color) +
ggplot2::geom_point(data = dta, ggplot2::aes(x = v2, y = delta1)) +
ggplot2::ylab(ylab) +
ggplot2::xlab(xlab)
} else {
gp <- ggplot2::ggplot(data = dta, ggplot2::aes(x = v2, y = delta1)) +
ggplot2::geom_line(color = color) +
ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "lower", ymax = "upper"),
alpha = 0.15,
fill = color) +
ggplot2::ylab(ylab) +
ggplot2::xlab(xlab)
}
# add distribution to the plot
if (rug == TRUE) {
gp <- gp + ggplot2::geom_rug(data = ep[var2],
ggplot2::aes_string(x = var2),
inherit.aes = FALSE,
sides = "b",
col = color,
alpha = 0.1)
}
# add horizontal line for 0
if (zeroline == TRUE) {
gp <- gp + ggplot2::geom_hline(yintercept = 0, linetype = "dashed")
}
return(gp)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.