Lab09W - Correlation & Regression

Loading Packages & Data

Load tidyverse to get access to ggplot2 and dplyr. Load ggpmisc to add regression equation. Load broom to convert model outputs to dataframes. Load gridExtra and ggExtra for additional ggplot features such as marginal distributions on scatterplot.

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}
packages <- c("tidyverse", "ggpmisc", "broom", "gridExtra", "ggExtra")
ipak(packages)

Load in the data for the exercise.
NOTE: You will need to change this path on your machine

ex8_data <- read.csv("/Users/Julian/GDrive/Misc/Classes/InterStats/Ex8_LabData.csv") %>% tbl_df()

Make light theme default (white backgronud) for ggplot

old <- theme_set(theme_light(base_size = 12))

Scatterplot

Create a scatterplot with CESD scores at Time 1 (CESD1) and Time 2 (CESD2). Begin the plot with ggplot(dataframe, aes(X variable, Y variable)). Then add the scatterplot by adding the geom_point() layer.

ggplot(ex8_data, aes(CESD1, CESD2)) + 
  geom_point()

Add a LOESS line (i.e., smooth / running average) by adding the geom_smooth() layer.

ggplot(ex8_data, aes(CESD1, CESD2)) + 
  geom_point() + 
  geom_smooth()

Correlation

Compute the correlation between CESD1 and CESD2 using the cor(Variable 1, Variable 2) function.

cor(ex8_data$CESD1, ex8_data$CESD2)

[1] 0.4845713

Test the correlation between CESD1 and CESD2 using the cor.test(Variable 1, Variable 2) function.

cor.test(ex8_data$CESD1, ex8_data$CESD2)

    Pearson's product-moment correlation

data:  ex8_data$CESD1 and ex8_data$CESD2
t = 7.5954, df = 188, p-value = 0.0000000000014
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3675733 0.5864698
sample estimates:
      cor 
0.4845713 

BONUS: Boostrapped Confidence Interval

This code is more complicated, but here’s the general sequence of events: set the random number generator to today’s date; specify the number of simulations (boostrap resamplings) we are conducting (in this case 1,000); boostrap our data; test the correlation on each replication; rename the variable name to “pearson_r”; arrange these estimates from smallest to largest; isolate the 2.5% and 97.5% percentiles; round these estimates to four decimals; then add a column to specify the lower and upper bound.

set.seed(10312016)
nSims <- 1000
ex8_data %>% 
  bootstrap(nSims) %>% 
  do(tidy(cor(.$CESD1, .$CESD2))) %>% 
  ungroup() %>% select(pearson_r = x) %>% 
  arrange(pearson_r) %>% 
  slice( c(nSims*.025, nSims*.975)) %>% 
  mutate_all(funs(round(.,4))) %>% 
  add_column(bound = c("lower", "upper"))

Regression

Now we will conduct a linear regression of CESD time 1 scores predicting CESD at time 2. Here we use the lm() function to create a linear model. The first argument takes a formula in model notation (Y variable ~ X variables; outcomes ~ predictors). The second argument takes a dataframe. We assign this function to a model object called “m1”. We then use summary() to retrieve key details about the linear model we just estimated.

NOTE: You may notice that the significance codes (e.g., ***) have a strange output in HTML files, but they will render more cleanly in your RStudio window.

m1 <- lm(CESD2 ~ CESD1, data=ex8_data)
summary(m1)

Call:
lm(formula = CESD2 ~ CESD1, data = ex8_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-21.931  -7.460  -1.998   6.900  40.830 

Coefficients:
            Estimate Std. Error t value        Pr(>|t|)    
(Intercept)  2.56758    1.64667   1.559           0.121    
CESD1        0.45031    0.05929   7.595 0.0000000000014 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.27 on 188 degrees of freedom
Multiple R-squared:  0.2348,    Adjusted R-squared:  0.2307 
F-statistic: 57.69 on 1 and 188 DF,  p-value: 0.0000000000014

Plug in estimates to specify equation for line of best fit. Copy and paste the values we obtained above. The estimate of (Intercept) corresponds to a, whereas the estimate of CESD1 corresponds to b. \[ \hat{Y} = 2.57 + 0.45x \]

Examine confidence interval of the estimates using the confint(model object) function.

confint(m1)

                 2.5 %    97.5 %
(Intercept) -0.6807409 5.8158967
CESD1        0.3333569 0.5672637

Manually add linear regression fit onto scatterplot. Start by retrieving the estimate of the intercept (i.e., the first coefficient in the model) using the coef(model object) function and specifying the first element [1]. Next retrieve the estimate of the slope using the same function but instead indexing the second element [2]. Now we add a layer to create a manual line where we specify our slope and intercept using geom_abline(intercept = intercept estimate, slope = slope estimate). Let’s also color it blue so we can distinguish it from the black data points.

m1_intercept <- coef(m1)[1]
m1_slope <- coef(m1)[2]
ggplot(ex8_data, aes(CESD1, CESD2)) + 
  geom_point() + 
  geom_abline(intercept = m1_intercept, slope = m1_slope, color="blue")

Automatically add linear regression fit (w/ 95% CI) onto scatterplot using geom_smooth() but this time we specify that the method argument will be a linear model (lm), rather than a LOESS line (default).

ggplot(ex8_data, aes(CESD1, CESD2)) + 
  geom_point() + 
  geom_smooth(method=lm)

Remember to add title and axes to your plots using xlab(), ylab(), and ggtitle().

ggplot(ex8_data, aes(CESD1, CESD2)) + 
  geom_point() + 
  geom_smooth(method=lm) + 
  xlab("CESD at Time 1") + 
  ylab("CESD at Time 2") + 
  ggtitle("CESD at Time 2 as a function of CESD at Time 1")

BONUS: Adding Equation

We can also superimpose the regression line using the stat_poly_eq() function.

ggplot(ex8_data, aes(CESD1, CESD2)) + 
  geom_point() + 
  geom_smooth(method=lm) +
  xlab("CESD at Time 1") + 
  ylab("CESD at Time 2") + 
  ggtitle("CESD at Time 2 as a function of CESD at Time 1") +
  stat_poly_eq(aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")),
               formula = y~x, parse = TRUE, label.x.npc = .5, color="blue")

BONUS: Marginal Histograms

We can also add histograms to convey marginal distributions

p1 <- ggplot(ex8_data, aes(CESD1, CESD2)) + 
  geom_point() + 
  geom_smooth(method=lm) +
  xlab("CESD at Time 1") + 
  ylab("CESD at Time 2") + 
  stat_poly_eq(aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")),
               formula = y~x, parse = TRUE, label.x.npc = .5, color="blue") + 
  
  # Remove gridlines 
  theme_bw(base_size = 16) +
  theme(axis.text.x     = element_text(size = 12),
        axis.title.y    = element_text(vjust = +1.5),
        panel.grid.major  = element_blank(),
        panel.grid.minor  = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.line  = element_line(colour = "black"))
ggMarginal(
  p = p1,
  type = 'histogram',
  bins = 20,
  margins = 'both',
  size = 2,
  col = 'black',
  fill="grey"
)