<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/how-to-build-a-logistic-regression-model-from-scratch-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://theautomatic.net/2018/10/02/how-to-build-a-logistic-regression-model-from-scratch-in-r/"> Open Source Automation</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><img alt="how to build a logistic regression model from scratch in r" src="https://i1.wp.com/theautomatic.net/wp-content/uploads/2018/09/logistic-regression-from-scratch-in-r.jpg?w=640" data-recalc-dims="1"></p>
<p>In a <a rel="nofollow" href="http://theautomatic.net/2017/12/11/vectorize-fuzzy-matching/" target="_blank">previous post</a>, we showed how using vectorization in R can vastly speed up fuzzy matching. Here, we will show you how to use R’s vectorization functionality to efficiently build a <a rel="nofollow" href="https://en.wikipedia.org/wiki/Logistic_regression" target="_blank">logistic regression</a> model. Now we could just use the <a rel="nofollow" href="https://topepo.github.io/caret/index.html" target="_blank">caret</a> or <a rel="nofollow" href="https://www.rdocumentation.org/packages/stats/versions/3.5.1?tap_a=5644-dce66f&tap_s=10907-287229" target="_blank">stats</a> packages to create a model, but building algorithms from scratch is a great way to develop a better understanding of how they work under the hood.</p>
<p>In writing the logistic regression algorithm from scratch, we will consider the following definitions and assumptions:</p>
<p><strong>x</strong> = A dxn matrix of <em>d</em> predictor variables, where each column x<sub>i</sub> represents the vector of predictors corresponding to one data point (with <em>n</em> such columns i.e. <em>n</em> data points)</p>
<p>d = The number of predictor variables (i.e. the number of dimensions)</p>
<p>n = The number of data points</p>
<p><strong>y</strong> = A vector of labels i.e. y<sub>i</sub> equals the label associated with x<sub>i</sub>; in our case we’ll assume y = 1 or -1</p>
<p><strong>Θ</strong> = The vector of coefficients, Θ<sub>1</sub>, Θ<sub>2</sub>…Θ<sub>d</sub> trained via gradient descent. These correspond to x<sub>1</sub>, x<sub>2</sub>…x<sub>d</sub></p>
<p>α = Step size, controls the rate of gradient descent</p>
<p>Logistic regression, being a binary classification algorithm, outputs a probability between 0 and 1 of a given data point being associated with a positive label. This probability is given by the equation below:</p>
<p><img alt="logistic regression probability formula" src="https://i2.wp.com/theautomatic.net/wp-content/uploads/2018/09/final_logistic_regression_probability_formula.png?w=640" data-recalc-dims="1"></p>
<p>Recall that <strong>Θ</strong>, x> refers to the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Dot_product" target="_blank">dot product</a> of <strong>Θ</strong> and <strong>x</strong>.</p>
<p>In order to calculate the above formula, we need to know the value of <strong>Θ</strong>. Logistic regression uses a method called <a rel="nofollow" href="https://en.wikipedia.org/wiki/Gradient_descent" target="_blank">gradient descent</a> to learn the value of <strong>Θ</strong>. There are a few variations of gradient descent, but the variation we will use here is based upon the following update formula for Θ<sub>j</sub>:</p>
<p><img alt="logistic regression update formula" src="https://i1.wp.com/theautomatic.net/wp-content/uploads/2018/09/logistic-regression-formula.png?w=640" data-recalc-dims="1"></p>
<p>This formula updates the j<sup>th</sup> element of the <strong>Θ</strong> vector. Logistic regression models run this gradient descent update of <strong>Θ</strong> until either 1) a maximum number of iterations has been reached or 2) the difference between the current update of <strong>Θ</strong> and the previous value is below a set threshold. To run this update of theta, we’re going to write the following function, which we’ll break down further along in the post. This function will update the entire <strong>Θ</strong> vector in one function call i.e. <em>all j elements</em> of <strong>Θ</strong>.</p>
<pre>

# function to update theta via gradient descent
update_theta &lt;- function(theta_arg,n,x,y,d,alpha)
{
  
  # calculate numerator of the derivative
  numerator &lt;- t(replicate(d , y)) * x
  
  # perform element-wise multiplication between theta and x,
  # prior to getting their dot product
  theta_x_prod &lt;- t(replicate(n,theta_arg)) * t(x)
  dotprod &lt;- rowSums(theta_x_prod)
  
  denominator &lt;- 1 + exp(-y * dotprod)
  
  # cast the denominator as a matrix
  denominator &lt;- t(replicate(d,denominator)) 
  
  # final step, get new theta result based off update formula
  theta_arg &lt;- theta_arg - alpha * rowSums(numerator / denominator) 
  
  
  return(theta_arg)
  
}

</pre>
<h2><strong>Simplifying the update formula</strong></h2>
<p>To simplify the update formula for Θ<sub>j</sub>, we need to calculate the <a rel="nofollow" href="http://mathworld.wolfram.com/Derivative.html" target="_blank">derivative</a> in the formula above.</p>
<p>Let’s suppose z = (y<sub>i</sub>)(<strong>Θ</strong>, x<sub>i</sub>>). We’ll abbreviate the summation of 1 to n by simply using Σ.</p>
<p>Then, calculating the derivative gives us the following:</p>
<p>Σ (1 / exp(1 + z)) * exp(z) * x<sub>i</sub>y<sub>i</sub></p>
<p>= Σ (exp(z) / exp(1 + z)) * x<sub>i</sub>y<sub>i</sub></p>
<p>Since exp(z) / (1 + exp(z)) is a known identity for 1 / (1 + exp(-z)), we can substitute above to get:</p>
<p>Σ 1 / (1 + exp(-z)) * x<sub>i</sub>y<sub>i</sub></p>
<p>= Σ x<sub>i</sub>y<sub>i</sub> / (1 + exp(-z))</p>
<p>Now, substituting (y<sub>i</sub>)(<strong>Θ</strong>, x<sub>i</sub>>) back for z:</p>
<p>= Σ x<sub>i</sub>y<sub>i</sub> / (1 + exp(-(y<sub>i</sub>)(<strong>Θ</strong>, x<sub>i</sub>>)))</p>
<p>Plugging this derivative result into the rest of the update formula, the below expression tells us how to update Θ<sub>j</sub>:</p>
<p>Θ<sub>j</sub> ← Θ<sub>j</sub> – αΣ x<sub>i</sub>y<sub>i</sub> / (1 + exp(-(y<sub>i</sub>)(<strong>Θ</strong>, x<sub>i</sub>>)))</p>
<p>To convert this math into R code, we’ll split up the formula above into three main steps:</p>
<li>Calculate the numerator: x<sub>i</sub>y<sub>i</sub></li>
<li>Calculate the denominator: (1 + exp(-(y<sub>i</sub>)(<strong>Θ</strong>, x<sub>i</sub>>)))</li>
<li>Plug the results from first two steps back into the formula above</li>
<p><br>
The idea to keep in mind throughout this post is that we’re not going to use a for loop to update each j<sup>th</sup> element of <strong>Θ</strong>. Instead, we’re going to use vectorization to update the entire <strong>Θ</strong> vector at once via element-wise matrix multiplication. This will vastly speed up the gradient descent implementation.</p>
<h2><strong>Step 1) Calculating the numerator</strong></h2>
<p>In the numerator of the derivative, we have the summation of 1 to n of y<sub>i</sub> times x<sub>i</sub>. Effectively, this means we have to calculate the following:</p>
<p><strong>y<sub>1</sub></strong><strong>x<sub>1</sub></strong> + <strong>y<sub>2</sub></strong><strong>x<sub>2</sub></strong>…+<strong>y<sub>n</sub></strong><strong>x<sub>n</sub></strong></p>
<p>This calculation needs to be done for all j elements in <strong>Θ</strong> to fully update the vector. So, we actually need to run the following calculations:</p>
<p><strong>y<sub>1</sub></strong><strong>x<sub>1,1</sub></strong> + <strong>y<sub>2</sub></strong><strong>x<sub>1,2</sub></strong>…+<strong>y<sub>n</sub></strong><strong>x<sub>1,n</sub></strong></p>
<p><strong>y<sub>1</sub></strong><strong>x<sub>2,1</sub></strong> + <strong>y<sub>2</sub></strong><strong>x<sub>2,2</sub></strong>…+<strong>y<sub>n</sub></strong><strong>x<sub>2,n</sub></strong></p>
<p>…</p>
<p><strong>y<sub>1</sub></strong><strong>x<sub>d,1</sub></strong> + <strong>y<sub>2</sub></strong><strong>x<sub>d,2</sub></strong>…+<strong>y<sub>n</sub></strong><strong>x<sub>d,n</sub></strong></p>
<p>Since <strong>y</strong> is a vector, or put another way, is an nx1 matrix, and <strong>x</strong> is a dxn matrix, where d is the number of predictor variables i.e. <strong>x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub>…x<sub>d</sub></strong>, and n is the number of labels (how many predicted values we have), then we can compute the above calculations by creating a dxn matrix where each row is a duplicate of <strong>y</strong>. This way we have d duplicates of <strong>y</strong>. Each ith element of <strong>x</strong> (i.e. x<sub>i</sub>) corresponds to the i<sup>th</sup> <em>column</em> of <strong>x</strong>. So x<sub>j,i</sub> refers to the element in the j<sup>th</sup> row and i<sup>th</sup> column of <strong>x</strong>.</p>
<p><img alt="update theta numerator of logistic regression" src="https://i0.wp.com/theautomatic.net/wp-content/uploads/2018/09/update-theta-numerator-logistic-regression-1.png?w=640" data-recalc-dims="1"></p>
<p>In the above expression, instead of doing traditional matrix multiplication, we are going to do element-wise multiplication i.e. the j<sup>th</sup>, i<sup>th</sup> element of the resultant matrix will equal the j<sup>th</sup>, i<sup>th</sup> elements of each matrix multiplied by each other. This can be done in one line of R code, like below:</p>
<pre>

# calculate numerator of the derivative of the loss function
numerator &lt;- t(replicate(d , y)) * x

</pre>
<p>Initially we create a matrix where each column is equal to y<sub>1</sub>, y<sub>2</sub>,…y<sub>n</sub>. This is <strong>replicate(d, y)</strong>. We then transpose this so that we can perform the element-wise multiplication described above. This element-wise multiplication gets us the following dxn matrix result:</p>
<p><img alt="" src="https://i0.wp.com/theautomatic.net/wp-content/uploads/2018/09/xy_product.png?w=640" data-recalc-dims="1"></p>
<p>Notice how the sum of each j<sup>th</sup> row corresponds to the each calculation above i.e. the sum of row 1 is Σx<sub>i</sub>y<sub>i</sub> for j = 1, the sum of row 2 is Σx<sub>i</sub>y<sub>i</sub> for j = 2…the sum of row d is Σx<sub>i</sub>y<sub>i</sub> for j = d. Rather then using slower for loops, this allows us to calculate the numerator of the derivative given above for each element (Θ<sub>j</sub>) in <strong>Θ</strong> using vectorization. We’ll postpone doing the summation until after we’ve calculated the denominator piece of the derivative.</p>
<h2><strong>Step 2) Calculating the denominator</strong></h2>
<p>To evaluate the denominator of our formula, we need to first calculate the dot product of <strong>Θ</strong> and <strong>x</strong>. We do this similarly to our numerator calculation:</p>
<pre>

theta_x_prod &lt;- t(replicate(n,theta_arg)) * t(x)
dotprod &lt;- rowSums(theta_x_prod)

</pre>
<p>The above code corresponds to the math below:</p>
<p><strong>t(replicate(n,theta_arg))</strong></p>
<p><img alt="theta transpose logistic regression" src="https://i2.wp.com/theautomatic.net/wp-content/uploads/2018/09/theta_transpose-1.png?w=640" data-recalc-dims="1"></p>
<p><strong>t(x)</strong></p>
<p><img alt="x transpose logisitic regression" src="https://i2.wp.com/theautomatic.net/wp-content/uploads/2018/09/x_transpose_logistic_regression.png?w=640" data-recalc-dims="1"></p>
<p>Thus, <strong>theta_x_prod</strong> equals:</p>
<p><img alt="x times theta logistic regression" src="https://i1.wp.com/theautomatic.net/wp-content/uploads/2018/09/x_times_theta_logistic_regression.png?w=640" data-recalc-dims="1"></p>
<p>The sum of each i<sup>th</sup> row is, by definition, the dot product of <strong>Θ</strong> and <strong>x<sub>i</sub></strong>. Thus, by calculating the sum of each row, we can get a vector like this:</p>
<p><strong>&lt;Θ, x<sub>1</sub>>, &lt;Θ, x<sub>2</sub>>, &lt;Θ, x<sub>3</sub>>, … &lt;Θ, x<sub>n</sub>></strong></p>
<p>The calculation above is handled in R by the <strong>rowSums(theta_x_prod)</strong> code above.</p>
<p>Next, we plug the <strong>dotprod</strong> result into our formula to calculate the denominator. We then use the familiar transpose / replicate idea to create a matrix where each j<sup>th</sup> row will be used in updating Θ<sub>j</sub>.</p>
<pre>

denominator &lt;- 1 + exp(-y * dotprod)
  
# cast the denominator as a matrix
denominator &lt;- t(replicate(d,denominator)) 
  
</pre>
<h2><strong>Step 3) Finishing the update formula</strong></h2>
<p>The last line of code in the function for updating theta is finishing the rest of the formula:</p>
<pre>

# final step, get new theta result based off update formula
theta_arg &lt;- theta_arg - alpha * rowSums(numerator / denominator) 

</pre>
<h2><strong>Finally…putting it all together</strong></h2>
<p>The next function we need will essentially call our update_theta function above until either the change in the updated versus previous theta value is below a threshold, or a maximum number of iterations has been reached.</p>
<pre>


# wrapper around update_theta to iteratively update theta until the maximum number of iterations is reached
get_final_theta &lt;- function(theta,x,y,threshold,alpha=0.9,max_iter = 100)
{
  
  n &lt;- ncol(x)
  d &lt;- length(theta)
  
  first_iteration &lt;- TRUE
  total_iterations &lt;- 0

  # while the number of iterations is below the input max allowable, 
  # update theta via gradient descent
  new_theta &lt;- theta
  while(total_iterations &lt;= max_iter)
  {
    
    first_iteration &lt;- FALSE
    
    # create copy of theta for comparison between new and old
    old_theta &lt;- new_theta
    new_theta &lt;- update_theta(new_theta,n,x,y,d,alpha)

    diff_theta &lt;- sqrt(sum((new_theta - old_theta)**2))

    if(diff_theta &lt; threshold) {break}
    
    # index the iteration number
    total_iterations &lt;- total_iterations + 1
    
  }
  
  # return the train theta parameter
  return(new_theta)
  
}

</pre>
<p>Lastly, we write a simple function to calculate the predicted probability from the logistic regression model given a final theta vector:</p>
<p><img alt="logistic regression probability formula" src="https://i2.wp.com/theautomatic.net/wp-content/uploads/2018/09/final_logistic_regression_probability_formula.png?w=640" data-recalc-dims="1"></p>
<pre>

# wrapper around the gradient descent algorithm implemented in the two functions above
# ~ outputs logistic regression results
lr &lt;- function(x, y, theta_arg)
{
    result &lt;- sapply(1:ncol(x_arg) , function(index) 1 / (1 + exp(sum(theta_arg * x_arg[,index]))))
    return(result)
  
  
}


</pre>
<p>Now, if we want to use our code to train a logistic regression model, we can do it like below, using a randomized dataset as an example.</p>
<pre>

# create initial values for theta
theta &lt;- rep(0.5, 100)

# initialize other needed parameters
alpha &lt;- 0.9
threshold &lt;- 1
max_iter &lt;- 100

# generate random matrix of predictors
train_x &lt;- lapply(1:100, function(elt) sample(1:10000, 100, replace = TRUE))
train_x &lt;- Reduce(cbind, x)

# set seed for reproducibility
set.seed(2000)

# generate random labels
train_y &lt;- sample(c(1, -1), 100, replace = TRUE)

# get trained theta
train_theta_result &lt;- get_final_theta(theta, train_x , y, 1, .9, max_iterations = 100)

# run trained model on dataset
train_result_predictions &lt;- lr(train_x, train_theta_result)


</pre>
<p>Then, if we want to use the trained logistic regression model on a new test data, we can use the <strong>train_theta_result</strong> vector above.</p>
<pre>

# set seed for reproducibility
set.seed(1234)

# generate randomized test dataset
test_x &lt;- lapply(1:100, function(elt) sample(1:10000, 100, replace = TRUE))
test_x &lt;- Reduce(cbind, x)

# get predictions on test dataset
test_result_predictions &lt;- lr(test_x, train_theta_result)

</pre>
<p>That’s it for this post! Please check out other R posts here: <a rel="nofollow" href="http://theautomatic.net/category/r/" target="_blank">http://theautomatic.net/category/r/</a></p>
<p>The post <a rel="nofollow" href="http://theautomatic.net/2018/10/02/how-to-build-a-logistic-regression-model-from-scratch-in-r/" target="_blank">How to build a logistic regression model from scratch in R</a> appeared first on <a rel="nofollow" href="http://theautomatic.net/" target="_blank">Open Source Automation</a>.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://theautomatic.net/2018/10/02/how-to-build-a-logistic-regression-model-from-scratch-in-r/"> Open Source Automation</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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/KQygKcPfs3g" height="1"></a></div>