<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/linear-regression-from-scratch-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://datascienceplus.com/linear-regression-from-scratch-in-r/"> DataScience+</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>One of the very first learning algorithms that you’ll encounter when studying data science and machine learning is least squares linear regression. Linear regression is one of the easiest learning algorithms to understand; it’s suitable for a wide array of problems, and is already implemented in many programming languages. Most users are familiar with the <code>lm()</code> function in R, which allows us to perform linear regression quickly and easily. But one drawback to the <code>lm()</code> function is that it takes care of the computations to obtain parameter estimates (and many diagnostic statistics, as well) on its own, leaving the user out of the equation. So, how does the <code>lm()</code> function obtain these parameter estimates?</p>
<p>In this post, I will outline the process from first principles in R. I will use only matrices, vectors, and matrix operations to obtain parameter estimates using the closed-form linear algebraic solution. After reading this post, you’ll see that it’s actually quite simple, and you’ll be able to replicate the process with your own data sets (though using the <code>lm()</code> function is of course much more efficient).</p>
<p>For this example, I’ll use the well-known “Boston” data set from the MASS library. For those who aren’t familiar with it, the Boston data set contains 14 economic, geographic, and demographic variables for 506 tracts in the city of Boston from the early 1990s. Our target variable will be the median home value for each tract — the column titled ‘medv.’ Let’s get an idea of what this data set looks like:</p>
<pre>library(MASS)
str(Boston)<em>
## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...</em></pre>
<p>To learn more about the definition of each variable, type <code>help(Boston)</code> into your R console.</p>
<p>Now we’re ready to start. Linear regression typically takes the form<br>
$$y=\beta X + \epsilon$$where ‘y’ is a vector of the response variable, ‘X’ is the matrix of our feature variables (sometimes called the ‘design’ matrix), and β is a vector of parameters that we want to estimate. \(\epsilon\) is the error term; it represents features that affect the response, but are not explicitly included in our model.<br>
The first thing we will need is a vector of our response variable, typically called ‘y’. In this case, our response variable is the median home value in thousands of US dollars: ‘medv’:</p>
<pre>y &lt;- Boston$medv</pre>
<p>Next, we’ll need our ‘X’ or ‘design’ matrix of the input features. We can do this with the <code>as.matrix()</code> function, grabbing all the variables except ‘medv’ from the Boston data frame. We’ll also need to add a column of ones to our X matrix for the intercept term.</p>
<pre># Matrix of feature variables from Boston
X &lt;- as.matrix(Boston[-ncol(Boston)])

# vector of ones with same length as rows in Boston
int &lt;- rep(1, length(y))

# Add intercept column to X
X &lt;- cbind(int, X)
</pre>
<p>Now that we have our response vector and our ‘X’ matrix, we can use them to solve for the parameters using the following closed form solution:<br>
$$\beta = (X^TX)^{-1}X^Ty$$The derivation of this equation is beyond the scope of this post. If you are interested in the derivation, a good place to start is <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Linear_regression">Wikipedia</a> or any good undergraduate textbook on linear regression.</p>
<p>We can implement this in R using our ‘X’ matrix and ‘y’ vector. The <code>%*%</code>  operator is simply matrix multiplication. The <code>t()</code> function takes the transpose of a matrix, and <code>solve()</code> calculates the inverse of any (invertible) matrix.</p>
<pre># Implement closed-form solution
betas &lt;- solve(t(X) %*% X) %*% t(X) %*% y

# Round for easier viewing
betas &lt;- round(betas, 2)
</pre>
<p>Now we have our parameter vector <em>β</em>. These are the linear relationships between the response variable ‘medv’ and each of the features. For example, we see that an increase of one unit in the number of rooms (rm) is associated with a $3,810 increase in home value. We can find the relationship between the response and any of the feature variables in the same manner, paying careful attention to the units in the data description. Notice that one of our features, ‘chas’, is a dummy variable which takes a value of 0 or 1 depending on whether or not the tract is adjacent to the Charles River. The coefficient of ‘chas’ tells us that homes in tracts adjacent to the Charles River (coded as 1) have a median price that is $2,690 higher than homes in tracts that do not border the river (coded as 0) when the other variables are held constant.</p>
<pre>print(betas)<em>
## [,1]
## int 36.46
## crim -0.11
## zn 0.05
## indus 0.02
## chas 2.69
## nox -17.77
## rm 3.81
## age 0.00
## dis -1.48
## rad 0.31
## tax -0.01
## ptratio -0.95
## black 0.01
## lstat -0.52</em></pre>
<p>Now let’s verify that these results are in fact correct. We want to compare our results to those produced by the <code>lm()</code>  function. Most of you will already know how to do this.</p>
<pre># Linear regression model
lm.mod &lt;- lm(medv ~ ., data=Boston)

# Round for easier viewing
lm.betas &lt;- round(lm.mod$coefficients, 2)

# Create data.frame of results
results &lt;- data.frame(our.results=betas, lm.results=lm.betas)

print(results)
<em>
## our.results lm.results
## int 36.46 36.46
## crim -0.11 -0.11
## zn 0.05 0.05
## indus 0.02 0.02
## chas 2.69 2.69
## nox -17.77 -17.77
## rm 3.81 3.81
## age 0.00 0.00
## dis -1.48 -1.48
## rad 0.31 0.31
## tax -0.01 -0.01
## ptratio -0.95 -0.95
## black 0.01 0.01
## lstat -0.52 -0.52</em></pre>
<p>The parameter estimates are identical! Now you know how these estimates are obtained. Keep in mind that you’re unlikely to favor implementing linear regression in this way over using <code>lm()</code> . The <code>lm()</code>  function is very quick, and requires very little code. Using it provides us with a number of diagnostic statistics, including \(R^2\), t-statistics, and the oft-maligned p-values, among others. Nevertheless, I wanted to provide a window into what <code>lm()</code>  does behind the scenes.</p>
<p>There are other ways to solve for the parameters in linear regression. For example, gradient descent can be used to obtain parameter estimates when the number of features is extremely large, a situation that can drastically slow solution time when using the closed-form method. Stochastic gradient descent is often used when both the number of features and the number of samples are very large. I hope to detail these in a future post.</p>
<p>I have made the code from this post available at my Github <a target="_blank" rel="nofollow" href="https://github.com/capt-calculator/linear-regression-from-scratch-r">here</a>.<br>
Cheers!</p>
<ul></ul>
<p><strong>
<p>Related Post</p>
<p>
<ol>
<li><a title="Permanent Link toR for Publication by Page Piccinini: Lesson 5 – Analysis of Variance (ANOVA)" target="_blank" rel="nofollow" href="https://datascienceplus.com/r-for-publication-by-page-piccinini-lesson-5-analysis-of-variance-anova/">R for Publication by Page Piccinini: Lesson 5 – Analysis of Variance (ANOVA)</a></li>
<li><a title="Permanent Link toR for Publication by Page Piccinini: Lesson 4 – Multiple Regression" target="_blank" rel="nofollow" href="https://datascienceplus.com/r-for-publication-by-page-piccinini-lesson-4-multiple-regression/">R for Publication by Page Piccinini: Lesson 4 – Multiple Regression</a></li>
<li><a title="Permanent Link toR for Publication by Page Piccinini: Lesson 3 – Logistic Regression" target="_blank" rel="nofollow" href="https://datascienceplus.com/r-for-publication-by-page-piccinini-lesson-3-logistic-regression/">R for Publication by Page Piccinini: Lesson 3 – Logistic Regression</a></li>
<li><a title="Permanent Link toR for Publication by Page Piccinini: Lesson 2 – Linear Regression" target="_blank" rel="nofollow" href="https://datascienceplus.com/r-for-publication-by-page-piccinini-lesson-2-linear-regression/">R for Publication by Page Piccinini: Lesson 2 – Linear Regression</a></li>
<li><a title="Permanent Link toHow to detect heteroscedasticity and rectify it?" target="_blank" rel="nofollow" href="https://datascienceplus.com/how-to-detect-heteroscedasticity-and-rectify-it/">How to detect heteroscedasticity and rectify it?</a></li>
</ol></strong>
<div>
<a target="_blank" rel="nofollow" href="http://feeds.feedburner.com/~ff/datascienceplus?a=bPPvIYXKMz4:Qzp4lTo3hjo:D7DqB2pKExk"><img src="http://feeds.feedburner.com/~ff/datascienceplus?i=bPPvIYXKMz4:Qzp4lTo3hjo:D7DqB2pKExk" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.feedburner.com/~ff/datascienceplus?a=bPPvIYXKMz4:Qzp4lTo3hjo:qj6IDK7rITs"><img src="http://feeds.feedburner.com/~ff/datascienceplus?d=qj6IDK7rITs" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.feedburner.com/~ff/datascienceplus?a=bPPvIYXKMz4:Qzp4lTo3hjo:yIl2AUoC8zA"><img src="http://feeds.feedburner.com/~ff/datascienceplus?d=yIl2AUoC8zA" border="0"></a>
</div>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://datascienceplus.com/linear-regression-from-scratch-in-r/"> DataScience+</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/eHh-sDhMuWw"></a></div>