<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/an-easy-way-to-accidentally-inflate-reported-r-squared-in-linear-regression-models/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.win-vector.com/blog/2017/06/an-easy-way-to-accidentally-inflate-reported-r-squared-in-linear-regression-models/"> R – Win-Vector Blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Here is an absolutely <em>horrible</em> way to confuse yourself and get an inflated reported <code>R-squared</code> on a simple linear regression model in <a rel="nofollow" href="https://www.r-project.org/" target="_blank"><code>R</code></a>.</p>
<p>We have written about this before, but we found a new twist on the problem (interactions with categorical variable encoding) which we would like to call out here.<span></span></p>
<p>First let’s set up our problem with a data set where the quantity to be predicted (<code>y</code>) has no real relation to the independent variable (<code>x</code>). We will first build our example data:</p>
<div>
<pre><code><span>library</span>(<span>&quot;sigr&quot;</span>)
<span>library</span>(<span>&quot;broom&quot;</span>)
<span>set.seed</span>(<span>23255</span>)

d &lt;-<span> </span><span>data.frame</span>(<span>y=</span> <span>runif</span>(<span>100</span>),
                <span>x=</span> <span>ifelse</span>(<span>runif</span>(<span>100</span>)>=<span>0.5</span>, <span>'a'</span>, <span>'b'</span>),
                <span>stringsAsFactors =</span> <span>FALSE</span>)</code></pre>
</div>
<p>Now let’s build a model and look at the summary statistics returned as part of the model fitting process:</p>
<div>
<pre><code>m1 &lt;-<span> </span><span>lm</span>(y~x, <span>data=</span>d)
<span>t</span>(broom::<span>glance</span>(m1))</code></pre>
</div>
<pre><code>##                        [,1]
## r.squared       0.002177326
## adj.r.squared  -0.008004538
## sigma           0.302851476
## statistic       0.213843593
## p.value         0.644796456
## df              2.000000000
## logLik        -21.432440763
## AIC            48.864881526
## BIC            56.680392084
## deviance        8.988463618
## df.residual    98.000000000</code></pre>
<div>
<pre><code>d$pred1 &lt;-<span> </span><span>predict</span>(m1, <span>newdata =</span> d)</code></pre>
</div>
<p>I <em>strongly</em> prefer to directly calculate the the model performance statistics off the predictions (it lets us easily compare different modeling methods), so let’s also do that also:</p>
<div>
<pre><code><span>cat</span>(<span>render</span>(sigr::<span>wrapFTest</span>(d, <span>'pred1'</span>, <span>'y'</span>),
           <span>format=</span><span>'markdown'</span>))</code></pre>
</div>
<p><strong>F Test</strong> summary: (<i>R<sup>2</sup></i>=0.0022, <em>F</em>(1,98)=0.21, <em>p</em>=n.s.).</p>
<p>So far so good. Let’s now remove the &quot;intercept term&quot; by adding the &quot;<code>0+</code>&quot; from the fitting command.</p>
<div>
<pre><code>m2 &lt;-<span> </span><span>lm</span>(y~<span>0</span>+x, <span>data=</span>d)
<span>t</span>(broom::<span>glance</span>(m2))</code></pre>
</div>
<pre><code>##                        [,1]
## r.squared      7.524811e-01
## adj.r.squared  7.474297e-01
## sigma          3.028515e-01
## statistic      1.489647e+02
## p.value        1.935559e-30
## df             2.000000e+00
## logLik        -2.143244e+01
## AIC            4.886488e+01
## BIC            5.668039e+01
## deviance       8.988464e+00
## df.residual    9.800000e+01</code></pre>
<div>
<pre><code>d$pred2 &lt;-<span> </span><span>predict</span>(m2, <span>newdata =</span> d)</code></pre>
</div>
<p>Uh oh. That <em>appeared</em> to vastly improve the reported <code>R-squared</code> and the significance (&quot;<code>p.value</code>&quot;)!</p>
<p>That does not make sense, anything <code>m2</code> can do <code>m1</code> can also do. In fact the two models make essentially identical predictions, which we confirm below:</p>
<div>
<pre><code><span>sum</span>((d$pred1 -<span> </span>d$y)^<span>2</span>)</code></pre>
</div>
<pre><code>## [1] 8.988464</code></pre>
<div>
<pre><code><span>sum</span>((d$pred2 -<span> </span>d$y)^<span>2</span>)</code></pre>
</div>
<pre><code>## [1] 8.988464</code></pre>
<div>
<pre><code><span>max</span>(<span>abs</span>(d$pred1 -<span> </span>d$pred2))</code></pre>
</div>
<pre><code>## [1] 4.440892e-16</code></pre>
<div>
<pre><code><span>head</span>(d)</code></pre>
</div>
<pre><code>##             y x     pred1     pred2
## 1 0.007509118 b 0.5098853 0.5098853
## 2 0.980353615 a 0.5380361 0.5380361
## 3 0.055880927 b 0.5098853 0.5098853
## 4 0.993814410 a 0.5380361 0.5380361
## 5 0.636617385 b 0.5098853 0.5098853
## 6 0.154032277 a 0.5380361 0.5380361</code></pre>
<p>Let’s double check the fit quality of the predictions.</p>
<div>
<pre><code><span>cat</span>(<span>render</span>(sigr::<span>wrapFTest</span>(d, <span>'pred2'</span>, <span>'y'</span>),
           <span>format=</span><span>'markdown'</span>))</code></pre>
</div>
<p><strong>F Test</strong> summary: (<i>R<sup>2</sup></i>=0.0022, <em>F</em>(1,98)=0.21, <em>p</em>=n.s.).</p>
<p>Ah. The prediction fit quality is identical to the first time (as one would expect). This is yet another reason to directly calculate model fit quality from the predictions: it isolates you from any foibles of how the modeling software does it.</p>
<p>The answer to our puzzles of &quot;what went wrong&quot; is something we have written about before <a rel="nofollow" href="http://www.win-vector.com/blog/2016/12/be-careful-evaluating-model-predictions/" target="_blank">here</a>.</p>
<p>Roughly what is going on is:</p>
<p>If the fit formula sent <code>lm()</code> has no intercept (triggered by the &quot;<code>0+</code>&quot;) notation then <code>summary.lm()</code> changes how it computes <code>r.squared</code> as follows (from <code>help(summary.lm)</code>):</p>
<pre>
r.squared   R^2, the ‘fraction of variance explained by the model’,
              R^2 = 1 - Sum(R[i]^2) / Sum((y[i]- y*)^2),
            where y* is the mean of y[i] if there is an intercept and zero otherwise.
</pre>
<p>This is pretty bad.</p>
<p>Then to add insult to injury the &quot;<code>0+</code>&quot; notation also changes how <code>R</code> encodes the categorical variable <code>x</code>.</p>
<p>Compare:</p>
<div>
<pre><code><span>summary</span>(m1)</code></pre>
</div>
<pre><code>## 
## Call:
## lm(formula = y ~ x, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53739 -0.23265 -0.02039  0.27247  0.47111 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.53804    0.04515  11.918   &lt;2e-16 ***
## xb          -0.02815    0.06088  -0.462    0.645    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3029 on 98 degrees of freedom
## Multiple R-squared:  0.002177,   Adjusted R-squared:  -0.008005 
## F-statistic: 0.2138 on 1 and 98 DF,  p-value: 0.6448</code></pre>
<div>
<pre><code><span>summary</span>(m2)</code></pre>
</div>
<pre><code>## 
## Call:
## lm(formula = y ~ 0 + x, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53739 -0.23265 -0.02039  0.27247  0.47111 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## xa  0.53804    0.04515   11.92   &lt;2e-16 ***
## xb  0.50989    0.04084   12.49   &lt;2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3029 on 98 degrees of freedom
## Multiple R-squared:  0.7525, Adjusted R-squared:  0.7474 
## F-statistic:   149 on 2 and 98 DF,  p-value: &lt; 2.2e-16</code></pre>
<p>Notice the second model directly encodes both levels of <code>x</code>. This means that if <code>m1</code> has <code>pred1 = a + b * (x=='b')</code> we can reproduce this model (intercept and all) as <code>m2</code>: <code>pred2 = a * (x=='a') + (a+b) * (x=='b')</code>. I.e., the invariant <code>(x=='a') + (x=='b') == 1</code> means <code>m2</code> can imitate the model with the intercept term.</p>
<p>The presumed (and I think weak) justification of <code>summary.lm()</code> switching the model quality assessment method is something along the lines that <code>mean(y)</code> may not be in the model’s concept space and this might lead to reporting negative <code>R-squared</code>. I don’t have any problem with negative <code>R-squared</code>, it can be taken to mean you did worse than the unconditional average. However, even if you accept the (no-warning) scoring method switch: that argument doesn’t apply here. <code>m2</code> can imitate having an intercept, so it isn’t unfair to check if it is better than using only the intercept.</p>
	<p>
<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/2017/06/an-easy-way-to-accidentally-inflate-reported-r-squared-in-linear-regression-models/"> 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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/sYX7A27lL3s" height="1"></a></div>