<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/predictability-in-network-models/"></div></div><div><a data-url="https://www.r-bloggers.com/predictability-in-network-models/" data-counturl="https://www.r-bloggers.com/predictability-in-network-models/" data-count="horizontal" data-text="Predictability in Network Models" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://jmbh.github.io//Predictability-in-network-models/"> Jonas Haslbeck - r</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Network models have become a popular way to abstract complex systems and gain insights into relational patterns among observed variables in <a target="_blank" rel="nofollow" href="http://www.sachaepskamp.com/files/NA/NetworkTakeover.pdf">almost any area of science</a>. The majority of these applications focuses on analyzing the structure of the network. However, if the network is not directly observed (Alice and Bob are friends) but <em>estimated</em> from data (there is a relation between smoking and cancer), we can analyze – in addition to the network structure – the predictability of the nodes in the network. That is, we would like to know: how well can an arbitrarily picked node in the network predicted by all remaining nodes in the network?</p>
<p>Predictability is interesting for several reasons:</p>
<ol>
<li>It gives us an idea of how <em>practically relevant</em> edges are: if node A is connected to many other nodes but these only explain, let’s say, only 1% of its variance, how interesting are the edges connected to A?</li>
<li>We get an indication of how to design an <em>intervention</em> in order to achieve a change in a certain set of nodes and we can estimate how efficient the intervention will be</li>
<li>It tells us to which extent different parts of the network are <em>self-determined or determined by other factors</em> that are not included in the network</li>
</ol>
<p>In this blogpost, we use the R-package <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/mgm/index.html">mgm</a> to estimate a network model and compute node wise predictability measures for a <a target="_blank" rel="nofollow" href="http://cpx.sagepub.com/content/3/6/836.short">dataset</a> on <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Posttraumatic_stress_disorder">Post Traumatic Stress Disorder (PTSD)</a> symptoms of <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/2008_Sichuan_earthquake">Chinese earthquake victims</a>. We visualize the network model and predictability using <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/qgraph/index.html">the qgraph package</a> and discuss how the combination of network model and node wise predictability can be used to design effective interventions on the symptom network.</p>
<h2>Load Data</h2>
<p>We load the data which the authors made freely available:</p>
<figure>
<pre><code data-lang="r"><span>data</span><span> </span><span>&lt;-</span><span> </span><span>read.csv</span><span>(</span><span>'http://psychosystems.org/wp-content/uploads/2014/10/Wenchuan.csv'</span><span>)</span><span>
</span><span>data</span><span> </span><span>&lt;-</span><span> </span><span>na.omit</span><span>(</span><span>data</span><span>)</span><span>
</span><span>p</span><span> </span><span>&lt;-</span><span> </span><span>ncol</span><span>(</span><span>data</span><span>)</span><span>
</span><span>dim</span><span>(</span><span>data</span><span>)</span><span>
</span><span>[</span><span>1</span><span>]</span><span> </span><span>344</span><span> </span><span>17</span></code></pre>
</figure>
<p>The datasets contains complete responses to 17 PTSD symptoms of 344 individuals. The answer categories for the intensity of symptoms ranges from 1 ‘not at all’ to 5 ‘extremely’. The exact wording of all symptoms is in the <a target="_blank" rel="nofollow" href="http://cpx.sagepub.com/content/3/6/836.short">paper of McNally and colleagues</a>.</p>
<h2>Estimate Network Model</h2>
<p>We estimate a <a target="_blank" rel="nofollow" href="http://www.jmlr.org/proceedings/papers/v33/yang14a.pdf">Mixed Graphical Model (MGM)</a>, where we treat all variables as continuous-Gaussian variables. Hence we set the type of all variables to <code>type = 'g'</code> and the number of categories for each variable to 1, which is the default for continuous variables <code>lev = 1</code>:</p>
<figure>
<pre><code data-lang="r"><span>install.packages</span><span>(</span><span>'mgm'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>mgm</span><span>)</span><span>
</span><span>fit_obj</span><span> </span><span>&lt;-</span><span> </span><span>mgmfit</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>data</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>rep</span><span>(</span><span>'g'</span><span>,</span><span> </span><span>p</span><span>),</span><span>
</span><span>lev</span><span> </span><span>=</span><span> </span><span>rep</span><span>(</span><span>1</span><span>,</span><span> </span><span>p</span><span>),</span><span>
</span><span>rule.reg</span><span> </span><span>=</span><span> </span><span>'OR'</span><span>)</span></code></pre>
</figure>
<p>For more info on how to estimate Mixed Graphical Models using the mgm package see <a target="_blank" rel="nofollow" href="http://jmbh.github.io/Estimation-of-mixed-graphical-models/">this previous post</a> or the <a target="_blank" rel="nofollow" href="https://arxiv.org/pdf/1510.06871v2.pdf">mgm paper</a>.</p>
<h2>Compute Predictability of Nodes</h2>
<p>After estimating the network model we are ready to compute the predictability for each node. Node wise predictability (or error) can be easily computed, because the graph is estimated by taking each node in turn and regressing all other nodes on it. As a measure for predictability we pick the propotion of explained variance, as it is straight forward to interpret: 0 means the node at hand is not explained at all by other nodes in the nentwork, 1 means perfect prediction. We centered all variables before estimation in order to remove any influence of the intercepts. For a detailed description of how to compute predictions and to choose predictability measures, <a target="_blank" rel="nofollow" href="https://arxiv.org/abs/1610.09108">check out this preprint</a>. In case there are additional variable types (e.g. categorical) in the network, we can choose an appropriate measure for these variables (e.g. % correct classification, see <code>?predict.mgm</code>).</p>
<figure>
<pre><code data-lang="r"><span>pred_obj</span><span> </span><span>&lt;-</span><span> </span><span>predict</span><span>(</span><span>fit_obj</span><span>,</span><span> </span><span>data</span><span>,</span><span> </span><span>error.continuous</span><span> </span><span>=</span><span> </span><span>'VarExpl'</span><span>)</span><span>
</span><span>></span><span> </span><span>pred_obj</span><span>$</span><span>error</span><span>
</span><span>Variable</span><span> </span><span>Error</span><span> </span><span>ErrorType</span><span>
</span><span>1</span><span> </span><span>intrusion</span><span> </span><span>0.583</span><span> </span><span>VarExpl</span><span>
</span><span>2</span><span> </span><span>dreams</span><span> </span><span>0.590</span><span> </span><span>VarExpl</span><span>
</span><span>3</span><span> </span><span>flash</span><span> </span><span>0.513</span><span> </span><span>VarExpl</span><span>
</span><span>4</span><span> </span><span>upset</span><span> </span><span>0.615</span><span> </span><span>VarExpl</span><span>
</span><span>5</span><span> </span><span>physior</span><span> </span><span>0.601</span><span> </span><span>VarExpl</span><span>
</span><span>6</span><span> </span><span>avoidth</span><span> </span><span>0.648</span><span> </span><span>VarExpl</span><span>
</span><span>7</span><span> </span><span>avoidact</span><span> </span><span>0.626</span><span> </span><span>VarExpl</span><span>
</span><span>8</span><span> </span><span>amnesia</span><span> </span><span>0.327</span><span> </span><span>VarExpl</span><span>
</span><span>9</span><span> </span><span>lossint</span><span> </span><span>0.419</span><span> </span><span>VarExpl</span><span>
</span><span>10</span><span> </span><span>distant</span><span> </span><span>0.450</span><span> </span><span>VarExpl</span><span>
</span><span>11</span><span> </span><span>numb</span><span> </span><span>0.333</span><span> </span><span>VarExpl</span><span>
</span><span>12</span><span> </span><span>future</span><span> </span><span>0.450</span><span> </span><span>VarExpl</span><span>
</span><span>13</span><span> </span><span>sleep</span><span> </span><span>0.531</span><span> </span><span>VarExpl</span><span>
</span><span>14</span><span> </span><span>anger</span><span> </span><span>0.483</span><span> </span><span>VarExpl</span><span>
</span><span>15</span><span> </span><span>concen</span><span> </span><span>0.604</span><span> </span><span>VarExpl</span><span>
</span><span>16</span><span> </span><span>hyper</span><span> </span><span>0.602</span><span> </span><span>VarExpl</span><span>
</span><span>17</span><span> </span><span>startle</span><span> </span><span>0.605</span><span> </span><span>VarExpl</span></code></pre>
</figure>
<p>We calculated the percentage of variance explained in each of the nodes in the network. Next, we visualize the estimated network and discuss its structure in relation to explained variance.</p>
<h2>Visualize Network &amp; Predictability</h2>
<p>We provide the estimated weighted adjacency matrix and the node wise predictability measures as arguments to <code>qgraph()</code> …</p>
<figure>
<pre><code data-lang="r"><span>install.packages</span><span>(</span><span>'qgraph'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>qgraph</span><span>)</span><span>
</span><span>jpeg</span><span>(</span><span>paste0</span><span>(</span><span>figDir</span><span>,</span><span> </span><span>'McNallyNetwork.jpg'</span><span>),</span><span> </span><span>width</span><span> </span><span>=</span><span> </span><span>1500</span><span>,</span><span> </span><span>height</span><span> </span><span>=</span><span> </span><span>1500</span><span>)</span><span>
</span><span>qgraph</span><span>(</span><span>fit_obj</span><span>$</span><span>wadj</span><span>,</span><span> </span><span># weighted adjacency matrix as input
</span><span> </span><span>layout</span><span> </span><span>=</span><span> </span><span>'spring'</span><span>,</span><span> </span><span>pie</span><span> </span><span>=</span><span> </span><span>pred_obj</span><span>$</span><span>error</span><span>$</span><span>Error</span><span>,</span><span> </span><span># provide errors as input
</span><span> </span><span>pieColor</span><span> </span><span>=</span><span> </span><span>rep</span><span>(</span><span>'#377EB8'</span><span>,</span><span>p</span><span>),</span><span>
</span><span>node.color</span><span> </span><span>=</span><span> </span><span>fit_obj</span><span>$</span><span>edgecolor</span><span>,</span><span>
</span><span>labels</span><span> </span><span>=</span><span> </span><span>colnames</span><span>(</span><span>data</span><span>))</span><span>
</span><span>dev.off</span><span>()</span></code></pre>
</figure>
<p>… and get the following network visualization:</p>
<p><img alt="center" src="https://i2.wp.com/jmbh.github.io/figs/2016-11-01-Predictability-in-network-models/McNellyNetwork.jpg?w=456" data-recalc-dims="1"></p>
<p><a target="_blank" rel="nofollow" href="http://jmbh.github.io/Predictability-in-network-models/">[Click here for the original post with larger figures]</a></p>
<p>Each variable is represented by a node and the edges correspond to partial correlations, because in this dataset the MGM consists only of conditional Gaussian variables. The green color of the edges indicates that all partial correlations in this graph are positive, and the edge-width is proportional to the absolute value of the partial correlation. The blue pie chart behind the node indicates the predictability measure for each node (more blue = higher predictability).</p>
<p>We see that intrusive memories, traumatic dreams and flashbacks cluster together. Also, we observe that avoidance of thoughts (avoidth) about trauma interacts with avoidance of acitivies reminiscent of the trauma (avoidact) and that hypervigilant (hyper) behavior is related to feeling easily startled (startle). But there are also less obvious interactions, for instance between anger and concentration problems.</p>
<p>Now, if we would like to reduce sleep problems, the network model suggests to intervene on the variables anger and startle. But what the network structure does not tell us is <em>how much</em> we could possibly change sleep through the variables anger and startle. The predictability measure gives us an answer to this question: 53.1%. If the goal was to intervene on amnesia, we see that all adjacent nodes in the network explain only 32.7% of its variance. In addition, we see that there are many small edges connected to amnesia, suggesting that it is hard to intervene on amnesia via other nodes in the symptom network. Thus, one would possibly try to find additional variables that are not included in the network that interact with amnesia or try to intervene on amnesia directly.</p>
<h2>Limitations!</h2>
<p>Of course, there are limitations to interpreting explained variance as predicted treatment outcome: first, we cannot know the causal direction of the edges, so any edge could point in one or both directions. However, if there is no edge, there is also no causal effect in any direction. Also, it is often reasonable to combine the network model with general knoweldge: for instance, it seems more likely that amnesia causes being upset than the other way around. Second, we estimated the model on cross-sectional data (each row is one person) and hence assume that all people are the same, which is an assumption that is always violated to some extent. To solve this problem we would need (many) repeated measurements of a single person, in order to estimate a model specific to that person. This also solves the first problem to some degree as we can use the direction of time as the direction of causality. One would then use models that predict all symptoms at time point t by all symptoms at an earlier time point, let’s say t-1. An example of such a model is the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Vector_autoregression">Vector Autoregressive (VAR) model</a>.</p>
<h2>Compare Within vs. Out of Sample Predictability</h2>
<p>So far we looked into how well we can predict nodes by all other nodes within our sample. But in most situations we are interested in the predictability of nodes in new, unseen data. In what follows, we compare the within sample predictability with the out of sample predictability.</p>
<p>We first split the data in two parts: a training part (60% of the data), which we use to estimate the network model and a test part, which we will use to compute predictability measures on:</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>1</span><span>)</span><span>
</span><span>ind</span><span> </span><span>&lt;-</span><span> </span><span>sample</span><span>(</span><span>c</span><span>(</span><span>TRUE</span><span>,</span><span>FALSE</span><span>),</span><span> </span><span>prob</span><span>=</span><span>c</span><span>(</span><span>.6</span><span>,</span><span> </span><span>.4</span><span>),</span><span> </span><span>size</span><span>=</span><span>nrow</span><span>(</span><span>data</span><span>),</span><span> </span><span>replace</span><span>=</span><span>T</span><span>)</span><span> </span><span>#</span><span> </span><span>divide</span><span> </span><span>data</span><span> </span><span>in</span><span> </span><span>2</span><span> </span><span>parts</span><span> </span><span>(</span><span>60</span><span>% training set, 40%</span><span> </span><span>test</span><span> </span><span>set</span><span>)</span></code></pre>
</figure>
<p>Next, we estimate the network only on the training data and compute the predictability measure both on the training data and the test data:</p>
<figure>
<pre><code data-lang="r"><span>fit_obj_cv</span><span> </span><span>&lt;-</span><span> </span><span>mgmfit</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>data</span><span>[</span><span>ind</span><span>,],</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>rep</span><span>(</span><span>'g'</span><span>,</span><span> </span><span>p</span><span>),</span><span>
</span><span>lev</span><span> </span><span>=</span><span> </span><span>rep</span><span>(</span><span>1</span><span>,</span><span> </span><span>p</span><span>),</span><span>
</span><span>rule.reg</span><span> </span><span>=</span><span> </span><span>'OR'</span><span>)</span><span>
</span><span>pred_obj_train</span><span> </span><span>&lt;-</span><span> </span><span>predict</span><span>(</span><span>fit_obj_cv</span><span>,</span><span> </span><span>data</span><span>[</span><span>ind</span><span>,],</span><span> </span><span>error.continuous</span><span> </span><span>=</span><span> </span><span>'VarExpl'</span><span>)</span><span> </span><span># Compute Preditions on training data 60%
</span><span>pred_obj_test</span><span> </span><span>&lt;-</span><span> </span><span>predict</span><span>(</span><span>fit_obj_cv</span><span>,</span><span> </span><span>data</span><span>[</span><span>!</span><span>ind</span><span>,],</span><span> </span><span>error.continuous</span><span> </span><span>=</span><span> </span><span>'VarExpl'</span><span>)</span><span>#</span><span> </span><span>Compute</span><span> </span><span>Predictions</span><span> </span><span>on</span><span> </span><span>test</span><span> </span><span>data</span><span> </span><span>40</span><span>%</span></code></pre>
</figure>
<p>We now look at the mean predictability over nodes for the training- and test dataset:</p>
<figure>
<pre><code data-lang="r"><span>mean</span><span>(</span><span>pred_obj_train</span><span>$</span><span>error</span><span>$</span><span>Error</span><span>)</span><span> </span><span># mean explained variance training data
</span><span>[</span><span>1</span><span>]</span><span> </span><span>0.5384118</span><span>
</span><span>mean</span><span>(</span><span>pred_obj_test</span><span>$</span><span>error</span><span>$</span><span>Error</span><span>)</span><span> </span><span># mean explained variance test data
</span><span>[</span><span>1</span><span>]</span><span> </span><span>0.4494118</span></code></pre>
</figure>
<p>As to be expected, the explained variance is higher in the training dataset. This is because we fit the model to structure that is specific to the training data and is not present in the population (noise). Note that both means are lower than the mean we would get by taking the mean of the explained variances above, because we used less observation to estimate the model and hence have less power to detect edges.</p>
<p>While the explained variance values are lower in the test set, there is a strong correlation between the explained variance of a node in the training- and the test set</p>
<figure>
<pre><code data-lang="r"><span>cor</span><span>(</span><span>pred_obj_train</span><span>$</span><span>error</span><span>$</span><span>Error</span><span>,</span><span> </span><span>pred_obj_test</span><span>$</span><span>error</span><span>$</span><span>Error</span><span>)</span><span>
</span><span>[</span><span>1</span><span>]</span><span> </span><span>0.8018155</span></code></pre>
</figure>
<p>which means that if a node has high explained variance in the training set, it tends to also have a high explained variance in the test set.</p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/predictability-in-network-models/"></div></div><div><a data-url="https://www.r-bloggers.com/predictability-in-network-models/" data-counturl="https://www.r-bloggers.com/predictability-in-network-models/" data-count="vertical" data-text="Predictability in Network Models" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://jmbh.github.io//Predictability-in-network-models/"> Jonas Haslbeck - r</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/urFEtV32Phg"></a></div>