<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/variables-can-synergize-even-in-a-linear-model/"></div></div><div><a data-url="https://www.r-bloggers.com/variables-can-synergize-even-in-a-linear-model/" data-counturl="https://www.r-bloggers.com/variables-can-synergize-even-in-a-linear-model/" data-count="horizontal" data-text="Variables can synergize, even in a linear model" 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://www.win-vector.com/blog/2016/09/variables-can-synergize-even-in-a-linear-model/"> R – Win-Vector Blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Introduction</h2>
<p>Suppose we have the task of predicting an outcome <code>y</code> given a number of variables <code>v1,..,vk</code>. We often want to “prune variables” or build models with fewer than all the variables. This can be to speed up modeling, decrease the cost of producing future data, improve robustness, improve explain-ability, even reduce over-fit, and improve the quality of the resulting model.</p>
<p>For some informative discussion on such issues please see the following:</p>
<ul>
<li><a target="_blank" href="http://www.win-vector.com/blog/2015/08/how-do-you-know-if-your-data-has-signal/">How Do You Know if Your Data Has Signal?</a></li>
<li><a target="_blank" href="http://www.win-vector.com/blog/2015/09/isyourmodelgoingtowork/">How do you know if your model is going to work?</a></li>
<li><a target="_blank" href="http://www.win-vector.com/blog/2016/08/variable-pruning-is-np-hard/">Variable pruning is NP hard</a></li>
</ul>
<p>In this article we are going to deliberately (and artificially) find and test one of the limits of the technique. We recommend simple variable pruning, but also think it is important to be aware of its limits.</p>
<p><span></span></p>
<p>To be truly effective in applied fields (such as data science) one often has to use (with care) methods that “happen to work” in addition to methods that “are known to always work” (or at least be aware, you are always competing against such); hence the interest in mere heuristic.</p>
<h2>The pruning heuristics</h2>
<p>Let <code>L(y;m;v1,...,vk)</code> denote the estimate loss (or badness of performance, so smaller is better) of a model for <code>y</code> fit using modeling method <code>m</code> and the variables <code>v1,...,vk</code>. Let <code>d(a;L(y;m;),L(y;m;a))</code> denote the portion of <code>L(y;m;)-L(y;m;a)</code> credited to the variable <code>a</code>. This could be the change in loss, something like <code>effectsize(a)</code>, or <code>-log(significance(a))</code>; in all cases <em>larger</em> is considered better.</p>
<p>For practical variable pruning (during predictive modeling) our intuition often implicitly relies on the following heuristic arguments.</p>
<ul>
<li><code>L(y;m;)</code> is monotone decreasing, we expect <code>L(y;m;v1,...,vk,a)</code> is no larger than <code>L(y;m;v1,...,vk,)</code>. Note this may be achievable “in sample” (or on training data), but is often false if <code>L(y;m;)</code> accounts for model complexity or is estimated on out of sample data (itself a good practice).</li>
<li>If <code>L(y;m;v1,...,vk,a)</code> is significantly lower than <code>L(y;m;v1,...,vk)</code> then we will be lucky enough have <code>d(a;L(y;m;),L(y;m;a))</code> not too small.</li>
<li>If <code>d(a;L(y;m;),L(y;m;a))</code> is not too small then we will be lucky enough to have <code>d(a;L(y;lm;),L(y;lm;a))</code> is non-negligible (where modeling method <code>lm</code> is one linear regression or logistic regression).</li>
</ul>
<p>Intuitively we are <em>hoping</em> variable utility has a roughly diminishing return structure and at least some non-vanishing fraction of a variable’s utility can be seen in simple linear or generalized linear models. Obviously this can not be true in general (interactions in decision trees being a well know situation where variable utility can increase in the presence of other variables, and there are many non-linear relations that escape detection by linear models).</p>
<p>However, if the above were true (or often nearly true) we could effectively prune variables by keeping only the set of variables <code>{ a | d(a;L(y;lm;),L(y;lm;a)) is non negligible}</code>. This is a (user controllable) heuristic built into our <a target="_blank" href="https://cran.r-project.org/package=vtreat"><code>vtreat</code></a> <a target="_blank" href="https://cran.r-project.org/">R</a> package and proves to be quite useful in practice.</p>
<p>I’ll repeat: we feel in real world data you can use the above heuristics to usefully prune variables. Complex models do eventually get into a regime of diminishing returns, and real world engineered useful variables usually (by design) have a hard time hiding. Also, remember data science is an empirical field- methods that happen to work will dominate (even if they do not apply in all cases).</p>
<h2>Counter-examples</h2>
<p>For every heuristic you should crisply know if it is true (and is in fact a theorem) or it is false (and has counter-examples). We stand behind the above heuristics, and will show their empirical worth in a follow-up article. Let’s take some time and show that they are not in fact laws.</p>
<p>We are going to show that per-variable coefficient significances and effect sizes are not monotone in that adding more variables can in fact improve them.</p>
<h3>First example</h3>
<p>First (using <a target="_blank" href="https://cran.r-project.org/">R</a>) we build a data frame where <code>y = a xor b</code>. This is a classic example of <code>y</code> being a function of two variable but not a <em>linear</em> function of them (at least over the real numbers, it is a linear relation over the field <a target="_blank" href="https://en.wikipedia.org/wiki/GF(2)">GF(2)</a>).</p>
<pre><code>d &lt;- data.frame(a=c(0,0,1,1),b=c(0,1,0,1))
d$y &lt;- as.numeric(d$a == d$b)</code></pre>
<p>We look at the (real) linear relations between <code>y</code> and <code>a</code>, <code>b</code>.</p>
<pre><code>summary(lm(y~a+b,data=d))</code></pre>
<pre><code>## 
## Call:
## lm(formula = y ~ a + b, data = d)
## 
## Residuals:
##    1    2    3    4 
##  0.5 -0.5 -0.5  0.5 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    0.500      0.866   0.577    0.667
## a              0.000      1.000   0.000    1.000
## b              0.000      1.000   0.000    1.000
## 
## Residual standard error: 1 on 1 degrees of freedom
## Multiple R-squared:  3.698e-32,  Adjusted R-squared:     -2 
## F-statistic: 1.849e-32 on 2 and 1 DF,  p-value: 1</code></pre>
<pre><code>anova(lm(y~a+b,data=d))</code></pre>
<pre><code>## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value Pr(>F)
## a          1      0       0       0      1
## b          1      0       0       0      1
## Residuals  1      1       1</code></pre>
<p>As we expect linear methods fail to find any evidence of a relation between <code>y</code> and <code>a</code>, <code>b</code>. This clearly violates our hoped for heuristics.</p>
<p>For details on reading these summaries we strongly recommend <a target="_blank" href="https://cran.r-project.org/doc/contrib/Faraway-PRA.pdf"><em>Practical Regression and Anova using R</em>, Julian J. Faraway, 2002</a>.</p>
<p>In this example the linear model fails to recognize <code>a</code> and <code>b</code> as useful variables (even though <code>y</code> is a function of <code>a</code> and <code>b</code>). From the linear model’s point of view variables are not improving each other (so that at least looks monotone), but it is largely because the linear model can not see the relation unless we were to add an interaction of <code>a</code> and <code>b</code> (denoted <code>a:b</code>).</p>
<h3>Second example</h3>
<p>Let us develop this example a bit more to get a more interesting counterexample.</p>
<p>Introduce new variables <code>u = a and b</code>, <code>v = a or b</code>. By the rules of logic we have <code>y == 1+u-v</code>, so there is a linear relation.</p>
<pre><code>d$u &lt;- as.numeric(d$a &amp; d$b)
d$v &lt;- as.numeric(d$a | d$b)
print(d)</code></pre>
<pre><code>##   a b y u v
## 1 0 0 1 0 0
## 2 0 1 0 0 1
## 3 1 0 0 0 1
## 4 1 1 1 1 1</code></pre>
<pre><code>print(all.equal(d$y,1+d$u-d$v))</code></pre>
<pre><code>## [1] TRUE</code></pre>
<p>We can now see the counter-example effect: together the variables work better than they did alone.</p>
<pre><code>summary(lm(y~u,data=d))</code></pre>
<pre><code>## 
## Call:
## lm(formula = y ~ u, data = d)
## 
## Residuals:
##          1          2          3          4 
##  6.667e-01 -3.333e-01 -3.333e-01 -1.388e-16 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.3333     0.3333       1    0.423
## u             0.6667     0.6667       1    0.423
## 
## Residual standard error: 0.5774 on 2 degrees of freedom
## Multiple R-squared:  0.3333, Adjusted R-squared:  5.551e-16 
## F-statistic:     1 on 1 and 2 DF,  p-value: 0.4226</code></pre>
<pre><code>anova(lm(y~u,data=d))</code></pre>
<pre><code>## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value Pr(>F)
## u          1 0.33333 0.33333       1 0.4226
## Residuals  2 0.66667 0.33333</code></pre>
<pre><code>summary(lm(y~v,data=d))</code></pre>
<pre><code>## 
## Call:
## lm(formula = y ~ v, data = d)
## 
## Residuals:
##          1          2          3          4 
##  5.551e-17 -3.333e-01 -3.333e-01  6.667e-01 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   1.0000     0.5774   1.732    0.225
## v            -0.6667     0.6667  -1.000    0.423
## 
## Residual standard error: 0.5774 on 2 degrees of freedom
## Multiple R-squared:  0.3333, Adjusted R-squared:      0 
## F-statistic:     1 on 1 and 2 DF,  p-value: 0.4226</code></pre>
<pre><code>anova(lm(y~v,data=d))</code></pre>
<pre><code>## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value Pr(>F)
## v          1 0.33333 0.33333       1 0.4226
## Residuals  2 0.66667 0.33333</code></pre>
<pre><code>summary(lm(y~u+v,data=d))</code></pre>
<pre><code>## Warning in summary.lm(lm(y ~ u + v, data = d)): essentially perfect fit:
## summary may be unreliable</code></pre>
<pre><code>## 
## Call:
## lm(formula = y ~ u + v, data = d)
## 
## Residuals:
##          1          2          3          4 
## -1.849e-32  7.850e-17 -7.850e-17  1.849e-32 
## 
## Coefficients:
##              Estimate Std. Error    t value Pr(>|t|)    
## (Intercept)  1.00e+00   1.11e-16  9.007e+15   &lt;2e-16 ***
## u            1.00e+00   1.36e-16  7.354e+15   &lt;2e-16 ***
## v           -1.00e+00   1.36e-16 -7.354e+15   &lt;2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.11e-16 on 1 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 4.056e+31 on 2 and 1 DF,  p-value: &lt; 2.2e-16</code></pre>
<pre><code>anova(lm(y~u+v,data=d))</code></pre>
<pre><code>## Warning in anova.lm(lm(y ~ u + v, data = d)): ANOVA F-tests on an
## essentially perfect fit are unreliable</code></pre>
<pre><code>## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq    F value    Pr(>F)    
## u          1 0.33333 0.33333 2.7043e+31 &lt; 2.2e-16 ***
## v          1 0.66667 0.66667 5.4086e+31 &lt; 2.2e-16 ***
## Residuals  1 0.00000 0.00000                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<p>In this example we see synergy instead of diminishing returns. Each variable becomes better in the presence of the other. This is on its own good, but indicates variable pruning is harder than one might expect- even for a linear model.</p>
<h3>Third example</h3>
<p>We can get around the above warnings by adding some rows to the data frame that don’t follow the designed relation. We can even draw rows from this frame to show the effect on a “more row independent looking” data frame.</p>
<pre><code>d0 &lt;- d
d0$y &lt;- 0
d1 &lt;- d
d1$y &lt;- 1
dG &lt;- rbind(d,d,d,d,d0,d1)
set.seed(23235)
dR &lt;- dG[sample.int(nrow(dG),100,replace=TRUE),,drop=FALSE]

summary(lm(y~u,data=dR))</code></pre>
<pre><code>## 
## Call:
## lm(formula = y ~ u, data = dR)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8148 -0.3425 -0.3425  0.3033  0.6575 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.34247    0.05355   6.396 5.47e-09 ***
## u            0.47235    0.10305   4.584 1.35e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4575 on 98 degrees of freedom
## Multiple R-squared:  0.1765, Adjusted R-squared:  0.1681 
## F-statistic: 21.01 on 1 and 98 DF,  p-value: 1.349e-05</code></pre>
<pre><code>anova(lm(y~u,data=dR))</code></pre>
<pre><code>## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## u          1  4.3976  4.3976   21.01 1.349e-05 ***
## Residuals 98 20.5124  0.2093                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<pre><code>summary(lm(y~v,data=dR))</code></pre>
<pre><code>## 
## Call:
## lm(formula = y ~ v, data = dR)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7619 -0.3924 -0.3924  0.6076  0.6076 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.7619     0.1049   7.263 9.12e-11 ***
## v            -0.3695     0.1180  -3.131   0.0023 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4807 on 98 degrees of freedom
## Multiple R-squared:  0.09093,    Adjusted R-squared:  0.08165 
## F-statistic: 9.802 on 1 and 98 DF,  p-value: 0.002297</code></pre>
<pre><code>anova(lm(y~v,data=dR))</code></pre>
<pre><code>## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## v          1  2.265 2.26503  9.8023 0.002297 **
## Residuals 98 22.645 0.23107                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<pre><code>summary(lm(y~u+v,data=dR))</code></pre>
<pre><code>## 
## Call:
## lm(formula = y ~ u + v, data = dR)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8148 -0.1731 -0.1731  0.1984  0.8269 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.76190    0.08674   8.784 5.65e-14 ***
## u            0.64174    0.09429   6.806 8.34e-10 ***
## v           -0.58883    0.10277  -5.729 1.13e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3975 on 97 degrees of freedom
## Multiple R-squared:  0.3847, Adjusted R-squared:  0.3721 
## F-statistic: 30.33 on 2 and 97 DF,  p-value: 5.875e-11</code></pre>
<pre><code>anova(lm(y~u+v,data=dR))</code></pre>
<pre><code>## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## u          1  4.3976  4.3976  27.833 8.047e-07 ***
## v          1  5.1865  5.1865  32.826 1.133e-07 ***
## Residuals 97 15.3259  0.1580                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<h3>Conclusion</h3>
<p>Consider the above counter example as <a target="_blank" href="https://en.wikipedia.org/wiki/Exception_that_proves_the_rule"><em>exceptio probat regulam in casibus non exceptis</em> (“the exception confirms the rule in cases not excepted”)</a>. Or roughly outlining the (hopefully labored and uncommon) structure needed to break the otherwise common and useful heuristics.</p>
<p>In later articles in this series we will show more about the structure of model quality and show the above heuristics actually working very well in practice (and adding a lot of value to projects).</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/variables-can-synergize-even-in-a-linear-model/"></div></div><div><a data-url="https://www.r-bloggers.com/variables-can-synergize-even-in-a-linear-model/" data-counturl="https://www.r-bloggers.com/variables-can-synergize-even-in-a-linear-model/" data-count="vertical" data-text="Variables can synergize, even in a linear model" 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://www.win-vector.com/blog/2016/09/variables-can-synergize-even-in-a-linear-model/"> R – Win-Vector Blog</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/TMQDlPNfFzw"></a></div>