Network analysis is a statistical methodology used to examine how variables within a system relate to one another. A network consists of nodes (representing variables) and edges (representing the relations between nodes). Networks can be either directed or undirected. In a directed network, edges indicate one-way relations, where one node is hypothesized to influence or predict another. In an undirected network, edges reflect bidirectional associations, such as co-occurence. Additionally, networks can be weighted or unweighted. In a weighted network, edges carry numerical values indicating the strength of the relation. In an unweighted network, edges simply indicate whether a relation is present or absent.

Depending on the modeling approach, edges may represent different types of relations. For example, one approach to building networks with continuous variables is called Gaussian Graphical Modeling, where we extract the partial correlations between variables to build the edges of our network. A partial correlation is the relation between two continuous variables while controlling for all other variables in the network. We can think of a partial correlation as reflecting how two variables change together, accounting for the influence of all other variables.

Another approach to building networks uses the Ising model (Brush, 1967; Van Borkulo et al., 2014). The Ising model is used specifically for binary variables, where edges \(\omega\) represent pairwise relations between nodes \(i\) and \(j\). An edge connecting two nodes is denoted as \(\omega_{ij}\), where \(\omega\) reflects the strength and direction (positive or negative) of the association between those nodes. In the Ising model, which uses binary nodes (e.g., 0 = inactive, 1 = active), a nonzero \(\omega_{ij}\) indicates the presence of an edge, whereas a value of zero indicates no direct connection between nodes. The pairwise relation models the log-odds of one node when another node is present and while controlling for all others. We can think of a pairwise relation as how the presence or absence of one node affects the likelihood of another node being present, while controlling for all other nodes in the network.

Depending on what our variables of interest are, we can examine if variables co-occur, predict each other (given we have longitudinal data), and how connected variables are within a system. A system would be the overarching construct of interest. For example, maybe we are interested in people’s preferences of ice cream flavor. The system would be flavor preferences, and the nodes in our network would be the different flavors, possibly chocolate, vanilla, or cookies n’ cream. A nonzero edge \(\omega_{ij}\) could indicate whether an individual who reports a preference for node \(i\), also reports a preference for node \(j\). In the context of clinical psychology research, a mental disorder can be viewed as a heterogeneous syndrome, or a system built from different symptoms with many relations that maintain the overall system.

Network temperature is a measure of how stable a network is. Lower values indicate a stable network with low variability and greater node alignment (i.e., how consistently nodes occur together), and higher values indicate a less stable network with higher variability and less node alignment.

The Ising model describes the probability of observing a particular configuration \(y\) of binary variables (e.g., symptoms being “on” or “off”, coded as 1 or 0) across a network of \(m\) nodes. The probability of observing a specific configuration 𝑦 is given by (Epskamp et al., 2018; Folk, 2025):

\[\Pr(Y=y)=\frac{exp(-\beta H(y;\tau,\Omega))}{Z(\tau,\Omega)}\]

Here we have:

In summary, this model allows us to examine how individual nodes via \(\tau\) and their pairwise interactions via \(\omega\) contribute to the likelihood of observing a particular node profile. In essence, the model quantifies how likely each node pattern is, based on the structure of the network and the parameters governing node activation and interactions.

The Hamiltonian function \(H(y;\tau,\Omega)\) captures this alignment and determines the likelihood of a particular configuration \(y\), or the likelihood of nodes occurring together or not occurring together. When we compute the probability of a given symptom pattern, we first compute its Hamiltonian (i.e., an alignment score where lower values reflect stronger alignment among nodes). We then multiply this by \(β\) and subsequently apply the exponential function and normalize using \(Z\).

If \(\beta\) is large, node patterns that align with the network structure are assigned much higher probabilities than misaligned ones—leading to a stable and peaked distribution. If \(\beta\) is small, many patterns are treated similarly, and the distribution is flatter and less stable. In this way, \(\beta\) directly controls how sharply the model distinguishes between likely and unlikely node patterns. To use the Ising model to construct a network, we must make the following assumptions:

  1. Our variables must be binary, meaning the variables must have two levels (e.g., yes or no, present or not present, 0/1).

  2. We assume that the dependencies of variables (i.e., the statistical associations between two variables) can be captured through pairwise interactions. This means we are only examining interactions between two variables at once, we would not examine a three way interaction or interaction between more variables.

  3. We assume that the relations between nodes are bidirectional, so the relation between A to B is the same as the relation from B to A.

  4. We assume nodes are conditionally independent given their neighbors, so if a node is not connected to another, their relation can be explained by the other nodes it is connected to.

  5. We assume that there are no unmeasured confounders, and that all variables of the system are included. Omitting variables may lead to spurious edges.

  6. We assume a homogeneous structure in group models. If we model a network among a group (for example, individuals with depression) we assume that the edges of the network are the same for everyone in the group.

This tutorial will examine two networks, one with binary depressive symptoms (i.e., measures of whether symptoms were present) and binary reasons for unmet mental health service need (i.e., measures of whether someone experienced a reason for unmet need). We will use the Ising model to construct networks and then examine how network temperature changes between male and female individuals in the depressive symptom network, and between sexual orientation groups in the reasons for unmet need network. We must fix the \(\beta\) parameter of a reference group for model identification to measure the relative temperature of other groups (Grimes et al., 2025), thus male and heterosexual individuals will have a temperature of one in their respective networks. We aim to answer the following two questions:

  1. How does the stability of depression as a system of symptoms differ from female individuals compared to male individuals?

  2. How does the stability of unmet need for mental health services as a system of reasons for unmet need differ among individuals with varying sexual orientations compared to heterosexual individuals?

Step 1: Load Data and Packages

library(psychonetrics)
library(tidyverse)
library(VIM)
library(purrr)
library(MLMusingR)
data <- read.csv("symptoms_nettemp.csv")
ndata <- read.csv("reasons_nettemp.csv")

The cross-sectional data used in this tutorial comes from the Collaborative Psychiatric Epidemiology Surveys and National Survey on Drug Use and Health and includes either full case analyses or already imputed data respectively.

Step 2: Prepare Data

#all nodes should be binary with 0/1 levels, 0 = not present, 1 = present
#Depressive Symptoms Network
#for sex, 1 = male, 2 = female
symptoms <- data %>%
  mutate(Sex = factor(Sex, levels = c(1,2), labels = c("male","female")),
         Sex = relevel(Sex, ref = "male"))

To prepare the data for Ising estimation, the dataframe should include only the binary node variables and the grouping variable. Once these are properly formatted, the model estimates the probability of every possible symptom configuration within each group, allowing us to compare how symptoms organize and interact across populations.

It’s essential that all node variables (in this case, symptoms) are binary. The Ising model assumes each node can only take on two values, typically coded as 0 (absent) or 1 (present). This binary coding is crucial because the model’s interaction term, \(\omega_ij y_i y_j\), relies on binary multiplication to capture whether two symptoms co-occur. For example, if both symptoms are present (1 × 1), the product is 1, signaling co-occurrence. But if either symptom is absent (0 × 1 or 0 × 0), the product is 0, indicating no co-occurrence. If node values are not binary (e.g., coded as 1 and 2), the product becomes uninterpretable in this context (e.g., 1 × 2 = 2), and the model can no longer accurately isolate meaningful pairwise relationships. As a result, the edge parameters (\(\omega\)) lose their psychological meaning, and the integrity of the model breaks down.

Step 3: Network Construction using the Ising Model

#nodes in the network
set.seed(123)
#vars argument just encodes which variables will be used from df, we technically don't need it since we only include variables of interest but doesn't hurt to have :D
vars <- names(symptoms)[1:9]

#building saturated model and running. This is an overfitted model with all relations between nodes before we adjust for significance. This is a dense model because it contains all edges
model1 <- Ising(symptoms, vars = vars, groups = "Sex") %>% runmodel

#equal networks, so now network structure between groups is equal. This means edges between nodes are equal across groups
model2 <- model1 %>% groupequal("omega") %>% runmodel

#equal thresholds, so now threshold/intercept structure between groups is equal, so network structure is equal and now external fields are also equal. External fields are the threshold structure of a network, so we are saying that the baseline probability of nodes occurring between groups is equal.
model3 <- model2 %>% groupequal("tau") %>% runmodel

#prune the saturated model and stepup to include only significant edges (finding a sparse model) -- we get rid of any non significant edges then add back in edges that significantly improve model fit
model3b <- model3 %>% prune(alpha = 0.05, adjust = "bonferroni") %>% stepup(mi = "mi_equal", alpha = 0.05, greedyadjust = "bonferroni")

#something to note: in Ising models beta must be fixed to 1 for a reference group in order to maintain model fit and estimate the relative differences of other groups. Allows us to model how temperature changes between groups or over time.

#below code <- to examine temp for diff groups
parameters(model3b) %>% filter(matrix == "beta")

final_model <- model3b

We build models that constrain the \(\Omega\) and \(\tau\) parameters when comparing groups before extracting the output from a final model to examine the differences in network temperature between groups. We are interested in how temperature changes between male and female individuals, but we do not want any observed differences in temperature to be due to differences in network structure or node thresholds. By constraining \(\Omega\) and \(\tau\) to be equal across groups, we can ensure that any group differences in network dynamics are attributable to differences in temperature, rather than changes to presence of edges or baseline node activity.

To ensure that any observed group differences in network temperature truly reflect differences in network dynamics, rather than differences in symptom prevalence or network structure, we adopt a constrained modeling approach. Specifically, we:

  1. Start with a saturated model that allows all node-to-node relationships to vary freely across groups.

  2. First, we constrain the edge weights (\(\Omega\)) to be equal across groups, fixing the structure of the network so that all groups share the same pattern of symptom connections.

  3. Next, we constrain the thresholds (\(\tau\)) to be equal, ensuring that the baseline likelihood of each symptom occurring is the same in all groups.

With both structure and baseline symptom activation held constant, we allow the temperature parameter (\(\beta\)) to vary freely across groups. This setup allows us to test whether groups differ in network rigidity—that is, whether their symptom systems favor a narrow set of co-occurring symptoms (low temperature) or are more variable and less predictable (high temperature).

To improve interpretability, we then apply a prune–stepup procedure:

This results in a final, sparse model that retains only the most statistically reliable symptom connections. We evaluate the fit of each model using standard metrics such as the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Root Mean Square Error of Approximation (RMSEA), and chi-squared model comparisons.

By holding structure and base rates constant, this approach lets us meaningfully interpret differences in temperature (\(\beta\)) as differences in how tightly or loosely symptoms cluster, shedding light on important variations in symptom dynamics across populations.

Step 4: Process output to plot temperature

#temperature estimates and 95% CIs from standard errors of the beta parameter
betas_est <- final_model@parameters$est[final_model@parameters$matrix == "beta"]
betas_se <- final_model@parameters$se[final_model@parameters$matrix == "beta"]
#pull the critical value from the standard normal dist for making CIs
z = qnorm(0.975)
upperCI <- betas_est + (z*betas_se)
lowerCI <- betas_est - (z*betas_se)
#set up temp plotting
sex <- c("Male","Female")
symptom_data <- data.frame(Sex = sex,
                           Temperature = 1/betas_est,
                           UpperCI = 1/upperCI,
                           LowerCI = 1/lowerCI)

To calculate network temperature, we extract the group-specific estimates of the \(\beta\) (inverse temperature) parameter as a matrix. Using the standard errors of these estimates, we compute 95% confidence intervals based on the normal approximation, using a critical value of 1.96. We then create a vector that matches the levels of our grouping variable (e.g., sex or sexual orientation). Finally, we construct a data frame that includes the grouping variable, the estimated temperature (calculated as the inverse of the \(\beta\) parameter), and the corresponding upper and lower bounds of the 95% confidence intervals, also inverted from the \(\beta\) estimates. This allows us to visualize differences in temperature to directly compare network stability across groups and determine whether any observed differences are statistically meaningful.

Step 5: Visualize Temperature

#plot depp symptoms temp
ggplot(symptom_data, aes(x = Sex, y = Temperature, group = 1, fill = Sex)) +
  geom_point(shape = 16) +
  geom_col(width = 0.6) +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.1) +
  scale_x_discrete() +
  labs(x = "Sex", y = "Temperature", title = "Temperature Difference") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, size = 12), legend.position='none') + 
  coord_cartesian(ylim = c(0.95, 1.1))

We then plot the network temperature using the data frame containing the estimated temperature and corresponding confidence intervals. The reference group (i.e., male individuals) has a temperature value of 1, with no confidence interval, because the \(\beta\) parameter for that group was fixed for model identification purposes.

Above is the bar graph showing the difference in temperature between female and male individuals in a network of depression symptoms. The results suggest that depression symptoms may be slightly more stable in female individuals than in male individuals; however, this conclusion is uncertain, as the confidence interval includes values above 1. A more stable network indicates less variability in symptom activation and a higher likelihood that symptoms co-occur or persist together for female individuals.

Below is another example examining network temperature differences in a network of reasons for unmet need between different sexual orientation groups using the same approach described above.

#Reasons network --- data is already imputed by NSDUH, only contains the first 10 reasons so we actually violate one of the assumptions (having all variables of a system present), it takes too long for my laptop to run
#only select variables of interest
reasons <- ndata %>%
  select(MHNTCOST, MHNTINSCV, MHNTENFCV, MHNTWHER, MHNTNOFND, MHNTPROBS, MHNTTIME, MHNTPRIV, MHNTPTHNK, MHNTHNDL, SEXIDENT22) %>%
  mutate(SEXIDENT22 = factor(SEXIDENT22, levels = c(1,2,3,4,5), labels = c("Heterosexual","Gay or Lesbian","Bisexual","Different Term","Unsure")),
         SEXIDENT22 = relevel(SEXIDENT22, ref = "Heterosexual")) %>%
  filter(!SEXIDENT22 %in% c(6,85,94,97,98)) %>%
  filter(!is.na(SEXIDENT22)) %>%
  filter(if_all(c(MHNTCOST, MHNTINSCV, MHNTENFCV, MHNTWHER, MHNTNOFND, MHNTPROBS, MHNTTIME, MHNTPRIV, MHNTPTHNK, MHNTHNDL), ~ !is.na(.)))
nmiss(reasons)
rvars <- names(reasons)[1:10]

#models for reasons (same process as above)
rmodel1 <- Ising(reasons, vars = rvars, groups = "SEXIDENT22") %>% runmodel
#Equal networks (omega = network structure, edges btw nodes equal across groups):
rmodel2 <- rmodel1 %>% groupequal("omega") %>% runmodel
#Equal thresholds (tau = threshold/intercept structure, omega still equal plus external fields equal):
rmodel3 <- rmodel2 %>% groupequal("tau") %>% runmodel
#Prune-stepup to find a sparse model:
rmodel3b <- rmodel3 %>% prune(alpha = 0.05, adjust = "bonferroni") %>% stepup(mi = "mi_equal", alpha = 0.05, greedyadjust = "bonferroni")

parameters(rmodel3b) %>% filter(matrix == "beta")

rfinal_model <- rmodel3b
#reasons 
rbetas_est <- rfinal_model@parameters$est[rfinal_model@parameters$matrix == "beta"]
rbetas_se <- rfinal_model@parameters$se[rfinal_model@parameters$matrix == "beta"]
rupperCI <- rbetas_est + (z*rbetas_se)
rlowerCI <- rbetas_est - (z*rbetas_se)
#set up temp plotting
SEXIDENT22 <- c("Heterosexual","Gay or Lesbian","Bisexual","Different Term","Unsure")
reasons_data <- data.frame(Sex_Ori = SEXIDENT22,
                           Temperature = 1/rbetas_est,
                           UpperCI = 1/rupperCI,
                           LowerCI = 1/rlowerCI)
#plot reasons temp
ggplot(reasons_data, aes(x = Sex_Ori, y = Temperature, group = 1, fill = Sex_Ori)) +
  geom_point(shape = 16) +
  geom_col(width = 0.6) +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.1) +
  scale_x_discrete() +
  labs(x = "Sexual Orientation", y = "Temperature", title = "Temperature Differences") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, size = 12), legend.position='none')+ 
  coord_cartesian(ylim = c(0.95, 1.525))

Above is the bar graph of the differences in temperature across sexual orientation groups in a network of reasons for unmet mental health need. We can conclude that individuals who are unsure of their sexual orientation, bisexual, or gay or lesbian individuals may have slightly less stable networks compared to heterosexual individuals, as indicated by higher temperature estimates. We can make the same conclusion for individuals who use a different term; however, we cannot be confident in our conclusion about individuals who use a different term as the confidence intervals include the value 1. A less stable network means there is greater variability in the presence of reasons within the network and that these reasons co-occur less consistently across individuals in that group. We should acknowledge that our reasons-for-unmet-need network violates one Ising-model assumption: that all relevant nodes be included. For this tutorial, we omitted additional reasons from the 2023 NSDUH because adding every node would substantially increase computational time.

References

Brush, S. G. (1967). History of the Lenz-Ising Model. Reviews of Modern Physics, 39(4), 883–893. https://doi.org/10.1103/RevModPhys.39.883

Epskamp, S., Maris, G., Waldorp, L. J., Borsboom, D., Irwing, P., Hughes, D. J., & Booth, T. (2018). Network Psychometrics. In The Wiley Handbook of Psychometric Testing (pp. 953–986). John Wiley & Sons. https://doi.org/10.1002/9781118489772.ch30

Finnemann, A., Borsboom, D., Epskamp, S., & Van Der Maas, H. L. J. (2021). The Theoretical and Statistical Ising Model: A Practical Guide in R. Psych, 3(4), 593–617. https://doi.org/10.3390/psych3040039

Folk, R. (2025). A new look at Ernst Ising’s thesis. The European Physical Journal B, 98(6). https://doi.org/10.1140/epjb/s10051-025-00954-x

Grimes, P. Z., Murray, A. L., Smith, K., Allegrini, A. G., Piazza, G. G., Larsson, H., Epskamp, S., Whalley, H. C., & Kwong, A. S. F. (2025). Network temperature as a metric of stability in depression symptoms across adolescence. Nature Mental Health, 3(5), 548–557. https://doi.org/10.1038/s44220-025-00415-5

Van Borkulo, C. D., Borsboom, D., Epskamp, S., Blanken, T. F., Boschloo, L., Schoevers, R. A., & Waldorp, L. J. (2014). A new method for constructing networks from binary data. Scientific Reports, 4(1), 5918. https://doi.org/10.1038/srep05918