R/extended_MPLE.R

Defines functions ex_dh ex_ddiagonal ex_dttriads ex_dctriads ex_drecip ex_dedgeweight ex_dout2star ex_din2star ex_net2xy ex_dtexp ex_pl ex_mple

# maximum pseudo-likelihood estimates
ex_mple <- function(GERGM_Object,
                    verbose = TRUE) {
  xy <- ex_net2xy(GERGM_Object)
  x <- xy$x
  y <- xy$y
  # why do we select this initialization
  est <- coef(lm(y ~ x - 1))
  ests <- NULL
  if (verbose) {
    ests <- optim(par = est,
                  ex_pl,
                  y = y,
                  x = x,
                  method = "BFGS",
                  hessian = TRUE,
                  control = list(fnscale = -1, trace = 6))
  } else {
    ests <- optim(par = est,
                  ex_pl,
                  y = y,
                  x = x,
                  method = "BFGS",
                  hessian = TRUE,
                  control = list(fnscale = -1, trace = 0))
  }
  return(ests)
}

# Log likelihood function calculations

# pseudolikelihood given theta#
ex_pl <- function(theta, y, x) {
  return(sum(log(ex_dtexp(y, x %*% theta))))
}

# The conditional density of each weight from a sample
ex_dtexp <- function(x, lambda) {
  den <- numeric(length(x))
  inds <- which(lambda != 0)
  den[inds] <- exp(x[inds] * lambda[inds]) /
    (1 / lambda[inds] * (exp(lambda[inds]) - 1))
  den[which(lambda == 0)] <- 1
  return(den)
}


# Convert an observed network to edge weight vectors x and y
ex_net2xy <- function(GERGM_Object) {
  net <- GERGM_Object@bounded.network
  y <- NULL
  x <- NULL
  nodes <- GERGM_Object@num_nodes
  if (GERGM_Object@include_diagonal) {
    for (i in 1:nodes) {
      for (j in (1:nodes)) {
        y <- c(y, net[i, j])
        x <- rbind(x, ex_dh(GERGM_Object, i, j))
      }
    }
  } else {
    for (i in 1:nodes) {
      for (j in (1:nodes)[-i]) {
        y <- c(y, net[i, j])
        x <- rbind(x, ex_dh(GERGM_Object, i, j))
      }
    }
  }

  return(list(y = y, x = x))
}


# din2star
ex_din2star <- function(i, j, net, node_set) {
  val <- 0
  if (i %in% node_set & j %in% node_set) {
    if (i == j) {
      others <- node_set[-j]
    } else {
      others <- node_set[-c(i, j)]
    }
    set <- cbind(others, j)
    val <- sum(net[set])
  }
  return(val)
}

# dout2star
ex_dout2star <- function(i, j, net, node_set) {
  val <- 0
  if (i %in% node_set & j %in% node_set) {
    if (i == j) {
      others <- node_set[-i]
    } else {
      others <- node_set[-c(i, j)]
    }
    set <- cbind(i, others)
    val <- sum((net[set]))
  }
  return(val)
}

# dedgeweight
ex_dedgeweight <- function(i, j, net, node_set) {
  val <- 0
  if (i %in% node_set & j %in% node_set) {
    if (i != j) {
      val <- 1
    }
  }
  return(val)
}

# drecip
ex_drecip <- function(i, j, net, node_set) {
  val <- 0
  if (i %in% node_set & j %in% node_set) {
    if (i != j) {
      val <- net[j, i]
    }
  }
  return(val)
}

# dctriads
ex_dctriads <- function(i, j, net, node_set) {
  val <- 0
  if (i %in% node_set & j %in% node_set) {
    if (i == j) {
      others <- node_set[-i]
    } else {
      others <- node_set[-c(i, j)]
    }
    triples <- cbind(i, j, others)
    val <- sum(net[triples[, c(2, 3)]] * net[triples[, c(3, 1)]])
  }
  return(val)
}

# dttriads
ex_dttriads <- function(i, j, net, node_set) {
  val <- 0
  if (i %in% node_set & j %in% node_set) {
    if (i == j) {
      others <- node_set[-i]
    } else {
      others <- node_set[-c(i, j)]
    }
    triples <- cbind(i, j, others)
    t2 <- sum(net[triples[, c(2, 3)]] * net[triples[, c(1, 3)]])
    t3 <- sum(net[triples[, c(3, 2)]] * net[triples[, c(3, 1)]])
    t4 <- sum(net[triples[, c(3, 2)]] * net[triples[, c(1, 3)]])
    val <- t2 + t3 + t4
  }
  return(val)
}

ex_ddiagonal <- function(i, j, net, node_set) {
  val <- 0
  if (i %in% node_set & j %in% node_set) {
    if (i == j) {
      val <- 1
    }
  }
  return(val)
}

# dh function weight w_{i,j} will be conditioned upon
# Calculate the marginal change in the network
ex_dh <- function(GERGM_Object, i, j, network = NULL) {
  # find out what stats we are using
  stats_to_use <- GERGM_Object@stats_to_use
  value <- rep(0, length(stats_to_use))

  # allow the user to pass in a network
  if (!is.null(network)) {
    net <- network
  } else {
    net <- GERGM_Object@bounded.network
  }

  # now loop through them to calculate d statistics
  for (k in 1:length(stats_to_use)) {
    # determine the node set
    node_set <- GERGM_Object@endogenous_statistic_node_sets[[k]]
    # if there was no reduced nodeset provided, then use all nodes
    if (length(node_set) == 1) {
      node_set <- 1:GERGM_Object@num_nodes
    }
    # now get the value depending on which statistic we are dealing with
    if (stats_to_use[k] == 1) {
      value[k] <- ex_dout2star(i, j, net, node_set)
    } else if (stats_to_use[k] == 2) {
      value[k] <- ex_din2star(i, j, net, node_set)
    } else if (stats_to_use[k] == 3) {
      value[k] <- ex_dctriads(i, j, net, node_set)
    } else if (stats_to_use[k] == 4) {
      value[k] <- ex_drecip(i, j, net, node_set)
    } else if (stats_to_use[k] == 5) {
      value[k] <- ex_dttriads(i, j, net, node_set)
    } else if (stats_to_use[k] == 6) {
      value[k] <- ex_dedgeweight(i, j, net, node_set)
    } else if (stats_to_use[k] == 7) {
      value[k] <- ex_ddiagonal(i, j, net, node_set)
    } else {
      stop("MPLE is not currently implemented for this statistic.")
    }
  }
  return(value)
}
matthewjdenny/GERGM documentation built on May 24, 2023, 1:28 a.m.