Correlations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Correlations\n",
    "\n",
    "We're now going to look at the relationships between two __numerical__ variables.  This will allow us to make this final leap from ANOVA to linear regression.  \n",
    "\n",
    "Correlations are a numerical representation of the strength of the relationship between two numerical variables - and in a way reflects the ability to predict the second variable given the value of the first variable.  \n",
    "\n",
    "<img src=\"images/horse.png\" width=\"300\" height=\"400\">\n",
    "\n",
    "\n",
    "\n",
    "The bivariate (relationship between __two__ variables) correlation tells us:\n",
    "- If the association exists\n",
    "- The strength of the association\n",
    "- The direction of the association\n",
    "\n",
    "Correlations specifically tell us about *__linear__* relationships between two variables.  They are represented by the lower-case letter $r$, and range from -1 to 1, where 0 is no correlation (or association) between the two variables, -1 as the strongest possible _negative_ correlation and 1 as the strongest possible _positive_ correlation.\n",
    "\n",
    "As you might have guessed, these $r$ values are related to our r-squared ($r^2$) values we've looked at previously.  $r$ is the \"coefficient of correlation,\" and when we square that value we get $r^2$ the \"coefficient of determination.\"  Keep in mind that while these values are related, the interpretations are different.\n",
    "\n",
    "So, what are linear relationships?  What is a positive vs. a negative correlation?  Let's look at some scatterplots!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## loading some libraries!\n",
    "library(tidyverse) ## all of our normal functions for working with data\n",
    "library(ggcorrplot) ## make pretty corrplot\n",
    "library(GGally) ## scatterplot matrix function\n",
    "library(gtrendsR) ## Google Trends API\n",
    "\n",
    "options(repr.plot.width=4, repr.plot.height=3) ## set options for plot size within the notebook -\n",
    "# this is only for jupyter notebooks, you can disregard this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## loading up our old friend the mpg dataset.  \n",
    "## remember this is a dataset that is \"built-in\" to R\n",
    "## this is not how you load data from elsewhere\n",
    "data(mpg) #load built-in dataset\n",
    "head(mpg) #peek at the first 6 observations\n",
    "summary(mpg) #look at summary of variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "?mpg ## obtain documentation for r functions or built-in datasets using ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we're going to plot hwy mpg vs. cty mpg.  What type of relationship would we expect between these variables?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg  %>% ggplot (aes(x = hwy, y = cty)) ## the variables we want to plot\n",
    "    + geom_point() ## the type of plot we want"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This graph shows an example of positive correlation - see how all of the dots \"line up\" and show a general trend that as the x variable increases, the y variable increases.  Because they are both increasing together, the correlation is positive.\n",
    "\n",
    "Remember the correlation value is an estimation of the _linear_ relationship between the two variables - and this is a clump of dots.  We can fit a \"best fit\" line to this graph to show that estimated linear relationship."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg %>% ggplot(aes(x=hwy, y=cty)) + ## the variables we want to plot\n",
    "  geom_point()+  ## the type of plot we want\n",
    "  geom_smooth(method=lm, se=FALSE) ## adding that \"best fit\" line using a linear model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is our \"best fit\" line.  These linear relationships are going to form the basis of Linear Regression that we will get into next week, but right now we're going to focus simply on the correlation between these two variables.  Note that the distance of the dots from the line affects the strength of the relationship - the closer to the line that the dots are clumped, the stronger the relationship.  \n",
    "\n",
    "Let's look at the relationship between hwy and another numerical variable - displ - instead.  displ is the engine displacement, in litres."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg %>% ggplot(aes(x=hwy, y=displ)) + \n",
    "  geom_point()+\n",
    "  geom_smooth(method=lm, se=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now what we see is an example of a negative correlation, because the line starts at the top left corner and goes down toward the bottom right.  So as hwy increases, displ decreases.  These dots seem to be more spread out from the line, so I would guess that the strength of the relationship between hwy and displ (the correlation between the variables) is lower than that between hwy and cty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg %>% ggplot(aes(x=hwy, y=cyl)) + \n",
    "  geom_point()+\n",
    "  geom_smooth(method=lm, se=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You need to make sure your variables are truely numeric, and not simply ordinal.  Because cylinder can only be 4, 5, 6, or 8, and no values in between those - we see that all of the values line up at those values of y.  The line indicates that there may be a general trend of negative association between these variables, but since cyl is ordinal it would not be best suited for a correlation analysis (however we can do an ANOVA using cyl as a categorical grouping variable, and hwy as the numerical DV).  There are non-parametric versions of the correlation coefficient that could also be used, but we will not cover in this course.\n",
    "\n",
    "There are various statistical tests that are most appropriate for certain types of data.  Once you have \"stocked\" your statistical \"toolbelt\" you can do an inference test on any combination of data and variable types, by using the one that is most appropriate for your data types and your RQ.\n",
    "\n",
    "We can also quickly inspect the correlations among all of the numerical variables in a dataset.  This will be more important when we get into modeling and visualize our data as a preparation for creating our linear models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr <- round(cor(mtcars), 1)  ## obtain all of the correlations within pairs of all the num vars\n",
    "ggcorrplot(corr) ## plot those correlations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This (above) is a heatmap of the strength and direction of the correlations between these variables where very blue is -1 and very read is +1.\n",
    "\n",
    "We can also create a scatterplot matrix where a number of variables are compared in scatterplots (like our examples above) in a grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=8, repr.plot.height=4) ## plot size options for Jupyter notebook ONLY\n",
    "pairs(iris[,1:4]) ## plot the correlations between pairs of variables in columns 1 through 4 of iris dataset\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use a categorical variable to color the dots in the grid by groups, to see how the relationship between the two numerical variables might be associated with an additional categorical variable....... more about this to come soon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=8, repr.plot.height=4)\n",
    "pairs(iris[,1:4], col=iris$Species) ## col= adds color based on the grouping variable specified"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also have situations when there is no correlation between the two variables: <BR>\n",
    "<img src=\"images/r2model.png\" width=\"500\" >\n",
    "    \n",
    "To summarize:\n",
    "<br>\n",
    "<img src=\"images/strength.PNG\" width=\"1000\" >\n",
    "\n",
    "There are also cases where variables have obvious associations, but they are not linear.  They have a __*correlation*__ of 0, but they are __*associated*__.  We have to be careful how we use the word correlation in writing our statistical results.\n",
    "\n",
    "<img src=\"images/curvy.PNG\" width=\"600\" >"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assumptions\n",
    "The assumptions are pretty basic, you need two variables, both numeric.  \n",
    "\n",
    "_IF_ you want to do significance testing, they would need to be normally distributed.\n",
    "\n",
    "## Significance Testing?\n",
    "Correlation coefficients by themselves are interpretable as the size of the relationship between two variables. However, there is also a significance test we can conduct on a correlation which will generate a t-score we can compare to the t-distribution to obtain a p-value.\n",
    "\n",
    "The hypotheses for this type of test is pretty basic - is the correlation coefficient significantly different from 0?  This tells us nothing about the relative strength, only if a significant effect exists (or not).\n",
    "\n",
    "#### Non-directional (two-tailed):\n",
    "$H_0: r = 0$ <BR>\n",
    "$H_A: r \\neq 0$ <BR>\n",
    "\n",
    "#### Directional (one-tailed):\n",
    "$H_0: r = 0$ <BR>\n",
    "$H_A: r > 0$ <BR>\n",
    "    __OR__ <BR>\n",
    " $H_A: r < 0$ <BR>   \n",
    "    \n",
    "## Reminder: Correlation does not equal Causation\n",
    "<img src=\"images/venti.jpg\" width=\"300\" >\n",
    "\n",
    "Correlation _could_ be evidence of potential causality, but:\n",
    "- there could be a third variable that is actually causing the effect (ice cream sales -> rise in crime)\n",
    "\n",
    "- we don't know which direction the effect occurs - does X predict Y or does Y predict X?\n",
    "\n",
    "## Calculating Correlation:\n",
    "\n",
    "## $r = \\frac{cov_{xy}}{s_xs_y} = \\frac{\\sum{(x-\\bar{x})(y-\\bar{y})}}{\\sqrt{\\sum{(x-\\bar{x})^2}\\sum{(y-\\bar{y})^2}}}$\n",
    "\n",
    "\n",
    "Let's look at some examples now. For fun, we're going to connect to the Google Trends API using the R package `gtrendsR` and get data about certain keyword searches over a period of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=10, repr.plot.height=4) ## set options for plot size within the notebook -\n",
    "# this is only for jupyter notebooks, you can disregard this.\n",
    "\n",
    "## call the google trends api and return hits info for keyword searches\n",
    "## hits are scaled to values between 0 - 100\n",
    "trends <- gtrends(c(\"statistics\", \"pugs\"), geo = \"US-MD\", time = \"2016-01-01 2019-09-30\", low_search_volume = T)\n",
    "plot(trends)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep in mind that to be correlated the variables __*do not*__ have to have the same magnitude - they just have to trend together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## extract the interest_over_time df from the gtrends object\n",
    "\n",
    "trend_time <- as_tibble(trends$interest_over_time)\n",
    "glimpse(trend_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## look at basic summary statistics\n",
    "\n",
    "trend_time %>%\n",
    "  group_by(keyword) %>%\n",
    "  summarise(mean(hits), median(hits), var(hits))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use this data we need to pivot it so that our \"long\" format is in \"wide\" format - so that we have hits for statistics and hits for pugs in their own columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trend_wide <- \n",
    "  trend_time %>%\n",
    "  spread(key = keyword, value = hits)\n",
    "glimpse(trend_wide)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at a scatterplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trend_wide  %>% \n",
    "    ggplot (aes(x = statistics, y = pugs))  + \n",
    "    geom_point() +\n",
    "    geom_smooth(method=lm, se=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, let's calculate the correlation between these two variables.  For this we will use the function (base R) cor().  The arguments for cor are simple just specify x and y (your two variables to compare)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cor(trend_wide$statistics, trend_wide$pugs) \n",
    "# correlation between hits for statistics and hits for pugs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This correlation is not as low as we may have expected.  It is a \"weak\" correlation, but is it significantly different from a zero correlation?  For that we need to use cor.test() which performs the hypothesis test as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cor.test(trend_wide$statistics, trend_wide$pugs) ## correlation with CI and t-test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So our p-value is below 0.05, so the correlation is statistically significant from zero, and therefore there is at least some correlation.  However, the confidence interval for the correlation is 0.09 to 0.36, so the true population correlation may be anywhere between 0.09 (very very low) or 0.36 (small/mediumish).\n",
    "\n",
    "__NOTE: the values in the output are subject to change as Google Trends samples from the overall data and only returns a small sample of their massive dataset. This is a real-life example of all of that sampling variance we've been talking about.__\n",
    "\n",
    "\n",
    "### Reporting a correlation\n",
    "How would we report this formally?\n",
    "\n",
    "Keyword searches in Google for \"statistics\" has a small correlation with keyword searches for \"pugs\" ($r = 0.24$) over the time period from 2016 to today within the state of Maryland.  While the correlation is significant (p < 0.001), the magnitude of the correlation is small."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R-squared\n",
    "Remember, our value of the proportion of variance explained is literally this correlation coefficient, $r$, squared - $r^2$.  Let's take a look:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cor(trend_wide$statistics, trend_wide$pugs)^2 ## calculate the correlation, and square it - ^2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpretation of this is that this % of the variance in the hits for statistics is explained in the hits for pugs.  So there are not many people out there like me that search for statistics and pugs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fun stuff related to correlations:\n",
    "- <a href=\"http://guessthecorrelation.com/\">Guess the Correlation game</a>\n",
    "- <a href=\"https://www.tylervigen.com/spurious-correlations\"> Spurious Correlations </a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}