inst/algorithm.md
In scNPM: Single-Cell Network Propensity Matching

scNPM algorithm: staged inference and BFDR-selected networks

The state includes latent expression, cluster memberships, mixture weights, component means, precisions, spike/slab allocations and observation parameters. Nonzero expression fixes its latent value to log X; missing latent values at zero observations are sampled using fixed-step HMC with Metropolis correction. Component labels are aligned to fixed initial mean anchors. After freezing, paired mode permutes whole matched pairs jointly; independent mode aligns each condition separately. Both conditions use the same 2K network references and any measured cell covariates. Condition membership (T or C) is observed throughout. By default, each candidate cluster pair has its own balanced ridge logistic regression; scaling is learned within each training fold. Global-score mode retains the pooled regression and scaling used in earlier versions.

Matching and edge feedback operate during adaptation. All clusters must be nonempty and the complete mapping must repeat for freeze_patience sweeps, no earlier than adapt_min. At adapt_max, the default consensus policy freezes the most frequent complete map among the last freeze_window valid states (default 100). Tied maps are ranked by maximum cost, then total cost, computed from the mean cost matrix over that window, then lexicographically. This is a budget-based decision: freeze_reason='budget_consensus' and matching_stable=FALSE distinguish it from a stable freeze. Invalid or empty-cluster states do not enter the window. With no valid state at all, the sampler still reports an error. freeze_policy='error' restores the strict stability requirement. Starting with the next sweep after freezing, adaptation feedback stops and the selected final prior is activated (paired by default in 1.1.0). With L=num_iter and S=num_save, retained sweeps are L-S+1 through L. A freeze at F leaves L-S-F final-model warmup sweeps, required to be at least final_burnin. The default minimum retained sample count for BFDR selection is 200.

Default budgets are 5000 total sweeps, 1000 retained sweeps, adapt_max=3000, freeze_patience=5 and at least 1000 post-freeze warmup sweeps. Short custom runs derive smaller default budgets satisfying the same retention constraints. Neither freezing rule by itself establishes MCMC convergence or good balance.

The diagnostic criterion is augmented likelihood plus a network-prior kernel minus the matching SMD sum. It is not a jointly optimized objective, and stochastic updates need not increase it.

# These settings are also recorded in fit$paired_prior and fit$provenance.
# Add them to the run_scNPM call above:
# network_prior = "paired",
# paired_discordance = 0.75,
# shared_edge_correlation = 0.995,
# ppi_reference = "control"

For each frozen pair and gene pair, S=(gamma_T,gamma_C) takes states 00 (neither), 10 (treatment only), 01 (control only), 11 (both). Let b=ppi_prior_low and d=paired_discordanceb(1-b). Base state factors are (1-b-d,d,d,b-d). For an edge with PPI weight pi, define R={pi/(1-pi)}/{b/(1-b)}. ppi_reference="control" multiplies the 01 and 11 factors by R. "treatment" multiplies 10 and 11. "both" multiplies 10/01 by R and 11 by R^2. Normalize the four factors before using them. The default control setting matches control-derived PPI experiments; use both when PPI evidence is condition-neutral. A missing PPI edge has R=1, not a forced absence.

The Gaussian coefficient kernels have covariance matrices

$$\Sigma_{00}=\operatorname{diag}(v_0^2,v_0^2),\quad \Sigma_{10}=\operatorname{diag}(v_1^2,v_0^2),\quad \Sigma_{01}=\operatorname{diag}(v_0^2,v_1^2),$$ $$\Sigma_{11}=v_1^2\begin{pmatrix}1&r\r&1\end{pmatrix}.$$

Here r=shared_edge_correlation; the default 0.995 strongly favors similar signed coefficients for shared edges. It is not a difference-rate parameter. All factors define one globally normalized prior with positive-definite precision constraints; the four factors are not unconditional state marginals. No gene modules, module learning, latent-expression Gibbs sampler or temperature ladder is included. The existing observation-model HMC remains unchanged. The paired graph kernel integrates one treatment/control coefficient pair when drawing its state and immediately restores both coefficients.

BFDR levels and tie handling are unchanged. A matched edge's p10/p01 now comes from this joint conditional, and group probabilities are p10+p11 or p01+p11. Matching freeze is not a convergence certificate. Performance from private development runs using other observation samplers is not release validation for this package, and no real-data precision or recall is guaranteed.

For each condition d and subtype k, use all H=G*(G-1)/2 unordered gene pairs. Let p[e] be the retained mean of the final-model conditional slab probability. It is recovered from the saved joint quantities as p10+p11 for treatment or p01+p11 for control. The marginal is independent of which counterpart is used; select_group_edges checks this numerical consistency. It uses the frozen-pair slice, as in the manuscript's group-network analysis, without resampling.

For E(c)={e:p[e]>=c}, estimate BFDR as mean(1-p[e]) over E(c). Sort probabilities in decreasing order and evaluate complete tie-block endpoints. Keep the largest prefix with estimated BFDR <= group_fdr_level (default 0.05). Its smallest selected probability is c; c=Inf denotes an empty set. This is an inclusive threshold, unlike the strict contrast threshold for differential edges. For fewer than min_draws retained draws, keep the curve but return an empty set.

For any reported set, conditional on the data, linearity of expectation gives E[FDP|X] = sum_{e in E} Pr(z[e]=0|X) / max(1,|E|). Independence across edges is not needed. With estimated probabilities, if their maximum absolute error is epsilon, the selected posterior expected FDP is at most alpha+epsilon. The bound concerns spike/slab allocations under the selected final posterior.

Every condition/subtype and gene-panel fit is a separate family. Select before restriction to shared genes; reuse the same saved mask for figures, exports and sensitivity analysis. A union across clusters is a summary of those selections and does not inherit a new per-union BFDR bound. Empirical allocation frequencies remain in edge_treat_est/edge_ctrl_est; group_fdr contains the Rao-Blackwell probabilities and BFDR masks. Default plotting requires these saved selections; selection="frequency" explicitly enables descriptive plots.

For each candidate (T_k,C_l), fit logistic regression to its cells using their known condition indicators. Equal total weight is assigned to T and C. The non-intercept ridge penalty is 0.001 times the total observation weight; the intercept receives a numerical penalty of 1e-8. A maximum of five folds is allocated by cell index, without consuming the sampler RNG. Each prediction uses a fit and predictor scaling from the other folds. For very small pairs where a training fold contains only one condition, an in-sample fit supplies that fold's predictions; such pairs should be inspected separately. The network features themselves are estimated from the full current state, so this is cross-fitting of the score regression, not the entire network fit.

The candidate cost is the absolute difference between mean predicted scores, divided by sqrt(max((var_T+var_C)/2,1e-12)). Each candidate thus uses a distinct conditional score. Scores across different pairs are not a single population propensity score. The known T/C labels are responses, not unknown cell types. matching_score='global' retains the original pooled propensity-score costs. Low assignment cost is relative to competing pairs; inspect the absolute SMDs and reference-feature balance even when the correspondence is stable.

At each unfrozen matching call, construct D[k,c] from the cluster means and sample variances of the current cell propensity scores. Minimize the ordered pair (max_k D[k,phi(k)], sum_k D[k,phi(k)]) over complete bijections. For K<=8, examine permutations in ascending order and update only for a strict improvement in one of the two objectives. K=4 examines 24 possibilities. For K>8, binary-search the sorted finite cost levels for a perfect matching, then use minimum-sum assignment restricted to that minimum bottleneck. Normalize allowed costs to [0,1] and use K+1 for forbidden entries, so a feasible allowed assignment always beats one containing a forbidden edge. Ascending solver traversal deterministically handles equal costs. Computation uses floating-point costs; no SMD upper bound or user-chosen near-tie tolerance is imposed.

If no all-finite permutation exists, every complete permutation has infinite maximum and total cost. Return the identity by the fixed tie rule. Missing SMDs from empty clusters are kept missing in diagnostics and treated as infinite for provisional assignment. Empty clusters prevent freezing. Finite but large SMDs are not discarded. Stable nonempty complete maps freeze after the configured patience interval; feedback is disabled in subsequent iterations and final warmup precedes retention. Label alignment still uses minimum squared distance to fixed anchors, which is a separate assignment problem from matching T to C clusters.

The optimizer's guarantee refers to each current cost matrix. A frozen map is not reoptimized at the final state. Its final maximum SMD need not remain optimal, nor must the maximum SMD decrease from one iteration to the next.

verbose=TRUE reports completed iterations, sampling stage, elapsed time, approximate ETA, freeze reason and retained count. The time interval defaults to 60 seconds, with additional reports at start, freezing, first retention and completion. Reports occur between sweeps; a long sweep can delay a report. progress_file optionally appends CSV records with run IDs, including when console messages are suppressed. No file is created by default. Reporting does not change the RNG state or retained draws.

Let e=(k,j,t), j<t. Using final draws, average the joint four-state conditional to estimate p10[e] and p01[e] for a frozen paired-prior counterpart. For independent priors, use within-draw conditional products. For alternative counterparts, use products of the paired-model conditional marginals within each draw; this does not change the prior or refit an alternative mapping. Set Delta[e]=p10[e]-p01[e], direction[e]=sign(Delta[e]), and ell[e]=1-max(p10[e],p01[e]). The direction asserts unequal spike/slab allocation states in the indicated groups, not gene absence or directed regulation.

For S(tau)={e:|Delta[e]|>tau}, define

BFDR(tau) = sum_{e in S(tau)} ell[e] / |S(tau)|.

All KG(G-1)/2 candidates are a single family, with at most one direction per candidate. For any data-determined nonempty set S, conditional on the data and the selected final posterior,

E[FDP(S)|X] = (1/|S|) * sum_{e in S} Pr(reported direction is false|X).

This is linearity of conditional expectation; edge independence is unnecessary. It applies to data-selected thresholds because the selection is fixed once X is conditioned on. It concerns the selected pairing and model allocation states; it does not provide posterior uncertainty for biological correspondence. Monte Carlo estimates, rather than exact probabilities, are used in code.

Evaluate every distinct absolute-contrast cutoff with exact score ties kept together. Among sets with estimated mean error <= alpha, take the largest. The average error can be nonmonotone in this contrast ordering. Do not stop at the first failure. The canonical strict cutoff is the next excluded score, or zero if all nonzero directions are included. Infinity means no discoveries. Below the retained-draw minimum, no set is selected even if a candidate curve passes alpha. Recompute the rule if defining a different inferential family.

The previous error-sorted prefix remains available as posterior_error. all_matchings replaces ell by the worst error over all control counterparts; it is a conservative sensitivity rule, not a posterior over mappings. Neither rule is claimed to give uniform frequentist FDR control. Empirical calibration requires repeated simulations with an explicitly specified truth.

Reference for the posterior expected error principle: Mueller, Parmigiani and Rice, FDR and Bayesian Multiple Comparisons Rules, https://web.ma.utexas.edu/users/pmueller/pap/MPR07.pdf .