log-likelihood

In the case of a single network, the model’s log-likelihood is given by

\[ \text{log}\left(\text{P}\left(\cdot\right)\right) = - \text{log}\left( % \sum_{\mathbf{y}'\in\mathbf{Y}}\text{exp}\left\{{\theta}^{\text{t}}\Delta g\left(\mathbf{y}'\right)\right\} % \right) \] In general, we can reduce the computational complexity of this calculations by looking at the isomorphic sufficient statistics, this is, group up elements based on the set of unique vectors of sufficient statistics:

\[ - \text{log}\left( % \sum_{\mathbf{s}\in\mathcal{S}}\mathbf{w}_\mathbf{s}\text{exp}\left\{{\theta}^{\text{t}}\mathbf{s}\right\} % \right) \]

Where \(\mathcal{S}\) is the support of the sufficient statistics under \(\mathcal{Y}\), \(\mathbf{s}\in\mathcal{S}\) is one of its realizations, and \(\mathbf{w}_\mathbf{s}\equiv|\left\{\mathbf{y}\in\mathcal{Y}: \Delta{}g\left(\mathbf{y}\right) = \mathbf{s}\right\}|\) is the number of networks in \(\mathcal{Y}\) which centered sufficient statistics equal \(\mathbf{s}\). Furthermore, we can write this in matrix form:

\[ - % \text{log}\left( % \mathbf{W}\times \text{exp}\left\{\mathbf{S}\times \theta\right\} % \right) \]

Where \(\mathbf{W}\equiv\{\mathbf{w}_\mathbf{s}\}_{\mathbf{s}\in\mathbf{S}}\) is a row vector of, length \(w\) and \(\mathbf{S}\) is a matrix of size \(w\times k\). The log-likelihood of a pooled model is simply the sum of the individual ones.

Gradient

The partial derivative of the log-likelihood with respect to \(j\)-th parameter equals:

\[ \frac{\delta}{\delta\theta_j}\text{log}\left(\text{P}\left(\cdot\right)\right) = % - \frac{ % \sum_{\mathbf{s}\in\mathcal{S}}\mathbf{w}_\mathbf{s}\mathbf{s}_j\text{exp}\left\{{\theta}^{\text{t}}\mathbf{s}\right\} % }{ % \sum_{\mathbf{s}\in\mathcal{S}}\mathbf{w}_\mathbf{s}\text{exp}\left\{{\theta}^{\text{t}}\mathbf{s}\right\} % }, \quad\forall j \]

We can also write this using matrix notation as follows:

\[ \nabla\text{log}\left(\text{P}\left(\cdot\right)\right) = % - % {\mathbf{S}}^{\text{t}}\times \left[ \mathbf{W}\circ \text{exp}\left\{\mathbf{S}\times \theta\right\} \right]/ % \lambda({\theta}) \]

Where \(\circ\) is the element-wise product and \(\lambda({\theta}) = \mathbf{W}\times \text{exp}\left\{\mathbf{S}\times \theta\right\}\).

Hessian

In the case of the hessian, each \((j, l)\) element of, \(\frac{\delta^2}{\delta\theta_k\delta\theta_u}\text{log}\left(\text{P}\left(\cdot\right)\right)\), can be computed as:

\[ \frac{% -\left(\sum_{\mathbf{s}'\in\mathcal{S}}\mathbf{s}_j'\mathbf{s}_l'\mathbf{w}_\mathbf{s}\text{exp}\left\{{\theta}^{\text{t}} \mathbf{s}\right\}\right) }{% \lambda(\theta)% } + \frac{% \left(\sum_{\mathbf{s}'\in\mathcal{S}}\mathbf{s}_j'\mathbf{w}_\mathbf{s}\text{exp}\left\{{\theta}^{\text{t}} \mathbf{s}\right\}\right) \left(\sum_{\mathbf{s}'\in\mathcal{S}}\mathbf{s}_l'\mathbf{w}_\mathbf{s}\text{exp}\left\{{\theta}^{\text{t}} \mathbf{s}\right\}\right) }{% \lambda(\theta)^2% } \] Where \(\mathbf{s}_j\) as the \(j\)-th element of the vector \(\mathbf{s}\). Once again, we can simplify this using matrix notation:

\[ \frac{% -\mathbf{W}\times\left[\mathbf{S}_j\circ\mathbf{S}_l \circ \text{exp}\left\{\mathbf{S}\times\theta\right\}\right]% }{% \lambda(\theta)% } + \frac{% \left(\mathbf{W}\times\left[\mathbf{S}_j \circ \text{exp}\left\{\mathbf{S}\times\theta\right\}\right]\right) \left(\mathbf{W}\times\left[\mathbf{S}_l \circ \text{exp}\left\{\mathbf{S}\times\theta\right\}\right]\right) }{% \lambda(\theta)^2% } \]

Where \(\mathbf{S}_j\) is the \(j\)-th column of the matrix \(\mathbf{S}\).

Limiting values

In the case that the MLE does not exists, i.e. \(\theta_k\to\pm\infty\), which occurs when the observed sufficient statistics are not in the interior of the space, for example, a fully connected network, the limit of the log-likelihood, the gradient, and hessian are finite. This is relevant as some special care needs to be taken when dealing with these cases.

While in general models in which the MLEs diverge may seem uninteresting, the fact that ERGMs describe a discrete and not continuous random variable makes this type of event easier to observed when compared to other families of distributions. Furthermore, in the case of methods such as bootstrapping, forward/backward selection, or other classes of algorithms that involve some level of automatization during the model fitting process, it is important to know how to correctly deal with the non-existence of MLEs. Fortunately, as we will show, the log-likelihood and its derivatives are well defined in such cases.

The main principle here is what happens in the argument of the exponential function that populates these functions. Depending on what is \(\mathbf{s}_k\), two cases follow:

1. The observe value of the \(k\)-th sufficient statistic is in the lower bound \(\mathbf{s}_k = \min_{\mathbf{s}'\in\mathcal{S}}\mathbf{s}_{k}'\) This happens, for example, when trying to fit a model with triangles where there are no observed triangles. In this case, with \(\mathbf{s}_k' \geq 0\) for all cases, the theoretical MLE for \(\theta_k\) goes to \(-\infty\), thus, in the limit, the expression \[ \lim_{\theta_k\to-\infty}\text{exp}\left\{{\theta}^{\text{t}}\mathbf{s}'\right\} = \left\{\begin{array}{ll} % \text{exp}\left\{\sum_{j\neq k}\mathbf{s}_j'\theta_j\right\} &\quad\text{if }\mathbf{s}_k' = 0 \\ 0 & \quad\text{if }\mathbf{s}_k' > 0 \end{array}\right. \]

2. The observe value of the \(k\)-th sufficient statistic is in the upper bound \(\mathbf{s}_k = \min_{\mathbf{s}'\in\mathcal{S}}\mathbf{s}_k'\) The most common case is when the sufficient statistic is satturated, for example, in a fully connected graph, where the MLE goes to infinity, \(\theta_k\to+\infty\). Similar to before, since \(\mathbf{s}_k' \leq 0\) for all cases, in the limit the previous expression is well defined: \[ \lim_{\theta_k\to+\infty}\text{exp}\left\{{\theta}^{\text{t}}\mathbf{s}'\right\} = \left\{\begin{array}{ll} % \text{exp}\left\{\sum_{j\neq k}\mathbf{s}_j'\theta_j\right\} &\quad\text{if }\mathbf{s}_k' = 0 \\ 0 & \quad\text{if }\mathbf{s}_k' < 0 \end{array}\right. \]

The two previous points can be interpreted as only graphs whose \(k\)-th sufficient statistic matches that of the observed data influence the model. Therefore, a key to compute the limiting values of the log-likelihood, gradient, and hessian will be on partitioning the summations over \(\mathbf{s}'\in\mathcal{S}\) into two different sets as follows:

\[\begin{align} \mathcal{S}_0 & \equiv \{\mathbf{s}\in\mathcal{S}: \mathbf{s}_u = 0, \forall u\in U\} \\ \mathcal{S}_1 & \equiv \mathcal{S}\setminus\mathcal{S}_0 \end{align}\]

Where \(U\) is the set of observed sufficient statistics that are on the boundary and thus have an MLE that diverges to \(\pm\infty\). In this partition \(\mathcal{S}_0\) contains all the realizations of sufficient statistics for which the statistics in \(U\) are equal to the observed ones. With this definition in hand, we will now show what is the asymptotic behavior of the log-likelihood, gradient, and hessian when one or more observed sufficient statistics are on not in the interior of the support.

library(ergmito)
data(fivenets)
(ans <- ergmito(fivenets ~ edges + triadcensus(15)))

## Warning: The observed statistics (target.statistics) are near or at the
## boundary of its support, i.e. the Maximum Likelihood Estimates maynot exist or
## be hard to be estimated. In particular, the statistic(s) "edges",
## "triadcensus.300".

## 
## ERGMito estimates
## note: A subset of the parameters estimates was replaced with +/-Inf,  and the Hessian matrix is not p.s.d.
## 
## Call:
## ergmito(model = fivenets ~ edges + triadcensus(15))
## 
##  Coefficients:
##           edges  triadcensus.300  
##         -0.6851             -Inf

# Approach using the limiting value
l <- with(ans$formulae, {
  sapply(1:nnets(ans), function(i) {
    # Preparing suff stats
    S <- t(t(stats_statmat[[i]]) - target_stats[i, ])
    W <- stats_weights[[i]]
    theta <- coef(ans)["edges"]
    
    # Index of members of S0
    s0idx <- which(S[,2] == 0)
    
    - log(sum(W[s0idx] * exp(S[s0idx,1] * theta)))
  })
})
sum(l)

## [1] -38.16338

# Which should be equivalent to the approach with large negative number
ans$formulae$loglik(c(coef(ans)["edges"], -1e100))

## [1] -38.16338

g <- with(ans$formulae, {
  sapply(1:nnets(ans), function(i) {
    # Preparing suff stats
    S <- t(t(stats_statmat[[i]]) - target_stats[i, ])
    W <- stats_weights[[i]]
    theta <- coef(ans)["edges"]
    
    # Index of members of S0
    s0idx <- which(S[,2] == 0)
    
    - sum(W[s0idx] * S[s0idx,1] * exp(S[s0idx,1] * theta))/
      sum(W[s0idx] * exp(S[s0idx,1] * theta))
  })
})
sum(g)

## [1] 0.002926372

# Which is equivalent to
ans$formulae$grad(c(coef(ans)["edges"], -1e100))

##             [,1]
## [1,] 0.002926372
## [2,] 0.000000000

h <- with(ans$formulae, {
  sapply(1:nnets(ans), function(i) {
    # Preparing suff stats
    S <- t(t(stats_statmat[[i]]) - target_stats[i, ])
    W <- stats_weights[[i]]
    theta <- coef(ans)["edges"]
    
    # Index of members of S0
    s0idx <- which(S[,2] == 0)
    
    # First part
    - sum(W[s0idx] * S[s0idx,1]^2 * exp(S[s0idx,1] * theta))/
      sum(W[s0idx] * exp(S[s0idx,1] * theta)) + 
      # Second bit
      sum(W[s0idx] * S[s0idx,1] * exp(S[s0idx,1] * theta)) ^ 2/
      sum(W[s0idx] * exp(S[s0idx,1] * theta))^2
  })
})
sum(h)

## [1] -12.97199

ans$formulae$hess(c(coef(ans)["edges"], -1e100))

##           [,1] [,2]
## [1,] -12.97199    0
## [2,]   0.00000    0