Regression in machine learning
Syed Zakir Abbas Zaidi
3/1/2022
In the discipline of machine learning, regression analysis is a key concept. It’s classified as supervised learning because the algorithm is taught both input and output labels. By estimating how one variable influences the other, it aids in the establishment of a link between the variables. Assume you’re in the market for a car and have decided that gas mileage will be a deciding factor in your purchase. How would you go about predicting the miles per gallon of some prospective rides? Because you know the car’s many characteristics (weight, horsepower, displacement, and so on), regression is a viable option. You can use regression techniques to identify the relationship between the MPG and the input data by plotting the average MPG of each automobile given its features. The regression function might be written as \(Y = f(X)\), with Y representing the MPG and X being the input features like as weight, displacement, horsepower, and so on. The desired function is \(f\), and this curve lets us determine if buying or not buying is helpful. Regression is the name for this technique. ## What is R? R is a statistical computing and graphics language and environment. It is a GNU project that is similar to the S language and environment established by John Chambers and colleagues at Bell Laboratories (previously AT&T, now Lucent Technologies). R can be thought of as a more advanced version of S. Although there are some significant differences, much of the code built for S works in R without modification. R is highly extendable and offers a wide range of statistical (linear and nonlinear modelling, classical statistical tests, time-series analysis, classification, clustering, etc.) and graphical tools. The S programming language is frequently used for statistical methods research, while R provides an Open-Source option for getting involved. One of R’s advantages is how simple it is to create well-designed publication-quality graphs, complete with mathematical symbols and calculations when needed. The defaults for small design choices in visuals have been carefully chosen, but the user retains complete control. R is accessible in source code form as Free Software under the provisions of the Free Software Foundation’s GNU General Public License. It compiles and operates on a wide range of UNIX and related systems (including FreeBSD and Linux), as well as Windows and MacOS. Only half of the job is done when it comes to creating a linear regression model. The model must correspond to the assumptions of linear regression in order to be used in practice, and there are 10 assumptions on which a linear regression model is based. These ten assumptions are: The regression model is linear in parameters The mean of residuals is zero Homoscedasticity of residuals or equal variance No autocorrelation of residuals The X variables and residuals are uncorrelated The number of observations must be greater than number of Xs The variability in X values is positive The regression model is correctly specified No perfect multicollinearity Normality of residuals ## Assumption 1 The regression model is linear in parameters An example of model equation that is linear in parameters Y = a + (β1X1) + (β2X2) Though, the X2 is raised to power 2, the equation is still linear in beta parameters. So, the assumption is satisfied in this case. ## Assumption 2 The mean of residuals is zero Check the mean of the residuals. If it zero (or very close), then this assumption is held true for that model. This is default unless you explicitly make amends, such as setting the intercept term to zero.
mod <- lm(dist ~ speed, data=cars)
mean(mod$residuals)## [1] 8.65974e-17Since the mean of residuals is approximately zero, this assumption holds true for this model. ## Assumption 3 Homoscedasticity of residuals or equal variance Once the regression model is built, set par (mfrow=c(2, 2)), then, plot the model using plot(lm.mod). This produces four plots. The top-left and bottom-left plots shows how the residuals vary as the fitted values increase.
par(mfrow=c(2,2)) # set 2 rows and 2 column plot layout
mod_1 <- lm(mpg ~ disp, data=mtcars) # linear model
plot(mod_1)
From the first plot (top-left), as the fitted values along x increase, the residuals decrease and then increase. This pattern is indicated by the red line, which should be approximately flat if the disturbances are homoscedastic. The plot on the bottom left also checks this and is more convenient as the disturbance term in Y axis is standardized. In this case, there is a definite pattern noticed. So, there is heteroscedasticity. Lets check this on a different model.
mod <- lm(dist ~ speed, data=cars[1:20, ]) # linear model
plot(mod)
Now, the points appear random, and the line looks pretty flat, with no increasing or decreasing trend. So, the condition of homoscedasticity can be accepted. Assumption 4 No autocorrelation of residuals This is applicable especially for time series data. Autocorrelation is the correlation of a time Series with lags of itself. When the residuals are autocorrelated, it means that the current value is dependent of the previous (historic) values and that there is a definite unexplained pattern in the Y variable that shows up in the disturbances. Below, are 3 ways you could check for autocorrelation of residuals. Using acf plot
# Method 1: Visualize with acf plot
library(ggplot2)
data(economics)
lmMod <- lm(pce ~ pop, data=economics)
acf(lmMod$residuals) # highly autocorrelated from the picture.
The X axis corresponds to the lags of the residual, increasing in steps of 1. The very first line (to the left) shows the correlation of residual with itself (Lag0), therefore, it will always be equal to 1. If the residuals were not autocorrelated, the correlation (Y-axis) from the immediate next line onwards will drop to a near zero value below the dashed blue line (significance level). Clearly, this is not the case here. So we can conclude that the residuals are autocorrelated. Using runs test
# Method 2: Runs test to test for randomness
library(lawstat)
runs.test(lmMod$residuals)##
## Runs Test - Two sided
##
## data: lmMod$residuals
## Standardized Runs Statistic = -23.812, p-value < 2.2e-16#=> Runs Test - Two sidedWith a p-value < 2.2e-16, we reject the null hypothesis that it is random. This means there is a definite pattern in the residuals. Using Durbin-Watson test.
# Method 3: Durbin-Watson test
lmtest::dwtest(lmMod)##
## Durbin-Watson test
##
## data: lmMod
## DW = 0.002159, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0#=> Durbin-Watson test#So, DW test also confirms our finding. To rectify it, add lag1 of residual as an X variable to the original model. This can be conveniently done using the slide function in Data Combine package.
library(DataCombine)
econ_data <- data.frame(economics, resid_mod1=lmMod$residuals)
econ_data_1 <- slide(econ_data, Var="resid_mod1", NewVar = "lag1", slideBy = -1)##
## Remember to put econ_data in time order before running.##
## Lagging resid_mod1 by 1 time units.econ_data_2 <- na.omit(econ_data_1)
lmMod2 <- lm(pce ~ pop + lag1, data=econ_data_2)Let’s check if the problem of autocorrelation of residuals is taken care of using this method.
acf(lmMod2$residuals)
Unlike the acf plot of lmMod, the correlation values drop below the dashed blue line from lag1 itself. So autocorrelation can’t be confirmed.
runs.test(lmMod2$residuals) #runs test##
## Runs Test - Two sided
##
## data: lmMod2$residuals
## Standardized Runs Statistic = 0.20913, p-value = 0.8343#Runs Test - Two sidedp-value = 0.3362. Can’t reject null hypothesis that it is random. With a p-value = 0.3362, we cannot reject the null hypothesis. Therefore, we can safely assume that residuals are not autocorrelated.
lmtest::dwtest(lmMod2)##
## Durbin-Watson test
##
## data: lmMod2
## DW = 2.0309, p-value = 0.6126
## alternative hypothesis: true autocorrelation is greater than 0#=> Durbin-Watson testWith a high p value of 0.667, we cannot reject the null hypothesis that true autocorrelation is zero. So the assumption that residuals should not be autocorrelated is satisfied by this model. If, even after adding lag1 as an X variable, does not satisfy the assumption of autocorrelation of residuals, you might want to try adding lag2, or be creative in making meaningful derived explanatory variables or interaction terms. This is more like art than an algorithm. Assumption 5 The X variables and residuals are uncorrelated How to check? Do a correlation test on the X variable and the residuals.
mod.lm <- lm(dist ~ speed, data=cars)
cor.test(cars$speed, mod.lm$residuals) # do correlation test ##
## Pearson's product-moment correlation
##
## data: cars$speed and mod.lm$residuals
## t = 5.583e-16, df = 48, p-value = 1
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2783477 0.2783477
## sample estimates:
## cor
## 8.058406e-17# Pearson's product-moment correlationp-value is high, so null hypothesis that true correlation is 0 can’t be rejected. Therefore, the assumption holds true for this model. ## Assumption 6 The number of observations must be greater than number of Xs This can be directly observed by looking at the data. ## Assumption 7 The variability in X values is positive This means the X values in a given sample must not all be the same (or even nearly the same).
var(cars$speed)## [1] 27.95918The variance in the X variable above is much larger than 0. So, this assumption is satisfied.
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:lawstat':
##
## levene.testmod2 <- lm(mpg ~ ., data=mtcars)Assumption 8
The regression model is correctly specified This means that if the Y and X variable has an inverse relationship, the model equation should be specified appropriately: Y=β1+β2∗(1/X) ## Assumption 9 No perfect multicollinearity There is no perfect linear relationship between explanatory variables. # How to check? Using Variance Inflation factor (VIF). ## What is VIF? VIF is a metric computed for every X variable that goes into a linear model. If the VIF of a variable is high, it means the information in that variable is already explained by other X variables present in the given model, which means, more redundant is that variable. So, lower the VIF (<2) the better. VIF for a X var is calculated as: VIF= 1 / (1−Rsq) where, Rsq is the Rsq term for the model with given X as response against all other Xs that went into the model as predictors. Practically, if two of the X′s have high correlation, they will likely have high VIFs. Generally, VIF for an X variable should be less than 4 in order to be accepted as not causing multi-collinearity. The cutoff is kept as low as 2, if you want to be strict about your X variables.
library(corrplot)## corrplot 0.92 loadedcorrplot(cor(mtcars[, -1]))
How to rectify? Two ways: 1. Either iteratively remove the X var with the highest VIF or, 2. See correlation between all variables and keep only one of all highly correlated pairs.
mod <- lm(mpg ~ cyl + gear, data=mtcars)
vif(mod)## cyl gear
## 1.320551 1.320551par(mfrow=c(2,2))
mod <- lm(dist ~ speed, data=cars)
plot(mod)
The convention is, the VIF should not go more than 4 for any of the X variables. That means we are not letting the RSq of any of the Xs (the model that was built with that X as a response variable and the remaining Xs are predictors) to go more than 75%. => 1/(1-0.75) => 1/0.25 => 4. ## Assumption 10 Normality of residuals The residuals should be normally distributed. If the maximum likelihood method (not OLS) is used to compute the estimates, this also implies the Y and the Xs are also normally distributed. This can be visually checked using the qqnorm() plot (top right plot).
par(mfrow=c(2,2)) # draw 4 plots in same window
mod <- lm(dist ~ speed, data=cars)
gvlma::gvlma(mod)##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Coefficients:
## (Intercept) speed
## -17.579 3.932
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma::gvlma(x = mod)
##
## Value p-value Decision
## Global Stat 15.801 0.003298 Assumptions NOT satisfied!
## Skewness 6.528 0.010621 Assumptions NOT satisfied!
## Kurtosis 1.661 0.197449 Assumptions acceptable.
## Link Function 2.329 0.126998 Assumptions acceptable.
## Heteroscedasticity 5.283 0.021530 Assumptions NOT satisfied!plot(mod) The qqnorm() plot in top-right evaluates this assumption. If points lie exactly on the line, it is perfectly normal distribution. However, some deviation is to be expected, particularly near the ends (note the upper right), but the deviations should be small, even lesser that they are here. Check Assumptions Automatically The gvlma() function from gvlma offers a way to check the important assumptions on a given linear model.
mod <- lm(dist ~ speed, data=cars[-c(23, 35, 49), ])
gvlma::gvlma(mod)##
## Call:
## lm(formula = dist ~ speed, data = cars[-c(23, 35, 49), ])
##
## Coefficients:
## (Intercept) speed
## -15.137 3.608
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma::gvlma(x = mod)
##
## Value p-value Decision
## Global Stat 7.5910 0.10776 Assumptions acceptable.
## Skewness 0.8129 0.36725 Assumptions acceptable.
## Kurtosis 0.2210 0.63831 Assumptions acceptable.
## Link Function 3.2239 0.07257 Assumptions acceptable.
## Heteroscedasticity 3.3332 0.06789 Assumptions acceptable.Three of the assumptions are not satisfied. This is probably because we have only 50 data points in the data and having even 2 or 3 outliers can impact the quality of the model. So the immediate approach to address this is to remove those outliers and re-build the model. Take a look at the diagnostic plot below to arrive at your own conclusion. From the above plot the data points: 23, 35 and 49 are marked as outliers. Lets remove them from the data and re-build the model.
Though the changes look minor, it is more closer to conforming with the assumptions. There is one more thing left to be explained. That is, the plot in the bottom right. It is the plot of standardized residuals against the leverage. Leverage is a measure of how much each data point influences the regression. The plot also contours values of Cook’s distance, which reflects how much the fitted values would change if a point was deleted. A point far from the centroid with a large residual can severely distort the regression. For a good regression model, the red smoothed line should stay close to the mid-line and no point should have a large cook’s distance (i.e. should not have too much influence on the model.)
influence.measures(mod)## Influence measures of
## lm(formula = dist ~ speed, data = cars[-c(23, 35, 49), ]) :
##
## dfb.1_ dfb.sped dffit cov.r cook.d hat inf
## 1 0.087848 -0.08003 0.08834 1.184 3.99e-03 0.1187 *
## 2 0.351238 -0.32000 0.35320 1.138 6.25e-02 0.1187 *
## 3 -0.145914 0.12652 -0.15010 1.114 1.14e-02 0.0735
## 4 0.285653 -0.24768 0.29384 1.075 4.31e-02 0.0735
## 5 0.047920 -0.04053 0.05012 1.113 1.28e-03 0.0615
## 6 -0.136783 0.11208 -0.14670 1.083 1.09e-02 0.0511
## 7 -0.047436 0.03725 -0.05287 1.089 1.43e-03 0.0422
## 8 0.081425 -0.06394 0.09076 1.083 4.19e-03 0.0422
## 9 0.212931 -0.16721 0.23734 1.031 2.80e-02 0.0422
## 10 -0.103835 0.07682 -0.12283 1.064 7.64e-03 0.0349
## 11 0.047151 -0.03488 0.05578 1.080 1.59e-03 0.0349
## 12 -0.163139 0.11031 -0.21176 1.008 2.22e-02 0.0292
## 13 -0.092988 0.06288 -0.12070 1.054 7.37e-03 0.0292
## 14 -0.047239 0.03194 -0.06132 1.071 1.92e-03 0.0292
## 15 -0.001863 0.00126 -0.00242 1.077 2.99e-06 0.0292
## 16 -0.052208 0.03031 -0.07843 1.061 3.13e-03 0.0250
## 17 0.020094 -0.01167 0.03019 1.071 4.66e-04 0.0250
## 18 0.020094 -0.01167 0.03019 1.071 4.66e-04 0.0250
## 19 0.130480 -0.07576 0.19602 1.003 1.90e-02 0.0250
## 20 -0.063700 0.02683 -0.12078 1.040 7.35e-03 0.0224
## 21 0.004170 -0.00176 0.00791 1.070 3.20e-05 0.0224
## 22 0.174775 -0.07362 0.33138 0.870 5.06e-02 0.0224
## 24 -0.087733 0.00892 -0.24379 0.948 2.86e-02 0.0213
## 25 -0.059046 0.00600 -0.16408 1.011 1.34e-02 0.0213
## 26 0.068553 -0.00697 0.19049 0.992 1.79e-02 0.0213
## 27 -0.023886 -0.02060 -0.13480 1.031 9.12e-03 0.0218
## 28 -0.005806 -0.00501 -0.03276 1.067 5.48e-04 0.0218
## 29 0.000274 -0.06235 -0.19077 1.002 1.80e-02 0.0238
## 30 0.000118 -0.02688 -0.08223 1.058 3.44e-03 0.0238
## 31 -0.000072 0.01639 0.05015 1.066 1.28e-03 0.0238
## 32 0.017941 -0.05284 -0.11169 1.054 6.32e-03 0.0274
## 33 -0.014171 0.04174 0.08822 1.062 3.95e-03 0.0274
## 34 -0.063461 0.18692 0.39507 0.848 7.09e-02 0.0274 *
## 36 0.081746 -0.16394 -0.27847 0.976 3.77e-02 0.0326
## 37 0.034106 -0.06840 -0.11619 1.062 6.84e-03 0.0326
## 38 -0.067855 0.13608 0.23115 1.007 2.64e-02 0.0326
## 39 0.182368 -0.30828 -0.45545 0.875 9.51e-02 0.0393
## 40 0.062731 -0.10604 -0.15666 1.060 1.24e-02 0.0393
## 41 0.034787 -0.05880 -0.08688 1.080 3.84e-03 0.0393
## 42 0.007121 -0.01204 -0.01778 1.088 1.62e-04 0.0393
## 43 -0.048260 0.08158 0.12053 1.071 7.37e-03 0.0393
## 44 -0.020499 0.02947 0.03716 1.108 7.06e-04 0.0573
## 45 0.200260 -0.27525 -0.33127 1.051 5.43e-02 0.0687
## 46 0.024652 -0.03277 -0.03811 1.138 7.42e-04 0.0816 *
## 47 -0.358515 0.47655 0.55420 0.979 1.46e-01 0.0816
## 48 -0.377456 0.50173 0.58348 0.964 1.60e-01 0.0816
## 50 -0.195430 0.25314 0.28687 1.118 4.14e-02 0.0961