Correlation, Scatterplots, and Regression Analysis
Daniel Lefevre (original by: Rachel Saidi)
May 25, 2020
Create a Scatterplot
In this example, look at US crime rates at the state level, in 2005, with rates per 100,000 population for crime types such as murder, robbery, and aggravated assault, as reported by the Census Bureau. There are 7 crime types in total. The dataset is clean to begin with.
library(readr)
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutcrime <- read_csv('http://datasets.flowingdata.com/crimeRatesByState2005.csv')##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## state = col_character(),
## murder = col_double(),
## forcible_rape = col_double(),
## robbery = col_double(),
## aggravated_assault = col_double(),
## burglary = col_double(),
## larceny_theft = col_double(),
## motor_vehicle_theft = col_double(),
## population = col_double()
## )# source: U.S. Census Bureau and Nathan YauCheck out the first few lines
## # A tibble: 6 x 9
## state murder forcible_rape robbery aggravated_assa… burglary larceny_theft
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 United S… 5.6 31.7 141. 291. 727. 2286.
## 2 Alabama 8.2 34.3 141. 248. 954. 2650
## 3 Alaska 4.8 81.1 80.9 465. 622. 2599.
## 4 Arizona 7.5 33.8 144. 327. 948. 2965.
## 5 Arkansas 6.7 42.9 91.1 387. 1085. 2711.
## 6 Californ… 6.9 26 176. 317. 693. 1916.
## # … with 2 more variables: motor_vehicle_theft <dbl>, population <dbl>Notice
The data has a column for the state and then the rest are rates for various crimes. Now make a quick scatterplot.
Map variables in the data onto the X and Y axes and change the axes labels and theme
The default gray theme of ggplot2 has a rather academic look. See here and here for how to use the theme option to customize individual elements of a chart. Use one of the ggplot2 built-in themes, and then customize the fonts.
# Change the theme
ggplot(crime, aes(x = burglary, y = murder)) +
xlab("Burglary rates in each state per 100,000") +
ylab("Murder rates in each state per 100,000") +
theme_classic(base_size = 12)
p1 <- ggplot(crime, aes(x = burglary, y = murder)) +
labs(title = "MURDERS VERSUS BURGLARIES IN US STATES PER 100,000",
caption = "Source: U.S. Census Bureau and Nathan Yau") +
xlab("Burglary rates in each state per 100,000") +
ylab ("Murder rates in each state per 100,000") +
theme_minimal(base_size = 12)
p1 + geom_point()
What is going on with the outlier?
The one point far higher than the rest represents Washington, D.C., which had a much higher murder rate of 35.4. The states with the next highest murder rate at that time were Louisiana and Maryland at 9.9 per 100,000.
Remove D.C. and US averages and replot:
crime2 <- crime[crime$state != "District of Columbia",]
crime2 <- crime2[crime2$state != "United States",]
p2 <- ggplot(crime2, aes(x = burglary, y = murder)) +
labs(title = "MURDERS VERSUS BURGLARIES IN US STATES PER 100,000",
caption = "Source: U.S. Census Bureau and Nathan Yau") +
xlab("Burglary rates in each state per 100,000") +
ylab ("Murder rates in each state per 100,000") +
theme_minimal(base_size = 12)
p2 + geom_point()
Now the scatterplot appears to show a correlation
Fix the axes to start at 0.
p3 <- p2 + xlim(250,1200)+ ylim(0,10)
p3 + geom_point()## Warning: Removed 1 rows containing missing values (geom_point).
Add a smoother in red with a confidence interval
p4 <- p3 + geom_point() + geom_smooth(color = "red")
p4## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
Add a linear regression with confidence interval
p5 <- p3 + geom_point() + geom_smooth(method='lm',formula=y~x)
p5## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
Add a title, make the line dashed, and remove the confidence interval band
The command se = FALSE takes away the CI band
p6 <- p3 + geom_point() + geom_smooth(method='lm',formula=y~x, se = FALSE, linetype= "dotdash", size = 0.3) +
ggtitle("BURGLARIES VERSUS MURDERS IN THE U.S.")
p6 ## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
What is the linear equation of that linear regression model?
In the form, y=mx + b, we use the command, lm(y~x), meaning, fit the predictor variable x into the model to predict y. Look at the values of (Intercept) and murder. The column, Estimate gives the value you need in your linear model. The column for Pr(>|t|) describes whether the predictor is useful to the model. The more asterisks, the more the variable contributes to the model.
cor(crime2$burglary, crime2$murder)## [1] 0.6231757fit1 <- lm(murder ~ burglary, data = crime2)
summary(fit1)##
## Call:
## lm(formula = murder ~ burglary, data = crime2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2924 -1.2156 -0.2142 1.1749 5.4978
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.395519 0.825748 0.479 0.634
## burglary 0.006247 0.001132 5.521 1.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.87 on 48 degrees of freedom
## Multiple R-squared: 0.3883, Adjusted R-squared: 0.3756
## F-statistic: 30.48 on 1 and 48 DF, p-value: 1.342e-06What does the output mean?
Cor stands for “correlation”. This is a value between (inclusively) -1 and 1. The correlation coefficient tells how strong or weak the correlation is. Values closer to +/- 1 are strong correlation (the sign is determined by the linear slope), values close to +/- 0.5 are weak correlation, and values close to zero have no correlation.
The model has the equation: murder = 0.0062(burglary) + 0.396
The slope may be interpreted in the following: For each additional burglary per 100,000, there is a predicted increase of 0.006 murders.
The p-value on the right of burglary has 3 asterisks which suggests it is a meaningful variable to explain the linear increase in murders. But we also need to look at the Adjusted R-Squared value. It states that about 38% of the variation in the observations may be explained by the model. In other words, 62% of the variation in the data is likely not explained by this model.
What about more variables?
Can a model with more predictors also be used? What would we be trying to predict?
Is there an easier way to compare multiple variables using a scatterplot matrix?
Check out the pairwise comparisons with density curves and correlation output
library(GGally)
ggpairs(crime2, columns = 2:8)
Another method: Use a correlation plot to explore the correlation among all variables
This correlation plot shows similar pairwise results as above, but in a heatmap of correlation values.
#install.packages("DataExplorer")
library(DataExplorer)
plot_correlation(crime)## Warning in dummify(data, maxcat = maxcat): Ignored all discrete features since
## `maxcat` set to 20 categories!
Collinearity
The key goal of multiple regression analysis is to isolate the relationship between EACH INDEPENDENT VARIABLE and the DEPENDENT VARIABLE. COLLINEARITY means variables are correlated and thus NOT INDEPENDENT. The more correlated the variables, the more difficult it is to change one variable without changing the other. This is important to keep in mind. The two different matrices gave slightly different correlation information. We are concerned with dependence of 2 or more variables.
Using library(GGally), the two variables with the highest correlation of 0.75 are robbery and murder. Using library(DataExplorer), murder and robbery are highly correlated at 0.92. Although different correlation values, they corroborate each other.
Now try to make a multiple regression model.
With multiple regression, there are several strategies for comparing variable inputs into a model. I will show you backward elimination. In backward elimination, start with all possible predictor variables with your response variable. In this case, we will use: burglary forcible_rape aggravated_assault larceny_theft motor_vehicle_theft Perform a model fit with all predictors.
Look at the p-value for each variable - if it is relatively small ( < 0.10), then it is likely contributing to the model.
Check out the residual plots. A good model will have a relatively straight horizontal red line across the scatterplot between residuals plotted with fitted values (see below for a good residuals plot). You can also look at the other plots (Normal QQ, Scale-Location, and Residuals vs Leverage), but for now we will focus on the residual vs. fitted plot. The more curved the red line, the more likely that a better model exists.
Look at the output for the Adjusted R-Squared value at the bottom of the output. The interpretation is:
__% (from the adjusted r-squared value) of the variation in the observations may be explained by this model. The higher the adjusted R-squared value, the better the model. We use the adjusted R-squared value because it compensates for more predictors mathematically increasing the normal R-squared value.
fit2 <- lm(murder ~ robbery + burglary + forcible_rape + aggravated_assault + larceny_theft + motor_vehicle_theft, data = crime2)
summary(fit2)##
## Call:
## lm(formula = murder ~ robbery + burglary + forcible_rape + aggravated_assault +
## larceny_theft + motor_vehicle_theft, data = crime2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7088 -0.7961 -0.0508 0.6757 3.4723
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9940985 1.0835208 0.917 0.364014
## robbery 0.0194331 0.0052193 3.723 0.000567 ***
## burglary 0.0041431 0.0013339 3.106 0.003352 **
## forcible_rape -0.0126884 0.0210395 -0.603 0.549627
## aggravated_assault 0.0045161 0.0023433 1.927 0.060576 .
## larceny_theft -0.0007841 0.0005622 -1.395 0.170246
## motor_vehicle_theft -0.0002426 0.0012751 -0.190 0.849982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.338 on 43 degrees of freedom
## Multiple R-squared: 0.7193, Adjusted R-squared: 0.6801
## F-statistic: 18.36 on 6 and 43 DF, p-value: 1.949e-10plot(fit2)
What are we really trying to predict?
If we are trying to predict murder rates, then we can see if any of the predictor variables contribute to this model. Note the adjusted R-squared value is 58.66% The only variable that does not appear to be as significant as the others is motor_vehicle_theft. So drop that and re-run the model.
Don’t forget, we found that murder and robbery are highly correlated, so it makes sense that robbery fits in the model.
fit3 <- lm(murder ~ robbery + burglary + forcible_rape + aggravated_assault + larceny_theft, data = crime2)
summary(fit3)##
## Call:
## lm(formula = murder ~ robbery + burglary + forcible_rape + aggravated_assault +
## larceny_theft, data = crime2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6923 -0.7545 -0.0751 0.6404 3.4836
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0611101 1.0134089 1.047 0.300785
## robbery 0.0189486 0.0045060 4.205 0.000126 ***
## burglary 0.0041189 0.0013131 3.137 0.003044 **
## forcible_rape -0.0134321 0.0204456 -0.657 0.514623
## aggravated_assault 0.0046152 0.0022596 2.042 0.047124 *
## larceny_theft -0.0008229 0.0005181 -1.588 0.119349
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.324 on 44 degrees of freedom
## Multiple R-squared: 0.719, Adjusted R-squared: 0.6871
## F-statistic: 22.52 on 5 and 44 DF, p-value: 3.917e-11plot(fit3)
The residuals plot looks worse once we drop motor_vehicle_theft, AND the adjusted R-squared value dropped to 68.7%. Maybe one last exploration is to change what we are predicting - select the simplest model (parsimonious) with burglary and aggravated assault on murders.
fit4 <- lm(murder ~ robbery + burglary + aggravated_assault , data = crime2)
summary(fit4)##
## Call:
## lm(formula = murder ~ robbery + burglary + aggravated_assault,
## data = crime2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6434 -0.7535 -0.0107 0.7229 3.7420
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.330470 0.601764 -0.549 0.5855
## robbery 0.021669 0.004075 5.318 3e-06 ***
## burglary 0.002732 0.001042 2.621 0.0118 *
## aggravated_assault 0.003570 0.002030 1.759 0.0853 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.338 on 46 degrees of freedom
## Multiple R-squared: 0.6998, Adjusted R-squared: 0.6802
## F-statistic: 35.74 on 3 and 46 DF, p-value: 4.451e-12plot(fit4)
Interesting!!
The adjusted R-squared went down to 68.02%. The residuals plot looks about the same.
The residuals plot shows observations 20, 18, and 46 have an effect on the residuals plot as well as all three have high leverage with Cooks distance.
Maryland is 20 Louisiana is 18 Virginia is 46
Back to simply murders and burglaries - bring in the state’s population as a size of the circle
p2 +
geom_point(aes(size = population), color = "red") + xlim(250,1200) + ylim(0,10) +
labs(title = "MURDERS VERSUS BURGLARIES IN US STATES PER 100,000",
caption = "Source: U.S. Census Bureau and Nathan Yau") +
xlab("Burglary rates in each state per 100,000") +
ylab ("Murder rates in each state per 100,000") +
theme_minimal(base_size = 12)## Warning: Removed 1 rows containing missing values (geom_point).
Finally, add some interactivity to the plot with plotly
p <- ggplot(crime2, aes(x = burglary, y = murder, size = population, text = paste("state:", state))) +
geom_point(alpha = 0.5, color = "red") + xlim(250,1200) + ylim(0,10) +
ggtitle("BURGLARIES VERSUS MURDERS IN THE U.S.", subtitle = "Sizes of circles are proportional to state populations") +
xlab("Burglary rates in each state per 100,000") +
ylab ("Murder rates in each state per 100,000") +
theme_minimal(base_size = 12)
p <- ggplotly(p)
pMake a series of charts from food stamps data
Now we will explore a series of other geom functions using the food stamps data.
Load the data, map variables onto the X and Y axes, and save chart template
# load data
food_stamps <- read_csv("food_stamps.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## year = col_double(),
## participants = col_double(),
## costs = col_double()
## )# save basic chart template
food_stamps_chart <- ggplot(food_stamps, aes(x = year, y = participants)) +
labs(title = "Food Stamps Participation Over the Years") +
xlab("Year") +
ylab("Participants (millions)") +
theme_minimal(base_size = 14)
food_stamps_chart
Make a line chart
food_stamps_chart +
geom_line()
Customize the line, and add a title
food_stamps_chart +
geom_line(size = 1.5, color = "red") +
ggtitle("Line chart")
Add a second layer to make a dot-and-line chart
food_stamps_chart +
geom_line() +
geom_point() +
ggtitle("Dot-and-line chart")
Make a column chart, then flip its coordinates to make a bar chart
# Make a column chart
food_stamps_chart +
geom_bar(stat = "identity") +
ggtitle("Column chart") +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
geom_bar works a little differently to the geoms we have considered previously. If you have not mapped data values to the Y axis with aes, its default behavior is to set the heights of the bars by counting the number of records for values along the X axis. If you have mapped a variable to the Y axis, and want the heights of the bars to represent values in the data, use you must use stat=“identity”.
coord_flip switches the X and Y axes.
# Make a bar chart
food_stamps_chart +
geom_bar(stat = "identity") +
ggtitle("Bar chart") +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()) +
coord_flip()
The difference between color and fill
For some geoms, notably geom_bar, you can set color for their outline as well as the interior of the shape.
When setting colors, color refers to the outline, fill to the interior of the shape.
# set color and fill
food_stamps_chart +
geom_bar(stat = "identity", color = "#888888", fill = "#CCCCCC", alpha = 0.5) +
ggtitle("Column chart")
Map color to the values of a continuous variable
# fill the bars according to values for the cost of the program
food_stamps_chart +
geom_bar(stat = "identity", color= "white", aes(fill = costs))
This code uses an aes mapping to color the bars according values for the costs of the program, in billions of dollars. ggplot2 recognizes that costs is a continuous variable, but its default sequential scheme applies more intense blues to lower values, which is counterintuitive.
Use a ColorBrewer sequential color palette
# use a colorbrewer gradient levels for intensity
food_stamps_chart +
geom_bar(stat = "identity", color = "#888888", aes(fill = costs)) +
scale_fill_gradient(name = "Cost\n($ billion", low = "#d1dee8", high = "#d92774")
scale_fill_distiller (and scale_color_distiller) work like scale_color_brewer, but set color gradients for ColorBrewer’s sequential and diverging color palettes; direction = 1 ensures that larger numbers are mapped to more intense colors (direction = -1 reverses the color mapping). Try changing the code I have: scale_fill_gradient() to scale_fill_distiller with different directions (1 or -1).
Notice also the in the title for the legend. This introduces a new line.
Control the position of the legend
This code uses the theme function to moves the legend from its default position to the right of the chart to use some empty space on the chart itself.
food_stamps_chart +
geom_bar(stat="identity", color = "#888888", aes(fill=costs)) +
scale_fill_gradient(name = "Cost\n($ billion", low = "#d1dee8", high = "#d92774") +
theme(legend.position=c(0.15,0.8))
The coordinates for the legend are given as a list: The first number sets the horizontal position, from left to right, on a scale from 0 to 1; the second number sets the vertical position, from bottom to top, again on a scale from 0 to 1.