This tutorial covers how to prepare for, complete, and report correlation in R. Before you start this tutorial, you should make sure you review the relevant lecture, as this tutorial assumes you already know what correlation is.

Setting Up

All you need is this tutorial and RStudio. Remember that you can easily switch between windows with the Alt + ↹ Tab (Windows) and ⌘ Command + ⇥ Tab (Mac OS) shortcuts.

Task 1

Create a new week_05 project and open it in RStudio. Then, create two new folders in the new week_5 folder: r_docs and data. Finally, open a new Markdown file and save it in the r_docs folder. Since we will be practicing reporting (writing up results), we will need Markdown later on. For the tasks, get into the habit of creating new code chunks as you go.

Task 2

Add and run the command to load the tidyverse package in the setup code chunk.

Task 2.1

If you don’t have it yet, install the GGally package IN THE CONSOLE. Do not add this command to a code chunk! Then, load the GGally package in your setup code chunk.

library(tidyverse)
library(GGally)

Task 3

As usual, we will need some data to work with. Since there’s so much interesting information in it, let’s keep looking at the gensex dataset we’ve been using so far. Use the link below to read in the data in a new code chunk.

Link: https://and.netlify.app/datasets/gensex_clean.csv

#Two options

#First, save the data at that link to the `data` folder
#Then load the data from that file location
gensex <- readr::read_csv("../data/gensex_clean.csv")

#Alternatively, load the data in from the link directly
gensex <- readr::read_csv("https://and.netlify.app/datasets/gensex_clean.csv")

Data Preparation

As usual, before beginning any analysis, we will want to take a look at our data to make sure it’s clean and correct before we proceed. Since we’re using our familiar gensex dataset that we’ve cleaned already, we should be pretty happy with this, but let’s have a look again at these variables. The main reason is that you can generate rs all day easily, but if you don’t know what the variables you are correlating actually mean, that doesn’t do you much good!

Codebook

For our practice today, let’s focus on the ratings of either romantic or sexual attraction, since we already looked at gender in the lecture.

Task 4

Overwrite the gensex object so that it only contains variables to do with sexual or romantic attraction. Use the codebook above to find each variable, and make sure you understand how each variable was measured.

gensex <- gensex %>% 
  dplyr::select(.data = ., #optional!
    Sexual_fewq_1, Sexual_gender_1, Sexual_strength_1,
    Romantic_freq_1, Romantic_gender_1, Romantic_strength_1
  )

gensex <- gensex %>% 
  dplyr::select(starts_with(c("Sexual", "Romantic")))

Correlation Matrices

Basic Matrix with cor()

To start, we can get a very basic correlation matrix very quickly with the function cor() from the stats package, which is automatically loaded. It isn’t very pretty, but it will get us a rough and ready correlation matrix quickly.

Task 5

Create a correlation matrix of all six variables.

gensex %>% 
  cor()

                    Sexual_strength_1 Sexual_fewq_1 Sexual_gender_1
Sexual_strength_1           1.0000000    0.62565969              NA
Sexual_fewq_1               0.6256597    1.00000000              NA
Sexual_gender_1                    NA            NA               1
Romantic_strength_1                NA            NA              NA
Romantic_freq_1             0.3430913    0.48449978              NA
Romantic_gender_1           0.0979998   -0.01233245              NA
                    Romantic_strength_1 Romantic_freq_1
Sexual_strength_1                    NA      0.34309127
Sexual_fewq_1                        NA      0.48449978
Sexual_gender_1                      NA              NA
Romantic_strength_1                   1              NA
Romantic_freq_1                      NA      1.00000000
Romantic_gender_1                    NA      0.07396634
                    Romantic_gender_1
Sexual_strength_1          0.09799980
Sexual_fewq_1             -0.01233245
Sexual_gender_1                    NA
Romantic_strength_1                NA
Romantic_freq_1            0.07396634
Romantic_gender_1          1.00000000

Notice that there are a lot of NA values in this table. Why do you think this is? What can we do about it?

Task 6

Look at the gensex data to figure out why the correlation matrix contains NAs. Then, produce a correlation matrix that doesn’t have any NAs in it.

gensex %>% summary()

 Sexual_strength_1 Sexual_fewq_1   Sexual_gender_1
 Min.   :1.000     Min.   :1.000   Min.   :1.00   
 1st Qu.:7.000     1st Qu.:5.000   1st Qu.:5.00   
 Median :8.000     Median :6.000   Median :8.00   
 Mean   :7.526     Mean   :6.154   Mean   :6.64   
 3rd Qu.:9.000     3rd Qu.:8.000   3rd Qu.:9.00   
 Max.   :9.000     Max.   :9.000   Max.   :9.00   
                                   NA's   :3      
 Romantic_strength_1 Romantic_freq_1 Romantic_gender_1
 Min.   :2.000       Min.   :1.000   Min.   :1.000    
 1st Qu.:7.000       1st Qu.:4.000   1st Qu.:5.000    
 Median :8.000       Median :6.000   Median :8.000    
 Mean   :7.776       Mean   :5.922   Mean   :6.631    
 3rd Qu.:9.000       3rd Qu.:7.000   3rd Qu.:9.000    
 Max.   :9.000       Max.   :9.000   Max.   :9.000    
 NA's   :2                                            

gensex <- gensex %>% 
  filter(!is.na(Sexual_gender_1) & !is.na(Romantic_strength_1))

gensex %>% 
  cor()

                    Sexual_strength_1 Sexual_fewq_1 Sexual_gender_1
Sexual_strength_1          1.00000000   0.601693729     0.081864835
Sexual_fewq_1              0.60169373   1.000000000     0.001622784
Sexual_gender_1            0.08186484   0.001622784     1.000000000
Romantic_strength_1        0.50969282   0.252900937     0.022522256
Romantic_freq_1            0.31539071   0.470851441     0.109874779
Romantic_gender_1          0.08464036  -0.025282400     0.870330126
                    Romantic_strength_1 Romantic_freq_1
Sexual_strength_1            0.50969282      0.31539071
Sexual_fewq_1                0.25290094      0.47085144
Sexual_gender_1              0.02252226      0.10987478
Romantic_strength_1          1.00000000      0.51362193
Romantic_freq_1              0.51362193      1.00000000
Romantic_gender_1            0.01380141      0.06792014
                    Romantic_gender_1
Sexual_strength_1          0.08464036
Sexual_fewq_1             -0.02528240
Sexual_gender_1            0.87033013
Romantic_strength_1        0.01380141
Romantic_freq_1            0.06792014
Romantic_gender_1          1.00000000

We can see that there’s a huge amount of variation in how strongly these variables are correlated. We can also see that diagonal line of 1s, which represent perfect correlations.

Task 7

In your Markdown document, write your own explanation of why there are absolutely perfect correlations in your matrix.

As we covered in the lecture, a correlation matrix like this is the same above the diagonal line of 1s and below it. So, more than of half this big wall of numbers is just repeated information! Let’s use a different function to get a bit more bang for our buck, as it were.

Snazzy Matrix with GGally

The GGally package contains a very useful function called ggscatmat(), which we can use to help us better understand the relationships between our variables. It outputs a correlation matrix like the one we saw above, but mixed in with some really useful plots.

Task 8

Put all of your variables into the GGally::ggscatmat() function and have a look at the output.

gensex %>% 
  ggscatmat()

gensex %>% 
  ggscatmat() +
  theme(aspect.ratio=.5)

This plot has a lot going on, but don’t worry! That’s just because there’s a lot of useful information here. Let’s focus on one bit at a time.

First, we can see that we have our six variable names listed along the top of the plot, and then the same six names along the side. Where those columns and rows intersect, some cells of the plot contain numbers, some contain dots, and some contain lines. Let’s look first at the lines, which are all on the diagonal.

We saw before that cor() gave us 1s on the diagonal, representing each variable correlated with itself. Here, instead of 1s (which aren’t very helpful or interesting), we actually have a histogram of each variable. From the shapes that mostly slope dramatically up to the right, we can see that low ratings are uncommon, and high ratings very common, for most of our variables. The exceptions are the two variables asking about frequency (of sexual and romantic attraction), for which the most common response is 6 - 7 on our 9-point scale.

Next, let’s look at the numbers. You might recognise these values from cor() above. These are the correlation coefficients for each pair of variables. To figure out which variables are involved in the correlation, go up to the column name and along to the row name.

Finally, we have the dots. These are scatterplots of each pair of variables. Again, to find which variables are in each plot, look at the column name and the row name. So, each unique pair of variables has a scatterplot and a correlation coefficient; and each individual variable has a histogram along the diagonal. Like I said - snazzy!

Mini-Quiz

Answer the following questions using the output from ggscatmat().

We now can easily investigate the correlation between any two variables in our dataset. However, in the lecture we also talked about reporting significance, degrees of freedom, and CIs, which we can’t get easily from this output. So, if we want to report the results of particular correlations, we should have a look at a different function, cor.test().

Pairwise Tests with cor.test()

Unlike cor(), cor.test() only tests the correlation between pairs of variables. If we want to look at lots of variables at once, we’d be better off with cor() (or, although we haven’t looked at it in this tutorial, Hmisc::rcorr()). However, if we want more information about a particular correlation, such as easily reportable degrees of freedom, p-values, and confidence intervals, cor.test() is the better choice.

Task 14

Open the help documentation for the cor.test() function.

?cor.test()

Under the heading Usage, we can see that the function has two methods: one is called “Default S3 method” and the other is mentions “class ‘formula’”. We could use either method, but we’ll use the formula version now - because it’s going to be extremely useful in the future!

Under the heading Arguments, the help documentation says about the formula argument: “a formula of the form ~ u + v, where each of u and v are numeric variables giving the data values for one sample.” So, we need to write the first argument as ~ variable_name_1 + variable_name_2. We then also have a data argument to specify the data that those variables come from.

Finally, we have two other important arguments with no default:

• alternative must be set to one of the three options "two.sided", "less", or "greater". This determines the type of alternative hypothesis. Because we are going to be using two-tailed tests, we will stick with "two-sided".
• method must be set to one of the three options "pearson", "kendall", or "spearman". Thus far, we have only been using Pearson correlations, so we’ll stick with that for now.
• ... are other options you can change. We’ll stick with the defaults for now.

Task 15

Use cor.test() to calculate the correlation between Sexual_strength_1 and Romantic_strength_1 as the output shows below.

gensex %>% 
  cor.test(~ Sexual_strength_1 + Romantic_strength_1, data = .,
           alternative = "two.sided", method = "pearson")

## Alternative without a pipe
cor.test(~ Sexual_strength_1 + Romantic_strength_1, data = gensex,
           alternative = "two.sided", method = "pearson")

    Pearson's product-moment correlation

data:  Sexual_strength_1 and Romantic_strength_1
t = 10.244, df = 299, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4208932 0.5888164
sample estimates:
      cor 
0.5096928 

When we run this function, the output is some reasonably easy to read information about the results of our test. We can see here that we have:

• The type of analysis we asked for (Pearson’s product-moment correlation) and the variables we used
• The value of r under “sample estimates: cor”
• The values of t, degrees of freedom, and p
• The 95% confidence interval around r

Task 16

Stop and look at the values that the output has given you. Based on what you have learned from the lecture, write down what you think this value of r means, and how you could interpret it for these variables.

Task 16.1

Create a scatterplot of these two variables to help you understand the correlation. Does this change your interpretation at all?

gensex %>% 
  ggplot(aes(Sexual_strength_1, Romantic_strength_1)) +
  geom_point(position = "jitter", alpha = .4) +
  scale_x_continuous(breaks = 0:9,
                     name = "Strength of Sexual Attraction")+
  scale_y_continuous(breaks = 0:9,
                     name = "Strength of Romantic Attraction")+
  theme_minimal()

Task 17

Save the output from cor.test() as a new object called sex_rom_cor, so we can use it again later.

sex_rom_cor <- gensex %>% 
  cor.test(~ Sexual_strength_1 + Romantic_strength_1, data = .,
           alternative = "two.sided", method = "pearson")

Our reporting in the tasks above could use some numbers, which we can get out of the output of cor.test(). Let’s look at how we can formally report the results of a correlation analysis.

Reporting Correlations

To finish our analysis, we have to tell people about it - so we should write up the results of our correlation test as in a report.

Results like this correlation test, that have a small number of values to report, are often reported in brackets in the text of your report. This takes the general form:

To get these values, let’s look again at the results that we already saved in the object sex_rom_cor.

Task 18

Call the name of the object sex_rom_cor to see the displayed output.

sex_rom_cor

    Pearson's product-moment correlation

data:  Sexual_strength_1 and Romantic_strength_1
t = 10.244, df = 299, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4208932 0.5888164
sample estimates:
      cor 
0.5096928 

This is the same output that we saw above when we ran the analysis. Although this is useful for us, the analysts, to understand the result, it’s not at all formatted well for reports.

Task 19

Use the general format above and the values to type out the result of the correlation analysis.

To complete our report, we should describe in words what this result means. We essentially include the statistical result as a citation to give evidence for our claim.

Task 20

Have a go writing out a report of this statistical analysis.

That’s looking pretty slick! This is exactly the sort of thing you would expect to see in a journal article.

If you’re happy with this and/or ready to be done, you skip down to the Recap to finish up - this is the end of the required bit of the tutorial. Well done!

If you’d like to learn how to get R to do the reporting for you, read on!

OPTIONAL: Inline Coding

There are two ways to do our reporting: typing out the results by hand, as we just did above, or using inline coding. For this module, you are only required to be able to type out the results by hand, as we did above.

However, inline coding is an extremely powerful feature of Markdown documents, but it requires a bit more coding skill. So, if you’re ready to give it a go, let’s get started!

Task 21

Use class() and str() to look at the class and structure of the sex_rom_cor object.

sex_rom_cor %>% class()

[1] "htest"

sex_rom_cor %>% str()

List of 9
 $ statistic  : Named num 10.2
  ..- attr(*, "names")= chr "t"
 $ parameter  : Named int 299
  ..- attr(*, "names")= chr "df"
 $ p.value    : num 2.64e-21
 $ estimate   : Named num 0.51
  ..- attr(*, "names")= chr "cor"
 $ null.value : Named num 0
  ..- attr(*, "names")= chr "correlation"
 $ alternative: chr "two.sided"
 $ method     : chr "Pearson's product-moment correlation"
 $ data.name  : chr "Sexual_strength_1 and Romantic_strength_1"
 $ conf.int   : num [1:2] 0.421 0.589
  ..- attr(*, "conf.level")= num 0.95
 - attr(*, "class")= chr "htest"

The first command tells us that this object, which contains the results from the function cor.test(), has the class htest. What does this mean? Objects of this class tend to have a specific structure, which we can see using str(). Now we can see that there is a lot more stored in here than it seems. If we want to get out specific values, we can simply subset them using $, as the output from str() suggests (notice the $s in the output). The output from cor.test() that we have here contains the following useful elements (among others):

• statistic: the calculated value of t
• parameter: the degrees of freedom for t
• p.value: the p-value
• estimate: the value of r
• conf.int: the values of the lower and upper bounds of the confidence interval for r

Task 22

Use subsetting to extract the value of r from the sex_rom_cor object.

sex_rom_cor$estimate

      cor 
0.5096928 

That’s very cool, but we can do a bit better…

Task 22.1

Use subsetting to extract the value of r from the sex_rom_cor object, then round it to two decimal places.

sex_rom_cor$estimate %>% 
  round(2)

 cor 
0.51 

That’s looking very promising! Now, if I want to report the correlation coefficient for this analysis as part of my writeup, I could use inline code to do this. Inline code is a small snippet of R code inside the text part of a Markdown document, no different from the code we use inside code chunks. The difference is, we can write lines of code that print out particular numbers - as we just did above - and then insert that code into the text where we want that number to appear. When we knit the document, R will run the inline code and replace it with whatever value the code produces. To do this, I just need to write `r some_r_code` in my Markdown.

If that’s clear as mud, it might help to see it in action. In the text of my Markdown document, I could write:

“Ratings of the strength of sexual and romantic attraction were correlated at r = `r sex_rom_cor$estimate %>% round(2)`”

When I knit the Markdown document, this will come out as:

“Ratings of the strength of sexual and romantic attraction were correlated at r = 0.51”

The code `r sex_rom_cor$estimate %>% round(2)` is evaluated by R to produce a number, and then that number appears in place of the code in the knitted document. Cool, eh?

Once we have the basic idea, we can get almost any value of interest from our sex_rom_cor object using subsetting. The only other thing you need to know is that if subsetting gives you multiple numbers, you can get out just one of them using square brackets with a number. For example, this code: object_name$variable_name[1] will give you the first observation in that variable.

Task 23

Use subsetting with square brackets to get out the lower and higher bounds of the confidence intervals from sex_rom_cor.

sex_rom_cor$conf.int[1]

[1] 0.4208932

sex_rom_cor$conf.int[2]

[1] 0.5888164

Task 24

Complete the report of this analysis below, replacing the square brackets with inline code that produces the correct number. Make sure you round all values to two decimal places, if necessary, except for p, which is rounded to three.

Strength of sexual attraction and strength of romantic attraction were significantly correlated, r([degrees of freedom]) = [value of r], 95% CI = [[lower bound], [upper bound]], p = [p-value].

Task 25

Knit your document to see your completed inline code!

Recap

Whew, well done! You’re welcome - and encouraged! - to keep working with the gensex data to practice what we’ve learned today.

In sum, we’ve covered:

• How to prepare and visualise data for correlation
• How to run the analysis, read the output, and store it in an object
• How to report and interpret the statistical results of these tests in real-world terms
  • Optionally, how to use the subsetting operator $ to report results with inline code

Remember, if you get stuck or you have questions, post them on Piazza, or bring them to StatsChats or to drop-ins.

Good job!

That’s all for today. See you soon!