Correlation and Regression
Manuela da Cruz Chadreque
1 de Abril de 2019
#install.packages("openintro")
library(openintro)## Warning: package 'openintro' was built under R version 3.4.4## Please visit openintro.org for free statistics materials##
## Attaching package: 'openintro'## The following objects are masked from 'package:datasets':
##
## cars, trees#install.packages("ggplot2")
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.4.4##
## Attaching package: 'ggplot2'## The following object is masked from 'package:openintro':
##
## diamonds#install.packages("knitr")
library(knitr)## Warning: package 'knitr' was built under R version 3.4.4#install.packages("dplyr")
library(dplyr)## Warning: package 'dplyr' was built under R version 3.4.4##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union#install.packages("plyr")
library(plyr)## Warning: package 'plyr' was built under R version 3.4.4## -------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## -------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarizeVisualizing two variables
Scatterplots
Scatterplots are the most common and effective tools for visualizing the relationship between two numeric variables.
The ncbirths dataset is a random sample of 1,000 cases taken from a larger dataset collected in 2004. Each case describes the birth of a single child born in North Carolina, along with various characteristics of the child (e.g. birth weight, length of gestation, etc.), the child’s mother (e.g. age, weight gained during pregnancy, smoking habits, etc.) and the child’s father (e.g. age). You can view the help file for these data by running ?ncbirths in the console.
# Scatterplot of weight vs. weeks
ggplot(data=ncbirths, aes(x=weeks, y=weight))+
geom_point()## Warning: Removed 2 rows containing missing values (geom_point).
Boxplots as discretized/conditioned scatterplots
If it is helpful, you can think of boxplots as scatterplots for which the variable on the x-axis has been discretized.
The cut() function takes two arguments: the continuous variable you want to discretize and the number of breaks that you want to make in that continuous variable in order to discretize it.
# Boxplot of weight vs. weeks
ggplot(data = ncbirths,
aes(x = cut(weeks, breaks =5), y = weight)) +
geom_boxplot()
###Creating scatterplots Creating scatterplots is simple and they are so useful that is it worthwhile to expose yourself to many examples. Over time, you will gain familiarity with the types of patterns that you see. You will begin to recognize how scatterplots can reveal the nature of the relationship between two variables.
The mammals dataset contains information about 39 different species of mammals, including their body weight, brain weight, gestation time, and a few other variables.
## Mammals scatterplot
ggplot(data=mammals, aes(x=BodyWt, y=BrainWt))+
geom_point()
The mlbBat10 dataset contains batting statistics for 1,199 Major League Baseball players during the 2010 season.
# Baseball player scatterplot
ggplot(data=mlbBat10, aes(x=OBP, y=SLG))+
geom_point()
The bdims dataset contains body girth and skeletal diameter measurements for 507 physically active individuals.
# Body dimensions scatterplot
ggplot(data=bdims, aes(x=hgt, y=wgt, colour=factor(sex)))+
geom_point()
The smoking dataset contains information on the smoking habits of 1,691 citizens of the United Kingdom.
# Smoking scatterplot
ggplot(data=smoking, aes(x=age, y=amtWeekdays))+
geom_point()## Warning: Removed 1270 rows containing missing values (geom_point).
Transformations
The relationship between two variables may not be linear. In these cases we can sometimes see strange and even inscrutable patterns in a scatterplot of the data. Sometimes there really is no meaningful relationship between the two variables. Other times, a careful transformation of one or both of the variables can reveal a clear relationship.
Recall the bizarre pattern that you saw in the scatterplot between brain weight and body weight among mammals in a previous exercise. Can we use transformations to clarify this relationship?
ggplot2 provides several different mechanisms for viewing transformed relationships. The coord_trans() function transforms the coordinates of the plot. Alternatively, the scale_x_log10() and scale_y_log10() functions perform a base-10 log transformation of each axis. Note the differences in the appearance of the axes.
# Scatterplot with coord_trans()
ggplot(data = mammals, aes(x = BodyWt, y = BrainWt)) +
geom_point() +
coord_trans(x = "log10", y = "log10")
# Scatterplot with scale_x_log10() and scale_y_log10()
ggplot(data = mammals, aes(x =BodyWt, y = BrainWt)) +
geom_point() +
scale_x_log10() + scale_y_log10()
Outliers
Identifying outliers
Identifying outliers In Chapter 5, we will discuss how outliers can affect the results of a linear regression model and how we can deal with them. For now, it is enough to simply identify them and note how the relationship between two variables may change as a result of removing outliers.
Recall that in the baseball example earlier in the chapter, most of the points were clustered in the lower left corner of the plot, making it difficult to see the general pattern of the majority of the data. This difficulty was caused by a few outlying players whose on-base percentages (OBPs) were exceptionally high. These values are present in our dataset only because these players had very few batting opportunities.
Both OBP and SLG are known as rate statistics, since they measure the frequency of certain events (as opposed to their count). In order to compare these rates sensibly, it makes sense to include only players with a reasonable number of opportunities, so that these observed rates have the chance to approach their long-run frequencies.
In Major League Baseball, batters qualify for the batting title only if they have 3.1 plate appearances per game. This translates into roughly 502 plate appearances in a 162-game season. The mlbBat10 dataset does not include plate appearances as a variable, but we can use at-bats (AB) – which constitute a subset of plate appearances – as a proxy.
# Scatterplot of SLG vs. OBP
mlbBat10 %>%
filter(AB>=200) %>%
ggplot(aes(x = OBP, y = SLG)) +
geom_point()
# Identify the outlying player
mlbBat10 %>%
filter(AB>=200, OBP<0.2)## name team position G AB R H 2B 3B HR RBI TB BB SO SB CS OBP
## 1 B Wood LAA 3B 81 226 20 33 2 0 4 14 47 6 71 1 0 0.174
## SLG AVG
## 1 0.208 0.146Correlation
Computing correlation
The cor(x, y) function will compute the Pearson product-moment correlation between variables, x and y. Since this quantity is symmetric with respect to x and y, it doesn’t matter in which order you put the variables.
At the same time, the cor() function is very conservative when it encounters missing data (e.g. NAs). The use argument allows you to override the default behavior of returning NA whenever any of the values encountered is NA. Setting the use argument to “pairwise.complete.obs” allows cor() to compute the correlation coefficient for those observations where the values of x and y are both not missing.
# Compute correlation
ncbirths %>%
summarize(r = cor(weight, mage))## r
## 1 0.05506589# Compute correlation for all non-missing pairs
ncbirths %>%
summarize(r = cor(weight, weeks, use = "pairwise.complete.obs"))## r
## 1 0.6701013Perception of correlation (2)
Estimating the value of the correlation coefficient between two quantities from their scatterplot can be tricky. Statisticians have shown that people’s perception of the strength of these relationships can be influenced by design choices like the x and y scales.
Nevertheless, with some practice your perception of correlation will improve. Toggle through the four scatterplots in the plotting window, each of which you’ve seen in a previous exercise. Jot down your best estimate of the value of the correlation coefficient between each pair of variables. Then, compare these values to the actual values you compute in this exercise.
# Correlation for all baseball players
mlbBat10 %>%
summarize(r = cor(OBP,SLG))## r
## 1 0.8145628# Correlation for all players with at least 200 ABs
mlbBat10 %>%
filter(AB>=200) %>%
summarize( r = cor(OBP, SLG))## r
## 1 0.6855364# Correlation of body dimensions
bdims %>%
group_by(sex) %>%
summarize(r = cor(hgt, wgt))## r
## 1 0.7173011# Correlation among mammals, with and without log
mammals %>%
summarize(r = cor(BodyWt, BrainWt),
r_log = cor(log(BodyWt), log(BrainWt)))## r r_log
## 1 0.9341638 0.9595748
Spurious correlation in random data Statisticians must always be skeptical of potentially spurious correlations. Human beings are very good at seeing patterns in data, sometimes when the patterns themselves are actually just random noise. To illustrate how easy it can be to fall into this trap, we will look for patterns in truly random data.
The noise dataset contains 20 sets of x and y variables drawn at random from a standard normal distribution. Each set, denoted as z, has 50 observations of x, y pairs. Do you see any pairs of variables that might be meaningfully correlated? Are all of the correlation coefficients close to zero?
noise<-read.csv("noise.csv", sep=";")
#Create faceted scatterplot
ggplot(data = noise, aes(x = x, y = y)) +
geom_point() +
facet_wrap(~ z)
The “best fit” line
The simple linear regression model for a numeric response as a function of a numeric explanatory variable can be visualized on the corresponding scatterplot by a straight line. This is a “best fit” line that cuts through the data in a way that minimizes the distance between the line and the data points.
We might consider linear regression to be a specific example of a larger class of smooth models. The geom_smooth() function allows you to draw such models over a scatterplot of the data itself. This technique is known as visualizing the model in the data space. The method argument to geom_smooth() allows you to specify what class of smooth model you want to see. Since we are exploring linear models, we’ll set this argument to the value “lm”.
Note that geom_smooth() also takes an se argument that controls the standard error, which we will ignore for now.
Uniqueness of least squares regression line
The least squares criterion implies that the slope of the regression line is unique. In practice, the slope is computed by R. In this exercise, you will experiment with trying to find the optimal value for the regression slope for weight as a function of height in the bdims dataset via trial-and-error.
To help, we’ve built a custom function for you called add_line(), which takes a single argument: the proposed slope coefficient.
# Scatterplot with regression line
ggplot(data = bdims, aes(x = hgt, y = wgt)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)
Fitting a linear model “by hand”
Recall the simple linear regression model: \(Y=b_0+b_1*X\)
Two facts enable you to compute the slope b1 and intercept b0 of a simple linear regression model from some basic summary statistics.
First, the slope can be defined as: \(b_1=r_x,y+S_x/S_y\)
bdims_summary<-read.csv("bdims_summary.csv", sep=";")
# Print bdims_summary
print(bdims_summary)## N r mean_hgt sd_hgt mean_wgt sd_wgt
## 1 507 0.7173011 171.1438 9.407205 69.14753 13.34576# Add slope and intercept
bdims_summary %>%
mutate(slope = r*(sd_wgt/sd_hgt),
intercept = mean_wgt-mean_hgt*slope)## N r mean_hgt sd_hgt mean_wgt sd_wgt slope intercept
## 1 507 0.7173011 171.1438 9.407205 69.14753 13.34576 1.017617 -105.0112Regression to the mean Regression to the mean is a concept attributed to Sir Francis Galton. The basic idea is that extreme random observations will tend to be less extreme upon a second trial. This is simply due to chance alone. While “regression to the mean” and “linear regression” are not the same thing, we will examine them together in this exercise.
One way to see the effects of regression to the mean is to compare the heights of parents to their children’s heights. While it is true that tall mothers and fathers tend to have tall children, those children tend to be less tall than their parents, relative to average. That is, fathers who are 3 inches taller than the average father tend to have children who may be taller than average, but by less than 3 inches.
The Galton_men and Galton_women datasets contain data originally collected by Galton himself in the 1880s on the heights of men and women, respectively, along with their parents’ heights.
Compare the slope of the regression line to the slope of the diagonal line. What does this tell you?
Galton_men<-read.csv("Galton_men.csv", sep=";")
Galton_women<-read.csv("Galton_women.csv", sep=";")
# Height of children vs. height of father
ggplot(data = Galton_men, aes(x = father, y = height)) +
geom_point() +
geom_abline(slope = 1, intercept = 0) +
geom_smooth(method = "lm", se = FALSE)
# Height of children vs. height of mother
ggplot(data = Galton_women, aes(x = mother, y = height)) +
geom_point() +
geom_abline(slope = 1, intercept = 0) +
geom_smooth(method = "lm", se = FALSE)
Fitting simple linear models
While the geom_smooth(method = “lm”) function is useful for drawing linear models on a scatterplot, it doesn’t actually return the characteristics of the model. As suggested by that syntax, however, the function that creates linear models is lm(). This function generally takes two arguments:
A formula that specifies the model
A data argument for the data frame that contains the data you want to use to fit the model
The lm() function return a model object having class “lm”. This object contains lots of information about your regression model, including the data used to fit the model, the specification of the model, the fitted values and residuals, etc.
# Linear model for weight as a function of height
lm(wgt ~ hgt, data = bdims)##
## Call:
## lm(formula = wgt ~ hgt, data = bdims)
##
## Coefficients:
## (Intercept) hgt
## -105.011 1.018# Linear model for SLG as a function of OBP
lm(SLG~OBP, data=mlbBat10)##
## Call:
## lm(formula = SLG ~ OBP, data = mlbBat10)
##
## Coefficients:
## (Intercept) OBP
## 0.009407 1.110323# Log-linear model for body weight as a function of brain weight
lm(log(BodyWt)~log(BrainWt), data=mammals)##
## Call:
## lm(formula = log(BodyWt) ~ log(BrainWt), data = mammals)
##
## Coefficients:
## (Intercept) log(BrainWt)
## -2.509 1.225The lm summary output
An “lm” object contains a host of information about the regression model that you fit. There are various ways of extracting different pieces of information.
The coef() function displays only the values of the coefficients. Conversely, the summary() function displays not only that information, but a bunch of other information, including the associated standard error and p-value for each coefficient, the R2, adjusted R2, and the residual standard error. The summary of an “lm” object in R is very similar to the output you would see in other statistical computing environments (e.g. Stata, SPSS, etc.)
mod<-lm(wgt~hgt, data=bdims)
# Show the coefficients
coef(mod)## (Intercept) hgt
## -105.011254 1.017617# Show the full output
summary(mod)##
## Call:
## lm(formula = wgt ~ hgt, data = bdims)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.743 -6.402 -1.231 5.059 41.103
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -105.01125 7.53941 -13.93 <2e-16 ***
## hgt 1.01762 0.04399 23.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.308 on 505 degrees of freedom
## Multiple R-squared: 0.5145, Adjusted R-squared: 0.5136
## F-statistic: 535.2 on 1 and 505 DF, p-value: < 2.2e-16Fitted values and residuals
Once you have fit a regression model, you are often interested in the fitted values \((y^i)\) and the residuals \((ei)\), where \(i\) indexes the observations. Recall that:
\(e_i=y_i ??? \hat{y_i}\)
The least squares fitting procedure guarantees that the mean of the residuals is zero (n.b., numerical instability may result in the computed values not being exactly zero). At the same time, the mean of the fitted values must equal the mean of the response variable.
In this exercise, we will confirm these two mathematical facts by accessing the fitted values and residuals with the \(fitted.values()\) and \(residuals()\) functions, respectively, for the following model:
\(mod <- lm(wgt ~ hgt, data = bdims)\)
# Mean of weights equal to mean of fitted values?
mean(bdims$wgt) == mean(fitted.values(mod))## [1] TRUE# Mean of the residuals
mean(residuals(mod))## [1] -1.266971e-15Tidying your linear model
As you fit a regression model, there are some quantities (e.g. \(R^2\)) that apply to the model as a whole, while others apply to each observation (e.g. \(y_i\)). If there are several of these per-observation quantities, it is sometimes convenient to attach them to the original data as new variables.
The augment() function from the broom package does exactly this. It takes a model object as an argument and returns a data frame that contains the data on which the model was fit, along with several quantities specific to the regression model, including the fitted values, residuals, leverage scores, and standardized residuals.
# Load broom
library(broom)## Warning: package 'broom' was built under R version 3.4.4mod<-lm(formula = wgt ~ hgt, data = bdims)
# Create bdims_tidy
bdims_tidy<-augment(mod)
# Glimpse the resulting data frame
glimpse(bdims_tidy)## Observations: 507
## Variables: 9
## $ wgt <dbl> 65.6, 71.8, 80.7, 72.6, 78.8, 74.8, 86.4, 78.4, 62....
## $ hgt <dbl> 174.0, 175.3, 193.5, 186.5, 187.2, 181.5, 184.0, 18...
## $ .fitted <dbl> 72.05406, 73.37697, 91.89759, 84.77427, 85.48661, 7...
## $ .se.fit <dbl> 0.4320546, 0.4520060, 1.0667332, 0.7919264, 0.81834...
## $ .resid <dbl> -6.4540648, -1.5769666, -11.1975919, -12.1742745, -...
## $ .hat <dbl> 0.002154570, 0.002358152, 0.013133942, 0.007238576,...
## $ .sigma <dbl> 9.312824, 9.317005, 9.303732, 9.301360, 9.312471, 9...
## $ .cooksd <dbl> 5.201807e-04, 3.400330e-05, 9.758463e-03, 6.282074e...
## $ .std.resid <dbl> -0.69413418, -0.16961994, -1.21098084, -1.31269063,...Making predictions
The fitted.values() function or the augment()-ed data frame provides us with the fitted values for the observations that were in the original data. However, once we have fit the model, we may want to compute expected values for observations that were not present in the data on which the model was fit. These types of predictions are called out-of-sample.
The ben data frame contains a height and weight observation for one person. The mod object contains the fitted model for weight as a function of height for the observations in the bdims dataset. We can use the predict() function to generate expected values for the weight of new individuals. We must pass the data frame of new observations through the newdata argument.
ben<-read.csv("ben.csv", sep=";")
# Predict the weight of ben
predict(mod, ben)## 1
## 81.00909Standard error of residuals
One way to assess strength of fit is to consider how far off the model is for a typical case. That is, for some observations, the fitted value will be very close to the actual value, while for others it will not. The magnitude of a typical residual can give us a sense of generally how close our estimates are.
However, recall that some of the residuals are positive, while others are negative. In fact, it is guaranteed by the least squares fitting procedure that the mean of the residuals is zero. Thus, it makes more sense to compute the square root of the mean squared residual, or root mean squared error (RMSE). R calls this quantity the residual standard error.
To make this estimate unbiased, you have to divide the sum of the squared residuals by the degrees of freedom in the model. Thus,
\(RMSE=\sqrt{ \frac{ \sum{ e^2_i} }{d.f}}=\sqrt{ \frac{ SSE }{d.f}}\)
# View summary of model
summary(mod)##
## Call:
## lm(formula = wgt ~ hgt, data = bdims)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.743 -6.402 -1.231 5.059 41.103
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -105.01125 7.53941 -13.93 <2e-16 ***
## hgt 1.01762 0.04399 23.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.308 on 505 degrees of freedom
## Multiple R-squared: 0.5145, Adjusted R-squared: 0.5136
## F-statistic: 535.2 on 1 and 505 DF, p-value: < 2.2e-16# Compute the mean of the residuals
mean(residuals(mod))## [1] -1.266971e-15# Compute RMSE
sqrt(sum(residuals(mod)^2) / df.residual(mod))## [1] 9.30804Assessing simple linear model fit
Recall that the coefficient of determination (\(R^2\)), can be computed as
\(R^2=1-\frac{SSE}{SST}=1-\frac{Var(e)}{Var(y)}\)
where e is the vector of residuals and y is the response variable. This gives us the interpretation of R2 as the percentage of the variability in the response that is explained by the model, since the residuals are the part of that variability that remains unexplained by the model.
bdims_tidy<-augment(mod)
# Compute R-squared
bdims_tidy %>%
summarize(var_y = var(wgt), var_e = var(.resid)) %>%
mutate(R_squared = 1-var_e/var_y)## var_y var_e R_squared
## 1 178.1094 86.46839 0.5145208Linear vs. average
The R2 gives us a numerical measurement of the strength of fit relative to a null model based on the average of the response variable:
\(\hat{y}_{null}=\bar{y}\) This model has an R2 of zero because SSE=SST. That is, since the fitted values (\(y_{null}\)) are all equal to the average (\(\bar{y}\)), the residual for each observation is the distance between that observation and the mean of the response. Since we can always fit the null model, it serves as a baseline against which all other models will be compared.
In the graphic, we visualize the residuals for the null model (mod_null at left) vs. the simple linear regression model (mod_hgt at right) with height as a single explanatory variable. Try to convince yourself that, if you squared the lengths of the grey arrows on the left and summed them up, you would get a larger value than if you performed the same operation on the grey arrows on the right.
Leverage
The leverage of an observation in a regression model is defined entirely in terms of the distance of that observation from the mean of the explanatory variable. That is, observations close to the mean of the explanatory variable have low leverage, while observations far from the mean of the explanatory variable have high leverage. Points of high leverage may or may not be influential.
The augment() function from the broom package will add the leverage scores (.hat) to a model data frame.
# Rank points of high leverage
mod %>%
augment() %>%
arrange(desc(.hat)) %>%
head()## # A tibble: 6 x 9
## wgt hgt .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 85.5 198. 96.6 1.26 -11.1 0.0182 9.30 0.0134 -1.20
## 2 90.9 197. 95.6 1.21 -4.66 0.0170 9.31 0.00221 -0.505
## 3 49.8 147. 44.8 1.13 5.02 0.0148 9.31 0.00221 0.543
## 4 80.7 194. 91.9 1.07 -11.2 0.0131 9.30 0.00976 -1.21
## 5 95.9 193 91.4 1.05 4.51 0.0126 9.32 0.00152 0.488
## 6 44.8 150. 47.1 1.04 -2.32 0.0124 9.32 0.000397 -0.251Influence
As noted previously, observations of high leverage may or may not be influential. The influence of an observation depends not only on its leverage, but also on the magnitude of its residual. Recall that while leverage only takes into account the explanatory variable (x), the residual depends on the response variable (y) and the fitted value (y^).
Influential points are likely to have high leverage and deviate from the general relationship between the two variables. We measure influence using Cook’s distance, which incorporates both the leverage and residual of each observation.
# Rank points of high leverage
mod %>%
augment() %>%
arrange(desc(.cooksd)) %>%
head()## # A tibble: 6 x 9
## wgt hgt .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 73.2 151. 48.8 0.974 24.4 0.0109 9.25 0.0386 2.64
## 2 116. 178. 75.9 0.507 40.5 0.00296 9.14 0.0282 4.36
## 3 104. 165. 63.0 0.491 41.1 0.00279 9.14 0.0273 4.42
## 4 109. 190. 88.8 0.946 19.8 0.0103 9.28 0.0238 2.13
## 5 67.3 152. 50.1 0.922 17.2 0.00982 9.29 0.0171 1.86
## 6 76.8 158. 55.3 0.729 21.5 0.00613 9.27 0.0166 2.32Removing outliers
Observations can be outliers for a number of different reasons. Statisticians must always be careful-and more importantly, transparent-when dealing with outliers. Sometimes, a better model fit can be achieved by simply removing outliers and re-fitting the model. However, one must have strong justification for doing this. A desire to have a higher R2 is not a good enough reason!
In the mlbBat10 data, the outlier with an OBP of 0.550 is Bobby Scales, an infielder who had four hits in 13 at-bats for the Chicago Cubs. Scales also walked seven times, resulting in his unusually high OBP. The justification for removing Scales here is weak. While his performance was unusual, there is nothing to suggest that it is not a valid data point, nor is there a good reason to think that somehow we will learn more about Major League Baseball players by excluding him.
Nevertheless, we can demonstrate how removing him will affect our model.
# Create nontrivial_players
nontrivial_players<- mlbBat10 %>%
filter(AB>=10, OBP<0.5)
# Fit model to new data
mod_cleaner<-lm(SLG~OBP, data=nontrivial_players)
# View model summary
summary(mod_cleaner)##
## Call:
## lm(formula = SLG ~ OBP, data = nontrivial_players)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31383 -0.04165 -0.00261 0.03992 0.35819
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.043326 0.009823 -4.411 1.18e-05 ***
## OBP 1.345816 0.033012 40.768 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07011 on 734 degrees of freedom
## Multiple R-squared: 0.6937, Adjusted R-squared: 0.6932
## F-statistic: 1662 on 1 and 734 DF, p-value: < 2.2e-16# Visualize new model
ggplot(data=nontrivial_players, aes(x=OBP,y=SLG))+
geom_point()+geom_smooth(method = "lm")
##High leverage points Not all points of high leverage are influential. While the high leverage observation corresponding to Bobby Scales in the previous exercise is influential, the three observations for players with OBP and SLG values of 0 are not influential.
This is because they happen to lie right near the regression anyway. Thus, while their extremely low OBP gives them the power to exert influence over the slope of the regression line, their low SLG prevents them from using it.
# Rank high leverage points
mod %>%
augment() %>%
arrange(desc(.hat),desc(.cooksd)) %>%
head()## # A tibble: 6 x 9
## wgt hgt .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 85.5 198. 96.6 1.26 -11.1 0.0182 9.30 0.0134 -1.20
## 2 90.9 197. 95.6 1.21 -4.66 0.0170 9.31 0.00221 -0.505
## 3 49.8 147. 44.8 1.13 5.02 0.0148 9.31 0.00221 0.543
## 4 80.7 194. 91.9 1.07 -11.2 0.0131 9.30 0.00976 -1.21
## 5 95.9 193 91.4 1.05 4.51 0.0126 9.32 0.00152 0.488
## 6 44.8 150. 47.1 1.04 -2.32 0.0124 9.32 0.000397 -0.251