R/qap_run.R

In ideanet: Integrating Data Exchange and Analysis for Networks ('ideanet')

#' Quadratic Assignment Procedure (\code{qap_run}).
#'
#' @description The \code{qap_run} function is a wrapper around \code{sna}'s Quadratic Assignment Procedure models \code{\link[sna:netlm]{sna::netlm}} and \code{\link[sna:netlogit]{sna::netlogit}}. It expects a networks objects containing dependent and independent variables of interest. It is required to use the output from \code{\link{qap_setup}}.
#'
#' @param net An \code{igraph} or \code{network} object.
#' @param dependent A string naming the dependent variable of interest. By default, the probability of a tie. Can also be the output of \code{\link{qap_setup}} using prefixes "same_", "diff_" or "abs_diff_".
#' @param variables A vector of strings naming the independent variables of interest. Must be the output of \code{\link{qap_setup}} using prefixes "same_", "diff_" and "abs_diff_", or suffixes "_ego" and "_alter".
#' @param directed A logical statement identifying if the network should be treated as directed. Defaults to \code{FALSE}.
#' @param family A string identifying the functional form. Options are \code{"linear"} and \code{"binomial"}. Defauts to \code{"linear"}.
#' @param reps A numeric value indicating the number of draws. Defaults to 500.
#' @return `qap_run` returns a list of elements that include:
#'
#' - \code{covs_df}, a data frame containing term labels, estimates, standard errors and p-values
#'
#' - \code{mods_df}, a data frame containing model-level information including the number of observations, AIC and BIC statistics.
#' @export
#'
#' @examples
#'
#' flor <- netwrite(nodelist = florentine_nodes,
#'                  node_id = "id",
#'                  i_elements = florentine_edges$source,
#'                  j_elements = florentine_edges$target,
#'                  type = florentine_edges$type,
#'                  directed = FALSE,
#'                  net_name = "florentine_graph")
#'
#' flor_setup <- qap_setup(flor$florentine_graph,
#'                         variables = c("total_degree"),
#'                         methods = c("difference"))
#'
#' flor_qap <- qap_run(flor_setup$graph,
#'                     variables = c("diff_total_degree"))
#'
#' # Inspect results
#' flor_qap$covs_df

qap_run <- function(net, dependent = NULL, variables, directed = FALSE, family = "linear", reps = 500) {

  tempvars <- c(dependent, variables)

  # remove isolates and convert to network
  if ("igraph" %in% class(net)) {
    iso <- which(igraph::degree(net) == 0)
    net <- igraph::delete_vertices(net, iso)

    # if graph is multiplex, have to merge it.
    if (igraph::any_multiple(net) == TRUE) {
      edge_attrs <- igraph::edge.attributes(net)
      edge_attr_comb <- lapply(edge_attrs, function(attr) {if (is.numeric(attr)) {"sum"} else {"random"}}) #sum for numeric, random for character.
      names(edge_attr_comb) <- names(edge_attrs)
      net <- igraph::simplify(net, edge.attr.comb = edge_attr_comb)}

    edges <- igraph::as_data_frame(net, what = "edges")
    nodes <- igraph::as_data_frame(net, what = "vertices")
    net <- network::network(x = edges, directed = directed)
  }

  # If there is an ego/alter call, create as many unique adjacency matrices as there are values.
  extra_mats <- list()
  if (directed == TRUE) {
    for (i in 1:length(variables)){
      name <- variables[[i]]
      if (grepl("^.*_ego$|^.*_alter$", name, fixed = FALSE)){
        varname <- gsub("^(.*)_.*$", "\\1", name)
        type <- gsub("^.*_(.*)", "\\1", name)
        unique_vals <- unique(nodes[[varname]])
        for (t in 1:length(unique_vals)) {
          mat <- as.matrix(net)
          mat[] <- 0
          varname_t <- paste0(name, "_", unique_vals[[t]])
          nodes[[varname_t]] <- ifelse(nodes[[varname]] == unique_vals[[t]], 1, 0)
          for (n in 1:nrow(nodes)) {
            if (type == "ego") {mat[rownames(mat) == nodes$name[[n]],] <- nodes[[varname_t]][[n]]}
            else if (type == "alter") {mat[,colnames(mat) == nodes$name[[n]]] <- nodes[[varname_t]][[n]]}

            extra_mats[[varname_t]] <- mat
          }
        }
      }
    }
  }

  # Create new values in case "both" is utilized
  for (i in 1:length(tempvars)){
    name <- tempvars[[i]]
    if (grepl("^both_", name, fixed = FALSE)) {
      varname <- gsub("^both_(.*)_.*", "\\1", name)
      varvalue <- gsub("^both_.*_(.*)", "\\1", name)
      nodes[[name]] <- ifelse(nodes[[varname]] == varvalue, 1, 0)
    }
  }

  # Get DV matrix
  if (is.null(dependent)) {
    dv <- as.matrix(net)
  } else if (grepl("^both_", dependent, fixed = FALSE)) {
    dv <- as.matrix(outer(nodes[[dependent]], nodes[[dependent]], "=="))
  } else if (grepl("^same_", dependent, fixed = FALSE)) {
    dv_name <- gsub("^same_", "", dependent)
    dv <- as.matrix(outer(nodes[[dv_name]], nodes[[dv_name]], "=="))
  } else if (grepl("^diff_", dependent, fixed = FALSE)) {
    dv_name <- gsub("^diff_", "", dependent)
    dv <- as.matrix(outer(nodes[[dv_name]], nodes[[dv_name]], "-"))
  } else if (grepl("^abs_diff_", dependent, fixed = FALSE)) {
    dv_name <- gsub("^abs_diff_", "", dependent)
    dv <- as.matrix(abs(outer(nodes[[dv_name]], nodes[[dv_name]], "-")))
  } else {
    stop(paste0(dependent, " is not an output of qap_setup() with prefix `same`,`diff`, `abs_diff` or `both`"))
  }

  # Get IV matrices
  ivs <- list()
  rem <- list() # list of empty matrices
  for (i in 1:length(variables)) {
    independent <- variables[[i]]
    if (grepl("same_", independent, fixed = FALSE)) {
      iv_name <- gsub("^same_", "", independent)
      iv <- as.matrix(outer(nodes[[iv_name]], nodes[[iv_name]], "=="))
    } else if (grepl("^both_", independent, fixed = FALSE)) {
      iv <- as.matrix(outer(nodes[[independent]], nodes[[independent]], "=="))
    } else if (grepl("^diff_", independent, fixed = FALSE)) {
      iv_name <- gsub("^diff_", "", independent)
      iv <- as.matrix(outer(nodes[[iv_name]], nodes[[iv_name]], "-"))
    } else if (grepl("^abs_diff_", independent, fixed = FALSE)) {
      iv_name <- gsub("^abs_diff_", "", independent)
      iv <- as.matrix(abs(outer(nodes[[iv_name]], nodes[[iv_name]], "-")))
    } else if (grepl("_ego$|_alter$", independent, fixed = FALSE)) {
      next
    } else {
      stop(paste0(independent, " is not an output of qap_setup() with prefix `same`,`diff`, `abs_diff` or `both`"))
    }

    if (sum(abs(iv), na.rm = TRUE) == 0) {# check if variables is empty
      warning(paste0("The variable ", independent, " is empty. It is excluded from the model."))
      rem[[i]] <- independent
    } else {ivs[[independent]] <- iv}
  }

  ivs <- ivs[lapply(ivs,length)>0]
  extra_mats <- extra_mats[-1] # holdout first element of categorical
  ivs <- c(ivs, extra_mats)

  # Run QAP
  if (length(ivs) != 0) {
    if (directed == TRUE) {mode = "digraph"} else {mode = "graph"}
    if (family == "binomial"){
      # if (all(dv %in% 0:1) == TRUE){
      res <- sna::netlogit(dv, ivs, reps = reps, mode = mode)
    } else if (family == "linear") {
      res <- sna::netlm(dv, ivs, reps = reps, mode = mode)
    } else {stop("Not an available family -- Try 'linear' or 'binomial'")}

    # Tidy results
    ivs_names <- names(ivs)
    rem <- unlist(rem[lapply(rem,length)>0])
    if (!is.null(rem)) {ivs_names <- ivs_names[-which(ivs_names %in% rem)]}
    ivs_names <- c("intercept", ivs_names)
    if (methods::is(res, "sna::netlogit")) {
      covs_df <- dplyr::tibble(covars = ivs_names, estimate = res$coefficients,
                               `exp(estimate)` = exp(res$coefficients), se = res$se,
                               pvalue = res$pgreqabs)
      mods_df <- dplyr::tibble(num_obs = res$n, aic = res$aic, bic = res$bic)
    } else {
      covs_df <- dplyr::tibble(covars = ivs_names, estimate = res$coefficients,
                               se = res$se, pvalue = res$pgreqabs)
      mods_df <- dplyr::tibble(num_obs = res$n, aic = res$aic, bic = res$bic)
    }

    output_list <- list(covs_df = covs_df, mods_df = mods_df)
    return(output_list)

  } else {stop("All IVs are empty. Try a different set of IVs")}
}