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-occurrence. 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(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 model is a dense network because it contains all relations between nodes (i.e., each node has an edge connected to every other node) before we adjust for significance.
model1 <- Ising(symptoms, vars = vars, groups = "Sex") %>% runmodel

#prune the saturated model to include only significant edges. We do this with the prune and stepup functions, which remove and add edges that significantly improve model fit, respectively. 
model1b <- model1 %>% prune(alpha = 0.05) %>% 
  stepup(mi = "mi_equal", alpha = 0.05)

#constrain omega parameter (i.e., network structure). This means that all groups' networks will have the same set of edges
model2 <- model1 %>% groupequal("omega") %>% runmodel

#prune-stepup to find a sparse model
model2b <- model2 %>% prune(alpha = 0.05) %>% 
  stepup(mi = "mi_equal", alpha = 0.05)

#constrain tau (i.e, network thresholds). This means that all groups' initial log odds of each node when no other nodes are present will be equal 
model3 <- model2 %>% groupequal("tau") %>% runmodel

#prune stepup to find a sparse model
model3b <- model3 %>% prune(alpha = 0.05) %>% 
  stepup(mi = "mi_equal", alpha = 0.05)

#constrain beta (i.e., inverse network temperature). This means that all groups' networks are scaled uniformly.
model4 <- model3 %>% groupequal("beta") %>% runmodel

#prune stepup to find a sparse model
model4b <- model4 %>% prune(alpha = 0.05) %>% 
  stepup(mi = "mi_equal", alpha = 0.05)

#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. This is why we will see the reference group (in this case male individuals) have a temperature of 1. 

#we can now examine the model fits and determine which model is best
psychonetrics::compare(
  `1. all parameters free (dense)` = model1,
  `2. all parameters free (sparse)` = model1b,
  `3. equal networks (dense)` = model2,
  `4. equal networks (sparse)` = model2b,
  `5. equal networks and thresholds (dense)` = model3,
  `6. equal networks and thresholds (sparse)` = model3b,
  `7. all parameters equal (dense)` = model4,
  `8. all parameters equal (sparse)` = model4b
) %>% arrange(BIC)
#the model with the best fit is model4b
final_model <- model4b

We build models that constrain the \(\Omega\), \(\tau\), and \(\beta\) parameters when comparing groups before extracting the output from a final model with the best fit using BIC as the criterion. This allows us to examine whether meaningful differences in network temperature (and structure and thresholds! :D) exist between groups. By constraining the \(\Omega\), \(\tau\), and \(\beta\) parameters to be equal across groups, we can investigate whether any group differences in network dynamics are attributable to singular parameters or differences in multiple system level dynamics.

This constraining procedure is well documented in previous work (Grimes et al., 2025). 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.

  4. Finally, we constrain the temperature (\(\beta\)) to be equal, ensuring that all groups’ distributions are scaled the same. This means that all configurations of symptoms have the same exact probability across groups.

If we were to keep networks with all edges (i.e., dense networks), we’d be assuming that all symptoms are conditionally dependent (i.e., each edge is conditioned by all other edges), which is unlikely in symptom networks. Therefore, we apply a prune–stepup procedure:

This results in a final, sparse model that retains the minimal edges needed to explain our data. 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.

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. Notice that the temperature for both groups is one. We actually don’t even have confidence intervals, because from the beginning our temperature was set to be equal across groups. This is because the model with the best fit (i.e., model4b) had constraints on all parameters. The results suggest that depression symptoms as a network do not meaningfully differ between female and male individuals. Of course, this is using cross-sectional data from a nationally representative dataset, and we could see differences unfold at other timescales (and across gender identity as opposed to sex).

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

#sparse network
rmodel1b <- rmodel1 %>% prune(alpha = 0.05) %>% 
  stepup(mi = "mi_equal", alpha = 0.05)

#Equal networks (omega = network structure, edges between nodes equal across groups):
rmodel2 <- rmodel1 %>% groupequal("omega") %>% runmodel

#sparse network
rmodel2b <- rmodel2 %>% prune(alpha = 0.05) %>% 
  stepup(mi = "mi_equal", alpha = 0.05)

#Equal thresholds (tau = threshold/intercepts/external fields, omega still equal plus tau equal now):
rmodel3 <- rmodel2 %>% groupequal("tau") %>% runmodel

#sparse network
rmodel3b <- rmodel3 %>% prune(alpha = 0.05) %>% 
  stepup(mi = "mi_equal", alpha = 0.05)

#equal temperature (beta = inverse temperature, both omega and tau are still equal)
rmodel4 <- rmodel3 %>% groupequal("beta") %>% runmodel

#sparse network
rmodel4b <- rmodel4 %>% prune(alpha = 0.05) %>% 
  stepup(mi = "mi_equal", alpha = 0.05)

#we can now examine the model fits and determine which model is best
psychonetrics::compare(
  `1. all parameters free (dense)` = rmodel1,
  `2. all parameters free (sparse)` = rmodel1b,
  `3. equal networks (dense)` = rmodel2,
  `4. equal networks (sparse)` = rmodel2b,
  `5. equal networks and thresholds (dense)` = rmodel3,
  `6. equal networks and thresholds (sparse)` = rmodel3b,
  `7. all parameters equal (dense)` = rmodel4,
  `8. all parameters equal (sparse)` = rmodel4b
) %>% arrange(BIC)

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 that 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