A Note on Ising Network Analysis with Missing Data

Abstract

The Ising model has become a popular psychometric model for analyzing item response data. The statistical inference of the Ising model is typically carried out via a pseudo-likelihood, as the standard likelihood approach suffers from a high computational cost when there are many variables (i.e., items). Unfortunately, the presence of missing values can hinder the use of pseudo-likelihood, and a listwise deletion approach for missing data treatment may introduce a substantial bias into the estimation and sometimes yield misleading interpretations. This paper proposes a conditional Bayesian framework for Ising network analysis with missing data, which integrates a pseudo-likelihood approach with iterative data imputation. An asymptotic theory is established for the method. Furthermore, a computationally efficient Pólya–Gamma data augmentation procedure is proposed to streamline the sampling of model parameters. The method’s performance is shown through simulations and a real-world application to data on major depressive and generalized anxiety disorders from the National Epidemiological Survey on Alcohol and Related Conditions (NESARC).

Recent years have witnessed the emergence of network psychometrics (Borsboom, 2008; Marsman and Rhemtulla, 2022; van der Maas et al., 2006), a family of statistical graphical models and related inference procedures, for analyzing and interpreting the dependence structure in psychometric data. These models embed psychometric items as nodes in an undirected or directed network (i.e., graph) and visualize their interrelationships through the network edges, which represent certain probabilistic conditional dependencies. Network psychometric methods concern the learning of the network structure. They have been developed under various settings, including undirected graphical models for cross-sectional data (Burger et al., 2022; Epskamp et al., 2018), directed networks for longitudinal data (Borsboom et al., 2021; Gile and Handcock, 2017; Ryan et al., 2022), and extended networks with latent variables for time series data or panel data (Epskamp, 2020). These methods have received wide applications in education (Koponen et al., 2019; Siew, 2020; Simon de Blas et al., 2021; Sweet et al., 2013; Willcox and Huang, 2017), psychology (Borsboom et al., 2021; Burgess and Hitch, 1999; Epskamp et al., 2018; Fried et al., 2017; Van Der Maas et al., 2017), and health sciences (Brunson and Laubenbacher, 2018; Kohler et al., 2022; Luke and Harris, 2007; Mkhitaryan et al., 2019).

Analyzing cross-sectional binary item response data with the Ising model (Ising, 1925) is common in network psychometric analysis. This analysis is typically performed based on a conditional likelihood (Besag, 1974) because the standard likelihood function is computationally infeasible when involving many variables. In this direction, Bayesian and frequentist methods have been developed, where sparsity-inducing priors or penalties are combined with the conditional likelihood for learning a sparse network structure (Epskamp and Fried, 2018; Li et al., 2019; Marsman et al., 2022; Mazumder and Hastie, 2012; Van Borkulo et al., 2014; Yuan and Lin, 2007). Besides, the Ising model is shown to be closely related to item response theory (IRT) models (Anderson and Yu, 2007; Holland, 1990). The log-multiplicative association models (Anderson and Yu, 2007), which are special cases of the Ising model, can be used as item response theory models and yield very similar results as IRT models. Furthermore, the Ising model and the conditional likelihood have been used for modeling the local dependence structure in locally dependent IRT models (Chen et al., 2018; Ip, 2002).

Due to its construction, the conditional likelihood does not naturally handle data with missing values, despite the omnipresence of missing data in psychometric data. To deal with missing values in an Ising network analysis, listwise deletion (Fried et al., 2020; Haslbeck and Fried, 2017) and single imputation (e.g., Huisman, 2009, Armour et al., 2017, Lin et al., 2020) are typically performed, which arguably may not be the best practice. In particular, it is well known that listwise deletion is statistically inefficient and requires the missing completely at random (MCAR) assumption (Little and Rubin, 2019) to ensure consistent estimation. Moreover, a naïve imputation procedure, such as mode imputation, likely introduces bias into parameter estimation. A sophisticated imputation procedure must be developed to ensure statistical validity and computational efficiency.

In this note, we propose an iterative procedure for learning an Ising network. The proposed procedure combines iterative imputation via full conditional specification (FCS; Liu et al., 2014, van Buuren, 2018), and Bayesian estimation of the Ising network. We show that the FCS leads to estimation consistency when the conditional models are chosen to take logistic forms. In terms of computation, we propose a joint Pólya–Gamma augmentation procedure by extending the Pólya–Gamma augmentation procedure for logistic regression (Polson et al., 2013). It allows us to efficiently sample parameters of the Ising model. Simulations are conducted to compare the proposed procedure with estimations based on the listwise deletion and single imputation. Finally, the proposed procedure and a complete-case analysis are applied to study the network of major depressive disorder (MDD) and generalized anxiety disorders (GAD) based on data from the National Epidemiological Survey on Alcohol and Related Conditions (NESARC; Grant et al., 2003). In this analysis, data missingness is mainly due to two screening items for GAD. That is, a respondent’s responses to the rest of the MDD items are missing if they answered “no” to both screening items. This missing mechanism is Missing at Random (MAR; Little & Rubin, 2019). The complete-case analysis of missing data caused by screening items is known to be problematic in the literature of network psychometrics (Borsboom et al., 2017; McBride et al., 2023). Our Bayesian estimate of the edge coefficient between the two screening items is negative based on the complete cases, which can be seen as a result of Berkson’s paradox (De Ron et al., 2021). In contrast, the proposed method makes use of all the observed data entries and obtains a positive estimate of this edge coefficient. An identifiability result about the Ising model under this special missing data setting in Appendix, the item content, and a simulation study mimicking this setting suggest that the result given by the proposed method is sensible. A GitHub repository for the proposed algorithm is available at https://github.com/slzhang-fd/IsingNetMissing-replication.

1. Proposed Method

1.1. Ising Model

Consider a respondent answering J binary items. Let $Y={(Y_1,\dots ,Y_J)}^{\top }\in {\{0,1\}}^J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}= (Y_1,\ldots , Y_J)^\top \in \{0,1\}^J$$\end{document} be a binary random vector representing the respondent’s responses. We say $Y$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}$$\end{document} follows an Ising model if its probability mass function satisfies

(1) $\begin{array}{c}P(Y=y\mid S)=\frac{1}{c\left(S\right)}exp\left(\frac{1}{2},y^{\top },S,y\right)=\frac{1}{c\left(S\right)}exp\left(\sum _{j=1}^J,s_{\mathrm{jj}},y_j,/,2,+,\sum _{j=1}^{J-1},\sum _{k=j+1}^J,s_{\mathrm{jk}},y_j,y_k\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P({{\textbf{Y}}}= {{\varvec{y}}}\mid {{\textbf{S}}}) = \frac{1}{c({{\textbf{S}}})}\exp \left[ \frac{1}{2}{{\varvec{y}}}^\top {{\textbf{S}}}{{\varvec{y}}}\right] =\frac{1}{c({{\textbf{S}}})}\exp \left[ \sum _{j=1}^J s_{jj}y_{j}/2+\sum _{j=1}^{J-1}\sum _{k=j+1}^Js_{jk}y_{j}y_k\right] , \end{aligned}$$\end{document}

where $S={\left(s_{\mathrm{ij}}\right)}_{J\times J}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}= (s_{ij})_{J\times J}$$\end{document} is a J by J symmetric matrix that contains parameters of the Ising model and

$\begin{array}{c}c\left(S\right)=\sum _{y\in {\{0,1\}}^J}exp\left(\sum _{j=1}^J,s_{\mathrm{jj}},y_j,/,2,+,\sum _{j=1}^{J-1},\sum _{k=j+1}^J,s_{\mathrm{jk}},y_j,y_k\right)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} c({{\textbf{S}}}) = \sum _{{{\varvec{y}}}\in \{0,1\}^J}\exp \left[ \sum _{j=1}^J s_{jj}y_{j}/2+\sum _{j=1}^{J-1}\sum _{k=j+1}^Js_{jk}y_{j}y_k\right] \end{aligned}$$\end{document}

is a normalizing constant. The parameter matrix $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} encodes a network with the J items being the nodes. More specifically, an edge is present between nodes i and j if and only if the corresponding entry $s_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{ij}$$\end{document} is nonzero. If an edge exists between nodes i and j, then $Y_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_i$$\end{document} and $Y_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_j$$\end{document} are conditionally dependent given the rest of the variables. Otherwise, the two variables are conditionally independent.

In Ising network analysis, the goal is to estimate the parameter matrix $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} . The standard likelihood function is computationally intensive when J is large, as it requires computing a normalizing constant $c\left(S\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c({{\textbf{S}}})$$\end{document} which involves a summation of all the $2^J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2^J$$\end{document} response patterns. To address this computational issue, Besag (1975) proposed a conditional likelihood which is obtained by aggregating the conditional distributions of $Y_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_j$$\end{document} given $Y_{-j}={(Y_1,\dots ,Y_{j-1},Y_{j+1},\dots ,Y_J)}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}_{-j} = (Y_1,\ldots , Y_{j-1}, Y_{j+1},\ldots , Y_J)^\top $$\end{document} , for $j=1,\dots ,J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots , J$$\end{document} , where the conditional distribution of $Y_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_j$$\end{document} given $Y_{-j}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}_{-j}$$\end{document} takes a logistic regression form. More precisely, the conditional likelihood with one observation $y$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{y}}}$$\end{document} is defined as

(2) $\begin{array}{c}{p}^{\ast }(y\mid S)=\prod _{j=1}^{J}p(y_j\mid y_{-j},S)=\prod _{j=1}^J\frac{exp\left((s_{\mathrm{jj}}/2+{\sum }_{k\ne j}{s}_{\mathrm{jk}}{y}_k),y_j\right)}{1+exp\left(s_{\mathrm{jj}},/,2,+,{\sum }_{k\ne j},s_{\mathrm{jk}},y_k\right)}.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p^*({{\varvec{y}}}\mid {{\textbf{S}}}) = \prod _{j=1}^J p(y_j\mid {{\varvec{y}}}_{-j},{{\textbf{S}}})= \prod _{j=1}^J\frac{\exp \left[ (s_{jj}/2 + \sum _{k\ne j} s_{jk}y_k)y_j\right] }{1+\exp \left( s_{jj}/2 + \sum _{k\ne j} s_{jk}y_k\right) }. \end{aligned}$$\end{document}

A disadvantage of the conditional likelihood is that it requires a fully observed dataset because missing values cannot be straightforwardly marginalized out from (2). In what follows, we discuss how missing data can be treated in the conditional likelihood.

1.2. Proposed Method

Consider a dataset with N observations. Let ${\Omega }_j\subset \{1,\dots ,N\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{j} \subset \{1,\ldots , N\}$$\end{document} be the subset of observations whose data on item j are missing. For each observation i and item j, $y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}$$\end{document} denotes the observed response if $i\notin {\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i \notin \Omega _j$$\end{document} , and otherwise, $y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}$$\end{document} is missing. Thus, the observed data include ${\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _j$$\end{document} and $y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}$$\end{document} , for $i\in \{1,\dots ,N\}\backslash {\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i\in \{1,\ldots , N\}{\setminus }\Omega _j$$\end{document} and $j=1,\dots ,J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j = 1,\ldots , J$$\end{document} .

The proposed procedure iterates between two steps—(1) imputing the missing values of $y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}$$\end{document} for $i\in {\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i \in \Omega _j$$\end{document} , $j=1,\dots ,J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots , J$$\end{document} , achieved via a full conditional specification, and (2) sampling the posterior distribution of $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} given the most recently imputed data. Let t be the current iteration number. Further, let $y_i^{(t-1)}={(y_{i1}^{(t-1)},\dots ,y_{\mathrm{iJ}}^{(t-1)})}^{\top },i=1,\dots ,N$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{y}}}^{(t-1)}_i = (y_{i1}^{(t-1)},\ldots , y_{iJ}^{(t-1)})^\top , i=1,\ldots , N$$\end{document} , be imputed data from the previous iteration, where $y_{\mathrm{ij}}^{(t-1)}=y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y^{(t-1)}_{ij} = y_{ij}$$\end{document} for $i\notin {\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i \notin \Omega _j$$\end{document} and $y_{\mathrm{ij}}^{(t-1)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y^{(t-1)}_{ij}$$\end{document} is imputed in the $(t-1)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(t-1)$$\end{document} th iteration for $i\in {\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i \in \Omega _j$$\end{document} . For the tth iteration, the imputation and sampling steps are described as follows.

Imputation We initialize the imputation in the tth iteration with the previously imputed data set $(y_{i1}^{(t-1)},\dots ,y_{\mathrm{iJ}}^{(t-1)})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y_{i1}^{(t-1)},\ldots ,y_{iJ}^{(t-1)})$$\end{document} , $i=1,\dots ,N$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,\ldots , N$$\end{document} . Then, we run a loop over all the items, $j=1,\dots ,J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j = 1,\ldots , J$$\end{document} . In step j of the loop, we impute $y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}$$\end{document} for all $i\in {\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i \in \Omega _j$$\end{document} , given the most recently imputed data $(y_{i1}^{\left(t\right)},\dots ,y_{i,j-1}^{\left(t\right)},y_{\mathrm{ij}}^{(t-1)},\dots ,y_{\mathrm{iJ}}^{(t-1)})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y_{i1}^{(t)},\ldots ,y_{i,j-1}^{(t)},y_{ij}^{(t-1)},\ldots ,y_{iJ}^{(t-1)})$$\end{document} , $i=1,\dots ,N$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i = 1,\ldots , N$$\end{document} . We then obtain an updated data set $(y_{i1}^{\left(t\right)},\dots ,y_{i,j}^{\left(t\right)},y_{i,j+1}^{(t-1)},\dots ,y_{\mathrm{iJ}}^{(t-1)})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y_{i1}^{(t)},\ldots ,y_{i,j}^{(t)},y_{i,j+1}^{(t-1)},\ldots ,y_{iJ}^{(t-1)})$$\end{document} by incorporating the newly imputed values for $y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}$$\end{document} .

The imputation of each variable j is based on the conditional distribution of $Y_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_j$$\end{document} given $Y_{-j}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}_{-j}$$\end{document} . Under the Ising model, this conditional distribution takes a logistic regression form. For computational reasons to be discussed in the sequel, we introduce an auxiliary parameter vector ${\beta }_j={({\beta }_{j1},\dots ,{\beta }_{\mathrm{jJ}})}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\beta }}_j = (\beta _{j1},\ldots , \beta _{jJ})^\top $$\end{document} as coefficients in the logistic regression, instead of directly using $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} from the previous iteration to sample the missing $y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}$$\end{document} s. Unlike the constraint of $s_{\mathrm{ij}}=s_{\mathrm{ji}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{ij} = s_{ji}$$\end{document} in the symmetric matrix $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} , no constraints are imposed on ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_j$$\end{document} , $j=1,\dots ,J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots , J$$\end{document} , which makes the sampling computationally efficient; see discussions in Sect. 1.4. The imputation of variable j consists of the following two steps:

1. 1. Sample auxiliary parameter vector ${\beta }_j^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\beta }}_j^{(t)}$$\end{document} from the posterior distribution

   (3) $\begin{array}{c}{p}^{(t,j)}\left({\beta }_j\right)\propto {\pi }_j\left({\beta }_j\right)\prod _{i=1}^N\frac{exp\left(({\beta }_{\mathrm{jj}}/2+{\sum }_{k\ne j}{\beta }_{\mathrm{jk}}{y}_{\mathrm{ik}}^{(t,j-1)}),y_{\mathrm{ij}}^{(t,j-1)}\right)}{1+exp({\beta }_{\mathrm{jj}}/2+{\sum }_{k\ne j}{\beta }_{\mathrm{jk}}{y}_{\mathrm{ik}}^{(t,j-1)})},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p^{(t,j)}({\varvec{\beta }}_j) \propto \pi _j({\varvec{\beta }}_j)\prod _{i=1}^N \frac{\exp \left[ (\beta _{jj}/2+\sum _{k\ne j}\beta _{jk}y_{ik}^{(t,j-1)})y_{ij}^{(t,j-1)}\right] }{1+\exp (\beta _{jj}/2+\sum _{k\ne j}\beta _{jk}y_{ik}^{(t,j-1)})}, \end{aligned}$$\end{document}
   where ${\pi }_j\left({\beta }_j\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi _j({\varvec{\beta }}_j)$$\end{document} is the prior distribution for the auxiliary parameters ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\beta }}_j$$\end{document} .

2. 2. Sample $y_{\mathrm{ij}}^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}^{(t)}$$\end{document} for each $i\in {\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i \in \Omega _j$$\end{document} from a Bernoulli distribution with success probability

   (4) $\begin{array}{c}\frac{exp({\beta }_{\mathrm{jj}}^{\left(t\right)}/2+{\sum }_{k\ne j}{\beta }_{\mathrm{jk}}^{\left(t\right)}{y}_{\mathrm{ik}}^{(t,j-1)})}{1+exp({\beta }_{\mathrm{jj}}^{\left(t\right)}/2+{\sum }_{k\ne j}{\beta }_{\mathrm{jk}}^{\left(t\right)}{y}_{\mathrm{ik}}^{(t,j-1)})}.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\exp (\beta _{jj}^{(t)}/2+\sum _{k\ne j}\beta _{jk}^{(t)}y_{ik}^{(t,j-1)})}{1+\exp (\beta _{jj}^{(t)}/2+\sum _{k\ne j}\beta _{jk}^{(t)}y_{ik}^{(t,j-1)})}. \end{aligned}$$\end{document}

After these two steps, we obtain $(y_{i1}^{\left(t\right)},\dots ,y_{i,j}^{\left(t\right)},y_{i,j+1}^{(t-1)},\dots ,y_{\mathrm{iJ}}^{(t-1)})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y_{i1}^{(t)},\ldots ,y_{i,j}^{(t)},y_{i,j+1}^{(t-1)},\ldots ,y_{iJ}^{(t-1)})$$\end{document} by incorporating the newly imputed values for $y_{\mathrm{ij}},i\in {\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}, i\in \Omega _j$$\end{document} . We emphasize that only the missing values are updated. For $i\notin {\Omega }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i \notin \Omega _j$$\end{document} , $y_{\mathrm{ij}}^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}^{(t)}$$\end{document} is always the observed value of $y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ij}$$\end{document} . After the loop over all the items, we have the imputed data set $(y_{i1}^{\left(t\right)},\dots ,y_{\mathrm{iJ}}^{\left(t\right)})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y_{i1}^{(t)},\ldots ,y_{iJ}^{(t)})$$\end{document} as the output from this imputation step.

Sampling $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} Given the most recently imputed data $y_i^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{y}}}_i^{(t)}$$\end{document} , $i=1,\dots ,N$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,\ldots , N$$\end{document} , update $S^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}^{(t)}$$\end{document} by sampling from the pseudo-posterior distribution

(5) $\begin{array}{c}p(S\mid y_1^{\left(t\right)},\dots ,y_N^{\left(t\right)})\propto \pi \left(S\right)\prod _{i=1}^N{p}^{\ast }(y_i^{\left(t\right)}\mid S),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p({{\textbf{S}}}\mid {{\varvec{y}}}_1^{(t)},\ldots ,{{\varvec{y}}}_N^{(t)})\propto \pi ({{\textbf{S}}})\prod _{i=1}^N p^*({{\varvec{y}}}_i^{(t)} \mid {{\textbf{S}}}), \end{aligned}$$\end{document}

where $\pi \left(S\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi ({{\textbf{S}}})$$\end{document} is the prior distribution for the Ising parameter matrix $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} and recall that ${\prod }_{i=1}^N{p}^{\ast }(y_i^{\left(t\right)}\mid S)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\prod _{i=1}^N p^*({{\varvec{y}}}_i^{(t)} \mid {{\textbf{S}}})$$\end{document} is the conditional likelihood.

Figure 1 Flowchart of the updating rule for the proposed method.

Algorithm 1 Ising Network Analysis with Iterative Imputation

Figure 1 visualizes the steps performed in the proposed method. Note that it is unnecessary to sample the parameter matrix $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} during the burn-in period and in every iteration after the burn-in period; thus, we employ a thinning step after the burn-in period. This is done to both decrease computational cost and reduce the auto-correlation in the imputed data. Moreover, we outline the proposed algorithm in Algorithm 1. The final estimate of $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} is obtained by averaging all the $S^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}^{(t)}$$\end{document} obtained after the burn-in period. The computational details, including the sampling of auxiliary parameters and Ising parameter matrix and discussions of the computational complexity, are given in Sect. 1.4.

We remark that our method imputes the missing variables one by one for each observation. This method is chosen because simultaneously imputing all the missing variables is typically computationally infeasible, especially when some observation units have many missing values. Simultaneous imputation requires evaluating the joint distribution of the missing variables given the observed ones, whose computational complexity grows exponentially with the number of missing values. In contrast, the proposed method is based on unidimensional conditional distributions, which is computationally more feasible. We also note that the proposed method has several variants that should also work well. These variants are discussed in Sect. 3.

1.3. Statistical Consistency

As our method is not a standard Bayesian inference procedure, we provide an asymptotic theory under the frequentist setting to justify its validity. In particular, we show that the $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} parameter sampled from the pseudo-posterior distribution converges to the true parameter $S_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}_0$$\end{document} , under the assumptions that the Ising model is correctly specified and the data are MAR.

Consider one observation with a complete data vector $Y={(Y_1,\dots ,Y_J)}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}= (Y_1,\ldots , Y_J)^\top $$\end{document} . Further, let $Z={(Z_1,\dots ,Z_J)}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Z}}}= (Z_1,\ldots , Z_J)^\top $$\end{document} be a vector of missing indicators, where $Z_{\mathrm{ij}}=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_{ij} = 1$$\end{document} if $Y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_{ij}$$\end{document} is observed and $Z_{\mathrm{ij}}=0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_{ij} = 0$$\end{document} otherwise. We further let $Y_{\mathrm{obs}}=\{Y_j:Z_j=1,j=1,\dots ,J\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}_{obs} = \{Y_j: Z_j = 1, j = 1,\ldots , J\}$$\end{document} and $Y_{\mathrm{mis}}=\{Y_j:Z_j=0,j=1,\dots ,J\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}_{mis} = \{Y_j: Z_j = 0, j = 1,\ldots , J\}$$\end{document} be the observed and missing entries of $Y$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}$$\end{document} , respectively. Consider the joint distribution of observable data $(Y_{\mathrm{obs}},Z)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({{\textbf{Y}}}_{obs}, {{\textbf{Z}}})$$\end{document} , taking the form

(6) $\begin{array}{c}P(Y_{\mathrm{obs}}=y_{\mathrm{obs}},Z=z\mid S,\varphi )=\sum _{y_j:z_j=0}\left(exp,\left(y^{\top },S,y,/,2\right),/,c\left(S\right)\right)q(z\mid y,\varphi ),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P({{\textbf{Y}}}_{obs} = {{\varvec{y}}}_{obs}, {{\textbf{Z}}}= {\textbf{z}}\mid {{\textbf{S}}}, \varvec{\phi }) = \sum _{y_j:z_j=0} \left( \exp \left( {{\varvec{y}}}^\top {{\textbf{S}}}{{\varvec{y}}}/2\right) /{c({{\textbf{S}}})} \right) q({\textbf{z}} \mid {{\varvec{y}}},\varvec{\phi }), \end{aligned}$$\end{document}

where $exp\left(y^{\top },S,y,/,2\right)/c\left(S\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\exp \left( {{\varvec{y}}}^\top {{\textbf{S}}}{{\varvec{y}}}/2\right) /{c({{\textbf{S}}})}$$\end{document} is the distribution of $Y=y$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}= {{\varvec{y}}}$$\end{document} under the Ising model, $q(z\mid y,\varphi )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q({\textbf{z}} \mid {{\varvec{y}}},\varvec{\phi })$$\end{document} denotes the conditional probability of $Z=z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Z}}}= {\textbf{z}}$$\end{document} given $Y=y$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{Y}}}= {{\varvec{y}}}$$\end{document} , and $\varphi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\phi }$$\end{document} denotes the unknown parameters of this distribution. The MAR assumption, also known as the ignorable missingness assumption, means that the conditional distribution $q(z\mid y,\varphi )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q({\textbf{z}} \mid {{\varvec{y}}},\varvec{\phi })$$\end{document} depends on $y$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{y}}}$$\end{document} only through the observed entries, i.e., $q(z\mid y,\varphi )=q(z\mid y_{\mathrm{obs}},\varphi )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q({\textbf{z}} \mid {{\varvec{y}}},\varvec{\phi }) = q(\textbf{z}\mid {{\varvec{y}}}_{obs},\varvec{\phi })$$\end{document} . In that case, (6) can be factorized as

(7) $\begin{array}{c}P(Y_{\mathrm{obs}}=y_{\mathrm{obs}},Z=z\mid S,\varphi )=q(z\mid y_{\mathrm{obs}},\varphi )\times \left(\sum _{y_j:z_j=0},exp,\left(y^{\top },S,y,/,2\right),/,c\left(S\right)\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P({{\textbf{Y}}}_{obs} = {{\varvec{y}}}_{obs}, {{\textbf{Z}}}= {\textbf{z}}\mid {{\textbf{S}}}, \varvec{\phi }) = q({\textbf{z}} \mid {{\varvec{y}}}_{obs},\varvec{\phi }) \times \left( \sum _{y_j:z_j=0} \exp \left( {{\varvec{y}}}^\top {{\textbf{S}}}{{\varvec{y}}}/2\right) /{c({{\textbf{S}}})} \right) . \end{aligned}$$\end{document}

Consequently, the inference of $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} does not depend on the unknown distribution $q(z\mid y,\varphi )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q({\textbf{z}} \mid {{\varvec{y}}},\varvec{\phi })$$\end{document} .

As shown in Liu et al. (2014), the MAR assumption, together with additional regularity conditions, ensures that the iterative imputation of the missing responses converges to the imputation distribution under a standard Bayesian procedure as the number of iterations and the sample size N go to infinity. A key to this convergence result is the compatibility of the conditional models in the imputation step—the logistic regression models are compatible with the Ising model as a joint distribution. The validity of the imputed samples further ensures the consistency of the estimated Ising parameter matrix. We summarize this result in Theorem 1.

Theorem 1

Assume the following assumptions hold: (1) The Markov chain for missing data, generated by the iterative imputation algorithm (Algorithm 1), is positive Harris recurrent and thus admits a unique stationary distribution; (2) the missing data process is ignorable; and (3) A regularity condition holds for prior distributions of Ising model parameters and auxiliary parameters, as detailed in the supplementary material. Let ${\pi }_N^{\ast }\left(S\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi _N^*({{\textbf{S}}})$$\end{document} be the posterior density of $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} implied by the stationary distribution of the proposed method. Given the true parameters $S_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}_0$$\end{document} for the Ising model and any $\epsilon >0,$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon >0,$$\end{document} we have ${\pi }_N^{\ast }\left(S\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi _N^*({{\textbf{S}}})$$\end{document} concentrates at $S_0,$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}_0,$$\end{document}

(8) $\begin{array}{c}{\int }_{B_{\epsilon }\left(S_0\right)}{\pi }_N^{\ast }\left(S\right)dS\to 1,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \int _{B_\varepsilon ({{\textbf{S}}}_0)}\pi _N^*({{\textbf{S}}})d{{\textbf{S}}}\rightarrow 1, \end{aligned}$$\end{document}

in probability as $N\to \infty$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty $$\end{document} . $B_{\epsilon }\left(S_0\right)=\{S:\Vert S-S_0\Vert <\epsilon \}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_\varepsilon ({{\textbf{S}}}_0)=\{{{\textbf{S}}}:\Vert {{\textbf{S}}}-{{\textbf{S}}}_0\Vert <\varepsilon \}$$\end{document} is the open ball of radius $\epsilon$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} at $S_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}_0$$\end{document} .

We provide intuitions about this consistency result. Suppose that the data are generated by an Ising model. The iterative imputation method ensures that the parameters of the logistic regressions are close to those implied by the true Ising model, and thus, the conditional distributions we use to impute the missing values are close to those under the true model. This further guarantees that the joint distribution of the imputed data given the observed ones is close to that under the true Ising model, and consequently, the Ising model parameters we learn from the imputed data are close to those of the true model.

1.4. Computational Details

In what follows, we discuss the specification of the prior distributions and the sampling of auxiliary parameters ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_j$$\end{document} and Ising model parameters $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} .

Sampling ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_j$$\end{document} We set independent mean-zero normal priors for entries of ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_j$$\end{document} . For the intercept parameter ${\beta }_{\mathrm{jj}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{jj}$$\end{document} , we use a weakly informative prior by setting the variance to 100. For the slope parameters ${\beta }_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{jk}$$\end{document} , $k\ne j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k\ne j$$\end{document} , we set a more informative prior by setting the variance to be 1, given that these parameters correspond to the off-diagonal entries of $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} , which are sparse and typically do not take extreme values. The sampling of the auxiliary parameters ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\beta }}_j$$\end{document} , following (3), is essentially a standard Bayesian logistic regression problem. We achieve it by a Markov chain Monte Carlo (MCMC) sampler called the Pólya–Gamma sampler (Polson et al., 2013).

To obtain ${\beta }_j^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }^{(t)}_j$$\end{document} , this sampler starts with ${\beta }_j^{(t-1)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }^{(t-1)}_j$$\end{document} from the previous step. It constructs an MCMC transition kernel by a data argumentation trick. More precisely, the following two steps are performed.

1. 1. Given ${\beta }^{(t-1)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }^{(t-1)}$$\end{document} , independently sample N augmentation variables, each from a Pólya–Gamma distribution (Barndorff-Nielsen et al., 1982).

2. 2. Given the N augmentation variables, sample ${\beta }^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }^{(t)}$$\end{document} from a J-variate normal distribution.

The details of these two steps are given in the supplementary material, including the forms of the Pólya–Gamma distributions and the mean and covariance matrix of the J-variate normal distribution. We choose the Pólya–Gamma sampler because it is very easy to construct and computationally efficient. It is much easier to implement than Metropolis–Hastings samplers which often need tuning to achieve good performance.

We comment on the computational complexity of the sampling of ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_j$$\end{document} . The sampling from the Pólya–Gamma distribution has a complexity O(NJ), and the sampling from the J-variate normal distribution has a complexity of $O\left(N{J}^2\right)+O\left(J^3\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(NJ^2)+O(J^3)$$\end{document} . Consequently, a loop of all the ${\beta }_j,j=1,\dots ,J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\beta }}_j, j=1,\ldots , J$$\end{document} , has a complexity of $O\left((N+J)J^3\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O((N+J)J^3)$$\end{document} .

Sampling $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} Since $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} is a symmetric matrix, we reparametrize it by vectorizing its off-diagonal entries (including the diagonal entries). Specifically, the reparameterization is done by half-vectorization, denoted by $\alpha =\text{vech}\left(S\right)={(s_{11},\dots ,s_{J1},s_{22},\dots ,s_{J2},\dots ,s_{\mathrm{JJ}})}^{\top }\in R^{J(J+1)/2}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}= \text {vech}({{\textbf{S}}}) = (s_{11},\ldots , s_{J1}, s_{22},\ldots , s_{J2},\ldots , s_{JJ})^\top \in {\mathbb {R}}^{J(J+1)/2}$$\end{document} . It is easy to see that $\text{vech}(·)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {vech}(\cdot )$$\end{document} is a one-to-one mapping between $R^{J(J+1)/2}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^{J(J+1)/2}$$\end{document} and $J\times J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J\times J$$\end{document} symmetric matrices. Therefore, we impose a prior distribution on $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}$$\end{document} and sample ${\alpha }^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}^{(t)}$$\end{document} in the tth iteration when $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} is sampled. Then we let $S^{\left(t\right)}={\text{vech}}^{-1}\left({\alpha }^{\left(t\right)}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}^{(t)} = \text {vech}^{-1}({\varvec{\alpha }}^{(t)})$$\end{document} .

Recall that a thinning step is performed, and $t_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0$$\end{document} is the gap between two samples of $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} . Let t be a multiple of $t_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0$$\end{document} and ${\alpha }^{(t-t_0)}=\text{vech}\left(S^{(t-t_0)}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}^{(t-t_0)} = \text {vech}({{\textbf{S}}}^{(t-t_0)})$$\end{document} be previous value of $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}$$\end{document} . The sampling of ${\alpha }^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}^{(t)}$$\end{document} is also achieved by a Pólya–Gamma sampler, which involves the following two steps similar to the sampling of ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_j$$\end{document} .

1. 1. Given ${\alpha }^{(t-t_0)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}^{(t-t_0)}$$\end{document} , independently sample NJ augmentation variables, each from a Pólya–Gamma distribution.

2. 2. Given the NJ augmentation variables, sample ${\alpha }^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}^{(t)}$$\end{document} from a $J(J+1)/2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J(J+1)/2$$\end{document} -variate normal distribution.

The details of these two steps are given in the supplementary material. We note that the computational complexity of sampling the NJ augmentation variables is $O\left(N{J}^2\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(NJ^2)$$\end{document} , and that of sampling ${\alpha }^{\left(t\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}^{(t)}$$\end{document} is $O\left(N{J}^5\right)+O\left(J^6\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(NJ^5)+O(J^6)$$\end{document} , resulting in an overall complexity $O\left((N+J)J^5\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O((N+J)J^5)$$\end{document} . Comparing the complexities of the imputation and sampling $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} steps, we notice that the latter is computationally much more intensive. This is the reason why we choose to impute data by introducing auxiliary parameters ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta _j}$$\end{document} s rather than using Ising network parameters $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} so that the iterative imputation mixes much faster in terms of the computation time. In addition, we only sample $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} every $t_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0$$\end{document} iterations for a reasonably large $t_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0$$\end{document} to avoid a high computational cost and also reduce the auto-correlation between the imputed data.

We remark that Marsman et al. (2022) considered a similar Ising network analysis problem based on fully observed data, in which they proposed a Bayesian inference approach based on a spike-and-slab prior to learning $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} . Their Bayesian inference is also based on a Pólya–Gamma sampler. However, they combined Gibbs sampling with a Pólya–Gamma sampler, updating one parameter in $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} at a time. This Gibbs scheme often mixes slower than the joint update of $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} as in the proposed method and, thus, is computationally less efficient. The proposed Pólya–Gamma sampler may be integrated into the framework of Marsman et al. (2022) to improve their computational efficiency.

2. Numerical Experiments

We illustrate the proposed method and show its power via simulation studies and a real-world data application. In Sect. 2.1, we conduct two simulation studies, evaluating the proposed method under two MAR scenarios, one of which involves missingness due to screening items. A further simulation study is carried out, applying our method to a 15-node Ising model governed by the MCAR mechanism. Detailed exposition of this study can be found in the supplementary materials.

2.1. Simulation

Study I We generate data from an Ising model with $J=6$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=6$$\end{document} variables. Missing values are generated under an MAR setting that is not MCAR. The proposed method is then compared with Bayesian inference based on (1) listwise deletion and (2) a single imputation, where the single imputation is based on the imputed data from the Tth iteration of Algorithm 1, recalling that $T_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_0$$\end{document} is the burn-in size.

We configure the true parameter matrix $S_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}_0$$\end{document} as follows. Since $S_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}_0$$\end{document} is a symmetric matrix, we only need to specify its upper triangular matrix and then the diagonal entries. For the upper triangular entries (i.e., $s_{\mathrm{jl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{jl}$$\end{document} , $j<l$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j<l$$\end{document} ), we randomly assign 50% of them to zero to introduce a moderately sparse setting. In addition, the nonzero parameters are then generated by sampling from a uniform distribution over the set $[-1,-0.4]\cup [0.4,1]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[-1, -0.4] \cup [0.4, 1]$$\end{document} . The intercept parameters $s_{\mathrm{jj}},j=1,\dots ,J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{jj},j=1,\ldots , J$$\end{document} are set to zero. The true parameter values are given in the supplementary material. Missing data are simulated by masking particular elements under an MAR mechanism. In particular, we have $z_{i6}=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_{i6}=1$$\end{document} , so that the sixth variable is always observed. We further allow the missingness probabilities of the first five variables (i.e., $z_{\mathrm{ij}}=0,j=1,\dots ,5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_{ij}=0,j=1,\ldots ,5$$\end{document} ) to depend on the values of $y_{i6}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{i6}$$\end{document} . The specific settings on $p(z_{\mathrm{ij}}=0\mid y_{i6}),j=1,\dots ,5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(z_{ij}=0\mid y_{i6}),j=1,\ldots , 5$$\end{document} are detailed in the supplementary material. Data are generated following the aforementioned Ising model and MAR mechanism for four different sample sizes, $N=1000,2000,4000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = 1000, 2000, 4000$$\end{document} , and 8000, respectively. For each sample size, 50 independent replications are created.

Three methods are compared—the proposed method, Bayesian inference with a single imputation, and Bayesian inference based on complete cases from listwise deletion. The Bayesian inference for complete data is performed by sampling parameters from the posterior implied by the pseudo-likelihood and a normal prior, which is a special case of the proposed method without iterative imputation steps. All these methods shared the same initial values $s_{\mathrm{jl}}^{\left(0\right)}\sim U(-0.1,0.1),1\le j\le l\le J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{jl}^{(0)}\sim U(-0.1,0.1),1\le j\le l\le J$$\end{document} . For our proposed method, we set the length of the Markov chain Monte Carlo (MCMC) iterations to $T=5000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T = 5000$$\end{document} and a burn-in size of $T_0=1000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_0 = 1000$$\end{document} , with a thinning parameter $k_0=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_0=10$$\end{document} . This setup leads to an effective total of 400 MCMC samples for the Ising parameter matrix $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\textbf{S}}}$$\end{document} . Notably, identical MCMC length and burn-in configuration are applied during parameters inference in the single imputation and complete-case analyses.

Figure 2 gives the plots for the mean squared errors (MSE) of the estimated edge parameters and intercept parameters under different sample sizes and for different methods. The MSE for each parameter $s_{\mathrm{jl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{jl}$$\end{document} is defined as

(9) $\begin{array}{c}\frac{1}{50}\sum _{k=1}^{50}{({\hat{s}}_{k,jl}-s_{0,jl})}^2.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{1}{50} \sum _{k=1}^{50}(\hat{s}_{k,jl}-s_{0,jl})^2. \end{aligned}$$\end{document}

Here, ${\hat{s}}_{k,jl}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{s}_{k,jl}$$\end{document} denotes the estimated value from the kth replication while $s_{0,jl}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{0,jl}$$\end{document} refers to the true value. Each box in panel (a) corresponds to the 15 edge parameters, and each box in panel (b) corresponds to the 6 intercept parameters. We notice that the listwise deletion procedure introduces biases into the edge and intercept estimation, resulting in the MSEs for certain parameters not decaying toward zero as the sample size grows. Additionally, both the proposed method and the single imputation method offer accurate parameter estimation, with MSEs decaying toward zero as the sample size increases. Notably, the proposed method is substantially more accurate than the single imputation method, suggesting that aggregating over multiple imputed datasets improves the estimation accuracy. Furthermore, for smaller sample sizes, the complete-case analysis seems to yield slightly more accurate estimates of the edge parameters than the single imputation method. Across four sample sizes, the median computational times for obtaining the results of the proposed method were 33, 50, 88, and 185 s, respectively.1

Figure 2 a Boxplots of MSEs for edge parameters $s_{\mathrm{jl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{jl}$$\end{document} . b Boxplots of MSEs for intercept parameters $s_{\mathrm{jj}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{jj}$$\end{document} .

Study II: Missing due to screening items Missingness due to screening items is commonly encountered in practice, posing challenges to the network analysis (Borsboom et al., 2017; McBride et al., 2023). This occurs, for example, in surveys where initial screening questions determine respondents’ eligibility or relevance for subsequent questions. Suppose respondents indicate a lack of relevant experience (i.e., their answers to the screening items are all negative). In that case, they are not prompted to answer certain follow-up questions, making the missingness of these responses depend on their answers to the screening questions and, thus, MAR. Our real data example in Sect. 2.2 involves two screening items, which results in a large proportion of missing data.

Table 1 MSEs and biases for edge parameters.

We consider a simulation setting involving two screening items to evaluate the proposed method’s performance under this setting. Similar to Study I, we consider a setting with six items, the first two of which are the screening items. The full data are generated under an Ising model, whose parameters are given in the supplementary material, where the corresponding network has six positive edges including one between the two screening items. The responses to the screening items are always set as observed for any observation. When an observation’s responses to the screening items are both zero, their responses to the rest of the four items are regarded as missing.

We consider a single sample size $N=8000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=8000$$\end{document} and generate 50 independent datasets. We apply the proposed method, the single imputation method, and the complete-case analysis. For each estimation procedure, we set the MCMC iterations $T=5000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T = 5000$$\end{document} , the burn-in size $T_0=1000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_0 = 1000$$\end{document} , and the thinning parameter $k_0=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_0=10$$\end{document} . These methods are compared in terms of MSEs and biases for parameter estimation.

Table 1 presents the result. For all the edge parameters except for $s_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{12}$$\end{document} , the three estimation methods work well, though the single imputation method is slightly less accurate, as indicated by its slightly larger MSEs. However, the complete-case estimate is substantially negatively biased for $s_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{12}$$\end{document} , the edge between two screening items. At the same time, the imputation-based methods are still accurate, with the proposed method having a smaller MSE than that of the single imputation method. This result confirms that running a complete-case analysis on data involving screening items is problematic while performing the imputation-based methods, especially the proposed method, yields valid results.

We provide discussions on this result. The negative bias for $s_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{12}$$\end{document} in the complete-case analysis is due to a selection bias, typically referred to as Berkson’s paradox (De Ron et al., 2021). The complete-case analysis excludes all the response vectors with negative responses to both screening items. Consequently, a positive response on one screening item strongly suggests a negative response on the other, regardless of the responses to the rest of the items. This results in a falsely negative conditional association between the two screening items. In fact, one can theoretically show that the frequentist estimate of $s_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{12}$$\end{document} based on the maximum pseudo-likelihood is negative infinity. The finite parameter estimate in Table 1 for $s_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{12}$$\end{document} is due to the shrinkage effect of the prior distribution that we impose. On the other hand, the proposed method makes use of the observed frequency of the (0, 0) response pattern for the two screening items, in addition to the frequencies of the fully observed response vectors. As shown by the identifiability result in Appendix, these frequencies are sufficient for identifying all the parameters of the Ising model.

2.2. A Real Data Application

We analyze the dataset for the 2001–2002 National Epidemiological Survey of Alcohol and Related Conditions (NESARC), which offers valuable insights into alcohol consumption and associated issues in the US population (Grant et al., 2003). The dataset consists of 43,093 civilian non-institutionalized individuals aged 18 and older. In this analysis, we focus on two specific sections of the survey that concern two highly prevalent mental health disorders—major depressive disorder (MDD) and generalized anxiety disorder (GAD). Because MDD and GAD have high symptom overlap (Hettema, 2008) and often co-occur (Hasin et al., 2005), it is important to perform a joint analysis of the symptoms of the two mental health disorders and study their separation. In particular, Blanco et al. (2014) performed factor analysis based on the same data and found that the two mental health disorders have distinct latent structures. We reanalyze the data, hoping to gain some insights from the network perspective of the two mental health disorders.

Table 2 Descriptions of MDD and GAD items and their missing rates.

Following Blanco et al. (2014), we consider data with nine items measuring MDD and six items measuring GAD. These items are designed according to the Diagnostic and Statistical Manual of Mental Disorders, Fourth Edition (DSM-IV) (American Psychiatric Association, 2000). These items ask the participants if they have recently experienced certain symptoms; see Table 2 for their short descriptions. After eliminating samples with entirely absent values across the 15 items, a total of 42,230 cases remain in the dataset. Note that any “Unknown” responses in the original data are converted into missing values. The dataset exhibits a significant degree of missingness, with only 2,412 complete cases for the 15 items, representing approximately 6% of the total cases. Specifically, the missing rate for each item is given in Table 2. Importantly, items D1 and D2 function as screening items and, thus, have a very low missing rate. The respondents did not need to answer items D3–D9 if the responses to D1 and D2 were “No” or “Unknown,” resulting in high missing rates for these items. This pattern suggests that the missing data in this study is not MCAR. The GAD items A1–A6 also have a screening item, which results in the high missing rates in A1 through A6. Following the treatment in Blanco et al. (2014), these screening items are not included in the current analysis.

Figure 3 Estimated network structure for MDD and GAD. a Complete-case analysis. b Proposed method.

We apply the proposed method and the complete-case analysis to the data. For each method, 10 MCMC chains with random starting values are used, each having 10,000 MCMC iterations and a burn-in size 5000. The Gelman–Rubin statistics are always below 1.018, confirming the satisfactory convergence of all 120 parameters for both methods. The estimated network structures for MDD and GAD items are presented in Fig. 3, where an edge is shown between two variables when the absolute value of the estimated parameter is greater than 0.5. We emphasize that this threshold is applied only for visualization purposes, rather than for edge selection. Consequently, the edges in Fig. 3 should only be interpreted as edges with large estimated parameters, rather than truly nonzero edges. The nine MDD items are shown as blue nodes at the bottom, and the six GAD items are shown as orange nodes at the top. The edges are colored blue and orange, which represent positive and negative parameter estimates, respectively. In addition, the line thickness of the edges indicates their magnitude. A clear difference between the two methods is the edge between D1 “depressed mood most of the day, nearly every day,” and D2 “markedly diminished interest or pleasure in all, or almost all, activities most of the day, nearly every day,” which are two screening questions in the survey that all the participants responded to. The estimated parameter for this edge has a large absolute value under each of the two methods, but the estimated parameter is negative in the complete-case analysis, while it is positive according to the proposed method. As revealed by the result of Study II in Sect. 2.1, the negative edge estimate of the edge between the screening items given by the complete-case analysis is spurious. Considering the content of these items, we believe that the estimate from the proposed method is more sensible. Other than this edge, the remaining structure of the two networks tends to be similar, but with some differences. In particular, we see that the complete-case analysis yields more edges than the proposed method; for example, the edges of A4–A5, A1–D5, D1–D6, D1–D7, D1–D8, and D8–D9 appear in the estimated network from the complete-case analysis but not in that of the proposed method, while only two edges, A3–A5 and D3–D4, are present in the network estimated by the proposed method but absent in the network from the complete-case analysis. We believe this is due to the higher estimation variance of the complete-case analysis caused by its relatively small sample size.

Finally, our analysis shows that the symptoms of each mental health disorder tend to densely connect with each other in the Ising network, while the symptoms are only loosely but positively connected between the two mental health disorders. The edges between the two mental health disorders identify the overlapping symptoms, including “D4: Insomnia or hypersomnia” and “A6: Sleep disturbance,” “A2: Easily fatigued” and “D6: Fatigue/loss of energy,” and “A3: Difficulty concentrating” and “D8: Diminished concentration.” These results suggest that MDD and GAD are two well-separated mental health disorders, despite their high symptom overlap and frequent co-occurrence. This result confirms the conclusion of Blanco et al. (2014) that GAD and MDD are closely related but different nosological entities.

3. Concluding remarks

In this paper, we propose a new method for Ising network analysis in the presence of missing data. The proposed method integrates iterative imputation into a Bayesian inference procedure based on conditional likelihood. An asymptotic theory is established that guarantees the consistency of the proposed estimator. Furthermore, a Pólya–Gamma machinery is proposed for the sampling of Ising model parameters, which yields efficient computation. The power of the proposed method is further shown via simulations and a real data application. An R package has been developed that will be made publicly available upon the acceptance of the paper.

The current work has several limitations that require future theoretical and methodological developments. First, this manuscript concentrates on parameter estimation for the Ising model in the presence of missing data. However, the problem of edge selection (Borsboom, 2022; Marsman et al., 2022; Noghrehchi et al., 2021; Roçkovóa, 2018) requires future investigation. There are several possible directions. One direction is to view it as a multiple testing problem and develop procedures that control certain familywise error rates or the false discovery rate for the selection of edges. To do so, one needs to develop a way to quantify the uncertainty for the proposed estimator. It is non-trivial, as the proposed method is not a standard Bayesian procedure, and we still lack a theoretical understanding of the asymptotic distribution of the proposed procedure. In particular, it is unclear whether the Bernstein–von Mises theorem that connects Bayesian and frequentist estimation holds under the current setting. Another direction is to view it as a model selection problem. In this direction, we can use sparsity-inducing priors to better explore the Ising network structure when it is sparse. We believe that the proposed method, including the iterative imputation and the Pólya–Gamma machinery, can be adapted when we replace the normal prior with the spike-and-slab prior considered in Marsman et al. (2022). This adaptation can be done by adding some Gibbs sampling steps. In addition, it is of interest to develop an information criterion that is computationally efficient while statistically consistent. This may be achieved by computing an information criterion, such as the Bayesian information criterion, for each imputed dataset and then aggregating them across multiple imputations. Finally, the proposed method has several variants that may be useful for problems of different scales. For problems of a relatively small scale (i.e., when J is small), we may perform data imputation using the sampled $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf{S}}$$\end{document} instead of using auxiliary parameters ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\beta }}_j$$\end{document} s. This choice will make the algorithm computationally more intensive, as the sampling of $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf{S}}$$\end{document} has a high computational complexity. On the other hand, it may make the estimator statistically more efficient as it avoids estimating the auxiliary parameters ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\beta }}_j$$\end{document} s, whose dimension is higher than $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf{S}}$$\end{document} . For very large-scale problems, one may estimate the Ising model parameters based only on the auxiliary parameters ${\beta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\beta }}_j$$\end{document} s. For example, we may estimate $s_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{ij}$$\end{document} by averaging the value of $({\beta }_{\mathrm{ij}}+{\beta }_{\mathrm{ji}})/2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\beta _{ij}+\beta _{ji})/2$$\end{document} over the iterations. This estimator is computationally more efficient than the proposed one, as it avoids sampling $S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf{S}}$$\end{document} given the imputed datasets. This estimator should still be consistent but may be statistically slightly less efficient than the proposed one.