{"cells":[{"metadata":{},"cell_type":"markdown","source":"I share my works on classifying Diabetic Retinopathy images using Fractal Dimensions (FD) and Topological Data Analysis (TDA), after repeatly failed to submit my works via Kaggle Kernel (time out error, out of resources, missing TDA packages [probably due to installation errors] etc...).  These features seem to discriminate eye fundus images with and without diabetic retinopathy very well, but they are not enough in further classifying images into various degree of severity. This findings are consistent with similar research papers below (not exhaustive list of papers in this topic), although these features are extracted differently as they did:\n- [1] Garside, K., Henderson, R., Makarenko, I., & Masoller, C. (2019). Topological data analysis of high resolution diabetic retinopathy images. PloS one, 14(5), e0217413. doi:10.1371/journal.pone.0217413\n- [2] Safitri, D. W., & Juniati, D. (2017, August). Classification of diabetic retinopathy using fractal dimension analysis of eye fundus image. In AIP conference proceedings (Vol. 1867, No. 1, p. 020011). AIP Publishing.\n\nMy approach divides the problems into two stages:\n- First, since healthy and non-healthy are very distinctive using FD and TDA features, a first model (model1) was built to remove high probability healthy cases.  By removing a large portion of sure healthy cases, it addresses the data imbalance problems at the same time;\n- Second, another model (model2) was built to classify the remaining cases into 5 classes (severity = 0, 1, 2, 3, 4).  But the accuracy of this model is not good, according to the CV results. See if we can try the testing dataset after the end of competition to see how well or bad it works in unseen data.\n\nSome more technical notes:\n- Only green channel is used, with images are standardised to width 1024 px.\n- GUDHI is much faster than Dionysus in extracting TDA features. It takes around 4 minutes to process an resized image.  Not sure if it is practical in real world application.\n- Under the standardised images with 1024px, usually around 200 hundred Homology components are found, much less than 200,000 - 300,000 components found in [1], but the TDA features are still useful in distinguishing healthy diabetic retinopathy images."},{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true},"cell_type":"code","source":"# Due to technical issues, and long processing time, I have extracted FD and TDA features from another platform. \n# For details of the feature extractions, please refer to the HTML file: /feature-extraction/Process_Images.html\n\n# read extracted feature datasets for both training and testing images\ntrain <- data.frame(readRDS( \"../input/feature-extraction/train_features.RDS\"))\ntrain.id <- train$id_code\ntrain.y <- factor(train$diagnosis)\n\n# Take first 10 PCA components from Persistence Homology Silhouettes\npca_obj <- prcomp(train[, 25:522], center = TRUE, scale. = TRUE, rank.=10)\nsummary(pca_obj)\ntrain.x <- cbind(train[, 3:24], pca_obj$x)\n\n# Prepare test data:\ntest <- data.frame(readRDS( \"../input/feature-extraction/test_features.RDS\"))\ntest.id <- test$id_code\ntest.x <- cbind(test[, c(3:24)], predict(pca_obj, newdata=test[,25:522]))\n\n# for normalisation of features\nlibrary('matrixStats')\nm <- colMeans(train.x)\ns <- colSds(as.matrix(train.x))\n\ntrain.x_ <- train.x\ntest.x_ <- test.x\n\nfor (i in 1:ncol(train.x)) {\n    train.x_[,i] <- (train.x_[,i]-m[i]) / s[i]\n    test.x_[,i] <- (test.x_[,i]-m[i]) / s[i]\n}\n\nsummary(train.x_)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Inspect the distribution of labels, i.e. severity of eye illnesses\n# The distribution is highly unbalanced, which posed problems for subsequent classification tasks.\naddmargins(table(train.y))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"library(umap)\nlibrary(ggplot2)\nlibrary(dplyr)\n\n# define a visualisation function\nvis_data <- function(inData, label, inTitle) {\n    # use umap to embed high dimension data into 2-D for visualisation\n    set.seed(345985)\n    umap_out <- umap(as.matrix(inData), n_components=2, metric='cosine', negative_sample_rate=30)\n    gdata <- data.frame('x1'=umap_out$layout[,1], 'x2'=umap_out$layout[,2], 'Severity'=as.factor(label))\n    # put label with higher values on top\n    gdata <- gdata %>% arrange(Severity)\n    ggplot(gdata, aes(x=x1, y=x2, color=Severity)) + geom_point(alpha=0.5)  + \n    ggtitle(inTitle) + scale_color_manual(values=c(\"green\", \"yellow\", \"orange\", \"red\", \"purple\"))\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We first inspect the distribution of severity using umap (a variant of t-SNE), with embedding using FD, TDA or FD+TDA features.  From the plots below, it shows that FD +TDA are best for discriminating healthy vs non-healthy cases.  We can first build a model to discriminate healthy vs non-healthy. It would address the data imbalance problem as well."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Visualise Distribution of Severity using Fractal Dimension Features\nvis_data(train.x[,2:22], train.y, 'Severity Distribution using Fractal Dimension Profile')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Visualise Distribution of Severity using Top 10 PCA components (accounted for 90% of variance) in Persistence Homology Silhouettes\nvis_data(scale(train.x_[, c(1, 23:32)]), train.y, 'Severity Distribution using Persistence Homology Summary Statistics')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Visualise Distribution of Severity using both Fractal Dimensions & Persistence Homology Silhouettes\nvis_data(scale(train.x_), train.y, 'Severity Distribution using both Fractal Dimensions & Persistence Homology Silhouettes')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Given the healthy cases are very different from non-healthy cases, we first use SVM to classify healthy (y=0) vs non-healthy (y!=0)\n# We only use Fractal Dimension as feature here.\ntrain.y.binary <- factor(ifelse(train.y==0, 'Y', 'N'))\naddmargins(table(train.y.binary))\n\nlibrary('caret')\n\n# use 2-fold CV, and repeat for 5 times\ntrctrl <- trainControl(method = \"repeatedcv\", number = 2, repeats=5, savePredictions = TRUE, classProbs =  TRUE)\n\n# grid search for optimal parameters\nsvmGrid <- expand.grid(sigma= 2^c(-4:-1), C= 2^c(-2:0))\n\nmodel1 <- train(x=train.x_, y=train.y.binary, method = \"svmRadial\", trControl=trctrl, tuneGrid = svmGrid)\n\n# show the model summary\nmodel1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"library(\"Metrics\")\n\npred1 <- as.numeric(predict(model1)) -1\n\nprint(paste0(\"Quadratic Weighted Kappa: \", round(ScoreQuadraticWeightedKappa(pred1, as.numeric(train.y.binary)-1), 3)))\naddmargins(table('Predict'=pred1, 'Actual'=as.numeric(train.y.binary)-1))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Investigating probability distribution, and get the cut off point to get high probability healthy case:\nprob1 <- predict(model1, type='prob')\nhist(prob1$Y, breaks=100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# define knn function, for getting number of healthy images from selected features\nauc <- function(prob, label, threshold) {\n    P <- sum(prob >=threshold)\n    FP <- sum(label !='Y' & prob>=threshold)\n    T <- sum(label=='Y')\n    TP <- sum(prob >= threshold & label =='Y')\n    \n    # Return number of predicted healthy case\n    print(paste0('Threshold Probability= ', threshold))\n    print(paste0('Number of healthy cases classified: ', P))\n    print(paste0('No. of false negative (misclassified healthy): ', FP))\n    print(paste0('% of false negative (misclassified healthy): ', round(FP/P, 3)))\n    print(paste0('% of recall (true healthy being identified): ', round(TP/T, 3)))\n}\n\nfor (p in seq(0.85, 0.99, by=0.01)) {\n    auc(prob1$Y, train.y.binary, p)\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# pick threshold at 0.92\nhealthy <- prob1$Y >= 0.92\n\n# Distribution of severity after removing high probability healthy cases:\naddmargins(table(train.y[!healthy]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# visualise data after removing high probability healthy case.  Those obvious healthy cases had been removed.\nvis_data(train.x_[!healthy,], train.y[!healthy], 'Severity Dist., after removing high probablity healthy case')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# filter for remaining cases\ntrain.x_.nonhealthy <-train.x_[!healthy,]\ntrain.y.nonhealthy <- train.y[!healthy]\n\n# compute weight to assign higher weight to minor classes\nwgt_summary <- length(train.y.nonhealthy)/table(train.y.nonhealthy) \nwgt <- c(wgt_summary[as.numeric(train.y.nonhealthy)])\nwgt_summary\nhead(wgt)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# SVM on the remaining case:\nlibrary(\"caret\")\n\ntrctrl <- trainControl(method = \"repeatedcv\", number = 5, repeats=5)\n\n# grid search for best parameters\nsvmGrid <- expand.grid(sigma= 2^c(-5:-3), C= 2^c(1:3))\n\nmodel2 <- train(x=train.x_.nonhealthy, y=train.y.nonhealthy, weight=wgt, method = \"svmRadial\", trControl=trctrl, tuneGrid = svmGrid)\n\n# see the summary\nmodel2\n\naddmargins(table('Predict'=predict(model2), 'Actual'=train.y.nonhealthy))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Stage 1, classify sure healthy case using model1\npred_test1 <- predict(model1, newdata=test.x_, type='prob')\n\n# list of high probility healthy case for test\ntest_healthy <- pred_test1$Y >=0.92\n\ntable(test_healthy)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Stage 2, classifying the remaining cases using model2\ntest.x_.nonhealthy <- test.x_[!test_healthy,]\n\npredict_test_nonhealthy <- as.numeric(as.character(predict(model2, newdata=test.x_.nonhealthy)))\n\n\ngdata <- rbind(data.frame('predicted'=0, test.x_[test_healthy,]),\n               data.frame('predicted'=predict_test_nonhealthy, test.x_[!test_healthy,]))\n\naddmargins(table(gdata$predicted))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vis_data(gdata[,-1], gdata$predicted, 'Predicted Severity Distribution in Testing Data')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# generate submission csv\nmy_submission <- data.frame('id_code' = c(test.id[test_healthy], test.id[!test_healthy]), 'diagnosis' = gdata$predicted)\nwrite.csv(my_submission, 'submission.csv', row.names=FALSE)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":1}