Network Temperature Informed Simulations of Multigroup Ising Models Using the NodeIdentifyR Algorithm

In silico interventions allow researchers to probe how changes in a system translate to different outcomes for individuals (typically clinical patients). For example, among individuals struggling with depression, an in silico intervention could involve simulating how activated symptom networks are after treating a specific symptom (e.g., reducing the occurrence of said symptom). The NodeIdentifyR Algorithm (NIRA) is an application of in silico interventions using network analysis, specifically for psychopathology research.

Network analysis is a statistical methodology that examines how a set of variables in a system are related to each other. Networks consist of nodes (i.e., variables) and edges (i.e., relations between variables). One application of network analysis involves examining undirected networks, where edges represent the co-occurrence between two nodes. Furthermore, edges within networks can be weighted (i.e., edges carry a numeric value denoting the strength of the relation between two nodes) or unweighted (i.e., edges indicate co-occurrence between two nodes).

NIRA utilizes network structures and thresholds estimated via the Ising model. The following explanation of the Ising model is taken from my previous tutorial detailing how to examine network temperature using a multigroup Ising framework. See https://rpubs.com/jpstand/1417075.

The Ising model (Brush, 1967; Van Borkulo et al., 2014) 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:

\(y = (y_1,y_2,...,y_m)\) is a vector of binary values representing the state (e.g., presence or absence) of each node. An example of a particular configuration where all 3 nodes are present is \(y = [1,1,1]\). In an Ising model nodes are binary and can take values of either 0 (indicating absence) or 1 (indicating presence).
\(\beta\) is the inverse temperature parameter, which scales the influence of the Hamiltonian (typically set to 1 in psychological applications). The \(\beta\) parameter, which sits outside the Hamiltonian and multiplies its entire value, is called the inverse temperature (Finnemann et al., 2021). The \(\beta\) parameter controls how strongly the network favors certain patterns of symptoms over others.
- A higher \(\beta\) (i.e., lower temperature because \(\beta\) is inverse temperature) means the network prefers some patterns (i.e., the network is more stable).
- A lower \(\beta\) means the network treats many patterns more equally.
\(Z(\tau,\omega)\) is the normalizing constant (i.e., the partition function), ensuring that the probabilities of all possible configurations of \(y\) sum to 1. Without it, the model would produce raw scores that are not on a common scale. The normalizing constant is calculated by summing the unscaled probabilities for every possible node configuration and is defined as \(Z(\tau,\Omega)=\sum_{y \in\{0,1\}^m}\mathit{exp}\!\bigl(-\beta\,H(y;\tau,\Omega)\bigr)\) where
- \(y \in \{0,1\}^m\) is the set of all length-\(m\) on/off patterns, where each \(y_i\) is 0 (absent) or 1 (present). We sum over this set because we need to collect the weight of every possible pattern before we normalize into probabilities.
- \(m\) is the total number of nodes in the network, so there are \(2^m\) possible patterns to consider.
- \(\beta\) is the inverse temperature, scaling how sharply weights favor low-energy patterns.
- \(H(y;\tau,\Omega)\) is the Hamiltonian (pattern’s energy), as defined above.
  This sum ensures the final probabilities form a valid probability distribution (Finnemann et al., 2021).
\(H(y; \tau, \Omega)\) is the Hamiltonian function, which quantifies the total energy (or misalignment) of the system. This function helps represent how likely nodes occur together and is defined as \(H(y; \tau, \Omega)= - \sum_{i=1}^m \tau_i \, y_i \;-\; \sum_{i=2}^m \sum_{j=1}^{i-1} \omega_{ij} \, y_i \, y_j\)
- \(\tau_i\) is the threshold (or bias) parameter for node \(i\), which captures the tendency of node \(i\) to be active (e.g., symptom present) in isolation.
  - A higher \(\tau_i\) means that node \(i\) is more likely to occur even if no other nodes are active.
- \(\omega_{ij}\) is the edge weight between nodes \(i\) and \(j\), reflecting the strength and direction of the statistical association between the two nodes.
  - Positive \(\omega_{ij}\): nodes \(i\) and \(j\) tend to co-occur
  - Negative \(\omega_{ij}\): the presence of one makes the other less likely
- The set of all threshold parameters is denoted \(\tau\), and the set of all edge weights is \(\Omega\).
- \(m\) denotes the total number of nodes/variables.
- The double summation term includes all unique node pairs \((i,j)\) to capture interactions (i.e., how likely two nodes are to co-occur).

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:

Our variables must be binary, meaning the variables must have two levels (e.g., yes or no, present or not present, 0/1).
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.
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.
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.
We assume that there are no unmeasured confounders, and that all variables of the system are included. Omitting variables may lead to spurious edges.
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.

NIRA utilizes the network structure (i.e., the \(\omega\) parameter) and network thresholds (i.e., the \(\tau\) parameter) of an estimated network to simulate how activated a network is (i.e., how many nodes are present) when the threshold of each node is adjusted. This allows researchers to then identify nodes of interest for interventions depending on the nodes with the largest shifts in total activation from baseline. NIRA involves alleviating simulations where nodes are adjusted to be less likely to occur, and aggravating simulations where nodes are adjusted to be more likely to occur. In each simulation type, a node will be perturbed by two standard deviations (Lunansky et al., 2022). That is, its odds will be decreased or increased by two standard deviations in alleviating or aggravating simulations, respectively.

NIRA uses the Metropolis-Hastings algorithm within the IsingSampler R package (Epskamp, 2025) to sample 5000 possible network states (Lunansky et al., 2022) from an estimated network’s structure and thresholds. NIRA begins by first sampling possible network states from the provided network structure and thresholds, with no node thresholds being adjusted yet. Node states from each iteration are stored in a matrix and returned via the IsingSampler function. After 5000 iterations, NIRA then begins perturbing nodes. Specifically, NIRA recalls the IsingSampler function for the node to be intervened upon after the initial 5000 samples of the baseline network. For each node, the threshold is altered by two standard deviations and then 5000 network states are sampled from the adjusted distribution. After another 5000 iterations, the IsingSampler function is called with the next node’s altered threshold and so on until all nodes in the network have been adjusted and had 5000 network states sampled. In total, NIRA results in \(5000 * (m +1)\) samples where \(m\) denotes the number of nodes included in the network.

Once all samples have been drawn, the sum scores of each sample is computed. That is, the sum of all node states (e.g., 0 or 1, -1 or 1) is computed for each sample, and afterwards, the mean of sums is computed for each node’s 5000 samples. This returns the average activated score of the network for the original network and each node’s simulated network. After all average activated scores are computed, the absolute difference is computed between the original network’s average activated score and each node’s simulated network’s score. The node with the largest absolute difference is then considered the most influential node to intervene on.

NIRA is traditionally used with a network estimated across one group of individuals. However, a common research question in clinical psychology is whether networks differ across groups of interest (e.g., across gender, sexual orientation, diagnosis). One popular modeling framework involves estimating multigroup Ising models with increasing parameter constraints. The following explanation of this modeling framework is also taken from my previous tutorial https://rpubs.com/jpstand/1417075.

Researchers 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:

Start with a saturated model that allows all node-to-node relationships to vary freely across groups.
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.
Next, we constrain the thresholds (\(\tau\)) to be equal, ensuring that the baseline likelihood of each symptom occurring is the same in all groups.
Finally, we constrain the inverse 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:

In the prune step, we remove edges that do not significantly improve model fit (p > .05).
In the step-up phase, we test whether re-adding any pruned edges significantly improves model fit. This helps us recover potentially meaningful relationships that were overly conservatively pruned.

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.

Because the inverse temperature parameter (\(\beta\)) can differ across groups, when using NIRA with networks estimated from a multigroup Ising framework, it’s important for the underlying IsingSampler function to account for differences in temperature. As stated earlier, the inverse temperature parameter scales the entire distribution produced by the Ising model, meaning that omitting these differences could lead to inaccurate results when simulating interventions across groups.

This tutorial will examine two networks and simulate alleviating and aggravating interventions for both; one with depressive symptoms across male and female individuals, and one with barriers to mental health care across sexual orientation groups (i.e., heterosexual individuals, sexual minority, individuals unsure of their sexual orientation). We will use the psychonetrics multigroup Ising modeling framework described by Grimes et al., (2025) and others to estimate networks of interest, and an adjusted version of the NIRA framework that accounts for differences in network temperature across groups to simulate alleviating and aggravating interventions (given differences are supported by estimated models).

First, we begin by estimating our network of depressive symptoms across male and female individuals. To summarize, the modeling framework below results in a best-fitting model with parameter constraints on the network structure, thresholds, and inverse temperature. Therefore, in our NIRA simulations, we only use one matrix, vector, and numeric for network structure, thresholds and inverse temperature, respectively.

Step 1: Estimate Networks

rm(list = ls ())
library(psychonetrics)
library(tidyverse)
library(VIM)
library(MLMusingR)
data <- read.csv("symptoms_nettemp.csv")

symptoms <- data %>%
  mutate(Sex = factor(Sex, levels = c(1,2), labels = c("male","female")),
         Sex = relevel(Sex, ref = "male"))
nmiss(symptoms)

#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

Now that we’ve identified our model of interest, we can use NIRA to identify which depressive sypmtoms may be most efficient to target in alleviating and aggravating simulations, respectively.

Step 2: Conduct NIRA Simulations

library(nodeIdentifyR)
library(IsingSampler)
dep_structure <- getmatrix(final_model, "omega")[[1]] #we only need to grab the structure from the first group because it's constrained across all groups. This same logic applies to the following model parameters
dep_thresholds <- as.vector(getmatrix(final_model, "tau")[[1]])
dep_inv_temp <- as.numeric(getmatrix(final_model, "beta")[[1]])

rownames(dep_structure) <- 1:9
colnames(dep_structure) <- 1:9

names(dep_thresholds) <- 1:9

#we adjust the simulateResponses function from NIRA to take beta as an argument
my_simulateResponses <- function (edge_weights, thresholds, beta, perturbation_type, amount_of_SDs_perturbation) 
{
    IsingSamples <- base::vector(mode = "list", length = base::ncol(x = edge_weights) + 
        1)
    base::names(x = IsingSamples) <- base::c(base::colnames(x = edge_weights), 
        "original")
    for (i in 1:base::length(x = IsingSamples)) {
        if (i %in% 1:(base::length(x = IsingSamples) - 1)) {
            perturbation <- thresholds
            if (perturbation_type == "aggravating") {
                perturbation[i] <- thresholds[i] + amount_of_SDs_perturbation * 
                  stats::sd(x = thresholds)
            }
            else if (perturbation_type == "alleviating") {
                perturbation[i] <- thresholds[i] - amount_of_SDs_perturbation * 
                  stats::sd(x = thresholds)
            }
            IsingModelState <- IsingSampler::IsingSampler(n = 5000, 
                graph = edge_weights, thresholds = perturbation, beta = beta, #add in beta argument
                responses = base::c(0L, 1L)) # you should adjust this to be the domain 
                                             # your variables are in (i.e., [0, 1] or [-1, 1])
            IsingSamples[[i]] <- IsingModelState
        }
        if (i == base::length(x = IsingSamples)) {
            IsingModelState <- IsingSampler::IsingSampler(n = 5000, 
                graph = edge_weights, thresholds = thresholds, beta = beta, #add in beta argument
                responses = base::c(0L, 1L)) # you should adjust this to be the domain 
                                             # your variables are in (i.e., [0, 1] or [-1, 1])
            IsingSamples[[i]] <- IsingModelState
        }
    }
    namesOrder <- base::c("original", base::names(x = IsingSamples)[1:base::ncol(x = edge_weights)])
    IsingSamples <- IsingSamples[namesOrder]
    return(IsingSamples)
}
environment(my_simulateResponses) <- asNamespace("nodeIdentifyR")

#now we can simulate our networks
#note that we only simulate once per simulation type because all model parameters were constrained across male and female individuals 

#we do not reset the seed because it was set from the earlier modeling framework
dep_all_sim <- my_simulateResponses(edge_weights = dep_structure,
                             thresholds = dep_thresholds,
                             beta = dep_inv_temp,
                             perturbation_type = "alleviating",
                             amount_of_SDs_perturbation = 2)

dep_agg_sim <- my_simulateResponses(edge_weights = dep_structure,
                             thresholds = dep_thresholds,
                             beta = dep_inv_temp,
                             perturbation_type = "aggravating",
                             amount_of_SDs_perturbation = 2)

We can now calculate the sum scores for alleviating and aggravating interventions and average them across all samples for each node. Afterwards, we can combine the plots of each intervention.

Step 3: Plot NIRA Scores for each symptom

dep_all_plot <- calculateSumScores(dep_all_sim) %>% 
  prepareDFforPlottingAndANOVA() %>% 
  plotSumScores(perturbation_type = "alleviating")

dep_all_plot_color <- dep_all_plot$sumScoresPlot + 
  labs(title = "Alleviating",
       x = "Symptom of which the threshold is altered") + 
  theme_bw(base_size = 16) +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title = element_text(size = 11),
        axis.title.x = element_text(size = 11),
        axis.title.y = element_text(size = 11),
        plot.title = element_text(hjust = .5, face = "bold"))

dep_all_plot_color$layers[[1]] <- geom_line(color = "#377EB8")
dep_all_plot_color$layers[[2]] <- geom_point(color = "#377EB8")
dep_all_plot_color$layers[[3]] <- geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), 
                                      width = 0.15, color = "#377EB8")

dep_agg_plot <- calculateSumScores(dep_agg_sim) %>% 
  prepareDFforPlottingAndANOVA() %>% 
  plotSumScores(perturbation_type = "aggravating")

dep_agg_plot_color <- dep_agg_plot$sumScoresPlot + 
  labs(title = "Aggravating",
       x = "Symptom of which the threshold is altered") + 
  theme_bw(base_size = 16) +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title = element_text(size = 11),
        axis.title.x = element_text(size = 11),
        axis.title.y = element_text(size = 11),
        plot.title = element_text(hjust = .5, face = "bold"))

dep_agg_plot_color$layers[[1]] <- geom_line(color = "#E41A1C")
dep_agg_plot_color$layers[[2]] <- geom_point(color = "#E41A1C")
dep_agg_plot_color$layers[[3]] <- geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), 
                                                width = 0.15, color = "#E41A1C")

dep_comb_plot <- dep_all_plot_color +
  labs(title = "Alleviating and Aggravating Comparison",
       x = "Symptom of which the threshold is altered", y = "Sum Score") + 
  geom_line(data = dep_agg_plot_color$data, color = "#E41A1C", linetype = "dashed") +
  geom_point(data = dep_agg_plot_color$data, color = "#E41A1C") +
  geom_errorbar(data = dep_agg_plot_color$data,
                aes(ymin = ciLower, ymax = ciUpper),
                width = 0.15, color = "#E41A1C") +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title=element_text(size = 15),
        axis.title.x=element_text(size = 15),
        axis.title.y =element_text(size = 15),
        plot.title = element_text(hjust = .5, face = "bold"))

dep_comb_plot

We see that for alleviating interventions (blue line) the best symptom to intervene on (i.e., focus on treating) is the fourth node, which is sleep disturbances. In the aggravating intervention (red-dashed line) the best symptom to intervene on (i.e., prevent from worsening) is the second node, which is anhedonia. However, the confidence intervals for many nodes across both alleviating and aggravating simulations overlap, meaning there may not be a difference in total network activation if we target sleep and anhedonia over, for example, concentration difficulties, in each simulation type, respectively. Let’s do another example with barriers to care as our nodes.

Below we estimate a network with ten barriers to care across three sexual orientation groups (i.e., heterosexual individuals, sexual minority individuals, individuals unsure of their sexual orientation) using the 2023 National Survey on Drug Use and Health (NSDUH). It is important to note that for the sake of computation time this network does not include all barriers to care measured in the 2023 NDSUH. This violates our modeling assumption that the network contains all variables of the system we’re measuring.

ndata <- read.csv("reasons_nettemp.csv")

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")),
         SEXORI = factor(case_when(SEXIDENT22 == "Heterosexual" ~ 1,
                                   SEXIDENT22 %in% c("Gay or Lesbian","Bisexual","Different Term") ~ 2,
                                   SEXIDENT22 == "Unsure" ~ 3),levels = c(1,2,3),labels = c("Heterosexual","Sexual Minority","Unsure")),
         SEXORI = relevel(SEXORI, 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 = "SEXORI") %>% 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)

bm <- rmodel3b

The best-fitting model in this case only has constraints on the network structure and thresholds, meaning that temperature was allowed to fluctuate across groups. Let’s plot the network temperature for each group to confirm.

rbetas_est <- bm@parameters$est[bm@parameters$matrix == "beta"]
rbetas_se <- bm@parameters$se[bm@parameters$matrix == "beta"]
z = qnorm(0.975)
rupperCI <- rbetas_est + (z*rbetas_se)
rlowerCI <- rbetas_est - (z*rbetas_se)
#set up temp plotting
SEXORI<- c("Heterosexual","Sexual Minority", "Unsure")
reasons_data <- data.frame(Sex_Ori = SEXORI,
                           Temperature = 1/rbetas_est,
                           UpperCI = 1/rupperCI,
                           LowerCI = 1/rlowerCI)
#plot barriers to care temp across sexual orientation
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))

Indeed, the temperature varies across groups. Therefore, when conducting the NIRA simulations we should account for differences in temperature across each group.

tau_mat <- getmatrix(bm, "tau")[[1]]
tau_vec <- drop(tau_mat)
names(tau_vec) <- 1:10
omega_mat <- getmatrix(bm,"omega")[[1]]
rownames(omega_mat) <- 1:10
colnames(omega_mat) <- 1:10
het_beta <- as.numeric(getmatrix(bm, "beta")[[1]])
sm_beta <- as.numeric(getmatrix(bm, "beta")[[2]])
unsure_beta <- as.numeric(getmatrix(bm, "beta")[[3]])

#alleviating simulation
het_all_sim <- my_simulateResponses(edge_weights = omega_mat,
                             thresholds = tau_vec,
                             beta = het_beta,
                             perturbation_type = "alleviating",
                             amount_of_SDs_perturbation = 2)
sm_all_sim <- my_simulateResponses(edge_weights = omega_mat,
                             thresholds = tau_vec,
                             beta = sm_beta,
                             perturbation_type = "alleviating",
                             amount_of_SDs_perturbation = 2)
unsure_all_sim <- my_simulateResponses(edge_weights = omega_mat,
                             thresholds = tau_vec,
                             beta = unsure_beta,
                             perturbation_type = "alleviating",
                             amount_of_SDs_perturbation = 2)
#aggravating simulation
het_agg_sim <- my_simulateResponses(edge_weights = omega_mat,
                             thresholds = tau_vec,
                             beta = het_beta,
                             perturbation_type = "aggravating",
                             amount_of_SDs_perturbation = 2)
sm_agg_sim <- my_simulateResponses(edge_weights = omega_mat,
                             thresholds = tau_vec,
                             beta = sm_beta,
                             perturbation_type = "aggravating",
                             amount_of_SDs_perturbation = 2)
unsure_agg_sim <- my_simulateResponses(edge_weights = omega_mat,
                             thresholds = tau_vec,
                             beta = unsure_beta,
                             perturbation_type = "aggravating",
                             amount_of_SDs_perturbation = 2)


het_all_plot <- calculateSumScores(het_all_sim) %>% 
  prepareDFforPlottingAndANOVA() %>% 
  plotSumScores(perturbation_type = "alleviating")

het_all_plot_color <- het_all_plot$sumScoresPlot + 
  labs(title = "Alleviating",
       x = "Symptom of which the threshold is altered") + 
  theme_bw(base_size = 16) +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title = element_text(size = 11),
        axis.title.x = element_text(size = 11),
        axis.title.y = element_text(size = 11),
        plot.title = element_text(hjust = .5, face = "bold"))

het_all_plot_color$layers[[1]] <- geom_line(color = "#377EB8")
het_all_plot_color$layers[[2]] <- geom_point(color = "#377EB8")
het_all_plot_color$layers[[3]] <- geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), 
                                                width = 0.15, color = "#377EB8")

het_agg_plot <- calculateSumScores(het_agg_sim) %>% 
  prepareDFforPlottingAndANOVA() %>% 
  plotSumScores(perturbation_type = "aggravating")

het_agg_plot_color <- het_agg_plot$sumScoresPlot + 
  labs(title = "Aggravating",
       x = "Symptom of which the threshold is altered") + 
  theme_bw(base_size = 16) +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title = element_text(size = 11),
        axis.title.x = element_text(size = 11),
        axis.title.y = element_text(size = 11),
        plot.title = element_text(hjust = .5, face = "bold"))

het_agg_plot_color$layers[[1]] <- geom_line(color = "#E41A1C")
het_agg_plot_color$layers[[2]] <- geom_point(color = "#E41A1C")
het_agg_plot_color$layers[[3]] <- geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), 
                                                width = 0.15, color = "#E41A1C")

het_comb_plot <- het_all_plot_color +
  labs(title = "Alleviating and Aggravating Comparison for Heterosexual Individuals",
       x = "Barrier of which the threshold is altered", y = "Sum Score") + 
  geom_line(data = het_agg_plot_color$data, color = "#E41A1C", linetype = "dashed") +
  geom_point(data = het_agg_plot_color$data, color = "#E41A1C") +
  geom_errorbar(data = het_agg_plot_color$data,
                aes(ymin = ciLower, ymax = ciUpper),
                width = 0.15, color = "#E41A1C") +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title=element_text(size = 15),
        axis.title.x=element_text(size = 15),
        axis.title.y =element_text(size = 15),
        plot.title = element_text(hjust = .5, face = "bold", size = 15))

sm_all_plot <- calculateSumScores(sm_all_sim) %>% 
  prepareDFforPlottingAndANOVA() %>% 
  plotSumScores(perturbation_type = "alleviating")

sm_all_plot_color <- sm_all_plot$sumScoresPlot + 
  labs(title = "Alleviating",
       x = "Barrier of which the threshold is altered") + 
  theme_bw(base_size = 16) +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title = element_text(size = 11),
        axis.title.x = element_text(size = 11),
        axis.title.y = element_text(size = 11),
        plot.title = element_text(hjust = .5, face = "bold"))

sm_all_plot_color$layers[[1]] <- geom_line(color = "#377EB8")
sm_all_plot_color$layers[[2]] <- geom_point(color = "#377EB8")
sm_all_plot_color$layers[[3]] <- geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), 
                                               width = 0.15, color = "#377EB8")

sm_agg_plot <- calculateSumScores(sm_agg_sim) %>% 
  prepareDFforPlottingAndANOVA() %>% 
  plotSumScores(perturbation_type = "aggravating")

sm_agg_plot_color <- sm_agg_plot$sumScoresPlot + 
  labs(title = "Aggravating",
       x = "Barrier of which the threshold is altered") + 
  theme_bw(base_size = 16) +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title = element_text(size = 11),
        axis.title.x = element_text(size = 11),
        axis.title.y = element_text(size = 11),
        plot.title = element_text(hjust = .5, face = "bold"))

sm_agg_plot_color$layers[[1]] <- geom_line(color = "#E41A1C")
sm_agg_plot_color$layers[[2]] <- geom_point(color = "#E41A1C")
sm_agg_plot_color$layers[[3]] <- geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), 
                                               width = 0.15, color = "#E41A1C")

sm_comb_plot <- sm_all_plot_color +
  labs(title = "Alleviating and Aggravating Comparison for Sexual Minority Individuals",
       x = "Barrier of which the threshold is altered", y = "Sum Score") + 
  geom_line(data = sm_agg_plot_color$data, color = "#E41A1C", linetype = "dashed") +
  geom_point(data = sm_agg_plot_color$data, color = "#E41A1C") +
  geom_errorbar(data = sm_agg_plot_color$data,
                aes(ymin = ciLower, ymax = ciUpper),
                width = 0.15, color = "#E41A1C") +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title=element_text(size = 15),
        axis.title.x=element_text(size = 15),
        axis.title.y =element_text(size = 15),
        plot.title = element_text(hjust = .5, face = "bold", size = 15))


unsure_all_plot <- calculateSumScores(unsure_all_sim) %>% 
  prepareDFforPlottingAndANOVA() %>% 
  plotSumScores(perturbation_type = "alleviating")

unsure_all_plot_color <- unsure_all_plot$sumScoresPlot + 
  labs(title = "Alleviating",
       x = "Barrier of which the threshold is altered") + 
  theme_bw(base_size = 16) +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title = element_text(size = 11),
        axis.title.x = element_text(size = 11),
        axis.title.y = element_text(size = 11),
        plot.title = element_text(hjust = .5, face = "bold"))

unsure_all_plot_color$layers[[1]] <- geom_line(color = "#377EB8")
unsure_all_plot_color$layers[[2]] <- geom_point(color = "#377EB8")
unsure_all_plot_color$layers[[3]] <- geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), 
                                                   width = 0.15, color = "#377EB8")

unsure_agg_plot <- calculateSumScores(unsure_agg_sim) %>% 
  prepareDFforPlottingAndANOVA() %>% 
  plotSumScores(perturbation_type = "aggravating")

unsure_agg_plot_color <- unsure_agg_plot$sumScoresPlot + 
  labs(title = "Aggravating",
       x = "Barrier of which the threshold is altered") + 
  theme_bw(base_size = 16) +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title = element_text(size = 11),
        axis.title.x = element_text(size = 11),
        axis.title.y = element_text(size = 11),
        plot.title = element_text(hjust = .5, face = "bold"))

unsure_agg_plot_color$layers[[1]] <- geom_line(color = "#E41A1C")
unsure_agg_plot_color$layers[[2]] <- geom_point(color = "#E41A1C")
unsure_agg_plot_color$layers[[3]] <- geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), 
                                                   width = 0.15, color = "#E41A1C")

unsure_comb_plot <- unsure_all_plot_color +
  labs(title = "Alleviating and Aggravating Comparison for Unsure Individuals",
       x = "Barrier of which the threshold is altered", y = "Sum Score") + 
  geom_line(data = unsure_agg_plot_color$data, color = "#E41A1C", linetype = "dashed") +
  geom_point(data = unsure_agg_plot_color$data, color = "#E41A1C") +
  geom_errorbar(data = unsure_agg_plot_color$data,
                aes(ymin = ciLower, ymax = ciUpper),
                width = 0.15, color = "#E41A1C") +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1), 
        axis.title=element_text(size = 15),
        axis.title.x=element_text(size = 15),
        axis.title.y =element_text(size = 15),
        plot.title = element_text(hjust = .5, face = "bold", size = 15))

het_comb_plot

sm_comb_plot

unsure_comb_plot

Here we see that accounting for differences in network temperature does lead to interesting results across sexual orientation groups. While the specific barriers to care of interest for alleviating and aggravating simulations remain the same across groups, the pattern of which barriers to care should be prioritized changes across groups. In the alleviating intervention (blue line), across all groups it seems that intervening on the first, fourth, and tenth nodes (i.e., cost barriers, not knowing where to find treatment, believing one can handle problems alone) will have an equal effect on total network activation across sexual orientation groups. However, in aggravating interventions (red-dashed line), targeting the third node (i.e., insurance not covering enough) appears to be the most effective barrier to prevent from becoming more rampant for all sexual orientation groups. Of course, this network omits many key barriers, so these interpretations do not hold any true value. However, we do see that accounting for network temperature differences is important when comparing differences in Ising model networks across groups.

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 (2025). IsingSampler: Sampling Methods and Distribution Functions for the Ising Model. R package version 0.2.4, https://github.com/SachaEpskamp/IsingSampler.

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

Lunansky, G., Naberman, J., Van Borkulo, C. D., Chen, C., Wang, L., & Borsboom, D. (2022). Intervening on psychopathology networks: Evaluating intervention targets through simulations. Methods, 204, 29–37. https://doi.org/10.1016/j.ymeth.2021.11.006

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

Network Temperature Informed Simulations of Multigroup Ising Models Using the NodeIdentifyR Algorithm

Philipp Standring

2026-04-17