1. Introduction

In this exercise, you will implement logistic regression and apply it to two different datasets. To get started with the exercise, you will need to download the starter code and unzip its contents to the directory where you wish to complete the exercise. You will be required to modify several parts of the code, by following the instructions in this assignment. Files needed for this practical include:


2. Logistic regression

In this part of the exercise, you will build a logistic regression model to predict whether a blood sample comes from a patient with malaria. Suppose that you are involved in a research project that aims to develop a novel biochemical test for the early diagnosis of the disease based on the levels of two different biomarkers. You have data from previous patients that you can use as a training set for logistic regression. For each training example, you have the levels of both biomarkers and the results of the biopsy that determines the presence of the parasite Plasmodium falciparum. Your task is to build a classiffication model that estimates a patient’s probability of having malaria based the levels of those two biomarkers.


2.1. Visualizing the data

Before starting to implement any learning algorithm, it is always good to visualize the data if possible. Here the axes are the two biomarker levels, and the positive (y = 1, malaria) and negative (y = 0, no) examples are shown with different markers.

data1 <- read.table("data1.txt", sep = ",")  # read comma-separated data

X <- as.matrix(data1[, 1:2])
y <- data1[, 3]

plot(X, xlab = "Biomarker 1", ylab = "Biomarker 2",
     col = ifelse(y == 1, "gold", "black"), 
     pch = ifelse(y == 1, 3, 1))
legend("topright", c("Malaria", "No"), 
       col = c("gold", "black"),
       pch = c(3, 1))


2.2. Implementation

2.2.1. Sigmoid function

Before you start with the actual cost function, recall that the logistic regression hypothesis is defined as:

\[ h_\theta(x) = g(\theta^TX); \] where function \(g\) is the sigmoid function. The sigmoid function is defined as:

\[ g(z) = \frac{1}{1+e^{-z}} \]

Your first task is to implement this function in sigmoid.R so it can be called by the rest of your program. When you are finished, try testing a few values by calling sigmoid(x) at the R command line. For large positive values of \(x\), the sigmoid should be close to 1, while for large negative values, the sigmoid should be close to 0. Evaluating sigmoid(0) should give you exactly 0.5. Your code should also work with vectors and matrices. For a matrix, your function should perform the sigmoid function on every element.


2.2.2. Compute cost and gradient

Now you will implement the cost function and gradient for logistic regression. Complete the code in costFunction.R and gradient.R to return the cost and gradient. Recall that the cost function in logistic regression is:

\[ J(\theta) = \frac{1}{2m}\sum\limits_{i=1}^m[-y^{(i)}log(h_{\theta}(x^{(i)})) - (1-y^{(i)})log(1 - h_{\theta}(x^{(i)}))] \]

and the gradient of the cost is a vector of the same length as \(\theta\) where the \(j^{th}\) element (for \(j = 0,1,...,n\)) is defined as follows:

\[ \frac{\partial J(\theta)}{\partial \theta_j} = \frac{1}{m}\sum\limits_{i=1}^m(h_{\theta}(x^{(i)}) - y^{(i)})x_j^{(i)} \]

Note that while this gradient looks identical to the linear regression gradient, the formula is actually different because linear and logistic regression have different definitions of \(h_\theta(x)\). Once you are done, the following code will call your costFunction using the initial \(\theta\) parameters. You should see that the cost is about 0.346 and grad[,2] about -12.

# Initialize some useful values 
m <- nrow(X) # number of training examples
n <- ncol(X) # number of features

# Add intercept term to X
X <- cbind(rep(1, m), X)
colnames(X) <- NULL

# Initialize fitting parameters
initial_theta <- rep(0, n+1)

# Compute and display initial cost and gradient 
source("sigmoid.R")
source("costFunction.R")
source("gradient.R")
cost <- costFunction(initial_theta, X, y)
grad <- gradient(initial_theta, X, y)


2.2.3. Learning parameters using optim

In the previous assignment, you found the optimal parameters of a linear regression model by implementing gradent descent. You wrote a cost function and calculated its gradient, then took a gradient descent step accordingly. This time, instead of taking gradient descent steps, you will use an R built-in function called optim.

R optim is an optimization solver that finds the minimum of an unconstrained function. Constraints in optimization often refer to constraints on the parameters, for example, constraints that bound the possible values \(\theta\) can take (e.g., \(\theta \le 1\)). Logistic regression does not have such constraints since \(\theta\) is allowed to take any real value. For logistic regression, you want to optimize the cost function \(J(\theta)\) with parameters \(\theta\). Particularly, you are going to use optim to find the best parameters \(\theta\) for the logistic regression cost function, given a fixed dataset (of \(X\) and \(y\) values). You will pass to optim the following inputs:

  • The initial values of the parameters we are trying to optimize.
  • Two functions that, when given the training set and a particular \(\theta\), compute i) the logistic regression cost and ii) the gradient with respect to \(\theta\) for the dataset (X, y)

We already have code written to call optim with the correct arguments below:

opt <- optim(initial_theta, fn = costFunction, gr = gradient, X, y,
             method = "BFGS", control = list(maxit = 1500))
theta <- opt$par
J <- opt$value

If you have completed the costFunction and gradient functions correctly, optim will converge on the right optimization parameters and return the final values of the cost and \(\theta\). Notice that by using optim, you did not have to write any loops yourself, or set a learning rate like you did for gradient descent. This is all done by optim; you only needed to provide functions calculating the cost and the gradient. Once optim completes, you should see that the cost is about 0.203. The final \(\theta\) value will then be used to plot the decision boundary on the training data. We also encourage you to look at the code in plotDecisionBoundary.R to see how to plot such a boundary using the \(\theta\) values.

source("plotDecisionBoundary.R")
plotDecisionBoundary(theta, X, y)

2.2.4. Evaluating logistic regression

After learning the parameters, you can use the model to predict whether a new patient will have malaria. For a patient with levels of biomarker 1 of 45 and an levels of biomarker 2 of 85, you should expect to see a probability of having the disease of 0.776.

sigmoid(t(theta) %*% c(1, 45, 85))

Another way to evaluate the quality of the parameters we have found is to see how well the learned model predicts on our training set. In this part, your task is to complete the code in predict.R. The predict function will produce “1” or “0” predictions given a dataset and a learnt parameter vector \(\theta\). After you have completed the code in predict.R, the following code will proceed to report the training accuracy of your classiffer by computing the percentage of examples it got correct. You should get an accouracy of 89%.

source("predict.R")
p <- predict(theta, X)

# Compute accuracy on our training set
acc <- mean(p == y)


3. Regularized logistic regression

In this part of the exercise, you will implement regularized logistic regression to predict whether a dog is affected by leishmaniosis. This parasitary disease, caused by a trypanosome called Leishmania can be detected using the results of two comercially available tests, which currently are not 100% trustable. Suppose you work in a veterninary clinic and you want to determine if a dog is affected by the disease and therefore should be treated. To help you make the decision, you have a dataset of test results and serologic diagnosis on past analyses, from which you can build a logistic regression model.


3.1. Visualizing the data

Here the axes are the two test scores, and the positive (y = 1, affected) and negative (y = 0, no) examples are shown with different markers.

data2 <- read.table("data2.txt", sep = ",")  # read comma-separated data

X <- as.matrix(data2[, 1:2])
y <- data2[, 3]

  plot(X, xlab = "Test 1", ylab = "Test 2",
     col = ifelse(y == 1, "gold", "black"), 
     pch = ifelse(y == 1, 3, 1))
legend("topright", c("Affected", "No"), 
       col = c("gold", "black"),
       pch = c(3, 1))

The plot shows that our dataset cannot be separated into positive and negative examples by a straight-line through the plot. Therefore, a straightforward application of logistic regression will not perform well on this dataset since logistic regression will only be able to find a linear decision boundary.


3.2. Feature mapping

One way to fit the data better is to create more features from each data point. In the provided function mapFeature.R, we will map the features into all polynomial terms of \(x_1\) and \(x_2\) up to the sixth power:

\[ mapFeature(x) = \begin{bmatrix} 1 \\ x_1 \\ x_2 \\ x_1^2 \\ x_1x_2 \\ x_2^2 \\ x_1^3 \\ \vdots \\ x_1^5x_2^5 \\ x_2^6 \end{bmatrix} \]

As a result of this mapping, our vector of two features (the scores on two diagnostic tests) has been transformed into a 28-dimensional vector. A logistic regression classifier trained on this higher-dimension feature vector will have a more complex decision boundary and will appear nonlinear when drawn in our 2-dimensional plot. While the feature mapping allows us to build a more expressive classifier, it also more susceptible to overfitting. In the next parts of the exercise, you will implement regularized logistic regression to fit the data and also see for yourself how regularization can help combat the overfitting problem.

# Add Polynomial Features
source("mapFeature.R")
X <- mapFeature(X[, 1], X[, 2])

m <- nrow(X) # no. of training examples
n <- ncol(X) - 1 # no. of features 

# Note that mapFeature also adds a column of ones for us, so the intercept
# term is already handled


3.3. Cost function and gradient

Now you will implement code to compute the cost function and gradient for regularized logistic regression. Complete the code in costFunctionReg.R and gradientReg.Rto return the cost and gradient, respectively. Recall that the regularized cost function in logistic regression is:

\[ J(\theta) = \frac{1}{2m}\sum\limits_{i=1}^m[-y^{(i)}log(h_{\theta}(x^{(i)})) - (1-y^{(i)})log(1 - h_{\theta}(x^{(i)}))] + \frac{\lambda}{2m}\sum\limits_{i=1}^n\theta_j^2 \]

Note that you should not regularize the parameter \(\theta_0\). In R, recall that indexing starts from 1, hence, you should not be regularizing the theta[1] parameter (which corresponds to \(\theta_0\)) in the code. The gradient of the cost function is a vector where the \(j^{th}\) element is defined as follows:

  • for \(j = 0\) \[ \frac{\partial J(\theta)}{\partial \theta_0} = \frac{1}{m}\sum\limits_{i=1}^m(h_{\theta}(x^{(i)}) - y^{(i)})x_j^{(i)} \]

  • for \(j \ge 1\)

\[ \frac{\partial J(\theta)}{\partial \theta_j} = \bigg( \frac{1}{m}\sum\limits_{i=1}^m(h_{\theta}(x^{(i)}) - y^{(i)})x_j^{(i)}\bigg) + \frac{\lambda}{m}\theta_j \]

Once you are done, the following code will call your costFunctionReg function using the initial value of \(\theta\) (initialized to all zeros). You should see that the cost is about 0.346.

# Initialize fitting parameters
initial_theta <- rep(0, n+1)

# Set regularization parameter lambda to 1
lambda <- 1

# Compute and display initial cost and gradient for regularized logistic
# regression
source("costFunctionReg.R")
source("gradientReg.R")
cost <- costFunctionReg(initial_theta, X, y, lambda)
grad <- gradientReg(initial_theta, X, y, lambda)


3.3.1. Learn parameters using optim

Similar to the previous part, you will use optim to learn the optimal parameters \(\theta\). If you have completed the cost and gradient for regularized logistic regression correctly, you should be able to step through the next code to learn the parameters.

opt <- optim(initial_theta, fn = costFunctionReg, gr = gradientReg, X, y, lambda, 
             method = "BFGS", control = list(maxit = 1500))
theta <- opt$par
J <- opt$value


3.4. Plotting the decision boundary

To help you visualize the model learned by this classifier, we have provided the function plotDecisionBoundary.R which plots the (nonlinear) decision boundary that separates the positive and negative examples. In plotDecisionBoundary.R, we plot the nonlinear decision boundary by computing the classifier’s predictions on an evenly spaced grid and then and draw a contour plot of where the predictions change from y = 0 to y = 1.

source("plotDecisionBoundary.R")
plotDecisionBoundary(theta, X, y)

Finally, try the following values of \(\lambda\): 0, 1, 10, 100. How does the decision boundary change when you vary \(\lambda\)? How does the training set accuracy vary? With a small \(\lambda\), you should find that the classifier gets almost every training example correct, but draws a very complicated boundary, thus overfitting the data. With a larger \(\lambda\), you should see a plot that shows an simpler decision boundary which still separates the positives and negatives fairly well. However, if \(\lambda\) is set to a too high a value, you will not get a good fit and the decision boundary will not follow the data so well, thus underfitting the data.