CRWCM

Maximum entropy model for the Conditionally Reciprocal Weighted Configuration Model (CRWCM).

The CRWCM (Di Vece, Pijpers & Garlaschelli (2023); RBCM+CRWCM in the NuMeTriS package) is the reciprocity-aware counterpart of the DCReM: a two-step null model for weighted, directed networks with continuous, positive weights that accounts for the different nature of the links weights sit on. The binary structure is supplied by an internally solved RBCM on the reciprocal degree sequences, giving the dyadic probabilities f⭢ᵢⱼ = p⭢ᵢⱼ (single link i→j) and f⭤ᵢⱼ = p⭤ᵢⱼ (reciprocated dyad). Conditional on the dyad state, the weights are exponential:

• a non-reciprocated link i→j carries a weight with rate θ⭢ᵢ + θ⭠ⱼ,
• a reciprocated pair carries two weights, wᵢⱼ with rate θ⭤ᵒᵢ + θ⭤ⁱⱼ and wⱼᵢ with rate θ⭤ᵒⱼ + θ⭤ⁱᵢ.

The 4N parameters θ = [θ⭢; θ⭠; θ⭤ᵒ; θ⭤ⁱ] constrain the four reciprocal strength sequences (s→, s←, s↔out, s↔in). The generalised (conditional) likelihood separates into two independent sub-problems — the non-reciprocated 2N system {θ⭢, θ⭠} and the reciprocated 2N system {θ⭤ᵒ, θ⭤ⁱ} — which are solved jointly as a single 4N problem for API uniformity (the Hessian is block-diagonal, so solvers benefit automatically).

Note: within a dyad the two weights wᵢⱼ and wⱼᵢ are correlated under the CRWCM (they are both non-zero only in the reciprocated state): Cov(wᵢⱼ, wⱼᵢ) ≠ 0, unlike for the DCReM. The layer-aware σₓ accounts for this.

CRWCM(G::T; kwargs...) where {T}
CRWCM(;d_out, d_in, d_rec, s_out, s_in, s_rec_out, s_rec_in, precision=Float64, kwargs...)

Constructor function for the CRWCM type.

By default the three reciprocal degree sequences (see nonreciprocated_outdegree, nonreciprocated_indegree, reciprocated_degree) and the four reciprocal strength sequences (see nonreciprocated_outstrength, nonreciprocated_instrength, reciprocated_outstrength, reciprocated_instrength) are computed from the weighted directed graph G. You can also pass all seven sequences directly as keyword arguments.

The CRWCM allows continuous, positive weights. Because the weights are strictly positive, a node has a zero strength in a channel iff its degree in that channel is zero; the constructor enforces this consistency (DomainError otherwise).

Examples

# generating a model from a weighted directed graph
julia> model = CRWCM(MaxEntropyGraphs.rhesus_macaques())
CRWCM{SimpleWeightedGraphs.SimpleWeightedDiGraph{Int64, Float64}, Float64} (16 vertices, 15 unique degree triples, 0.94 compression ratio)

solve_model!(m::CRWCM; kwargs...)

Compute the likelihood maximising parameters of the CRWCM model m. This is a two-step process: first the binary (conditional) RBCM layer is solved on the reciprocal degree sequences, then the weighted CRWCM layer is solved on the four reciprocal strength sequences.

By default both layers are computed with the fixed-point method. Channels whose strength constraint is zero have an undetermined parameter; after the solve those entries of θ are pinned to +Inf (an infinite rate, i.e. an exactly zero weight — consistent with the dyadic probability being zero).

Arguments (weighted CRWCM layer)

• method::Symbol: solution method, :fixedpoint (default) or any of :BFGS, :LBFGS and :Newton.
• initial::Symbol: initial guess, :strengths (default), :strengths_minor or :random.
• AD_method::Symbol: autodiff method, any of :AutoZygote, :AutoReverseDiff, :AutoForwardDiff and :AutoFiniteDiff (defaults to :AutoZygote).
• analytical_gradient::Bool: use the analytical gradient instead of autodiff (defaults to false).
• store_adjacency::Bool: cache the (binary) expected adjacency matrix m.Ĝ (defaults to false).

Arguments (binary conditional RBCM layer)

• method_conditional::Symbol: solution method for the binary layer (defaults to :fixedpoint).
• initial_conditional::Symbol: initial guess for the binary layer (defaults to :degrees).
• AD_method_conditional::Symbol: autodiff method for the binary layer (defaults to :AutoZygote).
• analytical_gradient_conditional::Bool: analytical gradient for the binary layer (defaults to false).

Common settings

• maxiters::Int: maximum number of iterations (defaults to 1000).
• verbose::Bool: show log messages (defaults to false).
• ftol::Union{Real, Nothing}: tolerance for the fixedpoint method (defaults to nothing, i.e. 1e-8), applied to the weighted layer and to the binary layer when either is solved with :fixedpoint (passing it when neither layer uses it warns). On the weighted layer it is a relative strength tolerance (see below); on the binary layer it bounds the fixed-point increment in parameter space.
• abstol, reltol: absolute/relative tolerances for the optimisation methods (default nothing).
• g_tol::Union{Number, Nothing}: gradient tolerance for the gradient-based methods (maps to Optim's g_abstol, default nothing). The gradient of the weighted layer is its constraint residual in strength units, but g_abstol is a stopping criterion rather than a guarantee: Optim can also stop on its function or parameter convergence checks and report success without the gradient ever reaching g_tol. Verify what was actually achieved with constraint_residual.

Examples

# default use
julia> model = CRWCM(MaxEntropyGraphs.rhesus_macaques());

julia> solve_model!(model);

initial_guess(m::CRWCM; method::Symbol=:strengths)

Compute an initial guess θ₀ = [θ⭢₀; θ⭠₀; θ⭤ᵒ₀; θ⭤ⁱ₀] for the maximum likelihood parameters of the weighted (CRWCM) layer of the model m.

The CRWCM parameters θ are the direct rate parameters, so every method returns non-negative, feasible values (dead channels get a zero placeholder; they are pinned to +Inf after the solve):

• :strengths (default): θᵢ = 𝟙[sᵢ > 0] / Σ_channels Σⱼ sⱼ per channel.
• :strengths_minor: θᵢ = 𝟙[sᵢ > 0] / (sᵢ + 1) per channel.
• :random: random values drawn from $U(0,1)$ (zeroed on dead channels).

rand(m::CRWCM; precomputed=false, rng=Random.default_rng())

Generate a random weighted directed graph from the CRWCM model m.

Each dyad (i,j) is drawn from its four-state RBCM distribution (single link i→j, single link j→i, reciprocated pair, absent); the realised links then get continuous exponential weights with the rates of the corresponding channel (a reciprocated pair gets two conditionally independent weights).

Arguments

• precomputed::Bool: not implemented for the CRWCM (the parameters are always used to generate the graph on the fly).
• rng::AbstractRNG: random number generator to use (defaults to Random.default_rng()).

Examples

julia> model = CRWCM(MaxEntropyGraphs.rhesus_macaques()); # generate a CRWCM model

julia> solve_model!(model); # compute the maximum likelihood parameters

julia> sample = rand(model); # sample a random weighted directed graph

julia> typeof(sample)
SimpleWeightedGraphs.SimpleWeightedDiGraph{Int64, Float64}

rand(m::CRWCM, n::Int; precomputed=false, rng=Random.default_rng())

Generate n random weighted directed graphs from the CRWCM model m. If multithreading is available, the graphs are generated in parallel; per-sample seeds are drawn from rng so the result is reproducible and independent of the thread schedule.

AIC(m::CRWCM)

Compute the Akaike Information Criterion (AIC) for the CRWCM model m. The parameters of the model must be computed beforehand. The (conditional) CRWCM has 4N parameters (four rate parameters $θ$ per node; the binary layer's dyadic probabilities are fixed inputs of the conditional likelihood).

See also AICc, L_CRWCM.

AICc(m::CRWCM)

Compute the corrected Akaike Information Criterion (AICc) for the CRWCM model m.

See also AIC, L_CRWCM.

BIC(m::CRWCM)

Compute the Bayesian Information Criterion (BIC) for the CRWCM model m.

See also AIC, L_CRWCM.

Return the number of vertices in the CRWCM network

L_CRWCM(θ::AbstractVector, s_out, s_in, s_rec_out, s_rec_in, nz_out, nz_in, nz_rec, x, y, z)

Compute the (generalised, conditional) log-likelihood of the CRWCM model given the weighted parameters θ = [θ⭢; θ⭠; θ⭤ᵒ; θ⭤ⁱ], the four reciprocal strength sequences, their non-zero index sets and the per-node binary (RBCM) fitnesses x, y, z, computing the dyadic probabilities on the fly.

The log-likelihood is

$\mathcal{G} = -\sum_i \left[ θ^{→}_i s^{→}_i + θ^{←}_i s^{←}_i + θ^{↔,o}_i s^{↔,out}_i + θ^{↔,i}_i s^{↔,in}_i \right] + \sum_{i≠j} \left[ f^{→}_{ij} \log(θ^{→}_i + θ^{←}_j) + f^{↔}_{ij} \log(θ^{↔,o}_i + θ^{↔,i}_j) \right]$

(Di Vece et al. (2023), Methods). It is block-separable into the non-reciprocated system {θ⭢, θ⭠} and the reciprocated system {θ⭤ᵒ, θ⭤ⁱ}. The log is domain-guarded (returns NaN outside θ + θ > 0) so the line search rejects infeasible steps. Channels whose strength constraint is zero have f = 0 in every pair term; they are excluded from the sums (their parameter is undetermined and pinned to +Inf after the solve).

This method is used when the dyadic probability matrices are not precomputed.

L_CRWCM(θ::AbstractVector, s_out, s_in, s_rec_out, s_rec_in, nz_out, nz_in, nz_rec, f⭢::AbstractMatrix, f⭤::AbstractMatrix)

Compute the log-likelihood of the CRWCM model, using the precomputed dyadic probability matrices f⭢ and f⭤. See also L_CRWCM(::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector, ::Vector, ::Vector, ::Vector, ::AbstractVector, ::AbstractVector, ::AbstractVector).

L_CRWCM(m::CRWCM)

Return the log-likelihood of the CRWCM model m based on the computed maximum likelihood parameters.

∇L_CRWCM!(∇L::AbstractVector, θ::AbstractVector, s_out, s_in, s_rec_out, s_rec_in, nz_out, nz_in, nz_rec, x, y, z)

Compute the gradient of the log-likelihood of the CRWCM model in a non-allocating manner, computing the dyadic probabilities on the fly.

∂L/∂θ⭢ᵢ = -s→ᵢ + Σⱼ≠ᵢ f⭢ᵢⱼ/(θ⭢ᵢ + θ⭠ⱼ), ∂L/∂θ⭠ᵢ = -s←ᵢ + Σⱼ≠ᵢ f⭢ⱼᵢ/(θ⭠ᵢ + θ⭢ⱼ) (note f⭠ᵢⱼ = f⭢ⱼᵢ), ∂L/∂θ⭤ᵒᵢ = -s↔outᵢ + Σⱼ≠ᵢ f⭤ᵢⱼ/(θ⭤ᵒᵢ + θ⭤ⁱⱼ) and ∂L/∂θ⭤ⁱᵢ = -s↔inᵢ + Σⱼ≠ᵢ f⭤ᵢⱼ/(θ⭤ⁱᵢ + θ⭤ᵒⱼ) (f⭤ is symmetric). Rows of dead channels (zero strength) are left at zero.

See also ∇L_CRWCM_minus!.

∇L_CRWCM_minus!(∇L::AbstractVector, θ::AbstractVector, s_out, s_in, s_rec_out, s_rec_in, nz_out, nz_in, nz_rec, x, y, z)

Compute minus the gradient of the log-likelihood of the CRWCM model (used for the minimisation carried out by Optimization.jl), computing the dyadic probabilities on the fly. Non-allocating.

See also ∇L_CRWCM!.

CRWCM_iter!(θ::AbstractVector, s_out, s_in, s_rec_out, s_rec_in, nz_out, nz_in, nz_rec, x, y, z, G::AbstractVector)

Compute the next fixed-point iteration for the CRWCM model, computing the dyadic probabilities on the fly. The pre-allocated buffer G is updated in place. The consistency equations mirror the CReM family, e.g. $θ^{→}_i = \left( \sum_{j≠i} f^{→}_{ij}/(1 + θ^{←}_j/θ^{→}_i) \right)/s^{→}_i$ for the four channels. Rows of dead channels (zero strength) are fixed at the current parameter value (residual zero).

set_xᵣ!(m::CRWCM)

Set the value of xᵣ to exp(-αᵣ[1:d_unique]) for the (binary layer of the) CRWCM model m.

set_yᵣ!(m::CRWCM)

Set the value of yᵣ to exp(-αᵣ[d_unique+1:2*d_unique]) for the (binary layer of the) CRWCM model m.

set_zᵣ!(m::CRWCM)

Set the value of zᵣ to exp(-αᵣ[2*d_unique+1:end]) for the (binary layer of the) CRWCM model m.

Ĝ(m::CRWCM)

Compute the expected (binary) adjacency matrix for the CRWCM model m from its RBCM layer ($⟨a_{ij}⟩ = p^{→}_{ij} + p^{↔}_{ij}$).

Note: under the RBCM layer the entries aᵢⱼ and aⱼᵢ are correlated within a dyad, so Ĝ alone does not characterise the dyadic joint distribution. The expected weights can be computed with Ŵ.

set_Ĝ!(m::CRWCM)

Set the expected (binary) adjacency matrix for the CRWCM model m.

Ŵ(m::CRWCM)

Compute the expected (unconditional) weighted adjacency matrix for the CRWCM model m, i.e. $⟨w_{ij}⟩ = f^{→}_{ij}/(θ^{→}_i + θ^{←}_j) + f^{↔}_{ij}/(θ^{↔,o}_i + θ^{↔,i}_j)$, so that sum(Ŵ(m), dims=2) ≈ s→ + s↔out (the total out-strength) and sum(Ŵ(m), dims=1) ≈ s← + s↔in.

set_Ŵ!(m::CRWCM)

Set the expected weighted adjacency matrix for the CRWCM model m.

σˣ(m::CRWCM)

Compute the standard deviation for the elements of the (binary) adjacency matrix for the CRWCM model m ($σ[a_{ij}] = \sqrt{⟨a_{ij}⟩(1-⟨a_{ij}⟩)}$, from the RBCM layer).

Note: the within-dyad correlation between aᵢⱼ and aⱼᵢ is captured separately and accounted for by σₓ. Read as "sigma star".

set_σ!(m::CRWCM)

Set the standard deviation for the elements of the (binary) adjacency matrix for the CRWCM model m.

σʷ(m::CRWCM)

Compute the standard deviation for the elements of the weighted adjacency matrix for the CRWCM model m. The (unconditional) weight wᵢⱼ is a three-component mixture (zero; exponential with rate θ⭢ᵢ + θ⭠ⱼ with probability f⭢ᵢⱼ; exponential with rate θ⭤ᵒᵢ + θ⭤ⁱⱼ with probability f⭤ᵢⱼ), so $⟨w_{ij}^2⟩ = 2f^{→}_{ij}/(θ^{→}_i + θ^{←}_j)^2 + 2f^{↔}_{ij}/(θ^{↔,o}_i + θ^{↔,i}_j)^2$ and $Var(w_{ij}) = ⟨w_{ij}^2⟩ - ⟨w_{ij}⟩^2$.

Note: within a dyad the weights wᵢⱼ and wⱼᵢ are correlated (see MaxEntropyGraphs._covʷ); σₓ(m, X; layer=:weighted) accounts for this.

set_σʷ!(m::CRWCM)

Set the standard deviation for the elements of the weighted adjacency matrix for the CRWCM model m.

_cov_dyads(m::CRWCM)

Compute the within-dyad covariance matrix of the binary layer of the CRWCM model m ($Cov(a_{ij}, a_{ji}) = p^{↔}_{ij} - ⟨a_{ij}⟩⟨a_{ji}⟩$, from the RBCM layer). Symmetric, zero diagonal.

_covʷ(m::CRWCM)

Compute the within-dyad covariance matrix of the weighted layer of the CRWCM model m. The weights wᵢⱼ and wⱼᵢ are both non-zero only in the reciprocated dyad state (where they are conditionally independent exponentials), so

$Cov(w_{ij}, w_{ji}) = \frac{f^{↔}_{ij}}{(θ^{↔,o}_i + θ^{↔,i}_j)(θ^{↔,o}_j + θ^{↔,i}_i)} - ⟨w_{ij}⟩⟨w_{ji}⟩$.

Symmetric, zero diagonal. This term is what distinguishes the CRWCM's weighted uncertainty from the DCReM's (where it vanishes).

precision(m::CRWCM)

Determine the compute precision of the CRWCM model m.

A(m::CRWCM, i::Int, j::Int)

Return the expected value of the (binary) adjacency matrix for the CRWCM model m at the node pair (i,j) ($⟨a_{ij}⟩ = p^{→}_{ij} + p^{↔}_{ij}$, from the RBCM layer).

❗ For performance reasons, the function does not check:

• if the node pair is valid.
• if the conditional parameters of the model have been computed.

p⭢(m::CRWCM, i::Int, j::Int)

Return the probability that the dyad (i,j) holds a single non-reciprocated link i→j under the binary (RBCM) layer of the CRWCM model m.

❗ For performance reasons, the function does not check if the node pair is valid or if the conditional parameters have been computed.

p⭠(m::CRWCM, i::Int, j::Int)

Return the probability that the dyad (i,j) holds a single non-reciprocated link j→i under the binary (RBCM) layer of the CRWCM model m.

p⭤(m::CRWCM, i::Int, j::Int)

Return the probability that the dyad (i,j) holds a reciprocated link pair under the binary (RBCM) layer of the CRWCM model m.

❗ For performance reasons, the function does not check if the node pair is valid or if the conditional parameters have been computed.

σₓ(m::CRWCM, X::Function; layer::Symbol=:binary, gradient_method::Symbol=:ReverseDiff)

Compute the standard deviation of metric X for the CRWCM model m via error propagation (the delta method of Squartini & Garlaschelli (2011), Eq. B.16, including the within-dyad covariance cross-terms).

Arguments

• layer::Symbol:
  • :binary (default): X is a function of the adjacency matrix; the gradient is evaluated at m.Ĝ, weighted by m.σ, plus the binary within-dyad covariance Cov(aᵢⱼ, aⱼᵢ) cross-term (requires set_Ĝ! and set_σ!).
  • :weighted: X is a function of the weight matrix; the gradient is evaluated at m.Ŵ, weighted by m.σʷ, plus the weighted within-dyad covariance Cov(wᵢⱼ, wⱼᵢ) cross-term (requires set_Ŵ! and set_σʷ!). This is the layer to use for weighted metrics such as the triadic fluxes.
• gradient_method::Symbol: :ForwardDiff, :ReverseDiff (default) or :Zygote.