Roadmap

Introduction

Networks Affect Outcomes

To study relationships among attributes in connected populations with interference, a unifying framework for attributes (\mathbf{X}, \mathbf{Y}) and connections \mathbf{Z} is needed

Networks Affect Outcomes

Advertisers target teenagers:

Networks Affect Outcomes

Teenagers have friends:

Networks Affect Outcomes

Teenagers have friends:

Networks Affect Outcomes

Network Interference

Predictors X_1, X_2 and outcomes Y_1, Y_2 of units 1, 2:

Regression Under Interference

Introduce comprehensive regression framework for studying relationships among attributes (\mathbf{X}, \mathbf{Y}) and connections \mathbf{Z}:

  • Scalable models, capturing heterogeneity along with dependence among (\mathbf{X}, \mathbf{Y}) and \mathbf{Z} in large populations.

  • Interpretable models, based on extensions of Generalized Linear Models (GLMs) to dependent (\mathbf{X}, \mathbf{Y}) and \mathbf{Z}.

  • Scalable methods, based on convex optimization of pseudolikelihoods, implemented by a divide and conquer approach using minorization-maximization (MM) methods.

  • Provable theoretical guarantees, based on a single observation of dependent random variables (\mathbf{X}, \mathbf{Y}) and \mathbf{Z}.

Motivation

  • Standard GLMs (lm, glm) assume independent units:
    • Attributes (X_i, Y_i) of unit i are unaffected by unit j
  • Problem: Violated in connected population because of spillover and contagion!
  • Solution: A regression framework for dependent attributes (\mathbf{X}, \mathbf{Y}) and network connections \mathbf{Z} implemented in R Package iglm

Publications

Software Paper

arXiv preprint (describing the iglm package):

Theory Paper

Journal of the American Statistical Association (2026):

R Package iglm

The Gap

  • Standard R functions (lm, glm) assume independent units
  • Ignoring network interference (spillover & contagion) can lead to bias & invalid inference
  • Available tools (e.g., ergm, ALAAMEE) do not handle joint (\mathbf{X}, \mathbf{Y}, \mathbf{Z}) models but only slices of it

The Solution

  • Extends GLMs: Familiar interface for regression under interference
  • Scalable Engine: Convex optimization using MM algorithms with C++ acceleration
  • Object-Oriented: Unified R6 class workflow for data, estimation & simulation

Advantages of R Package iglm

  • Comprehensive: Binary, count, and real-valued attributes \mathbf{X}, \mathbf{Y} & network connections \mathbf{Z}
  • Interpretable: Extension of standard GLMs with micro-behavioral foundation
  • Local Dependence: Neighborhoods \mathcal{N}_i to control dependencies
  • Scalable: Convex optimization via blockwise optimization
  • Guarantees: Error of O(N) weights decays at rate \sqrt{\log N / N}

Architecture & Setup

install.packages("iglm")
  • R6 class architecture (Chang, 2025):
    • iglm manages formula, sampler, and execution
    • iglm.data manages variables
  • Back-end implemented in C++ for high-performance computing

R6 Classes in R

Standard R (e.g., ergm)

Apply external functions to data/objects.

  • Syntax: Pass the object as an argument:

    # Fit the model:
fit <- ergm(nw ~ edges + mutual)
# Run external functions:
summary(fit)
sims <- simulate(fit, nsim = 10)

  • Copy-on-Modify: Modifying settings or data copies the entire object in memory

R6 Classes (e.g., iglm)

Data and functions live together in the object.

  • Syntax: Call functions/methods on the object:

    model <- iglm(nw ~ edges + mutual)
model$estimate()  # Fits in-place
model$simulate()  # Simulates in-place

  • In-Place Modification: Modifying the object in-place without copy overhead:

    model$set_control(max_it = 300)

Data Representation

Data Structure in iglm.data

  • Predictors \mathbf{X} = (X_i): binary, count, or real; fixed or random
  • Outcomes \mathbf{Y} = (Y_i): binary, count, or real
  • Connections \mathbf{Z} = (Z_{i,j}): binary Z_{i,j} \in \{0, 1\}, directed or undirected
  • Neighborhoods \mathcal{N}_i (exogenous): Overlap indicated by c_{i,j} = \mathbb{I}(\mathcal{N}_i \cap \mathcal{N}_j \neq \emptyset)

Running Example: Hate Speech on X

Is there a spillover effect from interacting with Republican politicians to employing offensive language in online communications?

  • Dataset: N = 495 state legislators across 6 states from Kim et al. (2022)
  • Outcome Y_i: Use of hate speech (1 = Yes, 0 = No)
  • Predictor X_i: Party affiliation (1 = Republican, 0 = Other) [fixed]
  • Network Z_{i,j}: Directed mention/repost network on X (Twitter)
  • Covariates: Gender (v_{i,1}), Race (v_{i,2}), and State (v_{i,3})

Running Example: Hate Speech on X

Is there a spillover effect from interacting with Republican politicians to employing offensive language in online communications?

Data Representation: iglm.data

  • Load the pre-packaged dataset directly:
data("state_twitter")
data.object <- state_twitter$iglm.data$clone()
  • Exogenous Covariates gender_attribute, white_attribute, match_gender, match_race, match_state are available via package
Show code 
# Extracting covariates from state_twitter object
match_gender     <- state_twitter$match_gender
match_race       <- state_twitter$match_race
match_state      <- state_twitter$match_state
white_attribute  <- state_twitter$white_attribute
gender_attribute <- state_twitter$gender_attribute
  • Import from raw files (vectors or matrices)
Show code 
z_network <- read.csv("z_network.csv")
x_attribute <- read.csv("x_attribute.csv")
y_attribute <- read.csv("y_attribute.csv")
neighborhood <- as.matrix(read.csv("neighborhood.csv"))

data.object <- iglm.data(
  x_attribute = x_attribute,  type_x = "binomial",
  y_attribute = y_attribute,  type_y = "binomial",
  z_network = z_network,      directed = TRUE,
  neighborhood = neighborhood
)

Descriptive Analysis: Data Object Summary

data.object
iglm.data object
  units                       : 495
  directed                    : TRUE
  edges (fixed = FALSE)       : 9218
  neighborhood edges          : 24398

Attribute summaries
  x_attribute (fixed = TRUE)  : binomial 1s=195, 0s=300, P(1)=0.394
  y_attribute                 : binomial 1s=204, 0s=291, P(1)=0.412

Descriptive Analysis: Visualizing the Network

  • The package provides a plotting function based on igraph
  • Node colors \rightarrow outcome variables \mathbf{X}
  • Node size \rightarrow predictor variables \mathbf{Y}
  • Grey lines \rightarrow directed connections \mathbf{Z}
  • Orange lines \rightarrow overlapping neighborhoods \mathcal{N}_i \cap \mathcal{N}_j \neq \emptyset
set.seed(123)
data.object$plot(
  edge.width = 0.5, 
  edge.arrow.size = 0.2, 
  legend_size = 0.5, 
  legend_size_n_levels = 2
)

Descriptive Analysis: Different Statistics

# Plot spillover degree distribution
data.object$spillover_degree_distribution()
  • Spillover degree distribution: The degree in the subnetwork of connections that can enable spillover (i.e., connections where c_{i,j} x_i y_j = 1 or c_{i,j} x_j y_i = 1)

# Plot outcome distribution
data.object$y_distribution()
  • Outcome distribution (y_distribution): The empirical distribution of the node-level outcomes Y (density or bar plot depending on the outcome type)

# Plot traditional degree distribution
data.object$degree_distribution()
  • Degree distribution: Distribution of number of outgoing and incoming connections

# Plot geodesic distances distribution
data.object$geodesic_distances_distribution()
  • Geodesic distance distribution: The distribution of shortest path distances (geodesic distances) between all pairs of connected actors in the symmetrized network

# Plot edgewise shared partner distribution
data.object$edgewise_shared_partner_distribution()
  • Edgewise shared partner distribution: The distribution of edgewise shared partner (ESP) counts (i.e., the number of mutual neighbors that two connected actors share in the network)

Framework

Joint Probability Model

\begin{aligned} f_{\theta}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \dfrac{1}{c(\theta)} &\left[ \prod_{i \in \mathcal{P}} a_{\mathbf{X}}(x_i) a_{\mathbf{Y}}(y_i) \exp(\theta_g^\top g_i(x_i, y_i)) \right] \\ \times &\left[ \prod_{(i,j) \in \mathcal{D}} a_{\mathbf{Z}}(z_{i,j}) \exp(\theta_h^\top h_{i,j}(x_i, x_j, y_i, y_j, \mathbf{z})) \right] \end{aligned}

  • \mathcal{P} is the set of units, \mathcal{D} is the set of dyads
  • a_{\mathbf{X}}, a_{\mathbf{Y}}, a_{\mathbf{Z}} are reference measures
  • g_i(\cdot) describes unit relationships
  • h_{i,j}(\cdot) describes network connections and spillover

Interpretable Models

  • If g_i and h_{i,j} are coordinatewise affine:
    • Conditionals of predictors X_i can be represented by GLMs
    • Conditionals of outcomes Y_i can be represented by GLMs
    • Conditionals of connections Z_{i,j} can be represented by GLMs
  • Key Advantage: Model results can be interpreted along the lines of standard GLMs

Illustrative Example: Unit Terms (g_i)

g_i(x_i, y_i) = \left( \begin{array}{c} y_i \\ x_i y_i \end{array} \right) \begin{array}{l} \longrightarrow \alpha_{\mathcal{Y}} \quad \text{(outcome intercept)} \\ \longrightarrow \beta_{\mathcal{X},\mathcal{Y}} \quad \text{(direct predictor effect)} \end{array}

Contribution to Linear Predictor \eta_i: \alpha_{\mathcal{Y}} + \beta_{\mathcal{X},\mathcal{Y}} x_i

Illustrative Example: Dyad Terms (h_{i,j})

h_{i,j} = \left(\begin{array}{c} z_{i,j} \\ -(1 - c_{i,j}) z_{i,j} \log N \\ c_{i,j} d_{i,j}(\mathbf{z}) z_{i,j} \\ c_{i,j} y_i x_j z_{j,i} \\ c_{i,j} y_i y_j z_{i,j} \end{array} \right) \begin{array}{lll} \longrightarrow & \theta_{\mathcal{O},i} + \theta_{\mathcal{I},j} & \text{(heterogeneity)} \\ \longrightarrow & \lambda & \text{(sparsity penalty)} \\ \longrightarrow & \gamma_{\mathcal{Z},\mathcal{Z}} & \text{(local transitivity)} \\ \longrightarrow & \color{#08696b}{\gamma_{\mathcal{X},\mathcal{Y},\mathcal{Z}}} & \color{#08696b}{\text{(spillover)}} \\ \longrightarrow & \color{#d9534f}{\gamma_{\mathcal{Y},\mathcal{Y},\mathcal{Z}}} & \color{#d9534f}{\text{(contagion)}} \end{array}

Contribution to Linear Predictor \eta_i: \color{#08696b}{\gamma_{\mathcal{X},\mathcal{Y},\mathcal{Z}} \sum_{j} c_{i,j} x_j z_{j,i}} + \color{#d9534f}{\gamma_{\mathcal{Y},\mathcal{Y},\mathcal{Z}} \sum_{j} c_{i,j} y_j (z_{i,j} + z_{j,i})}

Model Interpretation: Continuous Outcomes

  • Distribution: Y_i \mid (X, Y_{-i}, Z) \sim \mathcal{N}(\mu_i(\eta_i), \psi) with expectation: \mu_i(\eta_i) = \eta_i
  • Linear Predictor: \begin{aligned} \eta_i = \alpha_{\mathcal{Y}} + \beta_{\mathcal{X},\mathcal{Y}} x_i &+ \color{#08696b}{\gamma_{\mathcal{X},\mathcal{Y},\mathcal{Z}} \sum_{j} c_{i,j} x_j z_{j,i}} \\ &+ \color{#d9534f}{\gamma_{\mathcal{Y},\mathcal{Y},\mathcal{Z}} \sum_{j} c_{i,j} y_j (z_{i,j} + z_{j,i})} \end{aligned}
  • Extends Linear Regression (lm) to capture spillover and contagion
Show code 
# Data initialization with continuous outcome (type_y = "normal")
data.object <- iglm.data(
  y_attribute = y_attribute,
  type_y = "normal",
  scale_y = var(y_attribute)
)

# Or set on an existing data object:
data.object$set_type_y("normal")
data.object$set_scale_y(var(y_attribute))

Model Interpretation: Count Outcomes

  • Distribution: Y_i \mid (X, Y_{-i}, Z) \sim \text{Poisson}(\mu_i(\eta_i)) with expectation: \mu_i(\eta_i) = \exp(\eta_i)
  • Linear Predictor: \begin{aligned} \eta_i = \alpha_{\mathcal{Y}} + \beta_{\mathcal{X},\mathcal{Y}} x_i &+ \color{#08696b}{\gamma_{\mathcal{X},\mathcal{Y},\mathcal{Z}} \sum_{j} c_{i,j} x_j z_{j,i}} \\ &+ \color{#d9534f}{\gamma_{\mathcal{Y},\mathcal{Y},\mathcal{Z}} \sum_{j} c_{i,j} y_j (z_{i,j} + z_{j,i})} \end{aligned}
  • Extends Poisson Regression (glm(family = poisson)) to capture spillover and contagion
Show code 
# Data initialization with count outcome (type_y = "poisson")
data.object <- iglm.data(
  y_attribute = y_attribute,
  type_y = "poisson"
)

# Or set on an existing data object:
data.object$set_type_y("poisson")

Model Interpretation: Binary Outcomes

  • Distribution: Y_i \mid (X, Y_{-i}, Z) \sim \text{Bernoulli}(\mu_i(\eta_i)) with expectation: \mu_i(\eta_i) = \exp(\eta_i)/(1 + \exp(\eta_i))
  • Linear Predictor: \begin{aligned} \eta_i = \alpha_{\mathcal{Y}} + \beta_{\mathcal{X},\mathcal{Y}} x_i &+ \color{#08696b}{\gamma_{\mathcal{X},\mathcal{Y},\mathcal{Z}} \sum_{j} c_{i,j} x_j z_{j,i}} \\ &+ \color{#d9534f}{\gamma_{\mathcal{Y},\mathcal{Y},\mathcal{Z}} \sum_{j} c_{i,j} y_j (z_{i,j} + z_{j,i})} \end{aligned}
  • Extends Logistic Regression (glm(family = binomial)) to capture spillover and contagion
Show code 
# Data initialization with binary outcome (type_y = "binomial")
data.object <- iglm.data(
  y_attribute = y_attribute,
  type_y = "binomial"
)

# Or set on an existing data object:
data.object$set_type_y("binomial")

Local Dependence Assumption

  • Dependencies are localized using neighborhoods: h_{i,j}(x_i, x_j, y_i, y_j, \mathbf{z}) = h_{i,j}(z_{i,j}, z_{j,i}) \quad \text{for } \mathcal{N}_i \cap \mathcal{N}_j = \emptyset
  • Attributes of units with non-overlapping neighborhoods are conditionally independent.
  • Directed networks capture reciprocity via z_{j,i}.

Model Specification: Basics

data.object ~ term_1 + ...

  • Attribute Terms (g_i):
    • attribute_y, attribute_xy
  • Network Terms (h_{i,j}):
    • degrees,edges, transitive
  • Joint Terms (h_{i,j}):

Model Specification: Attribute Terms g_i

formula.iglm <- data.object ~
  attribute_y +                      # intercept for Y_i
  attribute_xy +                     # effect of X_i on Y_i
  cov_y(data = gender_attribute) +   # covariate effect: gender -> Y
  cov_y(data = white_attribute)      # covariate effect: race -> Y
  • attribute_y — baseline of the outcome Y_i: g_i(x_i,y_i) = y_i
  • attribute_xy — direct effect of predictor X_i on Y_i: g_i(x_i, y_i) = x_i y_i
  • cov_y(data = v) — additive effect of covariate v_i on Y_i: g_i(x_i, y_i) = v_i y_i

Model Specification: Network Terms h_{i,j}

formula.iglm <- data.object ~
  edges(mode = "alocal") +                      # baseline tie propensity (alocal)
  cov_z(data = match_gender, mode = "local") +  # gender homophily within N
  cov_z(data = match_race,   mode = "local") +  # race homophily within N
  cov_z(data = match_state,  mode = "local") +  # state homophily within N
  degrees                                       # in/out-degree node parameters
  • edges(mode = "alocal") — baseline for non-neighborhood dyads: h_{i,j}(x_i, x_j, y_i, y_j, \mathbf{z}) = (1 - c_{i,j})\, z_{i,j}
  • cov_z(data = m, mode = "local") — homophily within overlapping neighborhoods: h_{i,j}(x_i, x_j, y_i, y_j, \mathbf{z}) = c_{i,j}\, m_{i,j}\, z_{i,j}
  • degrees — \theta_{\mathcal{O},i} and \theta_{\mathcal{I},j} per node (depends if directed or undirected)

Model Specification: Network Terms h_{i,j}

formula.iglm <- data.object ~
  mutual(mode = "local") +                      # reciprocity: z_{i,j} z_{j,i} within N
  transitive                                    # transitivity: shared partners
  • mutual(mode = "local") — reciprocity within overlapping neighborhoods: h_{i,j}(x_i, x_j, y_i, y_j, \mathbf{z}) = \frac{1}{2}\, c_{i,j}\, z_{i,j}\, z_{j,i}
  • transitive — transitivity (local shared partners) within overlapping neighborhoods: h_{i,j}(x_i, x_j, y_i, y_j, \mathbf{z}) = c_{i,j}\, d_{i,j}(\mathbf{z})\, z_{i,j}, where d_{i,j}(\mathbf{z}) = \mathbb{I}\left(\sum_{k \neq i,j} c_{i,k} c_{j,k} z_{i,k} z_{k,j} \geq 1\right) is the indicator for at least one local shared partner.

Model Specification: Spillover Terms h_{i,j}(z_{i,j})

formula.iglm <- data.object ~
  spillover_yx(mode = "local") +   # X_j -> Y_i via z_{j,i}  (treatment spillover)
  spillover_yy(mode = "local")     # Y_j -> Y_i via z_{i,j} + z_{j,i}  (outcome contagion)
  • spillover_yx — treatment spillover from neighbor’s predictor (X_j) to unit’s outcome (Y_i): h_{i,j}(x_i, x_j, y_i, y_j, \mathbf{z}) = c_{i,j}\, y_i\, x_j\, z_{j,i}
  • spillover_yy — outcome contagion from neighbor’s outcome (Y_j) to unit’s outcome (Y_i): h_{i,j}(x_i, x_j, y_i, y_j, \mathbf{z}) = c_{i,j}\, y_i\, y_j\, z_{i,j}

Model Specification: Full Model Formula

formula.iglm <- data.object ~
  # -- Attribute terms --
  attribute_y + attribute_xy +
  cov_y(data = gender_attribute) + cov_y(data = white_attribute) +
  # -- Network terms --
  edges(mode = "alocal") +
  cov_z(data = match_gender, mode = "local") +
  cov_z(data = match_race,   mode = "local") +
  cov_z(data = match_state,  mode = "local") +
  degrees + mutual(mode = "local") + transitive +
  # -- Spillover/joint terms --
  spillover_yx(mode = "local") + spillover_yy(mode = "local")

Model Specification: Scaling and Localization

  • Localization restricts spillovers to overlapping neighborhoods (mode = "local")
  • Scaling controls spillover in large neighborhoods (important for count outcomes or high degree heterogeneity):

Unscaled (spillover_yy) \left( y_i y_j + y_j y_i \mathbb{I}_U(\mathbf{z}) \right) e_{i,j}^{(\mathtt{s})}

Scaled (spillover_yy_scaled) \left( \frac{y_i y_j}{\operatorname{deg}(i, \mathtt{s})} + \frac{y_j y_i}{\operatorname{deg}(j, \mathtt{s})} \mathbb{I}_U(\mathbf{z}) \right) e_{i,j}^{(\mathtt{s})}

Design Configurations

  • Fixed vs. random designs are specified via:
data.object$set_fix_x(TRUE)          # Fixed predictors
data.object$set_fix_z(TRUE)          # Fixed connections
data.object$set_fix_z_alocal(TRUE)   # Fixes non-overlapping network
  • Can also all be done when defining the data object with iglm.data()
iglm.data( ...,
  fix_x = TRUE,
  fix_z_alocal = TRUE)
  • Random designs allow generalization to superpopulations
  • Fixed designs limit inferences to observed states
    • iglm generalizes both glm and ergm

Estimation

Estimation: Pseudo-loglikelihood

Maximum likelihood is intractable due to normalizing constant

  • Maximize the pseudo-loglikelihood via an MM algorithm:

\ell(\theta; \mathbf{x}, \mathbf{y}, \mathbf{z}) = \sum_{i \in \mathcal{P}} \ell_{\mathbf{X},i}(\theta) + \sum_{i \in \mathcal{P}} \ell_{\mathbf{Y},i}(\theta) + \sum_{(i,j) \in \mathcal{D}} \ell_{\mathbf{Z},i,j}(\theta)

Optimization: Partitioned MM Algorithm

High-dimensional degree parameters (2N) require partitioned optimization:

  • Step 1 (MM): Update degree weights \theta_1 using a Minorization-Maximization algorithm (Quasi-Newton accelerated)
  • Step 2 (Newton-Raphson): Update parameters \theta_2

Estimation: Execution

# Initialize model from formula
model.iglm <- iglm(formula = formula.iglm)

# Customize optimizer and UQ parameters
model.iglm$set_control(control.iglm(
  max_it = 300,              
  tol = 1e-4,                
  offset_nonoverlap = 0,     
  var_method = "Mean-value" 
))

# Fit model parameters (runs in a few minutes)
model.iglm$estimate()
# Initialize model from formula
model.iglm <- iglm(formula = formula.iglm)

# Use pseudo Hessian standard errors
model.iglm$set_control(control.iglm(
  max_it = 100,              
  tol = 1e-3,                
  offset_nonoverlap = 0,     
  var_method = "Hessian" 
))

# Fit model parameters (runs in ~5 seconds)
model.iglm$estimate()
  • max_it & tol: Optimizer budget & stopping threshold
  • offset_nonoverlap: Penalty for sparse connections outside neighborhoods
  • var_method: Method for standard error calculation (MCMC-based mean-value vs. pseudo Hessian)

Estimation: Diagnostics

Plot convergence diagnostic paths via $plot(trace = TRUE) method:

model.iglm$plot(trace = TRUE)

Simulation

Simulation

Draw samples from the joint distribution of (X,Y,Z) under \theta via MCMC:

  • Simulating Outcomes and Predictor (Y \& X):
    • Gibbs sampler to cycle over each unit
  • Simulating Connections (Z):
    • Metropolis-Hastings or TNT sampler for local connections
    • Binomial sampling for non-overlapping connections

Simulation: Specification

# Define samplers for Y and Z only since X is fixed in this example
sampler.y.obj <- sampler.net.attr(n_proposals = data.object$n_actor * 10)
sampler.z.obj <- sampler.net.attr(n_proposals = nrow(data.object$overlap) * 10, tnt = TRUE)
sampler.obj <- sampler.iglm(
  n_burn_in = 1, n_simulation = 1, seed = 123,
  sampler_y = sampler.y.obj, sampler_z = sampler.z.obj
)
  • n_proposals: MCMC proposals per iteration
  • n_burn_in: Discarded initial iterations to reach stationary distribution
  • n_simulation: Number of simulated datasets to generate and store
  • seed: Random seed to ensure reproducibility

Simulation: Example

model.iglm$set_sampler(sampler.obj)
model.iglm$simulate() 
set.seed(123)
model.iglm$results$samples[[1]]$plot(
  edge.width = 0.5, 
  edge.arrow.size = 0.2, 
  legend_size = 2, 
  legend_size_n_levels = 2
)

Simulation: Applications

MCMC simulations are used for two main tasks:

  1. Uncertainty Quantification (UQ):
    • Approximation of standard errors \mathbb{V}_{\widehat{\theta}}(\widehat{\theta}) via simulations
  2. Model Assessment (Goodness-of-Fit):
    • Evaluates whether simulations from the fitted model replicate auxiliary network features
    • Via the $assess() method: Statistics from the Descriptive Analysis slide may be used

Goodness-of-Fit: Network Stats

# Other available functions include: degree_distribution, spillover_degree_distribution,
# dyadwise_shared_partner_distribution, y_distribution, and x_distribution
model.iglm$assess(formula = ~ edgewise_shared_partner_distribution +
                              geodesic_distances_distribution)

Goodness-of-Fit: Spillover Degrees

model.iglm$assess(formula = ~ spillover_degree_distribution)

Application

Application: Coefficient Table

Use the $summary() method on the fitted model to inspect results:

model.iglm$summary()
iglm object
----------------------------------------------------------------------------
Results: 

                                           Estimate    S.E. t-value Pr(>|t|)
attribute_y                                 -1.8592  0.3413   -5.45   <1e-04
attribute_xy                                -0.0518  0.2561   -0.20     0.84
cov_y(data = gender_attribute)              -0.1705  0.2627   -0.65     0.52
cov_y(data = white_attribute)                0.0807  0.2796    0.29     0.77
edges(mode = 'alocal')                      -1.0557  0.0035 -305.43   <1e-04
cov_z(data = match_gender, mode = 'local')   0.2224  0.0142   15.70   <1e-04
cov_z(data = match_race, mode = 'local')     0.6326  0.0101   62.75   <1e-04
cov_z(data = match_state, mode = 'local')    3.8671  0.0237  163.45   <1e-04
mutual(mode = 'local')                       2.5845  0.0454   56.97   <1e-04
transitive                                   0.6711  0.0138   48.76   <1e-04
spillover_yx(mode = 'local')                 0.0323  0.0394    0.82     0.41
spillover_yy(mode = 'local')                 0.0626  0.0130    4.81   <1e-04
                                              
attribute_y                                ***
attribute_xy                                  
cov_y(data = gender_attribute)                
cov_y(data = white_attribute)                 
edges(mode = 'alocal')                     ***
cov_z(data = match_gender, mode = 'local') ***
cov_z(data = match_race, mode = 'local')   ***
cov_z(data = match_state, mode = 'local')  ***
mutual(mode = 'local')                     ***
transitive                                 ***
spillover_yx(mode = 'local')                  
spillover_yy(mode = 'local')               ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Time for estimation: 11 mins

Degree Parameters:
  Outdegrees:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -11.2    -7.9    -6.9    -7.1    -6.2    -3.9 

  Indegrees:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -3.980  -1.108  -0.564  -0.616  -0.039   2.964

Interpretation: Conditional distribution of Y_i

Y_i \mid (X, Y_{-i}, Z) \sim \text{Bernoulli}(\mu_{\mathbf{Y},i}(\widehat\eta_{\mathbf{Y},i}))\\ \mu_{\mathbf{Y},i}(\widehat\eta_{\mathbf{Y},i}) = \dfrac{\exp(\widehat\eta_{\mathbf{Y},i})}{1 + \exp(\widehat\eta_{\mathbf{Y},i})}, Linear predictor for Hate Speech (Y_i \in \{0, 1\}): \begin{aligned} \widehat\eta_{\mathbf Y,i} &\coloneqq -1.859 \;-.052 \, x_{i} \;-.171\, v_{i,1} \;+.081\, v_{i,2} \\ &+ .032 \sum_{j \in \mathcal{P} \setminus \{i\}} c_{i,j}\, x_j\, z_{j,i} \;+ .063 \sum_{j \in \mathcal{P} \setminus \{i\}} c_{i,j}\, y_j\, (z_{i,j} + z_{j,i}) \end{aligned}

  • Contagion Effect: The outcome spillover coefficient (spillover_yy = .063) is positive and significant (p < .0001)

Interpretation: Conditional distribution of Z_{i,j}

Z_{i,j} \mid (X, Y, Z_{-i,j}) \sim \text{Bernoulli}(\mu_{\mathbf{Z},i,j}(\widehat\eta_{\mathbf{Z},i,j}))\\ \mu_{\mathbf{Z},i,j}(\widehat\eta_{\mathbf{Z},i,j}) = \dfrac{\exp(\widehat\eta_{\mathbf{Z},i,j})}{1 + \exp(\widehat\eta_{\mathbf{Z},i,j})},

Linear predictor for Z_{i,j} \in \{0, 1\}: \begin{aligned} \widehat{\eta}_{\mathbf{Z},i,j} &= \widehat{\theta}_{\mathcal{O},i} + \widehat{\theta}_{\mathcal{I},j} - 1.05(1 - c_{i,j}) \\ &\quad + \Big[ 3.867\, \text{state} + ... + \\ &\quad\,\,\,\,\quad 2.585\, \text{mutual} + .671\, \text{trans} + .032\, y_i x_j + .063\, y_i y_j \Big] c_{i,j} \end{aligned}

  • Strong state homophily, reciprocity, and transitivity

Model Comparison

Model Comparison

How can we compare different specifications of the model?

  • Compare full IGLM model containing spillover terms, reciprocity, and transitivity with a GLM ignoring interference 
formula.glm <- data.object ~
  attribute_y + attribute_xy +
  cov_y(data = white_attribute) + cov_y(data = gender_attribute) +
  edges(mode = "alocal") + edges(mode = "local") +
  cov_z(data = match_gender, mode = "local") +
  cov_z(data = match_race,   mode = "local") +
  cov_z(data = match_state,  mode = "local")

model.glm <- iglm(formula = formula.glm, name = "glm")
model.glm$estimate()

Predictive Performance: IGLM vs GLM

Show code 
pred.iglm <- model.iglm$predict(variant = "marginal")
pred.glm  <- model.glm$predict(variant = "marginal")
roc.iglm <- pROC::roc(pred.iglm$y$target, pred.iglm$y$prediction, quiet = TRUE)
roc.glm  <- pROC::roc(pred.glm$y$target, pred.glm$y$prediction, quiet = TRUE)

par(mar = c(4, 4, 1, 1))
plot(roc.iglm$specificities, roc.iglm$sensitivities, col = "#08696b",
     type = "l", lwd = 3, bty = "l", las = 1,
     xlab = "Specificity", ylab = "Sensitivity")
lines(roc.glm$specificities, roc.glm$sensitivities, col = "#e64173",
      type = "l", lwd = 3, bty = "l")
legend("bottomleft", legend = c("iglm", "glm"), lty = c(1,1), lwd = c(3,3),
       col = c("#08696b", "#e64173"), bty = "n")
  • IGLM (teal) outperforms GLM (red) in ROC-AUC for predicting hate speech (Y_i)
  • Network fit is also improved
  • Takeaway: Incorporating spillover and structural network effects improves predictions

Goodness-of-Fit: IGLM vs GLM

Show code 
# 1. Run simulations (plot = FALSE)
model.iglm$assess(
  formula = ~ edgewise_shared_partner_distribution + 
              spillover_degree_distribution, 
  plot = FALSE
)
model.glm$assess(
  formula = ~ edgewise_shared_partner_distribution + 
              spillover_degree_distribution, 
  plot = FALSE
)
# 2. Plot overlaid comparison
model.glm$results$plot(
  model_assessment = TRUE, 
  "iglm" = model.iglm$results$model_assessment
)

Comparison of IGLM and GLM via Goodness-of-Fit (GOF) simulations:

  1. Run $assess() on both models with plot = FALSE to compute simulations
  2. Call the plot method with model_assessment = TRUE, passing the comparison results

Custom Terms: Recipe

  • iglm uses a fast C++ backend.
  • Users can create their own custom terms:
create_userterms_skeleton()
  • Four Step Recipe:
    1. Decompose model terms g_i(\cdot) and h_{i,j}(\cdot).
    2. Specify change statistics \Delta_{\mathbf{X},i}, \Delta_{\mathbf{Y},i}, \Delta_{\mathbf{Z},i,j}.
    3. Write C++ code using XYZ_class methods and register.
    4. Write R wrapper InitIglmTerm.my_term.

Custom Terms: C++ Change Statistic

double my_stat_spillover(const XYZ_class &object, const int &unit_i, const int &unit_j,
                         const arma::mat &data, const double &type, const std::string &mode,
                         const bool &is_full_neighborhood) {
  double res = 0.0;
  if (mode == "x") { // x_i from 0 -> 1
    for (const int &k : object.adj_list_nb.at(unit_i)) res += object.y_attribute.get_val(k);
  } else if (mode == "y") { // y_i from 0 -> 1
    for (const int &k : object.adj_list_nb.at(unit_i)) res += object.x_attribute.get_val(k);  
  } else { // z_ij from 0 -> 1
    if (object.get_val_overlap(unit_i, unit_j)) {
      res = (object.x_attribute.get_val(unit_i) * object.y_attribute.get_val(unit_j)) +
            (object.x_attribute.get_val(unit_j) * object.y_attribute.get_val(unit_i));
    }
  }
  return res;
}
EFFECT_REGISTER("my_spillover", ::my_stat_spillover, "my_spillover", 0);

Custom Terms: R Wrapper Function

  • Define wrapper in the R folder of iglm.userterms:
InitIglmTerm.my_spillover <- function(data_object, arglist, ...) {
  arglist <- iglm:::check.IglmTerm(data_object, arglist, directed = FALSE)
  list(
    term_name = "my_spillover",
    coef_name = arglist$label
  )
}
  • Compile the generated package and load it:
library(iglm)
library(iglm.userterms)

Exercise

Exercise: Copenhagen Network Study

  • Dataset: N = 409 students tracked over 28 days (Sapiezynski, Stopczynski, Lassen, & Lehmann, 2019)
  • Outcome Y_i \in \mathbb{R}: Log-transformed phone call duration \log(t_i)
  • Predictor X_i: Gender (1 = Female, 0 = Male) [Fixed]
  • Network Z_{i,j}: Undirected friendship network
  • Neighborhood \mathcal{N}_i: Bluetooth proximity scans (\ge 24 hours)

Hints: Continuous Outcomes in iglm

Since Y_i is continuous (real-valued), we must use features not needed in the binary case:

  1. Outcome Scale Parameter: Before estimation, you must set the scale parameter (variance) on the data object using $set_scale_y():

    data("copenhagen")
copenhagen$set_scale_y(var(copenhagen$y_attribute))

  2. Model Assessment: To inspect the continuous outcome distribution fit, pass y_distribution to the $assess() method:

    # Note: assess() accepts: degree_distribution, spillover_degree_distribution,
# geodesic_distances_distribution, edgewise_shared_partner_distribution,
# dyadwise_shared_partner_distribution, y_distribution, x_distribution
model$assess(formula = ~ y_distribution + geodesic_distances_distribution)

Exercise: Setup & Data Preparation

  • Load the Copenhagen network study dataset
  • Compute and set the outcome scale parameter (variance of y_attribute) using $set_scale_y()
  • Check whether the distribution of y_attribute is set to "normal"
  • Load the Copenhagen network study dataset
  • Compute and set the outcome scale parameter (variance of y_attribute) using $set_scale_y()
  • Check whether the distribution of y_attribute is set to "normal"
data("copenhagen")
copenhagen$set_scale_y(var(copenhagen$y_attribute))
copenhagen
iglm.data object
  units                       : 409
  directed                    : FALSE
  edges (fixed = FALSE)       : 2517
  neighborhood edges          : 744

Attribute summaries
  x_attribute (fixed = TRUE)  : binomial 1s=90, 0s=319, P(1)=0.22
  y_attribute                 : normal mean=1.33, sd=1.78, scale= 3.16

Exercise: Model Specification

  • Formulate the joint model formula for continuous outcomes
  • Include the following terms: attribute_y, attribute_xy, degrees, edges(mode = "alocal"), transitive, spillover_xy, and spillover_yy
  • Formulate the joint model formula for continuous outcomes.
  • Include the following terms: attribute_y, attribute_xy, degrees, edges(mode = "alocal"), transitive, spillover_xy, and spillover_yy
formula.norm <- copenhagen ~
  attribute_y + attribute_xy + degrees +
  edges(mode = "alocal") + transitive +
  spillover_xy + spillover_yy

Exercise: Model Estimation

  • Initialize the model from the formula using iglm()
  • Estimate the model parameters using the $estimate() method for max. 300 iterations
  • Output the coefficient summary table using $summary()
  • Initialize the model from the formula using iglm()
  • Estimate the model parameters using the $estimate() method for max. 300 iterations
  • Output the coefficient summary table using $summary()
model.norm <- iglm(formula = formula.norm, control = control.iglm(max_it = 300))
model.norm$estimate()
model.norm$summary()
iglm object
--------------------------------------------------------------
Results: 

                       Estimate    S.E. t-value Pr(>|t|)    
attribute_y              1.0816  0.1078    10.0   <1e-04 ***
attribute_xy             0.4203  0.1781     2.4   0.0183 *  
edges(mode = 'alocal')  -1.2457  0.0013  -983.9   <1e-04 ***
transitive               2.8100  0.1653    17.0   <1e-04 ***
spillover_xy            -0.2447  0.1690    -1.4   0.1478    
spillover_yy             0.3416  0.1094     3.1   0.0018 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Time for estimation: 13 secs

Degree Parameters:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -4.51   -1.99   -1.51   -1.51   -1.00    0.77

Exercise: Model Assessment

  • Evaluate the goodness-of-fit via MCMC simulations
  • Include the continuous outcome distribution (y_distribution)
  • Evaluate the goodness-of-fit via MCMC simulations
  • Include the continuous outcome distribution (y_distribution)
model.norm$assess(formula = ~ y_distribution)

Exercise: Interpretation

  • Interpret the demographic, network structure, and treatment/outcome contagion effects
model.norm$summary()
iglm object
--------------------------------------------------------------
Results: 

                       Estimate    S.E. t-value Pr(>|t|)    
attribute_y              1.0816  0.1078    10.0   <1e-04 ***
attribute_xy             0.4203  0.1781     2.4   0.0183 *  
edges(mode = 'alocal')  -1.2457  0.0013  -983.9   <1e-04 ***
transitive               2.8100  0.1653    17.0   <1e-04 ***
spillover_xy            -0.2447  0.1690    -1.4   0.1478    
spillover_yy             0.3416  0.1094     3.1   0.0018 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Time for estimation: 13 secs

Degree Parameters:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -4.51   -1.99   -1.51   -1.51   -1.00    0.77
  • Outcome Contagion (spillover_yy): Significant (\widehat{\theta} = 0.327)
  • Gender Effect: Females call longer than men

Conclusion

What iglm Does

  • Joint regression for dependent outcomes \mathbf{Y}, predictors \mathbf{X}, and connections \mathbf{Z} under interference with support for binary, count, and continuous outcomes
  • Provides MCMC simulation for uncertainty quantification and goodness-of-fit
  • Fast C++ backend for simulation and optimization

How to Use & Contribute

Publications

Software Paper

arXiv preprint (describing the iglm package):

Theory Paper

Journal of the American Statistical Association (2026):

Literature

Chang, W. (2025). R6: Encapsulated classes with reference semantics. https://doi.org/10.32614/CRAN.package.R6
Fritz, C., & Schweinberger, M. (2025). Iglm: Regression under interference in connected populations. Retrieved from https://CRAN.R-project.org/package=iglm
Fritz, C., & Schweinberger, M. (2026). R Package iglm: Regression under Interference in Connected Populations. https://doi.org/10.48550/arXiv.2604.22791
Fritz, C., Schweinberger, M., Bhadra, S., & Hunter, D. R. (2025). A Regression Framework for Studying Relationships among Attributes under Network Interference. Journal of the American Statistical Association, (OnlineFirst), 1–22. https://doi.org/10.1080/01621459.2025.2565851
Kim, T., Nakka, N., Gopal, I., Desmarais, B. A., Mancinelli, A., Harden, J. J., … Boehmke, F. J. (2022). Attention to the COVID-19 pandemic on Twitter: Partisan differences among U.S. state legislators. Legislative Studies Quarterly, 47, 1023–1041. https://doi.org/10.1111/lsq.12367
Sapiezynski, P., Stopczynski, A., Lassen, D. D., & Lehmann, S. (2019). Interaction data from the Copenhagen Networks Study. Scientific Data, 6, 315.