report.Rmd

---
title: "Data Science Machine Learning 'Choose Your Own' Project - Indian Liver Patient Records"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE)
```

```{r load_packages}
if(!require(tidyverse)) install.packages("tidyverse", 
                                         repos = "http://cran.us.r-project.org")
if(!require(caret)) install.packages("caret", 
                                     repos = "http://cran.us.r-project.org")
if(!require(data.table)) install.packages("data.table", 
                                          repos = "http://cran.us.r-project.org")
if(!require(ggforce)) install.packages("ggforce", 
                                          repos = "http://cran.us.r-project.org")
if(!require(ggthemes)) install.packages("ggthemes", 
                                       repos = "http://cran.us.r-project.org")
if(!require(e1071)) install.packages("e1071", 
                                        repos = "http://cran.us.r-project.org")
if(!require(randomForest)) install.packages("randomForest", 
                                     repos = "http://cran.us.r-project.org")
if(!require(kernlab)) install.packages("kernlab", 
                                            repos = "http://cran.us.r-project.org")
if(!require(klaR)) install.packages("klaR", 
                                       repos = "http://cran.us.r-project.org")
if(!require(matrixStats)) install.packages("matrixStats", 
                                    repos = "http://cran.us.r-project.org")
if(!require(knitr)) install.packages("knitr", 
                                    repos = "http://cran.us.r-project.org")
if(!require(kableExtra)) install.packages("kableExtra", 
                                    repos = "http://cran.us.r-project.org")
if(!require(patchwork)) install.packages("patchwork", 
                                           repos = "http://cran.us.r-project.org")
```

# Overview
Being an important organ of the human body, diseases of the liver can be life-threatening. It is therefore highly important to diagnose and subsequently treat them as soon as possible. In clinical practice, liver function tests are used to this end. In liver function tests, a patient's blood is analyzed with regard to different biomarkers, e.g. the activity of enzymes which are specific to liver functions. Those biomarkers can be used as a proxy measure of the patient's liver's health.

The present data set contains the results of liver function tests of Indian patients and is available from the Machine Learning Repository of the University of California, Irvine (Dua and Graff, 2019). The data is classified into liver-patient ('affected' in the remainder of the report) and non-liver-patient ('unaffected' in the remainder of the report). 

The aim of this project is to use the available biomarker data to predict a person's affection status based on those same biomarkers. Firstly, standard methods of data exploration will be applied to become familiar with the data and visualize the data and the marker distribution. Subsequently, different machine learning algorithms will be trained on a part of the data and the best model will be used on the validation data set (split from the original data set in the first step of the analysis). 


# Data exploration (exploratory data analysis, EDA)
## The data set
The data set is contained in a comma-separated file (.csv), which is downloaded from the UCI repository. Since this file does not contain a header, the relevant column names are be taken from the description website (https://archive.ics.uci.edu/ml/datasets/ILPD+(Indian+Liver+Patient+Dataset)) and added to the resulting data frame. Here are the first three entries (note that the last column has been recoded as mentioned above):
```{r download_data}
all_data = fread("https://archive.ics.uci.edu/ml/machine-learning-databases/00225/Indian%20Liver%20Patient%20Dataset%20(ILPD).csv")
names(all_data) = c("Age", "Gender", "Total_bilirubin", "Direct_bilirubin", "Alkaline_phosphatase", "Alanine_aminotransferase", "Aspartate_aminotransferase", "Total_protein", "Albumin", "Albumin_globulin_ratio", "Affection_status")
all_data = all_data %>% mutate(Affection_status = ifelse(Affection_status == 1, "Affected", "Unaffected"))
kable(head(all_data, 3)) %>%
  row_spec(0, angle = 90)
```


In total, there are `r nrow(all_data)` observations of `r ncol(all_data)` variables in the data set. To make sure that the algorithms created later on will be valid, an important step is to see whether any data is missing. Using R's `is.na` function, it becomes apparent that information on the albumin-globulin-ratio is missing for four patients. These will be moved to a new data frame and analyzed separately in the end:

```{r missing_data}
missing_data = all_data %>% filter(is.na(Albumin_globulin_ratio))
all_data = all_data %>% filter(!is.na(Albumin_globulin_ratio))
kable(missing_data) %>% row_spec(0, angle = 90)
```


The new total number of observations in the data set is `r nrow(all_data)`.
The following table shows the distribution of the affection status in the data:

```{r aff_status}
all_data %>% dplyr::select(Affection_status) %>% table() %>% prop.table() %>% kable()
```


There is a slight imbalance in the data, with more affected than unaffected patient's available. However, this is still moderate and should not have consequences on the model generation.

To assess the accuracy of the different models, the data set is split into a `model` data set and a `validation` data set in a 9:1 ratio. This ensures that sufficient data is available for the training on the `model` data, while at the same time preventing overtraining. All further data visualization will be performed on the `model` data. The distribution of the affection status is similar between the two generated data sets:

`Model data`:
```{r split_validation}
# set.seed(1, sample.kind="Rounding")
set.seed(1)
test_index = createDataPartition(y = all_data$Affection_status, times = 1, p = 0.1, list = FALSE)
model_data = all_data[-test_index,]
validation_data = all_data[test_index,]
model_data %>% dplyr::select(Affection_status) %>% table() %>% prop.table() %>% kable
```

`Validation data`:
```{r validation_distribution}
validation_data %>% dplyr::select(Affection_status) %>% table() %>% prop.table() %>% kable
```

## Data visualization
The first variable is the age of the patients:

```{r age_plot, fig.width = 4.2, fig.height = 3}
model_data %>% ggplot(aes(x=Age)) + geom_histogram(bins=20) + ggtitle("Age distribution") + ylab("Count")
```

As can be seen, the age looks somewhat normally distributed. The youngest patient is `r min(model_data$Age)` years old. The age of the oldest patient is given as `r max(model_data$Age)` years. However, the data repository states that the age of every patient above the age of 90 was set to 90. The median age of the patients is `r median(model_data$Age)`. 

The second variable is the gender:

```{r gender_plot, fig.width = 4.2, fig.height = 3}
model_data %>% ggplot(aes(x=Gender)) + geom_bar() + ggtitle("Gender distribution") + ylab("Count")
```

There are approximately three times as many men as women in the data. As it has been shown that there are gender differences in liver disease prevalence (Buzzetti et al., 2017), this might introduce a bias into the data. 

Plotting both of the first two variables at once shows that the age distribution between the genders mimic each other fairly well:

```{r age_gender_plot, fig.width = 4.2, fig.height = 3}
model_data %>% ggplot(aes(x=Age, fill=Gender)) + geom_histogram(bins=20, alpha=0.5, position="identity") + ggtitle("Age distribution stratified by gender") + ylab("Count")
```

As the distribution of the affection status has been demonstrated in the tables above, we can instead complement this by the following visualization, which incorporates all the previous dimensions of information into one plot. In this plot, the points show the age and affection status stratified by gender (color), while the respective Tufte boxplots (minimalistic versions of standard boxplots; black) show this distribution irrespective of gender, i.e. age versus affection status:

```{r age_gender_aff_plot, fig.width = 5.6, fig.height = 4}
model_data %>% ggplot(aes(y=Age, x=Affection_status, color=Gender)) + geom_sina() + geom_tufteboxplot(colour="black") + ggtitle("Affection status stratified by age and gender")
model_data$Affection_status = as.factor(model_data$Affection_status)
model_data = model_data %>% mutate(Gender=ifelse(Gender=="Male", 1, 2))
```

This plot shows that neither age nor gender are perfectly matched between the two groups of affection status. This might be due to an imbalance in the data but might also reflect real differences in disease etiology and prevalence. 

Next up is the biomarker data. Present in the data set are:

* Total bilirubin (degradation product of hemoglobin)
* Direct bilirubin (sub-portion of the above)
* Alkaline phosphatase activity (ALP; enzyme)
* Alanine aminotransferase activity (ALT; enzyme)
* Aspartate aminotransferase activity (AST; enzyme)
* Total protein (protein content in the blood)
* Albumin (a protein family)
* Albumin-globulin-ratio (ratio between albumin and another protein family, globulin)

The following patchwork plot visualizes all the above mentioned data at once for an easier overview. Note that both bilirubin plot y-axes and all three enzymes plot y-axes have been `sqrt` transformed. 

```{r biomarker_data}
total_bilirubin = model_data %>% ggplot(aes(y=Total_bilirubin, x=Affection_status)) + geom_boxplot() + 
  scale_y_sqrt() + ylab("Total bilirubin (micromol/L)") + theme(axis.title.x = element_blank(), text = element_text(size=8))
direct_bilirubin = model_data %>% ggplot(aes(y=Direct_bilirubin, x=Affection_status)) + geom_boxplot() + 
  scale_y_sqrt() + ylab("Direct bilirubin (micromol/L)") + theme(axis.title.x = element_blank(), text = element_text(size=8))
alkaline = model_data %>% ggplot(aes(y=Alkaline_phosphatase, x=Affection_status)) + geom_boxplot() + 
  scale_y_sqrt() + ylab("ALP activity (U/L)") + theme(axis.title.x = element_blank(), text = element_text(size=8))
alanine = model_data %>% ggplot(aes(y=Alanine_aminotransferase, x=Affection_status)) + geom_boxplot() + 
  scale_y_sqrt() + ylab("ALT activity (U/L)") + theme(axis.title.x = element_blank(), text = element_text(size=8))
aspartate = model_data %>% ggplot(aes(y=Aspartate_aminotransferase, x=Affection_status)) + geom_boxplot() + 
  scale_y_sqrt() + ylab("AST activity (U/L)") + theme(axis.title.x = element_blank(), text = element_text(size=8))
total_protein = model_data %>% ggplot(aes(y=Total_protein, x=Affection_status)) + geom_boxplot()+ ylab("Total protein (g/dl)") + theme(axis.title.x = element_blank(), text = element_text(size=8))
albumin = model_data %>% ggplot(aes(y=Albumin, x=Affection_status)) + geom_boxplot() + ylab("Albumin (U/L)") + theme(axis.title.x = element_blank(), text = element_text(size=8))
ag_ratio = model_data %>% ggplot(aes(y=Albumin_globulin_ratio, x=Affection_status)) + geom_boxplot() + ylab("Albumin globulin ratio") + theme(axis.title.x = element_blank(), text = element_text(size=8))
(total_bilirubin + direct_bilirubin) / (alkaline + alanine + aspartate) / (total_protein + albumin + ag_ratio)
```


As can be seen from this plot combination, there are strong differences between affected and unaffected persons for the two bilirubin biomarkers and the three enzyme activity biomarkers. However, this difference is far more pronounced for the outliers in each category than for the median values. For total protein, albumin content and albumin-globulin-ratio, the two groups are very similar. This strongly indicates that bilirubin and the enzyme activities will be most useful for the model generation to distinguish between affected and unaffected persons. 

# Modelling and model assessment
To be able to generate different models and test them, while at the same time disregarding the `validation` data set, the `model` data set will again be split in a 9:1 ratio to obtain a `training` data set and a `test` data set.

```{r split_data}
# set.seed(1, sample.kind="Rounding")
set.seed(1)
test_index = createDataPartition(y = model_data$Affection_status, times = 1, p = 0.1, list = FALSE)
train_data = model_data[-test_index,]
test_data = model_data[test_index,]

train_data_x = train_data %>% dplyr::select(-Affection_status)
train_data_y = train_data %>% pull(Affection_status)
test_data_x = test_data %>% dplyr::select(-Affection_status)
test_data_y = test_data %>% pull(Affection_status)
```

As can be seen, the affection status distribution is still similar between the `training` data set and the `test` data set:

Training data set affection status distribution:
```{r train_distribution}
train_data %>% dplyr::select(Affection_status) %>% table() %>% prop.table() %>% kable
```

Test data set affection status distribution:
```{r test_distribution}
test_data %>% dplyr::select(Affection_status) %>% table() %>% prop.table() %>% kable
```

To create a baseline algorithm, we can use a model which (like tossing a coin) has a 50 % chance for predicting either 'affected' or 'unaffected':
```{r guess, echo=TRUE}
guess = function(data){
  n = nrow(data)
  prediction = sample(c("Affected", "Unaffected"), n, replace=TRUE)
  prediction = factor(prediction, levels=c("Affected", "Unaffected"))
  return(prediction)
}
```
To keep track of all model's performances, a separate data frame containing the obtained accuracies on the `test` data set is created:
```{r results_01}
# set.seed(2, sample.kind="Rounding")
set.seed(2)
results = data.frame(Model = "Guessing", Accuracy = mean(guess(test_data_x) == test_data_y))
kable(results)
```


A second option for a baseline algorithm is one which always predicts 'affected'. This will naturally have better performance then the guessing algorithm, due to the slightly skewed data:
```{r affected, echo=TRUE}
affected = function(data){
  n = nrow(data)
  prediction = rep("Affected", n)
  prediction = factor(prediction, levels=c("Affected", "Unaffected"))
  return(prediction)
}
```
``` {r results_02}
results = bind_rows(results, data.frame(Model = "All affected", Accuracy = mean(affected(test_data_x) == test_data_y)))
kable(results)
```
As expected, this algorithm is indeed superior to plain guessing. 


A final, relatively simple, algorithm, which is not based on machine learning itself is looking at the expected values for the present biomarkers and deciding on 'affected' status, once a certain number of biomarkers are outside of their given range. To do this, the following values could be used (Gowda et al., 2009):
```{r ref_values}
ref_values = data.frame(Parameter = c("Total bilirubin", "Direct bilirubin", "Alkaline phosphatase", "Alanine aminotransferase", "Aspartate aminotransferase"),
                        Unit = c("µmol/L", "µmol/L", "U/L", "U/L", "U/L"),
                        Reference_lower = c(2, 0, 41, 7, 0),
                        Reference_upper = c(21, 8, 133, 56, 35))
kable(ref_values)
```

The following algorithm can then be used to predict the affection status based on the number of biomarkers exceeding the threshold:
```{r clinical, echo=TRUE}
clinical = function(data=test_data_x, ref=ref_values, threshold=1) {
  data %>% mutate(ToBi = ifelse(Total_bilirubin < ref_values[1,3] |
                                  Total_bilirubin > ref_values[1,4], 1, 0)) %>%
    mutate(DiBi = ifelse(Direct_bilirubin < ref_values[2,3] |
                           Direct_bilirubin > ref_values[2,4], 1, 0)) %>%
    mutate(AlPh = ifelse(Alkaline_phosphatase < ref_values[3,3] |
                           Alkaline_phosphatase > ref_values[3,4], 1, 0)) %>%
    mutate(AlAm = ifelse(Alanine_aminotransferase < ref_values[4,3] |
                           Alanine_aminotransferase > ref_values[4,4], 1, 0)) %>%
    mutate(AsAm = ifelse(Aspartate_aminotransferase < ref_values[5,3] |
                           Aspartate_aminotransferase > ref_values[5,4], 1, 0)) %>%
    mutate(All = ToBi + DiBi + AlPh + AlAm + AsAm) %>%
    mutate(prediction = ifelse(All >= threshold, "Affected", "Unaffected")) %>%
    pull(prediction)
}
```

Setting this threshold to any value from one to five extends the accuracy table by five more values:
```{r results_03}
results = bind_rows(results, data.frame(Model = c("Clinical parameters (1)", "Clinical parameters (2)", "Clinical parameters (3)", "Clinical parameters (4)", "Clinical parameters (5)"), 
                                        Accuracy = c(mean(clinical(data=test_data_x, threshold=1) == test_data_y), mean(clinical(data=test_data_x, threshold=2) == test_data_y), mean(clinical(data=test_data_x, threshold=3) == test_data_y), mean(clinical(data=test_data_x, threshold=4) == test_data_y), mean(clinical(data=test_data_x, threshold=5) == test_data_y))))
kable(results)
```

Here, even the best-performing version of the algorithm (which predicts 'affected' once one of the used biomarkers exceeds its respective threshold) performs equal to predicting 'affected' for everyone (i.e. accuracy = proportion of affected patients).


To remedy this, 'true' machine learning algorithms will be used, which can determine more intricate relationships between the individual biomarkers (as well as age and gender). The `train` function of the `caret` package will be used for the generation of each model. Where applicable 50-fold cross-validation (9:1 ratio) is used.

First, five different machine learning models are trained. A k-nearest neighbor model (kNN); a random forest model (rf); a generalized linear model (glm); a support vector machine model (SVM); and a regularized discriminant analysis model (RDA):

```{r train_normal, echo=TRUE}
control = trainControl(method="cv", number=50, p=0.9)
set.seed(3)
fit_knn = train(train_data_x, train_data_y, method="knn", 
                tuneGrid=data.frame(k=seq(11, 31, 2)), 
                trControl=control)
set.seed(5)
fit_rf = train(train_data_x, train_data_y, method="rf", 
               tuneGrid=data.frame(mtry=seq(1, 10, 1)), 
               trControl=control)
set.seed(8)
fit_glm = train(train_data_x, train_data_y, method="glm")
set.seed(13)
fit_svm = train(train_data_x, train_data_y, method="svmRadialSigma", 
                tuneGrid=data.frame(expand.grid(C=c(2**seq(-5, 1, 2)), 
                                                sigma=c(2**seq(-5, 1, 2)))), 
                trControl=control)
set.seed(21)
fit_rda = train(train_data_x, train_data_y, method="rda")
```

For the k-nearest neighbors model, `k`-values between 11 and 31 were tested. A standard approach is to choose `k` as `sqrt(n)`, which in this case would be `r sqrt(nrow(train_data))`. Hence, sufficient values around this were chosen. For binary classification approaches, `k` should preferably be an uneven number. For the random forest model, `mtry` is tested for every possible value (based on the number of relevant variables). In a support vector machine, `C` is the penalty factor and `sigma` is the distance between the used training spheres. Both are usually rather small values. 

These models add upon the results data frame:
```{r results_04}
results = bind_rows(results, data.frame(Model= c("k-nearest neighbors", "Random forest", "Generalized linear model", "Support vector machine", "Regularized discriminant analysis"),
                    Accuracy = c(mean(predict(fit_knn, test_data_x) == test_data_y), mean(predict(fit_rf, test_data_x) == test_data_y), mean(predict(fit_glm, test_data_x) == test_data_y), mean(predict(fit_svm, test_data_x) == test_data_y), mean(predict(fit_rda, test_data_x) == test_data_y))))
kable(results)
```

While some of those models perform better than the approaches based on clinical parameters/threshold or simply assigning 'affected' to everyone, kNN performs worse and SVM and RDA show the exact same accuracy. This implies that those two models also predict 'affected' for every patient.


To further optimize the algorithms, dimension reduction via principal component analysis might help:

```{r pca}
pca = prcomp(train_data_x)
pca_importance = function(x) {
  vars = x$sdev^2
  vars = vars/sum(vars)
  rbind(`Standard deviation` = x$sdev, `Proportion of Variance` = vars, 
      `Cumulative Proportion` = cumsum(vars))
}
pca_summary = pca_importance(pca)
pca_summary = round(pca_summary, 2)
pca_summary = data.frame(pca_summary)
names(pca_summary) = c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10")
kable(pca_summary)
```
This shows that the first two principal components already account for more than 90 % of the variance in the data. Adding the third principal component increases the explained variance to more than 99 %.

Therefore, only three dimensions will be used for the models trained on dimension-reduced data; the remaining parameters are kept unchanged with regard to the previous model training:
```{r dimensions}
k = 3
x_train_dim = pca$x[, 1:k]
y_train_dim = train_data_y
x_test_dim = as.matrix(sweep(test_data_x, 2, colMeans(test_data_x))) %*% pca$rotation
x_test_dim = x_test_dim[, 1:k]
```
```{r train_dim_red, echo=TRUE}
set.seed(34)
fit_knn_dim = train(x_train_dim, y_train_dim, method="knn", 
                    tuneGrid=data.frame(k=seq(11, 31, 2)), 
                    trControl=control)
set.seed(55)
fit_rf_dim = train(x_train_dim, y_train_dim, method="rf", 
                   tuneGrid=data.frame(mtry=seq(1, 3, 1)), 
                   trControl=control)
set.seed(89)
fit_glm_dim = train(x_train_dim, y_train_dim, method="glm")
set.seed(144)
fit_svm_dim = train(x_train_dim, y_train_dim, method="svmRadialSigma", 
                    tuneGrid=data.frame(expand.grid(C=c(2**seq(-5, 1, 2)), 
                                                    sigma=c(2**seq(-5, 1, 2)))), 
                    trControl=control)
set.seed(233)
fit_rda_dim = train(x_train_dim, y_train_dim, method="rda")
```

This adds five more data points to the results data frame:
```{r results_05}
results = bind_rows(results, data.frame(Model = c("k-nearest neighbors (dim. red.)", "Random forest (dim. red.)", "Generalized linear model (dim. red.)", "Support vector machine (dim. red.)", "Regularized discriminant analysis (dim. red.)"),
                                        Accuracy = c(mean(predict(fit_knn_dim, x_test_dim) == test_data_y), mean(predict(fit_rf_dim, x_test_dim) == test_data_y), mean(predict(fit_glm_dim, x_test_dim) == test_data_y), mean(predict(fit_svm_dim, x_test_dim) == test_data_y), mean(predict(fit_rda_dim, x_test_dim) == test_data_y))))
kable(results)
```

Apparently, most of the five algorithms perform worse than or equal to the original model. Only the kNN model shows a slight increase in accuracy.

Since the performances so far vary widely, an ensemble model might leverage the strengths of each of the trained models. The following ensemble model is based on the models trained with the original data (before dimension reduction). The `threshold` parameter determines how many 'unaffected' predictions by the individual models are necessary for the ensemble model itself to predict 'unaffected':
```{r ensemble_normal, echo=TRUE}
ensemble = function(data=test_data_x, threshold=1){
  output = data.frame(knn = rep(NA, nrow(data)))
  output$knn = predict(fit_knn, data)
  output$rf = predict(fit_rf, data)
  output$glm = predict(fit_glm, data)
  output$svm = predict(fit_svm, data)
  output$rda = predict(fit_rda, data)
  output$count = rowCounts(as.matrix(output), value="Unaffected")
  output = output %>% mutate(prediction = ifelse(count>=threshold,
                                                 "Unaffected", "Affected"))
  output$prediction = factor(output$prediction, levels=c("Affected", "Unaffected"))
  return(output)
}
```


A similar model can be constructed for the models trained on the data after dimension reduction (based on the same underlying principle as described above):
```{r ensemble_dim_red, echo=TRUE}
ensemble_dim = function(data=x_test_dim, threshold=1){
  output = data.frame(knn = rep(NA, nrow(data)))
  output$knn = predict(fit_knn_dim, data)
  output$rf = predict(fit_rf_dim, data)
  output$glm = predict(fit_glm_dim, data)
  output$svm = predict(fit_svm_dim, data)
  output$rda = predict(fit_rda_dim, data)
  output$count = rowCounts(as.matrix(output), value="Unaffected")
  output = output %>% mutate(prediction = ifelse(count>=threshold,
                                                 "Unaffected", "Affected"))
  output$prediction = factor(output$prediction, levels=c("Affected", "Unaffected"))
  return(output)
}
```

Using those two ensemble models (with varying thresholds from one to five) further extends the results data frame to its final dimensions:

```{r results_06}
results = bind_rows(results, data.frame(Model = c("Ensemble (1)", "Ensemble (2)", "Ensemble (3)", "Ensemble (4)", "Ensemble (5)"), 
                                        Accuracy = c(mean(ensemble(data=test_data_x, threshold=1)$prediction == test_data_y), mean(ensemble(data=test_data_x, threshold=2)$prediction == test_data_y), mean(ensemble(data=test_data_x, threshold=3)$prediction == test_data_y), mean(ensemble(data=test_data_x, threshold=4)$prediction == test_data_y), mean(ensemble(data=test_data_x, threshold=5)$prediction == test_data_y))))
results = bind_rows(results, data.frame(Model = c("Ensemble (dim. red.) (1)", "Ensemble (dim. red.) (2)", "Ensemble (dim. red.) (3)", "Ensemble (dim. red.) (4)", "Ensemble (dim. red.) (5)"), 
                                        Accuracy = c(mean(ensemble_dim(data=x_test_dim, threshold=1)$prediction == test_data_y), mean(ensemble_dim(data=x_test_dim, threshold=2)$prediction == test_data_y), mean(ensemble_dim(data=x_test_dim, threshold=3)$prediction == test_data_y), mean(ensemble_dim(data=x_test_dim, threshold=4)$prediction == test_data_y), mean(ensemble_dim(data=x_test_dim, threshold=5)$prediction == test_data_y))))
kable(results)
```

These results show that the best performing model is still the GLM without dimension reduction, followed by kNN before dimension reduction and the ensemble before dimension reduction (using `threshold = 1`). The GLM also shows high values for further performance parameters:

```{r GLM_performance}
glm_performance = confusionMatrix(predict(fit_glm, test_data_x), test_data_y)$byClass
glm_performance = as.data.frame(glm_performance)
names(glm_performance) = c("Value")
kable(glm_performance)
```


GLM is therefore chosen as the final model and trained on the `model` data set:
```{r prep_model_data}
model_data_x = model_data %>% dplyr::select(-Affection_status)
model_data_y = model_data %>% pull(Affection_status)
```
```{r glm_train, echo=TRUE}
set.seed(377)
fit_glm_final = train(model_data_x, model_data_y, method="glm")
```

This final model could now be tested on the `validation` data. However, as a final check, the performance on the data set created from the observations with missing values of the albumin-globulin-ratio is tested. To prevent errors due to the `NA`s for this variable, the mean albumin-globulin-ratio of the `model` data will be used as a proxy value (`r mean(model_data$Albumin_globulin_ratio)`).

```{r prep_missing}
missing_data$Albumin_globulin_ratio = mean(model_data$Albumin_globulin_ratio)
missing_data$Affection_status = as.factor(missing_data$Affection_status)
missing_data = missing_data %>% mutate(Gender=ifelse(Gender=="Male", 1, 2))
missing_data_x = missing_data %>% dplyr::select(-Affection_status)
missing_data_y = missing_data %>% pull(Affection_status)
```

Here, GLM achieves an accuracy of `r mean(predict(fit_glm_final, missing_data_x) == missing_data_y)`, i.e. three out of four are predicted correctly.


Now for the performance of the GLM on the validation data:

```{r prep_validation}
validation_data$Affection_status = as.factor(validation_data$Affection_status)
validation_data = validation_data %>% mutate(Gender=ifelse(Gender=="Male", 1, 2))
validation_data_x = validation_data %>% dplyr::select(-Affection_status)
validation_data_y = validation_data %>% pull(Affection_status)
```

The accuracy of this model is `r mean(predict(fit_glm_final, validation_data_x) == validation_data_y)`. Further performance parameters can be seen in the following table:

```{r validation_check}
glm_performance = confusionMatrix(predict(fit_glm_final, validation_data_x), validation_data_y)$byClass
glm_performance = as.data.frame(glm_performance)
names(glm_performance) = c("Value")
kable(glm_performance)
```

While the accuracy of the model is close to 0.8, the very high sensitivity close to 1 in combination with the still low specificity indicates that predicting 'affected' still occurs very frequently, regardless of the actual outcome.


# Final results
Training the GLM on the `training` data and testing it on the `test` data led to an accuracy of about 0.77, thereby outperforming all other models and algorithms. The sensitivity was close to 1, with specificity being 0.27. Training this model on the whole `model` data (i.e. the combined `training` data and `test` data), made predictions on the `missing` data and validation on the `validation` data possible. 

For the `missing` data, the prediction accuracy was 0.75. However, since there were only a total of four observations in this data set, the meaningfulness of this is rather limited.

More importantly, the accuracy on the `validation` data showed a slight increase. Simultaneously, the specificity increased from 0.27 to below 0.35. 

The following table shows all obtained accuracies; note that those were obtained with different models, trained and tested on different fractions of the whole data set:
```{r results_final}
results = bind_rows(results, data.frame(Model = c("GLM on missing data", "GLM on validation data"), Accuracy = c(0.750, 0.797)))
results$Accuracy = round(results$Accuracy, 3)
kable(results)
```


# Conclusion

The aim of this project was to predict the affection status of liver patients based on different biomarkers. The available data (besides age and gender) included bilirubin and protein amounts in the patients' blood as well as the activity of three liver-specific enzymes. Based on the data visualization, the most striking differences could be observed for those enzymes as well as the bilirubin levels.

Training a number of machine learning models and ensembles on the data led to the emergence of a generalized linear model (GLM) as the best-performing model in terms of accuracy. This model also showed a very high sensitivity, although the specificity was rather low. Re-training this model on a larger fraction of the data set led to an acceptable performance on a subset of incomplete data, which had been filled with the mean of the data set (accuracy of 0.75). Finally, using the `validation` data showed that the model performed with an accuracy of 0.797. The sensitivity remained very high (0.976), while the specificity was slightly higher than before (0.351).

It became obvious that all models (including the final GLM) perform far better predictin the 'affected' status than predicting the 'unaffected' status. A possible reason for this is the slight imbalance in the data set. There were approximately two times as many affected patients as there were unaffected patients. Secondly, the visualization and comparison of the biomarker data implies that not every patient who is affected by a liver disease shows highly abnormal values for the biomarkers contained in this data set. It has been shown previously that patients with (chronic) liver disease may present without clinical symptons and values of liver function tests well within the accepted range (Ahmed et al., 2018), which is in line with these findings. Thus, affected patients which do not show significantly different levels of bilirubin, proteins or enzyme activity severely increase the difficulty of distinguishing between those. Furthermore, the affected patients with low levels of the biomarkers are likely to have lead the machine learning approaches to falsely flag unaffected patients as affected. 

In the future, these shortcomings could be remedied by different approaches. First and foremost would be an increase in the number of observations. Since machine learning is a data-driven approach, a larger number of available data often results in better results (as long as over-training is prevented). Another option, based on the present data set, would be either under-sampling of the data with the 'affected' flag or over-sampling of the data with the 'unaffected' flag. This way, the models might be able to better distinguish between those two outcomes. A third approach, which has shown promise in many areas, would be the application of deep learning. All those approached, or a combination thereof, would likely increase the prediction performance and therefore resolve the limitations of this project.


# References
Ahmed Z, Ahmed U, Walayat S, Ren J, Martin DK, Moole H, Koppe S, Yong S, Dhillon S. Liver function tests in identifying patients with liver disease. 2018. Clin Exp Gastroenterol. doi: 10.2147/CEG.S160537

Buzzetti E, Parikh PM, Gerussi A, Tsochatzis E. Gender differences in liver disease and the drug-dose gender gap. 2017. Pharmacol Res. doi: 10.1016/j.phrs.2017.03.014

Dua D, Graff C. UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. 2019. University of California, School of Information and Computer Science. doi: NA

Gowda S, Desai PB, Hull VV, Math AAK, Vernekar SN, Kulkarni SS.  A review on laboratory liver function tests. 2009. Pan Afr Med J. doi: 10.11604/pamj.25/11/2009.3.17.125

Hall P, Cash J. What is the Real Function of the Liver ‘Function’ Tests? 2012. Ulster Med J. doi: NA