<jupyter_start><jupyter_text>Correlation and Regression in R > describe relationships between two numerical quantities and characterize these relationships graphically.- author: <NAME>- toc: true- comments: true- categories: [correlation, regression]- image: images/corr-reg-shield.png  OverviewUltimately, data analysis is about understanding relationships among variables. Exploring data with multiple variables requires new, more complex tools, but enables a richer set of comparisons. We will describe relationships between two numerical quantities. We will characterize these relationships graphically, in the form of summary statistics, and through simple linear regression models.  Libraries<jupyter_code>library(ggplot2)
library(openintro)
library(dplyr)
library(tidyr)
library(readr)
library(broom)<jupyter_output>Warning message:
"package 'broom' was built under R version 3.6.3"<jupyter_text>Visualizing two variablesThe techniques for exploring bivariate relationships.  Visualizing bivariate relationships>  Bivariate  relationships- Both variables are numerical - Response variable   - a.k.a. y, dependent- Explanatory variable   - Something you think might be related to the response   - a.k.a. x, independent, predictor>  Graphical representations- Put response on vertical axis - Put explanatory on horizontal axis<jupyter_code>ggplot()+
    geom_line()+
    xlab("Independent (Explanatory Variable)") +
    ylab("Dependent (Response Variable)")
head(possum)<jupyter_output><empty_output><jupyter_text>ScatterplotScatterplots are the most common and effective tools for visualizing the relationship between two numeric variables.<jupyter_code>possum %>%
    ggplot(aes(x=tail_l, y=total_l)) +
    geom_point()
possum %>%
    ggplot(aes(x=tail_l, y=total_l)) +
    geom_point() +
    scale_x_continuous("Length of Possum Tail (cm)") +
    scale_y_continuous("Length of Possum Body (cm)")<jupyter_output><empty_output><jupyter_text>- Can think of boxplots as scatterplots......  - but with discretized explanatory variable- `cut()` function discretizes   - Choose appropriate number of "boxes"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.<jupyter_code>possum %>%
    ggplot(aes(x=cut(tail_l, breaks=5), y=total_l)) +
    geom_point()
possum %>%
    ggplot(aes(x=cut(tail_l, breaks=5), y=total_l)) +
    geom_boxplot()<jupyter_output><empty_output><jupyter_text>ncbirths<jupyter_code>head(ncbirths)<jupyter_output><empty_output><jupyter_text>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).<jupyter_code>ncbirths %>%
    ggplot(aes(x=weeks, y=weight)) +
    geom_point() +
    geom_smooth(method="lm", se=F)<jupyter_output>`geom_smooth()` using formula 'y ~ x'
Warning message:
"Removed 2 rows containing non-finite values (stat_smooth)."Warning message:
"Removed 2 rows containing missing values (geom_point)."<jupyter_text>Boxplots as discretized/conditioned scatterplots We you can think of boxplots as scatterplots for which the variable on the x-axis has been discretized.<jupyter_code>ncbirths %>%
    ggplot(aes(x=cut(weeks, breaks=5), y=weight)) +
    geom_boxplot()<jupyter_output><empty_output><jupyter_text>Characterizing bivariate relationships- Form(e.g.linear,quadratic,non-linear)- Direction(e.g.postive,negative)- Strength(how much scatter/noise?)- Outliers  Creating scatterplotswe will be using several datasets listed below. These data are available through the `openintro` package. Briefly:- The `mammals` dataset contains information about 39 different species of mammals, including their body weight, brain weight, gestation time, and a few other variables.- The `mlbBat10` dataset contains batting statistics for 1,199 Major League Baseball players during the 2010 season.- The `bdims` dataset contains body girth and skeletal diameter measurements for 507 physically active individuals.- The `smoking` dataset contains information on the smoking habits of 1,691 citizens of the United Kingdom.<jupyter_code>head(mammals)<jupyter_output><empty_output><jupyter_text>how the brain weight of a mammal varies as a function of its body weight.<jupyter_code>mammals %>%
    ggplot(aes(x=body_wt, y=brain_wt)) + 
    geom_point()
head(mlbbat10)<jupyter_output><empty_output><jupyter_text>how the slugging percentage (SLG) of a player varies as a function of his on-base percentage (OBP).<jupyter_code># Baseball player scatterplot
mlbbat10 %>%
  ggplot(aes(y=slg, x=obp)) +
  geom_point()

head(bdims)<jupyter_output><empty_output><jupyter_text>how a person's weight varies as a function of their height. We'll use color to separate by sex, which we'll need to coerce to a factor with `factor()`.<jupyter_code># Body dimensions scatterplot
bdims%>%
  ggplot(aes(x=hgt, y=wgt, color=factor(sex))) +
  geom_point()
head(smoking)<jupyter_output><empty_output><jupyter_text>how the amount that a person smokes on weekdays varies as a function of their age.<jupyter_code>smoking %>%
    ggplot(aes(x=age, y=amt_weekdays)) + 
    geom_point()
head(county)
table(county$median_edu)<jupyter_output><empty_output><jupyter_text>TransformationsThe 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.`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.<jupyter_code>mammals %>%
    ggplot(aes(x=body_wt, y=brain_wt)) +
    geom_point() +
    coord_trans(x="log10", y="log10")
mammals %>%
    ggplot(aes(x=body_wt, y=brain_wt)) +
    geom_point() +
    scale_x_log10() +
    scale_y_log10()<jupyter_output><empty_output><jupyter_text>Outliers<jupyter_code>mlbbat10 %>%
    ggplot(aes(x=stolen_base, y=home_run)) +
    geom_point()
mlbbat10 %>%
    ggplot(aes(x=stolen_base, y=home_run)) +
    geom_point(alpha=.5)<jupyter_output><empty_output><jupyter_text>>  Add some jitter<jupyter_code>mlbbat10 %>%
    ggplot(aes(x=stolen_base, y=home_run)) +
    geom_point(alpha=.5, position="jitter")<jupyter_output><empty_output><jupyter_text>>  Identify the outliers<jupyter_code>head(
    mlbbat10 %>%
        filter(stolen_base>60 | home_run>50) %>%
            select(name, team, position, stolen_base, home_run)
)<jupyter_output><empty_output><jupyter_text>Identifying outliersin the baseball dataset, 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.<jupyter_code># Filter for AB greater than or equal to 200
ab_gt_200 <- mlbbat10 %>%
  filter(at_bat>=200) 

# Scatterplot of SLG vs. OBP
ggplot(ab_gt_200, aes(x = obp, y = slg)) +
  geom_point()
# Identify the outlying player
ab_gt_200 %>%
  filter(obp<.2)<jupyter_output><empty_output><jupyter_text>CorrelationCorrelation as a means of quantifying bivariate relationships.  Quantifying the strength of bivariate relationships>  Correlation- Correlation coefcient between -1 and 1- Sign —> direction- Magnitude —> strength>  Non-linear correlation<jupyter_code>head(run09)
run09 %>%
    filter(div_place <= 10) %>%
        ggplot(aes(x=age, y=pace, color=gender)) +
        geom_point() +
        geom_smooth(se=F)<jupyter_output>`geom_smooth()` using method = 'loess' and formula 'y ~ x'<jupyter_text>>  Pearson product-moment correlation$$\begin{split}r(x,y)&=\frac{Cov(x,y)}{\sqrt{SXX.SYY}}\\&=\frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i - \bar{x})^2.(y_i - \bar{y})^2}}\end{split}$$Correlation coefficients cannot be greater than 1.  Computing correlationThe `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.<jupyter_code># Compute correlation
ncbirths %>%
  summarize(N = n(), r = cor(mage, weight))
# Compute correlation for all non-missing pairs
ncbirths %>%
  summarize(N = n(), r = cor(weeks, weight, use = "pairwise.complete.obs"))<jupyter_output><empty_output><jupyter_text>The Anscombe datasetIn 1973, <NAME> famously created four datasets with remarkably similar numerical properties, but obviously different graphic relationships. The `anscombe` dataset contains the `x` and `y` coordinates for these four datasets, along with a grouping variable, set, that distinguishes the quartet.<jupyter_code>head(anscombe)
head(anscombe2 <- anscombe %>%
    gather(x_name, x, x1:x4) %>%
        gather(y_name, y, y1:y4), 10)
anscombe %>%
    ggplot(aes(x=x1,y=y1)) +
    geom_point()
anscombe %>%
    ggplot(aes(x=x2,y=y2)) +
    geom_point()
anscombe %>%
    ggplot(aes(x=x3,y=y3)) +
    geom_point()
anscombe %>%
    ggplot(aes(x=x4,y=y4)) +
    geom_point()
table(anscombe2$x_name, anscombe2$y_name)
# Compute properties of Anscombe
anscombe2 %>%
  group_by(x_name) %>%
  summarize(
    N = n(), 
    mean_of_x = mean(x), 
    std_dev_of_x = sd(x), 
    mean_of_y = mean(y), 
    std_dev_of_y = sd(y), 
    correlation_between_x_and_y = cor(x,y)
  )
anscombe2 <- anscombe2 %>%
    filter((x_name=="x1"&y_name=="y1") | (x_name=="x2"&y_name=="y2") | (x_name=="x3"&y_name=="y3") | (x_name=="x4"&y_name=="y4"))
anscombe2 %>%
    ggplot(aes(x=x, y=y)) +
    geom_point() +
    facet_wrap(~x_name)
# Compute properties of Anscombe
anscombe2 %>%
  group_by(x_name) %>%
  summarize(
    N = n(), 
    mean_of_x = mean(x), 
    std_dev_of_x = sd(x), 
    mean_of_y = mean(y), 
    std_dev_of_y = sd(y), 
    correlation_between_x_and_y = cor(x,y)
  )<jupyter_output>`summarise()` ungrouping output (override with `.groups` argument)<jupyter_text>all of the measures are identical (ignoring rounding error) across the four different sets.  Perception of correlation (2)Estimating the value of the correlation coefficient between two quantities from their scatterplot can be tricky. [Statisticians have shown that](https://statistics.stanford.edu/sites/default/files/EFS%20NSF%20206.pdf) people's perception of the strength of these relationships can be influenced by design choices like the x and y scales.<jupyter_code># Run this and look at the plot
ggplot(data = mlbbat10, aes(x = obp, y = slg)) +
  geom_point()


# Correlation for all baseball players
mlbbat10 %>%
  summarize(N = n(), r = cor(obp, slg))
# Run this and look at the plot
mlbbat10 %>% 
    filter( at_bat> 200) %>%
    ggplot(aes(x = obp, y = slg)) + 
    geom_point()
# Correlation for all players with at least 200 ABs
mlbbat10 %>%
  filter(at_bat>=200) %>%
  summarize(N = n(), r = cor(obp, slg))
# Run this and look at the plot
ggplot(data = bdims, aes(x = hgt, y = wgt, color = factor(sex))) +
  geom_point() 


# Correlation of body dimensions
bdims %>%
  group_by(sex) %>%
  summarize(N = n(), r = cor(hgt, wgt))
# Run this and look at the plot
ggplot(data = mammals, aes(x = body_wt, y = brain_wt)) +
  geom_point() + scale_x_log10() + scale_y_log10()


# Correlation among mammals, with and without log
mammals %>%
  summarize(N = n(), 
            r = cor(body_wt, brain_wt), 
            r_log = cor(log(body_wt), log(brain_wt)))<jupyter_output><empty_output><jupyter_text>Spurious correlations  Spurious correlation in random dataStatisticians 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.  Simple linear regressionWith the notion of correlation under the belt, we'll now turn our attention to simple linear models.  Visualization of Linear Models Possums<jupyter_code>possum %>%
    ggplot(aes(y=total_l, x=tail_l)) +
    geom_point()<jupyter_output><empty_output><jupyter_text>Through the origin<jupyter_code>possum %>%
    ggplot(aes(y=total_l, x=tail_l)) +
    geom_point() +
    geom_abline(intercept=0, slope=2.5)<jupyter_output><empty_output><jupyter_text>Through the origin, better fit<jupyter_code>possum %>%
    ggplot(aes(y=total_l, x=tail_l)) +
    geom_point() +
    geom_abline(intercept=40, slope=1.3)<jupyter_output><empty_output><jupyter_text>The "best" fit lineThe 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 us 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 one 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"`.<jupyter_code>possum %>%
    ggplot(aes(y=total_l, x=tail_l)) +
    geom_point() +
    geom_smooth(method="lm")<jupyter_output>`geom_smooth()` using formula 'y ~ x'<jupyter_text>Ignore standard errors`geom_smooth()` also takes an se argument that controls the standard error, which we will ignore for now.<jupyter_code>possum %>%
    ggplot(aes(y=total_l, x=tail_l)) +
    geom_point() +
    geom_smooth(method="lm", se=F)<jupyter_output>`geom_smooth()` using formula 'y ~ x'<jupyter_text>a scatterplot of body weight as a function of height for all individuals in the `bdims` dataset with a simple linear model plotted over the data.<jupyter_code># Scatterplot with regression line
ggplot(data = bdims, aes(x = hgt, y = wgt)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)
add_line <- function (my_slope) {
    bdims_summary <- bdims %>%
        summarize(N = n(), r = cor(hgt, wgt), 
                  mean_hgt = mean(hgt), 
                  mean_wgt = mean(wgt),
                  sd_hgt = sd(hgt), 
                  sd_wgt = sd(wgt)) %>%
            mutate(true_slope = r * sd_wgt / sd_hgt, 
                   true_intercept = mean_wgt - true_slope * mean_hgt)
    p <- ggplot(data = bdims, aes(x = hgt, y = wgt)) + 
        geom_point() + 
        geom_point(data = bdims_summary,
                   aes(x = mean_hgt, y = mean_wgt),
                   color = "red", size = 3)
    my_data <- bdims_summary %>%
        mutate(my_slope = my_slope,
               my_intercept = mean_wgt - my_slope * mean_hgt)
    p + geom_abline(data = my_data, 
                    aes(intercept = my_intercept, slope = my_slope), color = "dodgerblue")
}<jupyter_output><empty_output><jupyter_text>Uniqueness of least squares regression lineThe least squares criterion implies that the slope of the regression line is unique. The slope is computed by R. We 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.<jupyter_code>add_line(1)<jupyter_output><empty_output><jupyter_text>Understanding Linear Models Generic statistical model$response = f(explanatory) + noise$ Generic linear model$response = intercept + (slope * explanatory) + noise$ Regression model$Y = \beta_{0} + \beta_{1}.X+\epsilon, \ \ \ \epsilon ~ N(0, \sigma_{\epsilon})$ Fitted values$\hat{Y} = \hat{\beta_0} + \hat{\beta_1} . X$ Residuals$e = Y -\hat{Y}$ Fitting procedure- Given $n$ observations of pairs $(x_i, y_i)$...- Find $\hat{\beta_0}$, $\hat{\beta_1}$  that minimize $\sum_{i=1}^{n} {e_i}^2$ Least squares- Easy, deterministic,  unique solution- Residuals sum to zero - Line must pass through ($\bar{x}$, $\bar{y}$)- Other criteria exist Key concepts- $\hat{Y}$ is expected value given corresponding $X$- $\hat{\beta}$ are estimates of true, unknown betas $\beta$- Residuals ($e$'s) are estimates of true, unknown epsilons $\epsilon$- "Error" may be misleading term—better: noise  Fitting a linear model "by hand"the simple linear regression model:$$Y = b_0 + b_1.X$$Two facts enable us to compute the slope $b_1$ and intercept $b_0$ of a simple linear regression model from some basic summary statistics.- First, the slope can be defined as:$$b_1 = r_{X, Y} . \frac{s_Y}{s_X}$$Where $r_{X, Y}$ is the correlation (`cor()`) of $X$ and $Y$ and $s_X$ and $s_Y$ represent the standard deviation (`sd()`) of $X$ and $Y$ respectively.- Second, the point $(\bar{x}, \bar{y})$ is always on the least squares regression line, where $\bar{x}$ and $\bar{y}$ denote the average of $x$ and $y$ respectively.<jupyter_code>(bdims_summary <- bdims %>%
        summarize(N = n(), r = cor(hgt, wgt), 
                  mean_hgt = mean(hgt), 
                  mean_wgt = mean(wgt),
                  sd_hgt = sd(hgt), 
                  sd_wgt = sd(wgt)) %>%
            mutate(true_slope = r * sd_wgt / sd_hgt, 
                   true_intercept = mean_wgt - true_slope * mean_hgt))<jupyter_output><empty_output><jupyter_text>Regression vs. regression to the mean Heredity- Galton's "regression to the mean"- **Thought experiment**: consider the heights of the children of NBAplayers Regression modeling- **"Regression"**: techniques for modeling a quantitative response- Types of regression models:  - Least squares  - Weighted  - Generalized  - Non parametric  - Ridge  - Bayesian. . .  Regression to the mean_Regression to the mean_ is a concept attributed to <NAME>. 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.  Interpreting regression modelsHow to interpret the coefficients in a regression model.  Interpretation of Regression Is that textbook overpriced?<jupyter_code>head(textbooks)<jupyter_output><empty_output><jupyter_text>Compared to the course number?<jupyter_code>textbooks %>%
    mutate(course_number = as.integer(gsub("\\D", "", course))) %>%
        ggplot(aes(x=course_number, y=ucla_new)) +
        geom_point()<jupyter_output><empty_output><jupyter_text>Compared to Amazon?<jupyter_code>names(textbooks)
textbooks %>%
    ggplot(aes(x=amaz_new, y=ucla_new)) +
    geom_point()
textbooks %>%
    ggplot(aes(x=amaz_new, y=ucla_new)) +
    geom_point() +
    geom_smooth(method="lm", se=F)<jupyter_output>`geom_smooth()` using formula 'y ~ x'<jupyter_text>Slope and intercept<jupyter_code>lm(ucla_new ~ amaz_new, data=textbooks)<jupyter_output><empty_output><jupyter_text>$$\widehat{ucla\_new} = 0.929 + {amaz\_new}*1.119$$ Units and scale<jupyter_code>textbooks %>%
    mutate(amaz_new_cents = amaz_new *100) %>%
        lm(ucla_new ~ amaz_new_cents, data=.)<jupyter_output><empty_output><jupyter_text>Fitting simple linear modelsWhile 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 modelThe `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.<jupyter_code># Linear model for weight as a function of height
lm(wgt ~ hgt, data = bdims)<jupyter_output><empty_output><jupyter_text>A person who is 170 cm tall is expected to weigh about 68 kg.<jupyter_code># Linear model for SLG as a function of OBP
lm(slg ~ obp, data=mlbbat10)


# Log-linear model for body weight as a function of brain weight
lm(log(body) ~ log(brain), data=mammals)<jupyter_output><empty_output><jupyter_text>linear model objectAn "lm" object contains a host of information about the regression model that we fit. There are various ways of extracting different pieces of information.<jupyter_code>class(mod<- lm(ucla_new ~ amaz_new, data=textbooks))
mod<jupyter_output><empty_output><jupyter_text>Fitted coefcientsThe `coef()` function displays only the values of the coefficients.<jupyter_code>coef(mod)<jupyter_output><empty_output><jupyter_text>SummaryThe `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 $R^2$, adjusted $R^2$, 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.)<jupyter_code>summary(mod)<jupyter_output><empty_output><jupyter_text>Fitted valuesOnce you have fit a regression model, you are often interested in the fitted values ($\hat{y}_i$) and the residuals ($e_i$) where $i$ indexes the observations. $$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.we will confirm these two mathematical facts by accessing the fitted values and residuals with the `fitted.values()` and `residuals()` functions, respectively<jupyter_code>head(fitted.values(mod), 10)
# Mean of weights equal to mean of fitted values?
head(mean(bdims$wgt) == fitted.values(mod), 10)<jupyter_output><empty_output><jupyter_text>Residuals<jupyter_code>head(residuals(mod), 10)
# Mean of the residuals
mean(residuals(mod))<jupyter_output><empty_output><jupyter_text>BroomAs we 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. $\hat{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.<jupyter_code>head(bdims_tidy <- augment(mod), 10)
glimpse(bdims_tidy)<jupyter_output>Rows: 73
Columns: 9
$ ucla_new   <dbl> 27.67, 40.59, 31.68, 16.00, 18.95, 14.95, 24.70, 19.50, ...
$ amaz_new   <dbl> 27.95, 31.14, 32.00, 11.52, 14.21, 10.17, 20.06, 16.66, ...
$ .fitted    <dbl> 34.44105, 38.26587, 39.29701, 14.74146, 17.96678, 13.122...
$ .se.fit    <dbl> 1.460257, 1.418184, 1.407411, 1.720649, 1.673789, 1.7446...
$ .resid     <dbl> -6.7710547, 2.3241308, -7.6170104, 1.2585386, 0.9832248,...
$ .hat       <dbl> 0.01944321, 0.01833896, 0.01806141, 0.02699567, 0.025545...
$ .sigma     <dbl> 10.51520, 10.54319, 10.50682, 10.54581, 10.54624, 10.544...
$ .cooksd    <dbl> 4.226829e-03, 4.686540e-04, 4.954865e-03, 2.059099e-04, ...
$ .std.resid <dbl> -0.65294223, 0.22399308, -0.73400193, 0.12183273, 0.0951...<jupyter_text>Using linear model  Examining residuals<jupyter_code>bdims_tidy %>%
    arrange(desc(.resid)) %>%
        head(10)<jupyter_output><empty_output><jupyter_text>Markup<jupyter_code>textbooks %>%
    filter(ucla_new == 197)<jupyter_output><empty_output><jupyter_text>Making predictionsThe `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.<jupyter_code>predict(mod, new_data <- data.frame(amaz_new=8.49))<jupyter_output><empty_output><jupyter_text>Visualize new observationsThe `geom_smooth()` function makes it easy to add a simple linear regression line to a scatterplot of the corresponding variables. And in fact, there are more complicated regression models that can be visualized in the data space with `geom_smooth()`. However, there may still be times when we will want to add regression lines to our scatterplot manually. To do this, we can use the `geom_abline()` function, which takes slope and intercept arguments. Naturally, we have to compute those values ahead of time, but we already saw how to do this (e.g. using `coef()`).<jupyter_code>isrs <- augment(mod, newdata=new_data)
textbooks %>%
    ggplot(aes(x=amaz_new, y=ucla_new)) +
    geom_point() +
    geom_smooth(method="lm") +
    geom_point(data=isrs,aes(y=.fitted), size=5, color="red")<jupyter_output>`geom_smooth()` using formula 'y ~ x'<jupyter_text>Model FitHow to assess the "fit" of a simple linear regression model.  Assessing Model Fit  How well does our textbook model t?<jupyter_code>textbooks %>%
    ggplot(aes(x=amaz_new, y=ucla_new)) +
    geom_point() +
    geom_smooth(method="lm", se=F)<jupyter_output>`geom_smooth()` using formula 'y ~ x'<jupyter_text>How well does our possum model t?<jupyter_code>possum %>%
    ggplot(aes(x=tail_l, y=total_l)) +
    geom_point() +
    geom_smooth(method="lm", se=F)<jupyter_output>`geom_smooth()` using formula 'y ~ x'<jupyter_text>SSEOne 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.<jupyter_code>(mod_possum <- lm(total_l ~ tail_l, data=possum)) %>%
    augment() %>%
        summarise(SSE=sum(.resid^2), SSE_also=(n()-1)*var(.resid))<jupyter_output><empty_output><jupyter_text>RMSEsome 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_ie_i^2}{d.f}} = \sqrt{\frac{SSE}{n-2}}$$You can recover the residuals from `mod` with `residuals()`, and the degrees of freedom with `df.residual()`.  Residual standard error (possums)<jupyter_code>summary(mod_possum)<jupyter_output><empty_output><jupyter_text>Residual standard error (textbooks)<jupyter_code>summary(mod)
# Compute the mean of the residuals
mean(residuals(mod))
# Compute RMSE
sqrt(sum(residuals(mod)^2) / df.residual(mod))<jupyter_output><empty_output><jupyter_text>Comparing model fits  Null (average) model- For all observations...$\hat{y} = \bar{y}$  SSE, null model<jupyter_code>(mod_null <- lm(total_l~1, data=possum) %>%
    augment(possum) %>%
        summarise(SST=sum(.resid^2)))<jupyter_output><empty_output><jupyter_text>SSE, our model<jupyter_code>augment(mod_possum) %>%
    summarise(SSE = sum(.resid^2))<jupyter_output><empty_output><jupyter_text>Coefcient of determination$$R^2 = 1 - \frac{SSE}{SST} = 1 - \frac{var(e)}{var(y)}$$$e$ is the vector of residuals and $y$ is the response variable. This gives us the interpretation of $R^2$ 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. Connection to correlation- For simple linear regression..$$r_{x,y}^2 = R^2$$<jupyter_code>summary(mod_possum)<jupyter_output><empty_output><jupyter_text>Over-reliance on R-squared> twitter: https://twitter.com/VictorOmondi197/status/1285197027429756928  Assessing simple linear model fit<jupyter_code>summary(mod_bdims <- lm(wgt ~ hgt, data=bdims))
# Compute R-squared
bdims_tidy <- augment(mod_bdims)
bdims_tidy %>%
  summarize(var_y = var(wgt), var_e = var(.resid)) %>%
  mutate(R_squared = 1 - (var_e/var_y))<jupyter_output><empty_output><jupyter_text>This means that 51.4% of the variability in weight is explained by height.  Linear vs. averageThe $R^2$ 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 $R^2$ of zero because $SSE = SST$. That is, since the fitted values ($\hat{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.  Unusual Points<jupyter_code>(regulars <- mlbbat10 %>%
    filter(at_bat > 400)
) %>%
    ggplot(aes(x=stolen_base, y=home_run)) +
    geom_point() +
    geom_smooth(method="lm", se=F)<jupyter_output>`geom_smooth()` using formula 'y ~ x'<jupyter_text>LeverageThe 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.$$h_i = \frac{1}{n} + \frac{(x_i - \bar{x})^2}{\sum_{i=1}^{n} (x_i - \bar{x})^2}$$<jupyter_code>mod_regs <- lm(home_run ~ stolen_base, data=regulars) 
mod_regs %>%
    augment() %>%
        arrange(desc(.hat)) %>%
            dplyr::select(home_run, stolen_base, .fitted, .resid, .hat) %>%
            head()<jupyter_output><empty_output><jupyter_text>Influence via Cook's distanceAs 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.while leverage only takes into account the explanatory variable ($x$) , the residual depends on the response variable ($y$) and the fitted value ($\hat{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.<jupyter_code>mod_regs <- lm(home_run ~ stolen_base, data=regulars) 
mod_regs %>%
    augment() %>%
        arrange(desc(.cooksd)) %>%
            dplyr::select(home_run, stolen_base, .fitted, .resid, .hat, .cooksd) %>%
            head()<jupyter_output><empty_output><jupyter_text>Dealing with OutliersObservations 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 $R^2$ is not a good enough reason!In the mlbBat10 data, the outlier with an OBP of 0.550 is <NAME>, 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.  Removing outliers<jupyter_code># Create nontrivial_players
nontrivial_players <- mlbbat10 %>%
  filter(at_bat >= 10, obp<.5)


# Fit model to new data
mod_cleaner = lm(slg ~ obp, data=nontrivial_players)

# View model summary
summary(mod_cleaner)


# Visualize new model
nontrivial_players %>%
  ggplot(aes(x=obp, y=slg))+
    geom_point() +
      geom_smooth(method="lm")<jupyter_output>`geom_smooth()` using formula 'y ~ x'<jupyter_text>High leverage pointsNot all points of high leverage are influential. While the high leverage observation corresponding to Bobby Scales 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.<jupyter_code># Rank high leverage points
mod_cleaner %>%
  augment() %>%
    arrange(desc(.hat)) %>%
      head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10417205785100205, 'fraction_numerical': 0.018603180079933293, 'mean_word_length': 4.934140931581641, 'pattern_counts': {'":': 0, '<': 264, '<?xml version=': 0, '>': 346, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '60', 'dir': 'jupyter-structured-clean-dedup', 'id': '477058', 'license': 'MIT', 'n_tokens_mistral': 11731, 'n_tokens_neox': 10685, 'n_words': 3755, 'path': 'Correlation and Regression in R.ipynb', 'repo_name': 'VictorOmondi1997/Correlation-and-Regression-in-R', 'size': '34814'}