R/utilities.R

In SSDALab/dnr: Simulate Dynamic Networks using Exponential Random Graph Models (ERGM) Family

##' General purpose regression engine for the methods bayesglm, glm and glmnet
##' @param XYdata matrix with X and Y columns. First column is named as y, other columns are X.
##' @param method string among ("glm", "glmnet", "bayesglm").
##' @param regIntercept Logical. Should intercept be included in the model? 
##' @param lambda for method "glmnet".
##' @param alpha for "glmnet"
##' @return list with elements: coef, se, lambda, fit (Coefficients, SE, lambda, if used, fit object.)
##' @author Abhirup

regEngine <- function(XYdata,
                      method = "bayesglm",
                      regIntercept = FALSE,
                      lambda = NA,
                      alpha = 1){
  if(method == "glmnet"){
    if(is.na(lambda)){
      mSelect <- glmnet::cv.glmnet(
        data.matrix(XYdata[, -1]),
        as.vector(XYdata[, 1]),
        family = "binomial",
        alpha = alpha)
      lambda <- mSelect$lambda.min
    }
    mFit <- glmnet::glmnet(data.matrix(XYdata[, -1]),
                           as.vector(XYdata[, 1]),
                           family = "binomial",
                           alpha = alpha,
                           lambda = lambda,
                           intercept = regIntercept)
    mCoef <- setNames(as.vector(mFit$beta),
                      colnames(XYdata[, -1]))
    mSE <- NA
    out <- list(coef = mCoef,
                se = mSE,
                lambda = lambda,
                fit = mFit)
  } else if (method == "glm") {
    if(regIntercept){
      mFit <- glm(y ~ ., data = XYdata,
                  family = binomial(link = "logit"))
      mSummery <- (summary(mFit))
      mCoef <- mFit$coefficients[-1]
      mSE <- numeric(length(mCoef))
      mSE[is.na(mCoef)] <- 0
      mSE[!is.na(mCoef)] <- mSummery$coefficients[-1 , 2]
      names(mSE) <- colnames(XYdata[, -1])
    } else{
      mFit <- glm(y ~ . -1, data = XYdata,
                  family = binomial(link = "logit"))
      mSummery <- (summary(mFit))
      mCoef <- mFit$coefficients
      mSE <- numeric(length(mCoef))
      mSE[is.na(mCoef)] <- 0
      mSE[!is.na(mCoef)] <- mSummery$coefficients[ , 2]
      names(mSE) <- colnames(XYdata[, -1])
    }
    out <- list(coef = mCoef,
                se = mSE,
                lambda = lambda,
                fit = mFit)
  } else if(method == "bayesglm"){
    if(regIntercept) {
      mFit <- arm::bayesglm(y ~., data = XYdata,
                            family = binomial(link = "logit"))
      mSummery <- (summary(mFit))
      mCoef <- mFit$coefficients[-1]
      mSE <- numeric(length(mCoef))
      mSE[is.na(mCoef)] <- 0
      mSE[!is.na(mCoef)] <- mSummery$coefficients[-1 , 2]
      names(mSE) <- colnames(XYdata[, -1])
    } else {
      mFit <- arm::bayesglm(y ~ . -1, data = XYdata,
                            family = binomial(link = "logit"))
      mSummery <- (summary(mFit))
      mCoef <- mFit$coefficients
      mSE <- numeric(length(mCoef))
      mSE[is.na(mCoef)] <- 0
      mSE[!is.na(mCoef)] <- mSummery$coefficients[ , 2]
      names(mSE) <- colnames(XYdata[, -1])
    }
    out <- list(coef = mCoef,
                se = mSE,
                lambda = lambda,
                fit = mFit)
  }
  return(out)
}




ungvectorize <- function(x,nvertex,gmode){
  if(gmode == "digraph"){
    out <- diag(nvertex)
    tmp <- c(out)
    tmp[tmp==0] <- x
    out <- matrix(tmp,nvertex)
    diag(out) <- 0
  } else{
    out <- diag(nvertex)
    tmp <- c(out)
    tmpid <- c(lower.tri(out))
    tmp[tmpid] <- x
    out <- matrix(tmp, nvertex)
    out[upper.tri(out)] <- t(out)[upper.tri(out)]
    diag(out) <- 0
  }
  out
}


##'binaryPlot
##' @title binaryPlot
##' @param x matrix
##' @param axlabs Binary, should the axis labels be shown.
##' @param ... title, xlabs, ylabs.
##' @description Plot for binary matrices, especially adjacency matrices.
##' @export

binaryPlot <- function(x, axlabs = TRUE,...){
  xmin <- min(x)
  xmax <- max(x)
  yLabels <- rownames(x)
  xLabels <- colnames(x)
  title <-c()
  if( length(list(...)) ){
    Lst <- list(...)
    if( !is.null(Lst$zlim) ){
      min <- Lst$zlim[1]
      max <- Lst$zlim[2]
    }
    if( !is.null(Lst$yLabels) ){
      yLabels <- c(Lst$yLabels)
    }
    if( !is.null(Lst$xLabels) ){
      xLabels <- c(Lst$xLabels)
    }
    if( !is.null(Lst$title) ){
      title <- Lst$title
    }
  }
  # check for null values
  if( is.null(xLabels) ){
    xLabels <- c(1:ncol(x))
  }
  if( is.null(yLabels) ){
    yLabels <- c(1:nrow(x))
  }
  # Reverse x
  reverse <- nrow(x) : 1
  yLabels <- yLabels[reverse]
  x <- x[reverse,]
  image(1:length(xLabels), 1:length(yLabels), t(x),
        col=c("#cccccc", "#525252"), xlab="",
        ylab="", axes=FALSE, zlim=c(xmin,xmax))
  if( !is.null(title) ){
    title(main=title)
  }
  if(axlabs){
    axis(BELOW<-1, at=1:length(xLabels), labels=xLabels, cex.axis=0.7)
    axis(LEFT <-2, at=1:length(yLabels), labels=yLabels, las= HORIZONTAL<-1,
         cex.axis=0.7)
  }
  grid(NULL, NULL, col = "#969696", lty = 6)
}


ilogit <- function(x) 1/(1+exp(-x))
thresh <- function(x, hi = 0.95, lo = 0.05) {
  x <- ifelse(x>hi, hi, x)
  x <- ifelse(x<lo, lo, x)
}

## supporting functions:
## network size (that allows NA)
network.size.1 <- function(x){
  if(!network::is.network(x)){
    if(is.na(x)) return(0)
  } else return(network::network.size(x))
}

## Remove networks that are of size 0.
rmNAnetworks <- function(netlist){
  netlens <- unlist(lapply(netlist, network.size.1))
  toremove <- which(netlens < 1)
  if(length(toremove) > 0) {
    return(netlist[-toremove])
  } else {
    return(netlist)
  }
}

network.vertex.names.1 <- function(x) {
  if(!network::is.network(x)){
    if(is.na(x)) return (NA)
  } else {
    return (network::network.vertex.names(x))
  }
}