Simple Regression

Regression can be done in n-dimension, but we will start with a two dimension (one input one output) example. When regression uses only one input to predict an output it is called simple regression. Once we have multiple inputs we call it multiple regression. If we are trying to fit a straight line or plane that best represents the relationship between input and output (or X and Y) then it is called linear regression. In 1 dimension this is just a point (for example the mean). In 2 dimensions this is a line, in 3 dimensions it becomes a plane and so on.

Equation

The simple linear regression equation is:

\[ y_i = w_0 + w_1x_i + err_i \]

where \(y_i\) stands for the output for input i, \(w_0\) is the intercept term, and \(w_1\) is the slope of the line. This is known as simple linear regression simple because you are just fitting a line. In calculus you must have seen the equation for a line - and this is kind of the same idea. The slightly complicated part is finding which line fits the set of points the “best”.

Measuring Error

To decide what the best line is, we have to come up with an idea for measuring how good one fit is from another. We can see in many cases which line seems best just by observation, but when two different lines look equally good we need a way to differentiate them. For linear regression the common error measurement is known as the residual sum of squares (RSS).

\[RSS = \frac{1}{n}\sum_{i=1}^{n}(y_i - [w_0 + w_1x_i])^2\]

This is a relative measure of distance from the real values to the fitted line. In the initial equation the noise (\(err_i\)) is the vector that makes up for difference between the actual point and the neighbouring point on the line. We have to square the value because we want a distance and so whether the line is above or below the real point should not affect the distance. In theory there are many options for this error (we can use an absolute value difference too just as an example), but we use this squared error because it has many properties. But how do we find which line has the lowest RSS?

Closed form solution

For linear regression, if we used a squared error there is a closed form solution we can calculate to create the output. Closed form means essentially we can create a set of math operations, that given certain conditions, we will always get our perfect line.

If we “matrixify” or “vectorize” the RSS equation above, we get \(argmin_w||y-Xw||^2\) where \(y_{n,1}\) is a vector representing the outputs, \(X_{n, 2}\) is a matrix with the first column containing 1’s, and second column containing X values, and \(w_{2,1}\) is a weight vector containing \(w_0, w_1\). In other words we are trying to pick the \(w\) vector such that the expression \(||y-Xw||^2\) is minimized. If we just expand the equations:

\[||y-Xw||^2 = (y-Xw)^T(y-Xw) = (y^Ty - 2w^TX^Ty + w^TX^TXw)\]

Since we want to minimize this, we can take the first derivative with respect to w and see when it equals 0, giving: \(0=-2X^Ty+2X^TXw\) which implies \(w=(X^TX)^{-1}(X^Ty)\). Easy way to pick the weight vector, as long as \((X^TX)^{-1}\) is invertible. However, assuming it is not invertible, or inverting it takes extra time, we need another method.

Gradient Descent

You know from calculus finding the minimum of things usually involves some sort of a derivative. If we take every possible line going through the same (similar) intercepts, they will converge at a minimum RSS. (Think about why there can’t be two lines which both have a minimum RSS).

The algorithm used to find the best fitting line (within some error) is known as gradient descent. Gradient descent is built into most regression features and so it is unlikely you will ever have to use it, but over here I will give an overview of how it works. There are many resources online that can explain one of the steps or theory if it is unclear.

The idea is that given a step size (n), a weight vector (w), and a tolerance (err), as long as we have not converged we update the w vector as follows:

\[w^{(t+1)} = w^{(t)} - n\nabla g(w^{(t)})\]

In the case of simple linear regression, we want:

\[min_{w_0,w_1}\sum_{i=1}^{N}(y_i-[w_0+w_1 x_i])^2\]

So we need the gradient of the RSS expression with respect to \(w_0\) and \(w_1\), which (you can check it yourself) evaluates to:

\[ \nabla RSS(w_0, w_1) = \begin{bmatrix} -2\sum_{i=1}^{N}[y_i-(w_0+w_1x_i)] \\ -2\sum_{i=1}^{N}[y_i-(w_0+w_1x_i)]x_i \end{bmatrix} \]

This idea can translate to more than 2 dimensions and in that case rather than finding a local min/max in a curve, it will involve finding the local min/max in a plane or hyperplane. By running this algorithm we will keep adjusting our line until we are not able to minimize the RSS further within some tolerance. This enforces that we will get a best-fitting line. As you can imagine, it also works for curves.

However, does having a low RSS gurantee that the line we predict will be useful?

Plotting

First we should plot variables. Here I will explore a few more specific plots versus the freestyle ones you used in Titanic.

First we will explore scatter plots (which can check for a trend), boxplots (which can identify outliers), and density plots (to see if the response variable does not deviate too much; is close to a normal distribution or bell curve).

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.6.1

library(cowplot)

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:ggplot2':
## 
##     ggsave

# Scatter Plots
pv1 <- ggplot(cars, aes(x=1:nrow(cars), y=sort(speed))) + geom_point()
pv2 <- ggplot(cars, aes(x=1:nrow(cars), y=sort(dist))) + geom_point()
pv1

pv1 + geom_point(x=1:nrow(cars)) + geom_smooth(method=lm)

pv2

pv2 + geom_smooth(method=lm)

# Box Plots
plot_grid(ggplot(cars, aes(y=cars$dist)) + geom_boxplot(), ggplot(cars, aes(y=cars$speed)) + geom_boxplot(), ncol=2)

# Density Plots
plot_grid(ggplot(cars, aes(x=cars$dist)) + geom_density(), ggplot(cars, aes(x=cars$speed)) + geom_density(), ncol=2)

Creating a linear model

In R, you can define and save a linear model (or any model) the same way we can save functions or variables (with the assignment operator, “<-”). Below, we will create a linear model for predicting distance using speed. The key arguments to lm are “formula” and data; formula is formatted as output/Y ~ input(s)/X.

linearModel <- lm(dist ~ speed, data=cars)

print(linearModel)

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Coefficients:
## (Intercept)        speed  
##     -17.579        3.932

help(lm)

## starting httpd help server ... done

6. QUESTION: In the model above, what does “Intercept” and “speed” stand for using the regression equation defined at the beginning of this notebook (\(w_0, w_1\))?

Summarize specific features

We can see more information abut the model we built below:

summary(linearModel)

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

We will discuss how to interpret this output in the next class.

Make predictions

Here is a method for using the linear model above to predict data for new points (or different speeds). To predict “dist” we need a dataframe that has a column named “speed” (as our formula requires a speed to predict a distance). Also note that the predictions initially are labelled “1,2,3…” and not by the input value, so below I rename them just to provide clarity.

# speeds to predict in vector form
speeds <- c(2,4,6,8,10)

speeds2 <- c(2,3,5,15,26,30)

# Build dataframe
b <- data.frame(speed = speeds2)

# Make predictions
preds <- predict(linearModel, b)
preds

##          1          2          3          4          5          6 
##  -9.714277  -5.781869   2.082949  41.407036  84.663533 100.393168

names(preds) <- as.character(speeds)
preds

##          2          4          6          8         10       <NA> 
##  -9.714277  -5.781869   2.082949  41.407036  84.663533 100.393168

Plot Regression Line

Finally, we will plot our regression line on the points we have from the data using the following code:

plot(cars$speed, cars$dist)
abline(linearModel)

Code for multiple linear regression

Code for it is as follows. We use a dataset of Swiss fertility measure and socio-economic indicators for each of 47 French-speaking provinces of Switzerland around 1888. There are 47 observations and 6 variables: Fertility, Agriculture, Examination, Education, Catholic, and Infant Mortality. Here we will build a model that predicts Education using all the other variables. The “.” represents “everything not including the output”.

data("swiss")
swissModel <- lm(Education ~ ., data=swiss)

This is the same structure as “lm” with additional arguments added to the “input” part of the formula.

Interpreting Results

We can interpret the results in more detail using the “summary” command.

summary(swissModel)

## 
## Call:
## lm(formula = Education ~ ., data = swiss)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.3949  -2.3716  -0.2856   2.8108  11.2985 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      32.74414    8.87888   3.688 0.000657 ***
## Fertility        -0.40851    0.08585  -4.758 2.43e-05 ***
## Agriculture      -0.16242    0.04488  -3.619 0.000804 ***
## Examination       0.41980    0.16339   2.569 0.013922 *  
## Catholic          0.10023    0.02150   4.663 3.29e-05 ***
## Infant.Mortality  0.20408    0.28390   0.719 0.476305    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.907 on 41 degrees of freedom
## Multiple R-squared:  0.7678, Adjusted R-squared:  0.7395 
## F-statistic: 27.12 on 5 and 41 DF,  p-value: 5.223e-12

Now is a good time to understand what the summary command for the model says:

• The first element (“Call”) tells you back the formula that was inputted.
• The residuals section shows how the residuals are distributed. Residuals are the difference between actual points and model prediction - and if the residuals are distributed symmetrically around 0 that is a good sign. Usually it is a good idea to plot residuals and check the distribution, but that feature may not always be available to you. Generally, if you see an even distribution around 0 and the relationship does not follow a pattern of dispersion it is a good sign.

plot(swissModel$residuals, main="Residuals Plot")
abline(h=0)

The coefficients (next row) have many components.

• The first one is the estimate of the value of \(w_0\), \(w_1\) etc. in the regression equation.

• Each estimate represents that with no other factors, education will be 32.744%… and with every unit increase in Fertility, the estimate can increase by ~.40, and so on for the other features.

• The standard error estimates the average amount that estimates vary from the actual response. This means (intuitevly) that we want the standard error to be very small compared to the estimate. Using that logic it appears that in the above Fertility is the only really good indicator, and Infant.Mortality is certainly a bad one.

• Coefficient t-value represents how far (how many standard deviations) the estimate is from 0. An answer closer to zero implies that this factor does not impact the output variable much. Therefore, for a useful predictor, we want this value to be as far from zero as possible (such as the intercept, fertility and even agriculture).

• Pr(>|t|) is the p-value which is used in hypothesis testing. To give a tldr version of hypothesis testing:

You are almost always testing whether or not the null hypothesis is true (meaning that whatever input you are testing has no value on the response) and if the p-value is larger than 0.05 you typically accept that this null hypothesis is true otherwise you can reject it, and rejecting it means the input is important. Typically a value will get rejected for a reason discussed above (i.e too much variation or too much an intercept too close to 0).

• In the last section, there are many more measures of “goodness of fit” out of which the most notable one is “Adjusted-R Squared”. The Adjusted-R squared is an estimate of how much of the models response variable is explained by an input (in this case 74.25) which is not bad, but this depends on the context too.

• Some people also look at the F-statistic which is measured based on how much larger it is than the number of variables inputted. Either way, the point is to assess how good the “fit” is.

Discarding Variables

One biproduct of all of these pieces of information is that we can now discard features that are taking up time. That way, when we are predicting we can input less fields and also run the model on less variables (all saving runtime).

There are two easy ways you can choose to discard variables:

• 1. Remove features with a high p-value.
• 2. Try all subsets of features and choose the combination with the highest R-squared.

If we remove the variables we see in our original model that are bad, we get the following final model:

swissModelFull <- lm(Education ~ Fertility + Agriculture + Examination + Catholic, data=swiss)

We see that we slightly increase the adjusted R-squared by removing a variable! This is because the adjusted R squared is calculated based on how well the input variables explain the actual resulting variable, with respect to how many input variables there are.

Categorical variables

In general, we treat categorical variables the same way we treat numeric variables. However, the regression reads each category as a separate impact. For example, let’s say we run a linear regression on an output variable Y using a matrix X that consists of one numeric variable and 4 categories (code below):

Y <- unlist(lapply(1:100, function(x) -x^3 + 2*x^2 + sin(x) + cos(x)))
Y <- as.vector(scale(Y)) + rnorm(50, 0, 0.5)

X <- data.frame(ind = 1:100, mFn = unlist(lapply(1:100, function(x) -x^3 + 2*x^2 + sin(x) + cos(x))), mC1 = as.factor(sample(1:4, 100, replace=TRUE)), mC2 = as.factor(sample(c("a","b", "c", "d"), 100, replace=TRUE)))
X$mFn <- as.vector(scale(X$mFn))

tmpData <- cbind(Y, X)

myTmpModel <- lm(Y ~ ., data=tmpData)

summary(myTmpModel)

## 
## Call:
## lm(formula = Y ~ ., data = tmpData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84351 -0.32991 -0.02175  0.27510  1.20962 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.222236   0.234956   0.946    0.347    
## ind         -0.000812   0.003908  -0.208    0.836    
## mFn          0.897455   0.113379   7.916 5.73e-12 ***
## mC12        -0.182761   0.136485  -1.339    0.184    
## mC13        -0.192612   0.139092  -1.385    0.170    
## mC14        -0.170273   0.130115  -1.309    0.194    
## mC2b        -0.131756   0.128094  -1.029    0.306    
## mC2c        -0.059502   0.122562  -0.485    0.628    
## mC2d        -0.168158   0.133305  -1.261    0.210    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.449 on 91 degrees of freedom
## Multiple R-squared:  0.8204, Adjusted R-squared:  0.8046 
## F-statistic: 51.95 on 8 and 91 DF,  p-value: < 2.2e-16

As we see, each individual category has been broken up into 3 separate variables. If we do not include any of these variables that represents the first category. Then each intercept is indicated as "adding category X impacts our estimate by __".

Summary

In summary, all of the above methods are used a lot in traditional experiments, and are a large part of statistical theory. It is usually ignored today, but still good to know about as in the future if you get funny results this can be a good way to verify what happened. Another very important part of these applications is having the ability to discard useless variables. In the real world regression is not used on the front end of most applications in a modern day as we have better methods to define relationships - but it is still an important method to learn and understand for the basics.

Set1

Chemical Manufacturing Process Data Info here: https://rdrr.io/rforge/AppliedPredictiveModeling/man/ChemicalManufacturingProcess.html

data("ChemicalManufacturingProcess")

Set2

Abalone Data Info here: https://archive.ics.uci.edu/ml/datasets/abalone Google “abalone” to see what they are

data("abalone")

Set3

Diamonds Data Info: https://ggplot2.tidyverse.org/reference/diamonds.html

data("diamonds")

Set4

Concrete Data Info: http://archive.ics.uci.edu/ml/datasets/concrete+compressive+strength

data("concrete")

Set 5

Boston Housing Data Info here: http://ugrad.stat.ubc.ca/R/library/mlbench/html/BostonHousing.html

data("BostonHousing2")

Set 6

Ozone Data Info here: https://www.rdocumentation.org/packages/mlbench/versions/2.1-1/topics/Ozone

data("Ozone")

Set 7

NYC Flights data

data("flights")

Set 8

Atmosphere data from NASA

data("atmos")

Set 9

Boston housing data

data("Boston")

Polynomial Regression

Add polynomial/quadratic terms to a regression. This is usually done by comparing each individual input to the output in a graph and then choosing a function that you think best generalizes to the input. One example was in a previous day. We will use the Boston Housing dataset and same comparison of variables.

modelSqBoston <- lm(medv ~ poly(lstat, 2, raw=TRUE), data=trainData)
modelPolyBoston <- lm(medv ~ poly(lstat, 6, raw = TRUE), data = trainData)

summary(modelSqBoston)

## 
## Call:
## lm(formula = medv ~ poly(lstat, 2, raw = TRUE), data = trainData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.8127  -3.9684  -0.5957   2.4612  25.2395 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 43.635690   0.978154   44.61   <2e-16 ***
## poly(lstat, 2, raw = TRUE)1 -2.414321   0.138008  -17.49   <2e-16 ***
## poly(lstat, 2, raw = TRUE)2  0.045510   0.004141   10.99   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.508 on 402 degrees of freedom
## Multiple R-squared:  0.6544, Adjusted R-squared:  0.6527 
## F-statistic: 380.6 on 2 and 402 DF,  p-value: < 2.2e-16

summary(modelPolyBoston)

## 
## Call:
## lm(formula = medv ~ poly(lstat, 6, raw = TRUE), data = trainData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.1645  -3.2842  -0.7063   2.1914  26.7562 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  7.314e+01  5.991e+00  12.208  < 2e-16 ***
## poly(lstat, 6, raw = TRUE)1 -1.486e+01  3.171e+00  -4.686 3.83e-06 ***
## poly(lstat, 6, raw = TRUE)2  1.844e+00  6.099e-01   3.024  0.00266 ** 
## poly(lstat, 6, raw = TRUE)3 -1.217e-01  5.544e-02  -2.194  0.02880 *  
## poly(lstat, 6, raw = TRUE)4  4.235e-03  2.560e-03   1.654  0.09893 .  
## poly(lstat, 6, raw = TRUE)5 -7.375e-05  5.787e-05  -1.274  0.20330    
## poly(lstat, 6, raw = TRUE)6  5.088e-07  5.069e-07   1.004  0.31608    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.193 on 398 degrees of freedom
## Multiple R-squared:  0.6957, Adjusted R-squared:  0.6912 
## F-statistic: 151.7 on 6 and 398 DF,  p-value: < 2.2e-16

ggplot(trainData, aes(lstat, medv) ) +
  geom_point() +
  stat_smooth(method = lm, formula = y ~ poly(x, 6, raw = TRUE))

predictions <- modelPolyBoston %>% predict(testData)

RMSE(predictions, testData$medv)

## [1] 5.323196

Spline Regression

Fit a smooth curve with a series of polynomial segments (essentially compiling polynomial segments called knots). We set the knots to be at 1/4, 1/2, 3/4.

knots <- quantile(trainData$lstat, p = c(0.25, 0.5, 0.75))
modelSpline <- lm (medv ~ bs(lstat, knots = knots), data = trainData)
# Make predictions
predictions <- modelSpline %>% predict(testData)
# Model performance
RMSE(predictions, testData$medv)

## [1] 5.32017

ggplot(trainData, aes(lstat, medv) ) +
  geom_point() +
  stat_smooth(method = lm, formula = y ~ splines::bs(x, df = 3))

Generalized additive models

Do spline regression but choose each “knot” in an automated way.

modelGAM <- gam(medv ~ s(lstat), data = trainData)
# Make predictions
predictions <- modelGAM %>% predict(testData)
# Model performance
RMSE(predictions, testData$medv)

## [1] 5.293

ggplot(trainData, aes(lstat, medv) ) +
  geom_point() +
  stat_smooth(method = gam, formula = y ~ s(x))

These are just some methods. To be fully honest I don’t do that much regression on the day-to-day, so it’s a good idea for you to explore these in depth during the scope of this mini-assignment you have.