R code Chapter 7
This document contains abridged sections from Discovering Statistics Using R and RStudio by Andy Field so there are some copyright considerations. You can use this material for teaching and non-profit activities but please do not meddle with it or claim it as your own work. See the full license terms at the bottom of the page.
Load the data
Remember to install the package for the book with library(discovr). After which you can load data into a tibble by executing:
name_of_tib <- discovr::name_of_data
For example, execute these lines to create the tibbles referred to in the chapter:
cats_tib <- discovr::roaming_cats
exam_tib <- discovr::exam_anxiety
liar_tib <- discovr::biggest_liar
If you want to read the file from the CSV and you have set up your project folder as I suggest in Chapter 1, then the general format of the code you would use is:
tibble_name <- here::here("../data/name_of_file.csv") %>%
readr::read_csv() %>%
dplyr::mutate(
...
code to convert character variables to factors
...
)
In which you’d replace tibble_name with the name you want to assign to the tibble and change name_of_file.csv to the name of the file that you’re trying to load. You can use mutate to convert categorical variables to factors. For example, for the exam_anxiety data you would load it by executing:
exam_tib <- here::here("data/exam_anxiety.csv") %>%
readr::read_csv() %>%
dplyr::mutate(
sex = forcats::as_factor(sex)
)
Bivariate correlations
Plot the data
GGally::ggscatmat(exam_tib, columns = c("exam_grade", "revise", "anxiety")) +
theme_minimal()
Pearson’s r
A single correlation:
exam_tib %>%
dplyr::select(exam_grade, revise) %>%
correlation::correlation()
## Parameter1 | Parameter2 | r | 95% CI | t | df | p | Method | n_Obs
## -------------------------------------------------------------------------------------
## exam_grade | revise | 0.40 | [0.22, 0.55] | 4.34 | 101 | < .001 | Pearson | 103
Let’s scale this up to include exam anxiety:
exam_cor <- exam_tib %>%
dplyr::select(exam_grade, revise, anxiety) %>%
correlation::correlation()
exam_cor
## Parameter1 | Parameter2 | r | 95% CI | t | df | p | Method | n_Obs
## ------------------------------------------------------------------------------------------
## exam_grade | revise | 0.40 | [ 0.22, 0.55] | 4.34 | 101 | < .001 | Pearson | 103
## exam_grade | anxiety | -0.44 | [-0.58, -0.27] | -4.94 | 101 | < .001 | Pearson | 103
## revise | anxiety | -0.71 | [-0.79, -0.60] | -10.11 | 101 | < .001 | Pearson | 103
Change the number of decimal places
exam_tib %>%
dplyr::select(exam_grade, revise, anxiety) %>%
correlation::correlation(digits = 3, ci_digits = 3)
## Parameter1 | Parameter2 | r | 95% CI | t | df | p | Method | n_Obs
## ----------------------------------------------------------------------------------------------
## exam_grade | revise | 0.397 | [ 0.220, 0.548] | 4.343 | 101 | < .001 | Pearson | 103
## exam_grade | anxiety | -0.441 | [-0.585, -0.271] | -4.938 | 101 | < .001 | Pearson | 103
## revise | anxiety | -0.709 | [-0.794, -0.598] | -10.111 | 101 | < .001 | Pearson | 103
Matrix format
exam_tib %>%
dplyr::select(exam_grade, revise, anxiety) %>%
correlation::correlation() %>%
summary()
## Parameter | anxiety | revise
## -------------------------------
## exam_grade | -0.44*** | 0.40***
## revise | -0.71*** |
Coefficient of determination
exam_cor <- exam_tib %>%
dplyr::select(exam_grade, revise, anxiety) %>%
correlation::correlation()
(exam_cor$r)^2
## [1] 0.1573873 0.1944752 0.5030345
Robust correlations
exam_tib %>%
dplyr::select(exam_grade, revise, anxiety) %>%
WRS2::winall()
## Call:
## WRS2::winall(x = .)
##
## Robust correlation matrix:
## exam_grade revise anxiety
## exam_grade 1.0000 0.3085 -0.3913
## revise 0.3085 1.0000 -0.6015
## anxiety -0.3913 -0.6015 1.0000
##
## p-values:
## exam_grade revise anxiety
## exam_grade NA 0.00183 7e-05
## revise 0.00183 NA 0e+00
## anxiety 0.00007 0.00000 NA
exam_tib %>%
dplyr::select(exam_grade, revise, anxiety) %>%
correlation::correlation(., method = "percentage")
## Parameter1 | Parameter2 | r | 95% CI | t | df | p | Method | n_Obs
## -------------------------------------------------------------------------------------------------
## exam_grade | revise | 0.34 | [ 0.15, 0.50] | 3.60 | 101 | < .001 | Percentage Bend | 103
## exam_grade | anxiety | -0.40 | [-0.55, -0.23] | -4.41 | 101 | < .001 | Percentage Bend | 103
## revise | anxiety | -0.61 | [-0.72, -0.47] | -7.66 | 101 | < .001 | Percentage Bend | 103
Spearman’s rho
Load the Biggest Liar data from the CSV file (assuming you’ve set up a project as described in chapter 1)
liar_tib <- readr::read_csv("../data/biggest_liar.csv") %>%
dplyr::mutate(
novice = forcats::as_factor(novice)
)
liar_tib %>%
dplyr::select(position, creativity) %>%
correlation::correlation(method = "spearman")
## Parameter1 | Parameter2 | rho | 95% CI | S | p | Method | n_Obs
## --------------------------------------------------------------------------------------
## position | creativity | -0.38 | [-0.56, -0.15] | 72123.51 | 0.002 | Spearman | 68
Kendall’s tau
liar_tib %>%
dplyr::select(position, creativity) %>%
correlation::correlation(method = "kendall")
## Parameter1 | Parameter2 | 95% CI | tau | z | p | Method | n_Obs
## ----------------------------------------------------------------------------------
## position | creativity | [-0.50, -0.07] | -0.30 | -3.24 | 0.001 | Kendall | 68
Bootstrapped confidence intervals
A quick aside, let’s look at the code for writing a function to print the mean:
print_mean <- function(variable){
mean <- sum(variable)/nrow(variable)
cat("Mean = ", mean)
}
exam_tib %>%
dplyr::select(exam_grade) %>%
print_mean()
## Mean = 56.57282
Or a sillier version:
print_mean <- function(harry_the_hungy_hippo){
mean <- sum(harry_the_hungy_hippo)/ nrow(harry_the_hungy_hippo)
cat("Harry the hungry hippo say that the mean = ", mean)
}
exam_tib %>%
dplyr::select(exam_grade) %>%
print_mean()
## Harry the hungry hippo say that the mean = 56.57282
Here’s a function to bootstrap Pearson’s r:
boot_r <- function(data, i){
cor(data[i, "exam_grade"], data[i, "revise"])
}
grade_revise_bs <- boot::boot(exam_tib, boot_r, R = 2000)
grade_revise_bs
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot::boot(data = exam_tib, statistic = boot_r, R = 2000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.3967207 -0.001476077 0.06983568
boot::boot.ci(grade_revise_bs)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot::boot.ci(boot.out = grade_revise_bs)
##
## Intervals :
## Level Normal Basic
## 95% ( 0.2613, 0.5351 ) ( 0.2673, 0.5424 )
##
## Level Percentile BCa
## 95% ( 0.2510, 0.5261 ) ( 0.2530, 0.5263 )
## Calculations and Intervals on Original Scale
A general bootstrap function
boot_r <- function(data, var1, var2, i){
cor(data[i, var1], data[i, var2])
}
grade_revise_bs <- boot::boot(exam_tib, boot_r, var1 = "exam_grade", var2 = "revise", R = 2000)
grade_revise_bs
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot::boot(data = exam_tib, statistic = boot_r, R = 2000, var1 = "exam_grade",
## var2 = "revise")
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.3967207 -0.001593843 0.06788372
boot::boot.ci(grade_revise_bs)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot::boot.ci(boot.out = grade_revise_bs)
##
## Intervals :
## Level Normal Basic
## 95% ( 0.2653, 0.5314 ) ( 0.2732, 0.5405 )
##
## Level Percentile BCa
## 95% ( 0.2529, 0.5202 ) ( 0.2518, 0.5198 )
## Calculations and Intervals on Original Scale
grade_anx_bs <- boot::boot(exam_tib, boot_r, var1 = "exam_grade", var2 = "anxiety", R = 2000)
grade_anx_bs
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot::boot(data = exam_tib, statistic = boot_r, R = 2000, var1 = "exam_grade",
## var2 = "anxiety")
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -0.4409934 0.002270827 0.0644994
boot::boot.ci(grade_anx_bs)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot::boot.ci(boot.out = grade_anx_bs)
##
## Intervals :
## Level Normal Basic
## 95% (-0.5697, -0.3168 ) (-0.5767, -0.3235 )
##
## Level Percentile BCa
## 95% (-0.5585, -0.3053 ) (-0.5609, -0.3088 )
## Calculations and Intervals on Original Scale
revise_anx_bs <- boot::boot(exam_tib, boot_r, var1 = "revise", var2 = "anxiety", R = 2000)
revise_anx_bs
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot::boot(data = exam_tib, statistic = boot_r, R = 2000, var1 = "revise",
## var2 = "anxiety")
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -0.7092493 0.001022998 0.108373
boot::boot.ci(revise_anx_bs)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot::boot.ci(boot.out = revise_anx_bs)
##
## Intervals :
## Level Normal Basic
## 95% (-0.9227, -0.4979 ) (-0.9525, -0.5483 )
##
## Level Percentile BCa
## 95% (-0.8702, -0.4660 ) (-0.8530, -0.3760 )
## Calculations and Intervals on Original Scale
Bootrapping Spearman’s rho and Kendall’s tau
boot_r <- function(data, var1, var2, method = "pearson", i){
cor(data[i, var1], data[i, var2], method = method)
}
liar_bs <- boot::boot(liar_tib, boot_r, var1 = "position", var2 = "creativity", method = "spearman", R = 2000)
boot::boot.ci(liar_bs)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot::boot.ci(boot.out = liar_bs)
##
## Intervals :
## Level Normal Basic
## 95% (-0.6204, -0.1376 ) (-0.6294, -0.1555 )
##
## Level Percentile BCa
## 95% (-0.5976, -0.1238 ) (-0.5818, -0.1042 )
## Calculations and Intervals on Original Scale
Point-biserial correlation
If you need to read the data in from the csv file:
cats_tib <- readr::read_csv("../data/roaming_cats.csv") %>%
dplyr::mutate(
sex = forcats::as_factor(sex)
)
cats_tib <- cats_tib %>%
dplyr::mutate(
sex_bin = ifelse(sex == "Male", 0, 1),
sex_bin_recode = ifelse(sex == "Male", 1, 0)
)
cats_tib %>%
dplyr::select(time, sex_bin) %>%
correlation::correlation()
## Parameter1 | Parameter2 | r | 95% CI | t | df | p | Method | n_Obs
## ---------------------------------------------------------------------------------------
## time | sex_bin | -0.38 | [-0.58, -0.14] | -3.11 | 58 | 0.003 | Pearson | 60
coefficient of determination
r_cat <- cats_tib %>%
dplyr::select(time, sex_bin) %>%
correlation::correlation()
(r_cat$r)^2
## [1] 0.1432276
cats_tib %>%
dplyr::select(time, sex_bin_recode) %>%
correlation::correlation()
## Parameter1 | Parameter2 | r | 95% CI | t | df | p | Method | n_Obs
## ---------------------------------------------------------------------------------------
## time | sex_bin_recode | 0.38 | [0.14, 0.58] | 3.11 | 58 | 0.003 | Pearson | 60
cats_tib %>%
dplyr::select(time, sex_bin) %>%
correlation::correlation(method = "biserial")
## Parameter1 | Parameter2 | rho | 95% CI | t | df | p | Method | n_Obs
## -----------------------------------------------------------------------------------------
## time | sex_bin | -0.47 | [-0.65, -0.25] | -4.03 | 58 | < .001 | Biserial | 60
Partial correlations
exam_tib %>%
dplyr::select(exam_grade, revise, anxiety) %>%
correlation::correlation(., partial = TRUE)
## Parameter1 | Parameter2 | r | 95% CI | t | df | p | Method | n_Obs
## -----------------------------------------------------------------------------------------
## exam_grade | revise | 0.13 | [-0.06, 0.32] | 1.35 | 101 | 0.182 | Pearson | 103
## exam_grade | anxiety | -0.25 | [-0.42, -0.06] | -2.56 | 101 | 0.024 | Pearson | 103
## revise | anxiety | -0.65 | [-0.75, -0.52] | -8.56 | 101 | < .001 | Pearson | 103
Comparing rs
The correlation between exam grade and anxiety for self-identifying males:
exam_tib %>%
dplyr::filter(sex == "Male") %>%
dplyr::select(exam_grade, anxiety) %>%
correlation::correlation()
## Parameter1 | Parameter2 | r | 95% CI | t | df | p | Method | n_Obs
## ----------------------------------------------------------------------------------------
## exam_grade | anxiety | -0.51 | [-0.68, -0.27] | -4.14 | 50 | < .001 | Pearson | 52
The correlation between exam grade and anxiety for self-identifying females:
exam_tib %>%
dplyr::filter(sex == "Female") %>%
dplyr::select(exam_grade, anxiety) %>%
correlation::correlation()
## Parameter1 | Parameter2 | r | 95% CI | t | df | p | Method | n_Obs
## ---------------------------------------------------------------------------------------
## exam_grade | anxiety | -0.38 | [-0.59, -0.12] | -2.89 | 49 | 0.006 | Pearson | 51
exam_male <- exam_tib %>%
dplyr::filter(sex == "Male")
exam_female <- exam_tib %>%
dplyr::filter(sex == "Female")
WRS2::twopcor(x1 = exam_male$exam_grade, y1 = exam_male$anxiety, x2 = exam_female$exam_grade, y2 = exam_female$anxiety)
## Call:
## WRS2::twopcor(x1 = exam_male$exam_grade, y1 = exam_male$anxiety,
## x2 = exam_female$exam_grade, y2 = exam_female$anxiety)
##
## First correlation coefficient: -0.5057
## Second correlation coefficient: -0.3814
## Confidence interval (difference): -0.4574 0.1315