<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/exploring-correlations-in-r-with-corrr/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://drsimonj.svbtle.com/exploring-correlations-in-r-with-corrr"> blogR</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><a rel="nofollow" href="https://twitter.com/drsimonj" target="_blank">@drsimonj</a> here to share a (sort of) readable version of <a rel="nofollow" href="https://www.meetup.com/en-AU/amst-R-dam/events/251102944/" target="_blank">my presentation at the amst-R-dam meetup</a> on 14 August, 2018: �Exploring correlations in R with corrr???.</p>
<p>Those who attended will know that I changed the topic of the talk, originally advertised as �R from academia to commerical business???. For anyone who�s interested, I gave that talk at useR! 2018 and, thanks to the R consortium, you can watch it <a rel="nofollow" href="https://www.youtube.com/embed/3eqJj7mj7lA" target="_blank">here</a>. I also gave a �Wrangling data in the Tidyverse??? tutorial that you can follow at <a rel="nofollow" href="https://www.youtube.com/embed/E-Vvg8uzcVM" target="_blank">Part 1</a> and <a rel="nofollow" href="https://www.youtube.com/embed/DwWH1mTerOc" target="_blank">Part 2</a>.</p>
<h1>
<a rel="nofollow" href="https://drsimonj.svbtle.com/exploring-correlations-in-r-with-corrr#the-story-of-corrr_1" target="_blank">�</a>The story of corrr</h1>
<p>Moving to corrr � the first package I ever created. It started when I was a postgrad student studying individual differences in decision making. My research data was responses to test batteries. My statistical bread and butter was regression-based techniques like multiple regression, path analysis, factor analysis (EFA and CFA), and structural equation modelling.</p>
<p>I spent a lot of time exploring correlation matrices to make model decisions, and diagnose poor fits or unexpected results! If you need proof, check out some of the correlations tables published in my academic papers like <a rel="nofollow" href="https://onlinelibrary.wiley.com/doi/abs/10.1002/bdm.1939" target="_blank">�Individual Differences in Decision Making Depend on Cognitive Abilities, Monitoring and Control???</a></p>
<h1>
<a rel="nofollow" href="https://drsimonj.svbtle.com/exploring-correlations-in-r-with-corrr#how-to-explore-correlations_1" target="_blank">�</a>How to explore correlations?</h1>
<p>To illustrate some of the challenges I was facing, let�s try explore some correlations with some very fancy data:</p>
<pre><code>d &lt;- mtcars
d$hp[3] &lt;- NA
head(d)
#>                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
#> Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
#> Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
#> Datsun 710        22.8   4  108  NA 3.85 2.320 18.61  1  1    4    1
#> Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
#> Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
#> Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1
</code></pre>
<p>We could be motivated by <a rel="nofollow" href="https://en.wikipedia.org/wiki/Multicollinearity" target="_blank">multicollinearity</a>:</p>
<pre><code>fit_1 &lt;- lm(mpg ~ hp,        data = d)
fit_2 &lt;- lm(mpg ~ hp + disp, data = d)
</code></pre>
<pre><code>summary(fit_1)
#> 
#> Call:
#> lm(formula = mpg ~ hp, data = d)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -5.7219 -2.2887 -0.8677  1.5815  8.1743 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 30.18457    1.69527  17.805  &lt; 2e-16 ***
#> hp          -0.06860    0.01039  -6.604 3.08e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.925 on 29 degrees of freedom
#>   (1 observation deleted due to missingness)
#> Multiple R-squared:  0.6006, Adjusted R-squared:  0.5868 
#> F-statistic: 43.61 on 1 and 29 DF,  p-value: 3.085e-07
</code></pre>
<pre><code>summary(fit_2)
#> 
#> Call:
#> lm(formula = mpg ~ hp + disp, data = d)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -4.9136 -2.3472 -0.7302  1.9798  6.7549 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 30.963202   1.373067  22.550  &lt; 2e-16 ***
#> hp          -0.024840   0.013480  -1.843 0.075991 .  
#> disp        -0.030992   0.007504  -4.130 0.000296 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.149 on 28 degrees of freedom
#>   (1 observation deleted due to missingness)
#> Multiple R-squared:  0.7518, Adjusted R-squared:  0.7341 
#> F-statistic: 42.41 on 2 and 28 DF,  p-value: 3.368e-09
</code></pre>
<p>Strange result. Let�s check the correlations between <code>mpg</code>, <code>hp</code>, and <code>disp</code> to try and diagnose this problem. It should be simple using the base R function, <code>cor()</code>. Right?</p>
<p>Err, what is with all the <code>NA</code>�s ?</p>
<pre><code>rs &lt;- cor(d)
rs
#>             mpg        cyl       disp hp        drat         wt
#> mpg   1.0000000 -0.8521620 -0.8475514 NA  0.68117191 -0.8676594
#> cyl  -0.8521620  1.0000000  0.9020329 NA -0.69993811  0.7824958
#> disp -0.8475514  0.9020329  1.0000000 NA -0.71021393  0.8879799
#> hp           NA         NA         NA  1          NA         NA
#> drat  0.6811719 -0.6999381 -0.7102139 NA  1.00000000 -0.7124406
#> wt   -0.8676594  0.7824958  0.8879799 NA -0.71244065  1.0000000
#> qsec  0.4186840 -0.5912421 -0.4336979 NA  0.09120476 -0.1747159
#> vs    0.6640389 -0.8108118 -0.7104159 NA  0.44027846 -0.5549157
#> am    0.5998324 -0.5226070 -0.5912270 NA  0.71271113 -0.6924953
#> gear  0.4802848 -0.4926866 -0.5555692 NA  0.69961013 -0.5832870
#> carb -0.5509251  0.5269883  0.3949769 NA -0.09078980  0.4276059
#>             qsec         vs          am       gear        carb
#> mpg   0.41868403  0.6640389  0.59983243  0.4802848 -0.55092507
#> cyl  -0.59124207 -0.8108118 -0.52260705 -0.4926866  0.52698829
#> disp -0.43369788 -0.7104159 -0.59122704 -0.5555692  0.39497686
#> hp            NA         NA          NA         NA          NA
#> drat  0.09120476  0.4402785  0.71271113  0.6996101 -0.09078980
#> wt   -0.17471588 -0.5549157 -0.69249526 -0.5832870  0.42760594
#> qsec  1.00000000  0.7445354 -0.22986086 -0.2126822 -0.65624923
#> vs    0.74453544  1.0000000  0.16834512  0.2060233 -0.56960714
#> am   -0.22986086  0.1683451  1.00000000  0.7940588  0.05753435
#> gear -0.21268223  0.2060233  0.79405876  1.0000000  0.27407284
#> carb -0.65624923 -0.5696071  0.05753435  0.2740728  1.00000000
</code></pre>
<p>Check the help page <a rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/stats/html/cor.html" target="_blank"><code>?cor</code></a>. Not so obvious. Default is <code>use = &quot;everything&quot;</code>, and buried down in the details:</p>
<blockquote>
<p>If use is �everything???, NAs will propagate conceptually, i.e., a resulting value will be NA whenever one of its contributing observations is NA.</p>
</blockquote>
<p>Have to handle missing values with <code>use</code>:</p>
<pre><code>rs &lt;- cor(d, use = &quot;pairwise.complete.obs&quot;)
rs
#>             mpg        cyl       disp         hp        drat         wt
#> mpg   1.0000000 -0.8521620 -0.8475514 -0.7749843  0.68117191 -0.8676594
#> cyl  -0.8521620  1.0000000  0.9020329  0.8298025 -0.69993811  0.7824958
#> disp -0.8475514  0.9020329  1.0000000  0.7860001 -0.71021393  0.8879799
#> hp   -0.7749843  0.8298025  0.7860001  1.0000000 -0.44258306  0.6505927
#> drat  0.6811719 -0.6999381 -0.7102139 -0.4425831  1.00000000 -0.7124406
#> wt   -0.8676594  0.7824958  0.8879799  0.6505927 -0.71244065  1.0000000
#> qsec  0.4186840 -0.5912421 -0.4336979 -0.7064809  0.09120476 -0.1747159
#> vs    0.6640389 -0.8108118 -0.7104159 -0.7162015  0.44027846 -0.5549157
#> am    0.5998324 -0.5226070 -0.5912270 -0.2196179  0.71271113 -0.6924953
#> gear  0.4802848 -0.4926866 -0.5555692 -0.1161968  0.69961013 -0.5832870
#> carb -0.5509251  0.5269883  0.3949769  0.7437843 -0.09078980  0.4276059
#>             qsec         vs          am       gear        carb
#> mpg   0.41868403  0.6640389  0.59983243  0.4802848 -0.55092507
#> cyl  -0.59124207 -0.8108118 -0.52260705 -0.4926866  0.52698829
#> disp -0.43369788 -0.7104159 -0.59122704 -0.5555692  0.39497686
#> hp   -0.70648093 -0.7162015 -0.21961793 -0.1161968  0.74378428
#> drat  0.09120476  0.4402785  0.71271113  0.6996101 -0.09078980
#> wt   -0.17471588 -0.5549157 -0.69249526 -0.5832870  0.42760594
#> qsec  1.00000000  0.7445354 -0.22986086 -0.2126822 -0.65624923
#> vs    0.74453544  1.0000000  0.16834512  0.2060233 -0.56960714
#> am   -0.22986086  0.1683451  1.00000000  0.7940588  0.05753435
#> gear -0.21268223  0.2060233  0.79405876  1.0000000  0.27407284
#> carb -0.65624923 -0.5696071  0.05753435  0.2740728  1.00000000
</code></pre>
<p>Can we focus on subset with dplyr? Nope.</p>
<pre><code>dplyr::select(rs, mpg, hp, disp)
#> Error in UseMethod(&quot;select_&quot;): no applicable method for 'select_' applied to an object of class &quot;c('matrix', 'double', 'numeric')&quot;
</code></pre>
<p>Riiiiiight! It�s a matrix and dplyr is for data frames.</p>
<pre><code>class(rs)
#> [1] &quot;matrix&quot;
</code></pre>
<p>So we can use square brackets with matrices? Or not�</p>
<pre><code>vars &lt;- c(&quot;mpg&quot;, &quot;hp&quot;, &quot;disp&quot;)
rs[rownames(rs) %in% vars]
#>  [1]  1.0000000 -0.8475514 -0.7749843 -0.8521620  0.9020329  0.8298025
#>  [7] -0.8475514  1.0000000  0.7860001 -0.7749843  0.7860001  1.0000000
#> [13]  0.6811719 -0.7102139 -0.4425831 -0.8676594  0.8879799  0.6505927
#> [19]  0.4186840 -0.4336979 -0.7064809  0.6640389 -0.7104159 -0.7162015
#> [25]  0.5998324 -0.5912270 -0.2196179  0.4802848 -0.5555692 -0.1161968
#> [31] -0.5509251  0.3949769  0.7437843
</code></pre>
<p>Mm, square brackets can take on different functions with matrices. Without a comma, it�s treated like a vector. With a comma, we can separately specify the dimensions.</p>
<pre><code>vars &lt;- c(&quot;mpg&quot;, &quot;hp&quot;, &quot;disp&quot;)
rs[rownames(rs) %in% vars, colnames(rs) %in% vars]
#>             mpg       disp         hp
#> mpg   1.0000000 -0.8475514 -0.7749843
#> disp -0.8475514  1.0000000  0.7860001
#> hp   -0.7749843  0.7860001  1.0000000
</code></pre>
<p>Aha! High correlation between input variables (multicollinearity).</p>
<p>But seriously, this syntax is pretty ugly.</p>
<pre><code>vars &lt;- c(&quot;mpg&quot;, &quot;hp&quot;, &quot;disp&quot;)
rs[rownames(rs) %in% vars, colnames(rs) %in% vars]
</code></pre>
<p>We diagnosed our multicollinearity problem. What if we want to something a bit more complex like exploring clustering of variables in high dimensional space? Could use exploratory factor analysis.</p>
<pre><code>factanal(na.omit(d), factors = 2)
#> 
#> Call:
#> factanal(x = na.omit(d), factors = 2)
#> 
#> Uniquenesses:
#>   mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb 
#> 0.155 0.072 0.101 0.145 0.295 0.174 0.154 0.264 0.160 0.245 0.387 
#> 
#> Loadings:
#>      Factor1 Factor2
#> mpg   0.694  -0.603 
#> cyl  -0.629   0.729 
#> disp -0.730   0.605 
#> hp   -0.335   0.862 
#> drat  0.810  -0.221 
#> wt   -0.809   0.414 
#> qsec -0.162  -0.905 
#> vs    0.280  -0.811 
#> am    0.910   0.110 
#> gear  0.859   0.133 
#> carb          0.781 
#> 
#>                Factor1 Factor2
#> SS loadings      4.506   4.342
#> Proportion Var   0.410   0.395
#> Cumulative Var   0.410   0.804
#> 
#> Test of the hypothesis that 2 factors are sufficient.
#> The chi square statistic is 68.49 on 34 degrees of freedom.
#> The p-value is 0.000414
</code></pre>
<pre><code>factanal(na.omit(d), factors = 5)
#> 
#> Call:
#> factanal(x = na.omit(d), factors = 5)
#> 
#> Uniquenesses:
#>   mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb 
#> 0.112 0.044 0.012 0.066 0.275 0.005 0.091 0.128 0.136 0.143 0.045 
#> 
#> Loadings:
#>      Factor1 Factor2 Factor3 Factor4 Factor5
#> mpg   0.647  -0.441  -0.483  -0.203         
#> cyl  -0.604   0.694   0.293   0.128         
#> disp -0.645   0.543   0.202   0.451   0.182 
#> hp   -0.264   0.659   0.544   0.205   0.305 
#> drat  0.806  -0.250                         
#> wt   -0.725   0.226   0.420   0.488         
#> qsec -0.184  -0.890  -0.269                 
#> vs    0.242  -0.848  -0.209  -0.155   0.163 
#> am    0.916   0.118                         
#> gear  0.892           0.237                 
#> carb  0.119   0.458   0.853                 
#> 
#>                Factor1 Factor2 Factor3 Factor4 Factor5
#> SS loadings      4.186   3.255   1.739   0.582   0.180
#> Proportion Var   0.381   0.296   0.158   0.053   0.016
#> Cumulative Var   0.381   0.676   0.835   0.887   0.904
#> 
#> Test of the hypothesis that 5 factors are sufficient.
#> The chi square statistic is 2.17 on 10 degrees of freedom.
#> The p-value is 0.995
</code></pre>
<p>So many questions! I�d much rather explore the correlations.</p>
<p>Let�s try to find all variables with a correlation greater than 0.90. Why doesn�t this work?!</p>
<pre><code>col_has_over_90 &lt;- apply(rs, 2, function(x) any(x > .9))
rs[, col_has_over_90]
#>             mpg        cyl       disp         hp        drat         wt
#> mpg   1.0000000 -0.8521620 -0.8475514 -0.7749843  0.68117191 -0.8676594
#> cyl  -0.8521620  1.0000000  0.9020329  0.8298025 -0.69993811  0.7824958
#> disp -0.8475514  0.9020329  1.0000000  0.7860001 -0.71021393  0.8879799
#> hp   -0.7749843  0.8298025  0.7860001  1.0000000 -0.44258306  0.6505927
#> drat  0.6811719 -0.6999381 -0.7102139 -0.4425831  1.00000000 -0.7124406
#> wt   -0.8676594  0.7824958  0.8879799  0.6505927 -0.71244065  1.0000000
#> qsec  0.4186840 -0.5912421 -0.4336979 -0.7064809  0.09120476 -0.1747159
#> vs    0.6640389 -0.8108118 -0.7104159 -0.7162015  0.44027846 -0.5549157
#> am    0.5998324 -0.5226070 -0.5912270 -0.2196179  0.71271113 -0.6924953
#> gear  0.4802848 -0.4926866 -0.5555692 -0.1161968  0.69961013 -0.5832870
#> carb -0.5509251  0.5269883  0.3949769  0.7437843 -0.09078980  0.4276059
#>             qsec         vs          am       gear        carb
#> mpg   0.41868403  0.6640389  0.59983243  0.4802848 -0.55092507
#> cyl  -0.59124207 -0.8108118 -0.52260705 -0.4926866  0.52698829
#> disp -0.43369788 -0.7104159 -0.59122704 -0.5555692  0.39497686
#> hp   -0.70648093 -0.7162015 -0.21961793 -0.1161968  0.74378428
#> drat  0.09120476  0.4402785  0.71271113  0.6996101 -0.09078980
#> wt   -0.17471588 -0.5549157 -0.69249526 -0.5832870  0.42760594
#> qsec  1.00000000  0.7445354 -0.22986086 -0.2126822 -0.65624923
#> vs    0.74453544  1.0000000  0.16834512  0.2060233 -0.56960714
#> am   -0.22986086  0.1683451  1.00000000  0.7940588  0.05753435
#> gear -0.21268223  0.2060233  0.79405876  1.0000000  0.27407284
#> carb -0.65624923 -0.5696071  0.05753435  0.2740728  1.00000000
</code></pre>
<p>The diagonal is 1. All cols have a value greater than .90!</p>
<p>Exclude diagonal:</p>
<pre><code>diag(rs) &lt;- NA
col_has_over_90 &lt;- apply(rs, 2, function(x) any(x > .9, na.rm = TRUE))
rs[, col_has_over_90]
#>             cyl       disp
#> mpg  -0.8521620 -0.8475514
#> cyl          NA  0.9020329
#> disp  0.9020329         NA
#> hp    0.8298025  0.7860001
#> drat -0.6999381 -0.7102139
#> wt    0.7824958  0.8879799
#> qsec -0.5912421 -0.4336979
#> vs   -0.8108118 -0.7104159
#> am   -0.5226070 -0.5912270
#> gear -0.4926866 -0.5555692
#> carb  0.5269883  0.3949769
</code></pre>
<p>Again, this syntax is pretty gross. Imagine showing this to a beginner and asking them to write down as much as they remember. Probably not much would be my guess.</p>
<p>What about vizualising correlations? I�d suggest giving up at this point.</p>
<p><a rel="nofollow" href="https://i1.wp.com/svbtleusercontent.com/vCAdmmx28nbzoKYcGVrh210xspap.gif?ssl=1" target="_blank"><img alt="depressed.gif" src="https://i2.wp.com/svbtleusercontent.com/vCAdmmx28nbzoKYcGVrh210xspap_small.gif?w=456&ssl=1" data-recalc-dims="1"></a></p>
<h2>
<a rel="nofollow" href="https://drsimonj.svbtle.com/exploring-correlations-in-r-with-corrr#exploring-data-with-the-tidyverse_2" target="_blank">�</a>Exploring data with the tidyverse</h2>
<p>Remember me as postgrad? I�d discovered the tidyverse and really liked it, because <em>exploring</em> data with the tidyverse is easy.</p>
<pre><code>library(tidyverse)
d %>% 
  select(mpg:drat) %>% 
  gather() %>% 
  ggplot(aes(value)) +
    geom_histogram() +
    facet_wrap(~key, scales = &quot;free&quot;)
</code></pre>
<p><a rel="nofollow" href="https://i2.wp.com/svbtleusercontent.com/v8a6SMH3qEUWDAnK2fbtr60xspap.png?ssl=1" target="_blank"><img alt="unnamed-chunk-17-1.png" src="https://i2.wp.com/svbtleusercontent.com/v8a6SMH3qEUWDAnK2fbtr60xspap_small.png?w=456&ssl=1" data-recalc-dims="1"></a></p>
<p>Can�t we have this for correlations? I don�t want to do any crazy mathematical operations or statistical tests. I just want to quickly explore the value. It�s not a big ask.</p>
<p>Good news! This is why I developed corrr as a tidyverse-style package for exploring correlations in R.</p>
<h1>
<a rel="nofollow" href="https://drsimonj.svbtle.com/exploring-correlations-in-r-with-corrr#a-hrefhttpgithubcomdrsimonjcorrrcorrra_1" target="_blank">�</a><a rel="nofollow" href="http://github.com/drsimonj/corrr/" target="_blank">corrr</a><br>
</h1>
<p><a rel="nofollow" href="https://i0.wp.com/svbtleusercontent.com/oDU6kAGErMQCpLiuDJoTRX0xspap.png?ssl=1" target="_blank"><img alt="corrr_logo.png" src="https://i1.wp.com/svbtleusercontent.com/oDU6kAGErMQCpLiuDJoTRX0xspap_small.png?w=456&ssl=1" data-recalc-dims="1"></a></p>
<p>Here�s a quick example to get a feel for the syntax:</p>
<pre><code>library(corrr)
d %>% 
  correlate() %>% 
  focus(mpg:drat, mirror = TRUE) %>% 
  network_plot()
#> 
#> Correlation method: 'pearson'
#> Missing treated using: 'pairwise.complete.obs'
</code></pre>
<p><a rel="nofollow" href="https://i0.wp.com/svbtleusercontent.com/u5EHjmYYpN2zSdm5i3ZTri0xspap.png?ssl=1" target="_blank"><img alt="unnamed-chunk-18-1.png" src="https://i0.wp.com/svbtleusercontent.com/u5EHjmYYpN2zSdm5i3ZTri0xspap_small.png?w=456&ssl=1" data-recalc-dims="1"></a></p>
<p>The first objective of corrr was to write a new function that uses <code>cor()</code> but:</p>
<ul>
<li> Handles missing by default with <code>use = &quot;pairwise.complete.obs&quot;</code>
</li>
<li> Stops diagonal from getting in the way by setting it to <code>NA</code>
</li>
<li> Makes it possible to use tidyverse tools by returning a <code>data.frame</code> instead of a <code>matrix</code>
</li>
</ul>
<p>Now, meet <code>correlate()</code></p>
<pre><code>rs &lt;- correlate(d)
#> 
#> Correlation method: 'pearson'
#> Missing treated using: 'pairwise.complete.obs'
rs
#> # A tibble: 11 x 12
#>    rowname     mpg     cyl    disp      hp     drat      wt     qsec
#>    </code></pre>