A note on large-scale logistic prediction: using an approximate graphical model to deal with collinearity and missing data

Large-scale prediction problems are often plagued by correlated predictor variables and missing observations. We consider prediction settings in which logistic regression models are used and propose a novel approach to make accurate predictions even when predictor variables are highly correlated and only partly observed. Our approach comprises three steps: first, to overcome the collinearity issue, we propose to model the joint distribution of the outcome variable and the predictor variables using the Ising network model. Second, to render the application of Ising networks feasible, we use a latent variable representation to apply a low-rank approximation to the network’s connectivity matrix. Finally, we propose an approximation to the latent variable distribution that is used in the representation to handle missing observations. We demonstrate our approach with numerical illustrations.


Introduction
Most large-scale or big data applications involve conditional models that utilize covariates to make predictions about a variable of interest.For instance, Google needs to predict which links to websites will be most advantageous based on millions of previous clicks, and Netflix needs to predict movie preferences based on millions of previous viewings and rankings.In these applications the interest is not in explaining why the connections between websites or movies exist, but in predicting which website will be most often requested or which movie will be preferred by an individual user.We will focus on the prediction problem where both the outcome variable and the covariates are binary, and the logistic regression model is an appropriate statistical model.
As is the case with all regression models, we observe that the logistic regression model is developed for situations where the covariates are independent and completely observed.However, a different situation is usually observed in largescale applications, where covariates are typically correlated.As a consequence of the correlations between covariates, i.e., collinearity, the obtained set of coefficients is no longer unique.This can be seen in, for instance, the coordinate descent algorithm (Hastie et al. 2015), where each covariate is treated separately.For two equivalent covariates, any solution with a linear combination of the two (normalized) coefficients is correct, even when regularization is applied.This is certainly an issue for the identification of relevant covariates, i.e., variable selection.In a particular sample, one of the collinear covariates will have a slightly larger coefficient and, therefore, ends up in the solution, while the other does not.But in another sample it could be the other way around.This means that variable selection with collinear covariates is unreliable.We will illustrate that collinearity is not a problem for prediction.
Another issue is that in most large-scale applications covariates are only partially observed (e.g., Rubin 1976;Rousseeuw 2016).For estimation and variable selection, it is then pertinent to know in which way the data came to be missing.For instance, data could be missing completely at random, which means that there is no connection between the missing observations and the data generating process.But data could also be missing precisely because of the data-generating process.For instance, a response to the question''do you drink more on average than others in your circle of friends'' will be missing if a negative response is observed to the question ''do you take alcohol''.In such cases, conditional on taking alcohol, the two covariates will be correlated, and so the missing observations cannot be ignored.We illustrate that predictions based on partially observed data can still be accurate, even when the process that generates missing data cannot be ignored in a statistical sense.
Our goal in this paper was to introduce a novel approach to make accurate predictions with the logistic regression model when covariates are highly correlated and only partly observed.Our approach comprises three steps: First, we propose to model the joint distribution of the outcome variable and the predictor variables with the Ising model.In the Ising model the correlations between the observed variables are explicitly modeled, which overcomes the collinearity issue.Second, we use recent results that relate Ising networks to latent variable models to render the application of Ising models computationally tractable.Specifically, we use a low-rank approximation to the network's connectivity matrix, which is opportune when variables are highly correlated (i.e., collinearity).Finally, we propose to approximate the latent variable distribution in the representation of the Ising model, which results in a model-based approximation to the full Ising model that is able to handle missing observations.Numerical illustrations are used to demonstrate different features of our approach.

Step I: the Ising model to overcome collinearity
For prediction, we are interested in the conditional distribution Pðx i j x ni Þ, where x i is an element of x 2 fÀ1; þ 1g p and x ni is the vector x excluding the ith element.Even though we are only interested in the predictive distribution Pðx i j x ni Þ, our observations are the realizations of a multivariate random variable x.The multivariate distribution PðxÞ that is consistent with the logistic regression model is the Ising network model (Lenz 1920;Ising 1925), where R 2 R pÂp is a symmetric matrix of pairwise interaction parameters r ij , and l 2 R p a vector of main effects.Observe that the correlations between elements in x are no problem for the Ising model as their interactions r ij are explicitly modeled.That is, there is no collinearity issue when estimating the full Ising model PðxÞ.
From the joint distribution we can then obtain the correct full-conditional Pðx i j x ni Þ, and it is easily seen that the full-conditional distribution that is obtained from the Ising model is a logistic regression model: Hence, we can overcome the collinearity issues with logistic regression by estimating the full Ising model.However, estimating an Ising model proves to be far more complex than estimating a logistic regression model.A first problem is the number of parameters that needs to be estimated for the Ising model.Whereas the number of parameters is linear in the number of covariates p for the logistic regression model, it is quadratic in p for the Ising model.A second problem is that the density of the Ising model is computationally intractable, except for small or heavily constrained networks.This computational burden is due entirely to the model's normalizing constant: which is the sum over all 2 p possible realizations of x.Thus, even though we have resolved the collinearity issue with the Ising model, we have also increased the number of parameters with an order of magnitude and need to deal with estimating a model that is computationally intractable. 3 Step II: low-rank approximations for computational tractability A latent variable representation of the Ising model, in combination with a low-rank approximation to the full-connectivity matrix R, is the two crucial ingredients to render large-scale applications of the Ising model entirely tractable.The latent variable representation of the Ising model was introduced by Kac (1968).Specifically, Kac showed that every eigenvector of the connectivity matrix R generates a latent variable, such that the manifest random variables x are independent given the full set of latent variables g.Since the diagonal elements from the connectivity matrix R are not identifiable from the data, we decompose it as where K is a diagonal matrix consisting of the eigenvalues of the original connectivity matrix, and the translation by c serves to ensure that all eigenvalues are positive, i.e., ensuring that UU T is positive (semi-)definite and at the same time preserve the off-diagonal elements of R. The latent variable representation of Kac then follows immediately from a clever use of the Gaussian identity: This latent variable representation has been further developed by Emch and Knops (1970) and has been independently rediscovered many times (e.g., Olkin and Tate 1961;Besag 1974;McCullagh 1994;Anderson and Vermunt 2000).
It was recently shown (Marsman et al. 2015;Epskamp et al. 2017) that the associated conditional distribution Pðx j gÞ is the multidimensional item response theory (MIRT) model (Reckase 2009) with u i being the ith column of U T .MIRT models are frequently used in psychological and educational measurement (Ackerman et al. 2003;Ackerman 1996), where the observed variables correspond to item responses on some test or questionnaire, and the latent variables relate to the trait or abilities being assessed (Borsboom and Molenaar 2015).Importantly, this representation inspired a fulldata-information estimation procedure that avoids having to compute the Ising model's intractable normalizing constant (Marsman et al. 2015).
A second crucial ingredient is the low-rank approach that Marsman et al. (2015) proposed to approximate the full connectivity matrix, such that the number of parameters becomes linear in p. Their low-rank approach makes use of the Eckart and Young Theorem (1936), which states that in a least squares sense the best rank-r approximation to the full connectivity matrix R is one in which all but the r largest eigenvalues are equated to zero.Low-rank approximations have become increasingly popular in prediction problems since their crucial role in winning the Netflix price competition (Koren et al. 2009;Bell and Koren 2007;Bell et al. 2010) and it has been part of Google's system ever since the very first implementation of the pageRank algorithm (Page et al. 1999;Brin and Page 2012).Most important for our present endeavors, however, is that a low-rank approximation to the full connectivity matrix is expected to make accurate predictions when predictor variables are highly correlated. 4 Step III: an approximate latent variable distribution for missing data An important feature of IRT models is that they are closed under marginalization.That is, because the manifest variables x i are independent given the latent variable g, we find that the marginal, is again an IRT model.That the IRT model is closed under marginalization makes it a valuable tool for applications with data that are subject to missing data (Eggen 2004), as one can simply marginalize over the missing observations.In contrast, the Ising model is not closed under marginalization.That is, in general we find that is itself not an Ising model.Unfortunately, the marginalization property of the IRT model does not transfer to the latent variable expression of the Ising model, which is entirely due to the latent variable distribution f ðgÞ.To train the Ising model in the face of incomplete data, we, therefore, either have to omit any incomplete cases from the analysis, use imputation techniques to artificially complete the observed data (Rubin 1987), or make use of approximate models that allow for missing data.We take the latter approach and show that we can make reliable predictions using a low-rank approximate model, even in the presence of correlated and missing data.
A key difference between regular applications of IRT models and the latent variable representation of the Ising model is in the prior (or population) distribution of the latent variables that are used.The distribution of latent variables is typically assumed to be multivariate normal, but in the representation of the Ising model it is the mixture: where to each of 2 p possible realizations of x we have a multivariate normal posterior distribution with mean U T x and variance 2I p .Even though this latent variable model is computationally intractable, it often takes a simple form in each of its dimensions.Specifically, it closely resembles either a single normal distribution with a mean of zero or a mixture of two normal distributions with their respective means placed symmetrically about zero.That the latent variable distribution tends to have either a single mode or has two modes can be seen by inspecting the derivative of logðf ðgÞÞ with respect to g (for ease of presentation assuming a single dimension); which shows that the latent variable distribution has modes (and minima) at the fixed points hðgÞ ¼ g.A plot of hðgÞ against g is shown in Fig. 1, together with the line g ¼ g, which rotates with respect to hðgÞ as a function of the variables in the model.Note that the line can cross the curve hðgÞ either once (at zero) such that there is a single mode, or three times (as depicted here) such that there are two modes (symmetric about zero) and a local minimum at zero.The latent variable distribution f ðgÞ can thus be closely approximated by a small mixture of normal distributions.
There are two important consequences of replacing the latent variable distribution of the Ising model with some other latent variable distribution, say gðgÞ.A first consequence is that the marginal distribution of the observed variables is, in general, not analytically available and requires numeric procedures to compute.A second consequence follows from the indeterminancy of the parameters U, l (and g) in the MIRT model.By replacing the latent variable distribution f ðgÞ with a distribution gðgÞ, the parameters are placed on a different scale.This Fig. 1 A plot of g versus hðgÞ illustrating that f ðgÞ can have either one or two modes.A local minimum and two local maxima are found on the intersection with the straight solid line indeterminancy does not affect the marginal distribution of observable variables (i.e., predictions), although it does affect parameter recovery.
The validity of our approximate model rests on how well f ðgÞ is approximated by gðgÞ.To assess the validity of this approach, we can make use of recent advances in plausible value methodology and explicitly consider whether or not the true latent variable distribution is equal (or similar) to gðgÞ.Plausible values are draws from the posterior distribution of the latent variables g (Mislevy 1991;von Davier et al. 2009) and are commonly used in large-scale educational surveys to accommodate researchers in the field that are not able to estimate the complex IRT models used for these surveys.Recently, it was shown that the marginal distribution of plausible values is a consistent estimator of the true latent variable distribution f ðgÞ (Marsman et al. 2016), meaning that one can assess the validity of using a single multivariate normal distribution by inspecting the (marginal) distribution of plausible values.

Numerical illustrations
Below we demonstrate the different aspects of our theory using three broad illustrations.The first illustration aims to showcase that the IRT model is able to make accurate predictions with correlated variables.The second illustration aims to showcase that the IRT model is able to accurately predict both observed and missing observations when the missing data mechanism is ignorable and the data are completely missing at random (MCAR; see Appendix A).We end this illustration with a comparison with logistic regression.The third illustration is used to demonstrate that there is a limit to the IRT model's capacity to accurately predict observed and missing data points when the missing data mechanism is nonignorable (data not missing at random; NMAR).We consider several situations in which we vary the effect of the missing data mechanism on the observed correlations and mix settings with MCAR and MNAR in the training and testing phase.The main results from our illustrations are shown in Fig. 2. We first discuss the methods and models that are used.

Generating correlated binary data
To generate correlated binary variables we use the Ising model, with a rank ten connectivity matrix R that is based on the following eigenvalues: :00; 0:80; 0:65; 0:30; 0:25; 0:20; 0:16; 0:11; 0:06; 0:01: The value of s modifies the strengths of pairwise correlations of variables in the network.Figure 3 illustrates that the more extreme correlations (i.e., AE1) occur for smaller values of s.In the three illustrations we will use the values s ¼ 0:5 and s ¼ 1:0.The Ising model's parameters Q and l will be sampled uniformly between À0:1=p and 0.1 / p for the p-variable networks (Q is made orthogonal), where scaling by p À1 ensured similar dynamics for the different sized networks.
Data were generated from the Ising model using a Gibbs sampler (Geman and Geman 1984) applied to the joint distribution f ðx;gÞ of the latent variables g and the data x.In each iteration of the Gibbs sampler, we sample from the full-conditional posterior distribution f ðg j xÞ of the latent variables, and the full-conditional distribution Pðx j gÞ of the data.Both full-conditional distributions are easy to sample from; the posterior distribution of the latent variable f ðg j xÞ is a multivariate Fig. 2 The main results for applications of a two-dimensional IRT model to data generated from a s ¼ 1:0 network.The y-axis shows the drop in accuracy for the test set predictions as compared to accuracy of the true model.That is, c true À c test if the complete data are used, or c true À c ðoÞ test when there are missing data.Predictions from the true model always make use of the complete data.The x-axis shows the different situations; collinearity (c.f.Table 1), ignorable missingness (denoted MCAR; c.f. Table 2), nonignorable missingness with moderate effects (denoted MNAR1; c.f. Table 5) and nonignorable missingness with severe effects (denoted MNAR2; c.f. Table 6) normal distribution with mean vector U T x (where as before) and a variance-covariance matrix 2I 10 , and Pðx j gÞ is a ten-dimensional IRT model.

Estimating the MIRT model
In each simulation we train the two-dimensional IRT model We assume a simple bivariate normal distribution for the latent variables; g $ N ð0;2I 2 Þ.We take a Bayesian approach to estimate the IRT model, which requires us to formulate prior distributions for U and l.We will use logistic prior distributions with location 0 and scale 1 for both U and l.A Gibbs sampler is used to produce samples from the joint multivariate posterior distribution of the parameters and latent variables, i.e., gðU; l; g j xÞ.Whereas the full-conditional distribution of the latent variables in the Ising model representation is a completely tractable normal posterior distribution, i.e., g j x; U; l $ N ðU T x; 2I 2 Þ, this is not the case when we replace the associated latent variable distribution.Specifically, when we utilize a normal prior distribution for the latent variables, we observe that the full-conditional posterior distribution is the intractable: Â gðgÞ: Similarly, we find that the full-conditional distributions of the ''item'' parameters U and l are also intractable: where gðu ij Þ and gðl i Þ are the logistic prior distributions, and v indexes the n observations.The problem of sampling from these full-conditional distributions has been addressed in several places (e.g., Patz and Junker 1999a, b;Maris and Maris 2002).We use a Metropolis approach (Metropolis et al. 1953;Hastings 1970;Tierney 1994) that was specifically designed to handle full-conditional distributions of this form (Marsman et al. 2015(Marsman et al. , 2017)).

Calculating prediction accuracy
The prediction accuracy can be calculated by 0-1 loss or Bayes risk.Computationally this has the advantage of being easy to compute, but it is also tied to convex alternatives like logistic loss that are asymptotically equivalent to 0-1 loss (see, e.g., Bartlett et al. 2006).In our simulations we use 0-1 prediction accuracy, which is defined as where 1 is the indicator function and we predict x vi using x Ã vi .Observe that cðx i ; x Ã i Þ is the ratio of correct predictions (true positives and true negatives) out of the n predictions that are made.
Since each of the p variables could be used as a dependent variable in Logistic regression-there are p full-conditionals Pðx i j x ni Þ-we calculate the prediction accuracy for each of the p variables and then average them.We furthermore repeat each procedure five times and average the results.

The two stages of our prediction procedure
Our prediction procedure consists of two stages, a training stage and a testing stage.This ends the training stage.Observe that the missing data were not used to estimate the parameters h train and g train .This ends the testing stage.Observe that testing data were not used to estimate the IRT parameters h train and that only the latent variables g test were estimated on the observed part from the testing data.

Illustration I: Collinearity
Table 1 reveals the prediction accuracy of a two-dimensional IRT model applied to data generated from a s ¼ 1:0 network, and data generated from a s ¼ 0:5 network.Evidently, the two-dimensional IRT model provides accurate predictions that crossvalidate well.As expected, the prediction accuracy is an increasing function of the observed sample correlations (c.f.Fig. 3).Furthermore, Table 1 shows that the approach is largely insensitive to variations in n and p, ensuring that the procedure scales when more observations n and/or more variables p become available.This scalability is important for situations where the number of observations n becomes too large, and one has to use a selection of the available observations, with the estimated IRT model cross-validating well in such applications.
The prediction accuracy of the IRT model is similar to the prediction accuracy of the true model.This indicates a good fit of the IRT model to the network data, which is somewhat surprising since we expect a bimodal or mixture distribution for the latent variables in a network of such correlated variables x i .This is indeed the case.Figure 4 shows the marginal distribution of the latent variables in the first dimension, i.e., the marginal distributions of plausible values.It is clear that the marginal distribution of plausible values is bimodal and diverges from the unimodal population distribution that we have used.The fact that our predictions are still on par irrespective of the fit of the latent variable distribution f ðgÞ shows that the prediction procedure is robust against misspecification of the latent variable model.

Illustration II: Ignorable missing data
Table 2 reveals the prediction accuracy of a two-dimensional IRT model applied to data generated from a s ¼ 1:0 network, with either 50 or 90% of the data missing completely at random (MCAR; see Appendix A).Similarly, in Table 3 we report the prediction accuracy of a two-dimensional IRT model applied to data generated from a s ¼ 0:5 network, with either 50 or 90% of the data MCAR.We report both the accuracy in predicting the observed data c The results that are reported in Tables 2 and 3 indicate that the IRT model efficiently operates when large portions of the data are MCAR.This is particularly evident when we compare the results in Tables 2 and 3 that are based on incomplete data with the results in Table 1 that are based on complete data.The prediction accuracy of the IRT model also compares favorably to the accuracy that is obtained from the true model that is also based on the complete data.Note also that the IRT model is entirely capable of making accurate predictions about the missing data, even with only 10% observations left to train the model.

Logistic regression
It is instructive to compare the MIRT model's performance with that of logistic regression.The ideal situation for this comparison would have more observations than variables (n [ p) so that there is no need for regularization in estimating the logistic regression model.We, therefore, use n ¼ 1000 observations on p ¼ 100 variables.
The results that are reported in Tables 2 and 3 are based on data with missing observations.To allow a meaningful comparison between logistic regression and MIRT when some observations are missing, we use multiple imputation to complete the datasets.Unfortunately, we now encounter a serious complication.To impute  the missing observations, we need to formulate an a priori distribution for the missing observations (c.f.Ibrahim et al. 2005).The most straightforward solution is to specify a joint distribution for a p þ 1 dimensional vector of binary variables x.That is, we assume that the variables are dependent on each other and inform about each other's missing values.Even though this is a straightforward strategy, it will boil down to an a priori distribution for the missing observations that is at least as complex as the computationally intractable Ising model.
To overcome this complication, we assume that the variables are a priori independent.Specifically, we assume for each missing observation that the prior probability that its value is þ1 is equal to some number p.We use the value p ¼ 0:5 since we have no a priori preference for a particular value of the missing observation.Denote the logistic regression model as with y 2 fÀ1; þ 1g the dependent variable and x 2 fÀ1; þ 1g pÀ1 a vector of covariates.The posterior distribution for a missing observation x i is easily computed: Observe that the distribution Pðx i j y; x ni Þ favors values of x i that minimize jð2y À 1Þ À EðyÞj ¼ jð2y À 1Þ À Pðy j xÞj.
We use the Gibbs sampler to estimate the logistic regression model and use logistic prior distributions with location 0 and scale 1 for the model's parameters a and b.Our imputation strategy expands the Gibbs sequence by two distinct steps.In the first step, missing values for the dependent variable are drawn from the predictive distribution Pðy j xÞ (i.e., the logistic regression model).In the second step, missing values for each of the p covariates are drawn from their respective posterior distributions Pðx i j y; x ni Þ.After these two steps the data are complete and we can simulate the model's parameters a and b from their full-conditional posterior distributions as if all data had been observed.
We generate 300 datasets from the s ¼ 1:0 network and 300 datasets from the s ¼ 0:5 network and randomly remove 0, 50 or 90 of the observations.The prediction accuracy for logistic regression applied to these datasets are reported in Table 4 and reveals three important results.The first result is that the MIRT model performs better on the completely observed data than the logistic regression model.This is likely due to collinearity, as the relative performance of logistic regression deteriorates with increasing correlations.For instance, when compared to the true model's prediction accuracy we observe an 8% accuracy drop for the s ¼ 1:0 network data and a 15% accuracy drop for the s ¼ 0:5 network data.
The second important result is that the logistic regression model's prediction accuracy on observed data is much improved when missing observations are introduced.In fact, logistic regression outperforms both the MIRT model and the true generating model on the remaining-observed-data. (This improvement was only seen in the dependent variable, not the covariates.)This striking increase in accuracy is due to the way that we impute the missing values.The imputation distribution Pðx i j y; x ni Þ tends to favor values that make the observed outcomes y more likely and minimize jð2y À 1Þ À Pðy j xÞj.As a result, prediction accuracy increases when more observations are missing and the model over-fits the remaining observed data.
The final important result that we observe from Table 4 is the poor prediction accuracy on the missing observations.Compared to the MIRT model the accuracy of predicting missing values drops approximately 25-35%.This striking difference between the accuracy on the missing data and on the observed data is a clear illustration of the poor cross-validation that follows from over-fitting on the observed data points.This is particularly problematic when one aims to predict nonobserved data points, e.g., classification of future preferences: whereas one believes to be doing quite a good job based on predictions of the observed data, one unknowingly is doing a very poor job in predicting non-observed data.

Illustration III: Nonignorable missing data
From the results that are reported in Tables 2 and 3 we have learned that the twodimensional IRT model provides accurate predictions when applied to correlated data where some of the observations are MCAR.Since there is no additional difficulty for the case when the data are MAR instead of MCAR, we consider here the situation where the IRT model is applied to data where the missing data mechanism is nonignorable, i.e., not missing at random (NMAR; see Appendix A).We compare situations where either the training data, the testing data, or both the training data and the testing data have 50% data NMAR or 50% data MCAR.
We use several mechanisms to produce nonignorable missing data patterns, with the missing data mechanism explicitly depending on the missing observations and/ or the parameters of the observed data model.In Appendix B we describe two procedures, one of which produces missing data patterns that have a moderate effect on observed correlations and parameter estimates-moderate nonignorability-and one which produces missing data patterns that have a severe effect on observed correlations and parameter estimates-severe nonignorability.The effect of the two described missing data mechanisms on the observed correlations is shown in Fig. 5a for moderate nonignorability, and Fig. 5c for severe nonignorability.In Fig. 5b, d we plot the corresponding estimates of u 1 , the first column of U, for the two situations described by Fig. 5a, c, respectively, against the estimates that are obtained from data where the missing observations are MCAR.Clearly, ignoring the missing data mechanism produces bias to the estimates of u 1 , especially for the severe nonignorability case.

estimates (MAR) estimates (nonignorable)
Fig. 5 The result of nonignorable missing data on observed correlations and parameter estimates.The two panels on the left show the distribution of the observed correlations from data generated from a s ¼ 1:0 network (black line), together with the observed correlations based on data subject to moderate nonignorability (top-left panel; gray line) and data subject to severe nonignorability (bottom-left panel; gray line).The two panels on the right show the estimates of u 1 , the first column of U, based on data that are MCAR against the estimates based on data subject to moderate nonignorability (top-right panel) and data subject to severe nonignorability (bottom-right panel) Table 5 reveals the prediction accuracy of the two-dimensional IRT model applied to data generated from a s ¼ 1:0 network, with three distinct situations in which either the training data, the test data, or both data sets have missing data patterns that have a moderate effect on observed correlations,i.e., moderate nonignorability.When we compare our predictions with the predictions that are made with the true model Pðx i j x ni Þ based on completely observed data, we see that the prediction accuracy drops approximately: 3-4% when some of the training data are NMAR and some of the testing data are NMAR, 4% when some training data are NMAR but testing data only have data MCAR, and 1% when the training data has data MCAR but the testing data have data NMAR.The 1% drop in accuracy that is found when the training data have some observations MCAR is within the range found for the complete data in Table 1 and data with observations MCAR in Table 2.This suggests that training the prediction model on data for which the missing data are missing at random does not threaten prediction accuracy.The 3-4% drop in accuracy that is found when training the model on data with nonignorable missing data patterns is slightly higher than the range found for the complete data and data with observations MCAR.Predictions about the missing testing data are, on average, more accurate than that of the observed testing data, whereas the opposite is found for the training data.
Table 6 reveals the prediction accuracy in the same situations as in Table 5, but where the nonignorable missing data mechanism has a severe effect on the observed correlations.When compared to predictions that are made by the true model Pðx i j x ni Þ based on complete data, we observe that the prediction accuracy drops approximately: 7-8% when some of the training data are NMAR and some of the testing data are NMAR, 8-9% when some training data are NMAR but testing data only have data MCAR, and 1% when the training data have data MCAR but the testing data have data NMAR.Thus, in comparison to the results reported in Table 5 for the missing data having a moderate effect on observed correlations, we have that the drop in prediction accuracy is roughly doubled when training the model on data

Discussion
We have illustrated that the combination of the Ising model and its latent variable approximation can be used to overcome collinearity and missing data issues in prediction applications of the logistic regression model.The prediction accuracy of the latent variable model on the observed data compares favorably to the true model used in Illustrations I-III (e.g., Fig. 2).The latent variable model was also able to accurately predict the non-observed data points in Illustration II (e.g., Tables 2, 3) and compares favorably to our illustration of the logistic regression model (e.g., Table 4).The model has its limits, which was clearly demonstrated in Illustration III using nonignorable missing data mechanisms.The prediction accuracy deteriorates when missing data mechanisms have a strong effect on the data that is used to train the model (e.g., Tables 5, 6).However, even with the poor quality of the data that were used to train the model in Illustration III, the prediction accuracy of the latent variable model compares favorably to our application of the logistic regression model (e.g., compare Tables 4 and 6).Therefore, we believe that our approach and the associated latent variable model are superior to the logistic regression model in prediction settings with correlated covariates and/or missing observations.We have considered a specific prediction setting with only binary random variables for our approach to overcome collinearity and missing data.Observe, however, that the two primary ideas that form our approach are entirely general.The first idea is to model the joint distribution of dependent and independent variables when the variables are correlated, e.g., Ising networks models, Gaussian graphical models (Lauritzen 1996), or mixtures thereof (Olkin and Tate 1961;Lauritzen and  Wermuth 1989).A specific feature of these particular models is that they specifically model the between variables and thus overcome collinearity issues in conditional regression models.However, the application of these models will also increase the number of parameters that need to be estimated.Our second idea is to approximate the full graphical model with a low-rank latent variable model, e.g., an IRT model, a factor model or mixtures thereof.This has two important benefits: it reduces the number of parameters that need to be estimated and introduces an elegant way of handling missing data.We believe that these two general ideas will inspire new avenues of future research and furthermore offer practical solutions to issues that are widespread in large-scale applications.
The latent variables g in this paper are used to summarize the observed data in order to make predictions.However, the latent variable, in combination with the model's parameters h ¼ fU;lg, can also be used to inform about the structure of the prediction problem.For instance, in regular applications of IRT models to educational tests (e.g., an end of primary school test), the model informs about the dimensionality of the test (e.g., separates a mathematics and language dimension) and informs how different aspects of the test covary across the ability spectrum.In much the same way we may study the latent variable model in prediction settings, which might improve the model and its predictions.For example, in psychological data we nearly always observe a pattern of positive correlations (in contrast to Fig. 3), which means that the entries in the first eigenvector tend to have the same sign.Without much loss of accuracy we can then replace the p unknown values in u 1 with a single (positive) number.This would significantly reduce the number of parameters that we need to estimate, and we obtain sparse models that can make accurate predictions.

Severe nonignorability
Severely nonignorable missing data were created as follows: First, we cycled through the consecutive column-pairs of the generated matrix (i.e., (1, 2), (2, 3),...,ðp À 1; pÞ, (p, When the correlation between columns i and j was positive we either removed all (1, 1) or ðÀ1; À1Þ observations, and when the correlation was negative we removed either all ð1; À1Þ or ðÀ1; 1Þ observations.We then cycled through the n rows of the generated matrix.For each row, we randomly generated a value y 2 f0;1;2;3g.We removed 50% of the þ1 observations when y equaled 0, 50% of the À1 observations when y equaled 1, removed 25% of the þ1 and 25% of the À1 responses if y equaled 2, and did nothing when y equaled 3.This created approximately 50% missing observations, and observations were either randomly omitted or placed back to make it an exact 50%.

Fig. 3
Fig. 3 Distribution of 499,500 correlations-Pearson's-/-in a p ¼ 1000 variable network based on n ¼ 1000 observations for s 2 f0:5;1:0;1:5;2:0g.With the value s ¼ 2:0 most of the observed correlations are close to zero, and with the value s ¼ 0:5 most of the observed correlations are near the extremes 5.4.1 The training stage comprises six steps (1) Generate training data x train from the Ising model.(2) Generate new predictions x Ã from the Ising model and compute c true ¼ cðx train ;x Ã Þ. (3) Split the data into an observed part x ðOÞ train and a 0%, 50% or 90% missing part x ðMÞ train .The missing part x ðMÞ train will only be used to evaluate predictions.(4) Use the Gibbs sampler to estimate the IRT parameters h train ¼ fU;lg and the latent variables g train using the observed training data x ðOÞ train .(5) Generate new predictions from the IRT model on the observed part,x ðOÞ Ã $ Pðx j g train ;h train Þ; Generate new predictions from the IRT model on the missing part,x ðMÞ Ã $ Pðx j g train ;h train Þ; Since the true model cannot be used with missing observations, we evaluate the predictions from the true model using the completely observed test data: c true ¼ cðx train ;x Ã Þ.

Fig. 4
Fig.4The marginal distribution of the n ¼ 1000 plausible values based on p ¼ 100 variables generated from a s ¼ 0:5 network

Table 5
Prediction accuracy of the two-dimensional IRT model applied to s ¼ 1:0 network data with moderate nonignorable missingness When the training data have data NMAR, we see that predicting the missing testing data is about as difficult as predicting the observed testing data.However, predicting the missing training data is clearly more difficult than predicting the observed training data.

Table 6
Prediction accuracy of the two-dimensional IRT model applied to s ¼ 1:0 network data with severe nonignorable missingness