{"cells":[{"metadata":{"_uuid":"4e1957454eaaa91410df84c153a660897a0bd95c"},"cell_type":"markdown","source":"This is a kernel for the 4th DSC challenge - unsupervised ML. The goal is to play around with clustering and/or PCA on the digits dataset. The idea is to do as much as possible without the labels. I'm going to use R as I've never used R before."},{"metadata":{"_uuid":"09dd689df9922d611a4ea363631de18207ca49af"},"cell_type":"markdown","source":"First read in the training data."},{"metadata":{"_uuid":"011bd4cfdd51fe6aa1b73bff008d86236f15383a","_execution_state":"idle","trusted":true},"cell_type":"code","source":"## Importing packages\nlibrary(tidyverse) # metapackage with lots of helpful functions\nlibrary(readr) # package to read in tabular data\nlibrary(caret) # package with regression and classification models\n\ntrain <- read_csv(\"../input/train.csv\")\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"077d4ae090652410f7952fdf1b1ea0a132233e4b"},"cell_type":"markdown","source":"Have a look at the training data..."},{"metadata":{"trusted":true,"_uuid":"5173ab97cc383bd7d749c2438353f2d100b3822d"},"cell_type":"code","source":"dim(train)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9523e0640a0717d9795b85455598d07e48285155"},"cell_type":"markdown","source":"Extract the labels now..."},{"metadata":{"trusted":true,"_uuid":"9001331386059871875c29efafb6fd3cb9bfd38f"},"cell_type":"code","source":"labels <- train[1]\ndim(labels)\nsummary(labels)\n\ntrain_predictors <- train[-1]\ndim(train_predictors)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"67c6ab581cb81bd0b2391232adfcb6ad3d31f04e"},"cell_type":"markdown","source":"The data has 42000 rows and 785 columns. \nLet's have alook at the digits by drawing them."},{"metadata":{"trusted":true,"_uuid":"3cbb6c03e774081dd5c17f27295e3f623738fe23"},"cell_type":"code","source":"sample <- sample(1:nrow(train),50)\nvar <- t(train[sample,-1])\nvar_matrix <- lapply(1:50, function(x) matrix(var[,x], ncol = 28))\noriginal_par <- par(no.readonly = T)\npar(mfrow = c(5,10), mar =  c(.1,.1,.1,.1))\n\nfor(i in 1:50) {\n  for(j in 1:28) {\n    var_matrix[[i]][j,] <- rev(var_matrix[[i]][j,])\n  }\n  image(var_matrix[[i]],col = grey.colors(255),axes = F)\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c66465893e5705bb93c84fb3cf63f495befaf38c"},"cell_type":"code","source":"par(original_par)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4f1e9a8fd3f9a2b732506d8b5fc1498fcbb54727"},"cell_type":"markdown","source":"Check for predictors with near zero variance and remove them."},{"metadata":{"trusted":true,"_uuid":"91141bb8b4eba643be0c0296acbddd4e1e05c6f0"},"cell_type":"code","source":"near_zero <- nearZeroVar(train[,-1], saveMetrics = T, freqCut = 10000/1, uniqueCut = 1/7)\nsum(near_zero$zeroVar)\nsum(near_zero$nzv)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9b52f6d6a947faf78a3c38ebe992f84030652b79"},"cell_type":"markdown","source":"There are 208 near-zero variance predictor, of which are 76 zero-variance predictors. Exclude these from the training data as they will not add much."},{"metadata":{"trusted":true,"_uuid":"72b132fb0b5cdd7b5a7fe9eb011de72943f80615"},"cell_type":"code","source":"cut_var <- rownames(near_zero[near_zero$nzv==TRUE,])\nvar <- setdiff(names(train_predictors),cut_var)\ntrain_no_nzv <- train_predictors[,var]\ndim(train_no_nzv)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5a2bc4010477bac95365481d253b46019965f1ed"},"cell_type":"markdown","source":"The training data is now 42000 x 577. The data has many dimensions and so we can try to use PCA to find the most important components. I will normalise the data first by the maximum value (255) and then get the covariance matrix."},{"metadata":{"trusted":true,"_uuid":"ab374c87bc3945f8667751857a65e7e1e703378f","scrolled":false},"cell_type":"code","source":"# Normalise the data to the max value (255)\ntrain_no_nzv <- train_no_nzv/255\n\n# Get the covariance matrix of the predictors\ncovtrain <- cov(train_no_nzv)\n\ndim(train_no_nzv)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c836add812b02aa5eaf388ea15e99890334c9dc6"},"cell_type":"markdown","source":"Try applying PCA to the covariance matrix... \nFind the components which explain the maximum variance - more explained variance = more information. Find the variance explained by each component by dividing the variance by the sum of the total variance."},{"metadata":{"trusted":true,"_uuid":"2d59c8f885a0ace340c85a3fae91e62e878f77ec"},"cell_type":"code","source":"# run the pca\ntrain_pc <- prcomp(covtrain)\n\n# get the variance \nvar_ex <- train_pc$sdev^2/sum(train_pc$sdev^2)\n\n# get the cumulative sum of the variances\nvar_cum <- cumsum(var_ex)\n\n# put the result into a data frame\nresult <- data.frame(num = 1:length(train_pc$sdev), varex = var_ex, varcum = var_cum)\n\n# plot number of components vs cumulative sum of variance\nplot(result$num,result$varcum,type=\"b\",xlim=c(0,100),\n     main=\"Variance of the data explained by the top 100 components\",\n     xlab=\"Number of components\",ylab=\"Variance explained\")\nabline(v = 25, lty = 2)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"876f644442b159f32ac59ef272de07f5370eb6bd"},"cell_type":"markdown","source":"It looks like most of the variance of the dataset is explained by the first 25 components. \n\nUseful measures of prcomp:\nCenter - mean of variables\nScale - standard deviation of variables\nRotation - principle component loading. Each column of rotation matrix contains the principle component loading vector.\nx -  the principal component score vectors.\n"},{"metadata":{"trusted":true,"_uuid":"4ea938f4df7f0dd7761db06e52900534397e21ff"},"cell_type":"code","source":"train_score <- as.matrix(train_no_nzv) %*% train_pc$rotation[,1:25]\nlabel <- as.factor(train[[1]])\ntrain_score_label <- cbind(label,as.data.frame(train_score))\n\ncolors <- rainbow(length(unique(train_score_label$label)))\nnames(colors) <- unique(train_score_label$label)\nplot(train_score_label$PC1,train_score_label$PC2,type=\"n\",main=\"First two principal components\")\ntext(train_score_label$PC1,train_score_label$PC2,label=train_score_label$label,col=colors[train_score_label$label])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"66f9a79c360e1891a53fe85ec0664fcdfc8df2d1"},"cell_type":"markdown","source":"It looks like just the first two components already do a good job of separating out the digits."},{"metadata":{"_uuid":"4e38508ca4a6fa860a85d5e6727b3205c865475f"},"cell_type":"markdown","source":"But what if I didn't have the labels? I will try using k-means clustering to see if I can separate out the first 25 principal components into 9 clusters."},{"metadata":{"trusted":true,"_uuid":"804eb9ad66640bfbb8ddbb4b4ab9a5d1d4e000f3","scrolled":true},"cell_type":"code","source":"# Initial cluster assignments are random, set seed to ensure reproducibility\nset.seed(20)\ntrain_score_cluster <- as.data.frame(train_score)\nk <- kmeans(train_score_cluster, 9, nstart=20, iter.max=1000)\nsummary(k)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0b4bcb430983b719cd89c67b98096cf14a35ba84"},"cell_type":"markdown","source":"Let's have a look at the same plot as above but this time using labels from clustering, not the actual labels..."},{"metadata":{"trusted":true,"_uuid":"ffb572e9482f622d79846b9c9fbc2fd449f9653d"},"cell_type":"code","source":"plot(train_score_cluster$PC1,train_score_cluster$PC2,type=\"n\",main=\"First 2 PCs after clustering\")\ntext(train_score_cluster$PC1,train_score_cluster$PC2,label=k$cluster,col=colors[k$cluster])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a46d0ab0e40abf479ee9c0e2f908dc4b119fd8c6"},"cell_type":"markdown","source":"It looks very similar to the plot with labels (ignoring the actual numbers). This seems like the k-means clustering has done a fairly good job of separating out the different types of digit."}],"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}