1. Introduction

In this practical, you will solve regression and classification problems using real datasets. During the hands-on, you will find a series of short questions (Q) and tasks (T) that you will have to complete. Remember that they can appear as possible questions in the weekly questionaires.

Required libraries

Before starting, install and load the required libraries:

library(ggplot2)
library(GGally)
library(caret)
library(MASS)
library(pROC)

2. Datasets

Prostate cancer data

This dataset comes from a study that examined the correlation between the level of prostate specific antigen (PSA) and a number of clinical measurements in men who were about to receive a radical prostatectomy (Stamey et al., 1989). It contains data from 97 patients and 9 clinical measurements: log of cancer volume (lcavol), log prostate weight (lweight), age, log of benign prostatic hyperplasia amount (bph), seminal vesicle invasion (svi), log capsular penetration (lcp), Gleason score (gleason), percent of Gleason scores 4 or 5 (pgg45) and log of PSA (lpsa).

Diabetes data

This dataset is a corrected version of the original one from the National Institute of Diabetes and Digestive and Kidney Diseases. It contains data from 768 female patients and 9 clinical variables: number of times pregnant (pregnant), plasma glucose concentration (glucose), diastolic blood pressure (pressure), triceps skin fold thickness (thickness), 2-Hour serum insulin (insulin), body mass index (mass), diabetes pedigree function (pedigree), age (age) and diabetes diagnosis (diabetes)

Diabetes multi-class data

This dataset (Reaven and Miller, 1979) presents three different groups of patients (1: overt diabetes,2: chemical diabetes, 3: normal) and 5 clinical variables measured on them: relative weight (relwt), fasting plasma glucose (glufast), area under the plasma glucose curve for the oral glucose tolerance test (glutest), area under the plasma insulin curve for the oral glucose tolerance test (insulin) and steady state plasma glucose response (steady). It contains data from 144 patients.

3. Regression: predict prostate specific antigen levels

Load the prostate dataset and inspect its content:

prostate <- read.delim("prostate.tsv", header = TRUE)  # read tab-separated file

Q1: In this dataset, which is the value of n and p? Which R functions did you use to calculate it?

When you analyze a new dataset, it is often useful to visualize it:

ggpairs(data = prostate, progress = FALSE)

As you can see, some correlations between the predictors and log(PSA) are large.

Before training a linear regression model, you need to split the data in a training and a test set. You will assign 70% of your data to the training set.

set.seed(123) # To ensure we can reproduce the result below

training.examples <- sample(1:nrow(prostate), round(nrow(prostate)*0.7)) # random selection of indices 

training.set <- prostate[training.examples, ] 
test.set <- prostate[-training.examples, ]

This step can be also done very easily using the function createDataPartition in the caret package:

set.seed(123)
training.examples <- createDataPartition(prostate$lpsa, p = 0.7, list = FALSE)

training.set <- prostate[training.examples, ]
test.set <- prostate[-training.examples, ]

Q2: Why do you need to split the data into training and test sets?

Now, fit a simple linear regression model using lpsa as response variable and lcavol as predictor, using the lm function:

model <- lm(lpsa ~ lcavol, data = training.set)
model
## 
## Call:
## lm(formula = lpsa ~ lcavol, data = training.set)
## 
## Coefficients:
## (Intercept)       lcavol  
##      1.5909       0.6581

Our hypothesis function, \(h(x)\), can be written as follows: \(h(lcavol) = 1.5909 + 0.6581lcavol\). Using this formula, for each value of lcavol, you can obtain a prediction of lpsa.

pred <- predict(model)

T1: Obtain the hypothesis function corresponding to regress lpsa on svi.

In this case it is straightforward to represent your data:

ggplot(training.set, aes(x = lcavol, y = lpsa)) +
  geom_point() +
  geom_smooth(formula = y ~ x, method = "lm", se = FALSE) +
  geom_segment(aes(xend = lcavol, yend = pred), lty = 2) +
  labs(x = "log(cancer volume)", y = "log(PSA)") +
  theme_bw() # Nicer plot

Q3: What do the blue line and the black dotted lines represent?

However, note that in most real-case scenarios you will not be able to visualize your model in 2D.

Now, update your model to include the rest of predictors:

model <- lm (lpsa ~ ., data = training.set)
# The formula specifies the equation y ~ x
# To include all the variables in the training set different from y, we can simply use the dot (.)

model
## 
## Call:
## lm(formula = lpsa ~ ., data = training.set)
## 
## Coefficients:
## (Intercept)       lcavol      lweight          age         lbph          svi  
##    1.146141     0.532557     0.534911    -0.024256     0.091608     0.755856  
##         lcp      gleason        pgg45  
##   -0.132859    -0.013903     0.005646

For a given predictor, the coefficient (\(\beta\)) can be interpreted as the average effect on the predicted response of an increase of one unit in the predictor, holding all other predictors fixed. For example, an increase of one unit in the log(cancer volume) would lead (holding the other predictors fixed) to an average increase in the predicted log(PSA) of approximately 0.53.

Q4: How much would change the the log(PSA) if the age increases in 2 years (holding other predictors fixed)?

If the assumptions of the linear regression model hold, for a given the predictor, we can use the t-statistic to evaluate whether or not there is significant association between the predictor and the outcome variable, that is whether the \(\beta\) coefficient of the predictor is significantly different from zero:

summary(model)$coef
##                 Estimate Std. Error     t value     Pr(>|t|)
## (Intercept)  1.146140971 1.52184733  0.75312481 4.542718e-01
## lcavol       0.532557052 0.11305978  4.71040230 1.472033e-05
## lweight      0.534911103 0.22657126  2.36089567 2.144233e-02
## age         -0.024255613 0.01320648 -1.83664414 7.113589e-02
## lbph         0.091608270 0.06952440  1.31764194 1.925510e-01
## svi          0.755856230 0.30049893  2.51533747 1.454194e-02
## lcp         -0.132858742 0.11802396 -1.12569292 2.647050e-01
## gleason     -0.013902816 0.18484239 -0.07521443 9.402902e-01
## pgg45        0.005646487 0.00508524  1.11036783 2.711984e-01

From the output above, the coefficients table shows the \(\beta\) coefficient estimates and their significance levels. Columns are:

Although here we are mostly interested in making predictions, you can review the model interpretation here.

Let’s use the model that we have trained to predict the log(PSA) in the test set dataset:

pred.test <- predict(model, newdata = test.set)

T2: Check out ?predict. What is its output in this case? What changed with respect to when we predicted using the training set above?

Evaluate the predicition using the \(MSE_{test} = \frac{1}{n_{test}}\sum\limits_{i=1}^{n_{test}}(h(x_{test}^{(i)})-y_{test}^{(i)})^2\):

MSE <- sum((pred.test - test.set$lpsa)^2)/nrow(test.set)
MSE
## [1] 0.4815295
# Other commonly used metrics include the RMSE
RMSE(pred = pred.test, obs = test.set$lpsa) # RMSE = sqrt(MSE) 
## [1] 0.6939233
# Or the R2, which represents the proportion of variation in the response variable that can be explained by predictor variables. It ranges from 0 to 1.
R2(pred = pred.test, obs = test.set$lpsa) 
## [1] 0.7772803

T3: Compute the RMSE in the training set. You expect it to be larger or smaller than in the test set? Why?

4. Classification: classify diabetes patients

Load the diabetes dataset and inspect its content. Plot the data.

diabetes <- read.csv("diabetes.csv", header = TRUE)  # read comma-separated file
ggpairs(data = diabetes, progress = FALSE)

Q5: How would you identify the variables that change the most between the two diabetes groups (pos and neg) by looking at this plot?

Split the data randomly into training and test sets:

set.seed(123)
training.examples <- createDataPartition(diabetes$diabetes, p = 0.8, list = FALSE)

training.set <- diabetes[training.examples, ]
test.set <- diabetes[-training.examples, ]

Use simple logistic regression to predict the probability of class membership (diabetes pos or neg) from the glucose predictor variable:

model <- glm(diabetes ~ glucose, data = training.set, family = "binomial")

model
## 
## Call:  glm(formula = diabetes ~ glucose, family = "binomial", data = training.set)
## 
## Coefficients:
## (Intercept)      glucose  
##    -6.15882      0.04327  
## 
## Degrees of Freedom: 313 Total (i.e. Null);  312 Residual
## Null Deviance:       398.8 
## Residual Deviance: 305.7     AIC: 309.7

The output above shows the estimate of the regression \(\beta\) coefficients and their significance levels. The intercept (\(\beta_0\)) is -6.16 and the coefficient of glucose variable is 0.043.

The logistic equation can be written as \(h(glucose) = \frac{1}{1 + e^{-(-6.16 + 0.043glucose)}}\)

Using this formula, for each new glucose plasma concentration value, you can predict the probability of the individuals being diabetes-positive.

Q6: Which is the equation of the decision boundary?

Predictions can be easily obtained using the predict function. Select the option type = "response" to directly obtain the probabilities.

probs <- predict(model, type = "response")

In this case it is straightforward to represent your data:

ggplot(data = training.set, aes(glucose, ifelse(diabetes == "pos", 1, 0))) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "glm", formula = y ~ x, method.args = list(family = "binomial"), se = F) +
  labs(
    x = "Plasma Glucose Concentration",
    y = "Probability of having diabetes"
    ) +
  theme_bw()

Q7: Approximately, by looking at the plot, which value of glucose corresponds to the threshold above which a patient is classified as diabetes-positive?

Let’s include now all the predictors in the model:

model <- glm(diabetes ~., data = training.set, family = "binomial")
summary(model)$coef
##                  Estimate  Std. Error     z value     Pr(>|z|)
## (Intercept) -1.053400e+01 1.439679266 -7.31690975 2.537464e-13
## pregnant     1.005031e-01 0.061266974  1.64041157 1.009196e-01
## glucose      3.709621e-02 0.006486093  5.71934633 1.069346e-08
## pressure    -3.875933e-04 0.013826185 -0.02803328 9.776356e-01
## triceps      1.417771e-02 0.019981885  0.70952823 4.779967e-01
## insulin      5.939876e-04 0.001508231  0.39383055 6.937061e-01
## mass         7.997447e-02 0.031798907  2.51500698 1.190300e-02
## pedigree     1.329149e+00 0.482291020  2.75590704 5.852963e-03
## age          2.718224e-02 0.020199295  1.34570257 1.783985e-01

From the output above, the coefficients table shows the \(\beta\) coefficient estimates and their significance levels, analogously to the case of linear regression.The coefficient estimate for the variable glucose is 0.037, which is positive. This means that an increase in the glucose levels is associated with increase in the probability of having diabetes. However the coefficient for the variable pressure is negative. This means that an increase in diastolic blood pressure is associated with a decreased probability of having diabetes.

Make predictions in the test dataset to evaluate the performance of your logistic regression model. First, obtain the probabilities of being diabetes-positive using the trained logistic regression model. Then, assign the individuals to with probability \(\ge\) 0.5 to the diabetes-positive class.

probs <- predict(model, newdata = test.set, type = "response")
pred <- as.factor(ifelse(probs >= 0.5, "pos", "neg"))

# Obtain confusion matrix for the test set
confusionMatrix(data = pred, reference = test.set$diabetes, positive = "pos")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction neg pos
##        neg  44  11
##        pos   8  15
##                                          
##                Accuracy : 0.7564         
##                  95% CI : (0.646, 0.8465)
##     No Information Rate : 0.6667         
##     P-Value [Acc > NIR] : 0.05651        
##                                          
##                   Kappa : 0.4356         
##                                          
##  Mcnemar's Test P-Value : 0.64636        
##                                          
##             Sensitivity : 0.5769         
##             Specificity : 0.8462         
##          Pos Pred Value : 0.6522         
##          Neg Pred Value : 0.8000         
##              Prevalence : 0.3333         
##          Detection Rate : 0.1923         
##    Detection Prevalence : 0.2949         
##       Balanced Accuracy : 0.7115         
##                                          
##        'Positive' Class : pos            
##

The accuracy of our classification model is 0.75.

Q8: What does this mean? And having a sensitivity of 0.58 and a specificity of 0.84? Hint: Check out ?caret::confusionMatrix

T4: Compute the missclassification rate.

Lets plot the corresponding ROC curve:

roc <- roc(test.set$diabetes, probs, smooth = TRUE)
df <- data.frame(sens = roc$sensitivities, spec = roc$specificities)
ggplot(data = df, aes(1-spec, sens))+
  geom_line() +
  geom_abline(intercept = 0, slope = 1, col = "red", lty = 2) +
  theme_bw()

# AUC
auc(roc)
## Area under the curve: 0.8401

Q9: What does the red dotted line represent?

T5: Compute the ROC AUC in the training set. You expect it to be larger or smaller than in the test set? Why?

Multi-class classification

Often it is common to have more than two classes for the response variable. For example, in our second diabetes dataset, we have three different classes (1: overt diabetes,2: chemical diabetes, 3: normal). In this case, we can make predictions using the one vs all strategy:

# Load and inspect the data
diabetes2 <- read.csv("diabetes2.csv", header = TRUE)  # read comma-separated file

# Remove the 'id' column (it is just an individual identifier) 
diabetes2$id <- NULL

# Split the data into training and test sets
set.seed(123)
training.examples <- createDataPartition(diabetes2$group, p = 0.7, list = FALSE)
training.set  <- diabetes2[training.examples, ]
test.set <- diabetes2[-training.examples, ]
    
# Empty vector to store the results
res <- c()

# Use a for loop to train 3 logistic regression models
for (i in 1:3){
    
    # Make a copy of the training set
    ts <- training.set
    Ts <- test.set
  
    # Convert levels != i to 0 and levels == i to 1
    ts$group[ts$group != i] <- 0
    ts$group[ts$group == i] <- 1
    Ts$group[Ts$group != i] <- 0
    Ts$group[Ts$group == i] <- 1

    # For simplicity, let's build a model using only one predictor: insulin
    model <- glm(group ~ insulin, data = Ts, family = "binomial")
    
    # Obtain probabilities
    probs  <- predict(model, newdata = Ts, type = "response")    

    # Store the result
    res <- cbind(res, probs)  

}

colnames(res) <- c("1", "2", "3")

# The prediction is the hypothesis that provides the maximum probability
prediction <- apply(res, 1, function(x){names(which.max(x))})

# Obtain confusion matrix for the test set
confusionMatrix(as.factor(prediction), as.factor(test.set$group))
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  1  2  3
##          1  5  1  0
##          2  2  7  2
##          3  0  5 20
## 
## Overall Statistics
##                                           
##                Accuracy : 0.7619          
##                  95% CI : (0.6055, 0.8795)
##     No Information Rate : 0.5238          
##     P-Value [Acc > NIR] : 0.00133         
##                                           
##                   Kappa : 0.5918          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3
## Sensitivity            0.7143   0.5385   0.9091
## Specificity            0.9714   0.8621   0.7500
## Pos Pred Value         0.8333   0.6364   0.8000
## Neg Pred Value         0.9444   0.8065   0.8824
## Prevalence             0.1667   0.3095   0.5238
## Detection Rate         0.1190   0.1667   0.4762
## Detection Prevalence   0.1429   0.2619   0.5952
## Balanced Accuracy      0.8429   0.7003   0.8295

Q10: How do you obtain your prediction in this case? Why do you have three different values for sensitivity and specificity? Which is the disease group more difficult to classify? Why?

Make sure you understand well the code above. Note that there are classification algorithms other than logistic regression capable to perform multi-class classification directly, as you will see in the next practical.

Remember to answer the short questions (Q) and complete the tasks (T). They can appear as possible questions in the weekly questionaires.