A. Distance matrix and clustering

Extract all the unique keywords across photos (or subset of photos):

words <- unlist(photos[,3:12])
words <- words[!duplicated(words)] #list of unique keywords

Next, we convert photos into a binary format and name it wordscore, with each row representing a photo, and each column representing a keyword. “1” is added if the word is present. We then convert wordscore into a sparse matrix. This will help reduce the load on the computer’s RAM, especially if the photo dataset is very large.

#parallel loop:
wordscore <- foreach(i = 1:length(photos[,1]), .combine=rbind) %dopar% {
  vec <- vector(mode = "integer",length = length(words))
  a <- match(photos[i,3:12], words)
  vec[a] <- 1
  
  cat(paste0(" row ", i), file=paste0("admin/log_wordscore.txt"), append=TRUE) #The loop's progress will be printed in this file
  vec
}
colnames(wordscore) <- words
rownames(wordscore) <- NULL
wordscore <- wordscore[,!is.na(colnames(wordscore))] #remove 'NA' keyword if present

library(Matrix) 
wordscore <- Matrix(wordscore, sparse = TRUE)  #convert to sparseMatrix to save memory

In the binary format, wordscore can now be converted into a distance matrix. To have a fair assessment of the similarity (and thus the distance) between two photos, we need to take into account if they have the same number of keywords generated. The Jaccard Index is used in the calculation, where the number of common keywords is divided by the total number of unique keywords between two photos.

To start with, we find out how many keywords each photo has (up to ten), and save the results as the vector lengword:

narmlength <- function(x){10-sum(is.na(x))} #create function
lengword <- apply(photos[,3:12], 1, narmlength) #apply function

Next, the similarity between each photo and all other photos is calculated manually in a loop, based on the Jaccard Index. Since most photos do not share keywords, the similarity value will tend to be “0” (less strain on computer’s RAM). The similarity matrix (loop output) is then converted into a distance matrix, and subsequently converted into a ‘dist’ object.

simimat <- foreach(i = 1:length(wordscore[,1]), .packages = "Matrix", .combine=cbind) %dopar% {
  
  ws <- wordscore[,which(wordscore[i,] == 1)] #for each photo, find the other photos (rows) with its keywords (cols)
  
  simi <- round(apply(as.matrix(ws),1,sum, na.rm=T)/(lengword+lengword[i]),2) #Jaccard index
  
  simi[1:i] <- 0 #only fill half the matrix
  simi[i] <- 1
  
  cat(paste("row",i), file=paste0("admin/log_simimat.txt"), append=TRUE)
  
  simi
}

colnames(simimat) <- NULL
rm(wordscore, lengword)


#convert similarity to distance
distmat <- 1-simimat
rm(simimat)

#Convert to a 'dist' object
dm <- as.dist(distmat)

Finally, we perform hierarchical clustering of photos, using Ward’s distance:

require(fastcluster)
require(graphics)

cluz <- fastcluster::hclust(dm, "ward.D2")

Go to ‘B. How many clusters?’ if the number of photo categories has not been determined.

B. How many clusters?

This section runs as a separate analysis from the final results. Note that the following script may take a long time to run if you have a large dataset.

In this analysis, we measure the average difference between within- and between-cluster variation, across different clustering scenarios. Thus, a higher value suggests distinct clusters that more ‘different’ from each other (i.e. greater variation/distance between clusters). As the number of clusters (k) increases, this value is expected to decrease. We plot these values, and use the L-Method to find the ‘knee’ of the evaluation graph. More information about the L-Method can be found at:

Salvador, S. & Chan, P. Determining the Number of Clusters / Segments in Hierarchical Clustering / Segmentation Algorithms. in 16th IEEE International Conference on Tools with Artificial Intelligence 576-584 (IEEE, 2004). doi:10.1109/ICTAI.2004.50

First, decide up to how many clusters (k) to test for. In this example, we test k from 2 to 20, and save it as the vector scenarios (19 scenarios):

scenarios <- numeric(length(2:20))

Create a function to measure the difference between within- and between-cluster variation across all photos. Run the function for different k values in scenarios.

differ <- function(dist, gr, pos){
  #dist is a single photo's vector of distances with all others
  #gr is the vector output of grp membership across all photos
  #pos is the position of the single photo in length(distmat[1,])
  
  gr2 <-numeric(length(gr)) #vector of "0"s
  gr2[gr==gr[pos]] <-1 #Which photos are in same cluster as the photo of interest?
  gr3 <- 1-tapply(dist,gr2, mean) #2 values generated: (1) mean distance compared to photos in other clusters, & (2) compared to photos within same cluster. Minus values from one to convert to similarity value.
  gr3[2]-gr3[1] #within-cluster minus between-cluster similarity (larger value means clusters are very different)
}


#Run function for different scenarios (numbers of clusters):
for(i in 2:(length(scenarios)+1)){
  grp <- cutree(cluz, k=i)  #cutree returns vector of grp memberships across all photos
  
  cat("\n<< Working on scenario k =", i, "/", (length(scenarios)+1), ">>\n")
  
  alldiffer <-numeric(length(distmat[,1])) #vector of "0"s"
  
  for(j in 1:length(distmat[,1])){ #run function for each photo (across rows)
    alldiffer[j]<- differ(distmat[j,], grp, j )
    
    cat("<row", j, "/", length(distmat[,1]), "photos>")
    
  }
  scenarios[i-1]<-mean(alldiffer) #find out the mean difference for each scenario (k)
  
  cat("\n<< Scenario k =", i, "COMPLETE >>")
}


#Create dataframe
scenarios <- cbind.data.frame(seq(2,(length(scenarios)+1),1), 1-scenarios) #convert to distance
colnames(scenarios) <- c("k", "distance")

Do note that the small number of photos in our example produces a relatively straight curve. To help with visualisation, we can also calculate the marginal change in the distance:

for(i in 2:length(scenarios$distance)){
  scenarios[i,3] <- scenarios$distance[i-1]-scenarios$distance[i]
}

colnames(scenarios) <- c("k","distance","marginalDelta")

Here are plots of the results across different clustering scenarios:

Differences between within- and between-cluster variation, across different clustering scenarios

Since such plots may not always allow us to visually determine the appropriate number of photo clusters, we can also use the L-Method as described in Salvador and Chan (2004). To do so, we plot possible pairs of best-fit lines to the curve, and calculate the total root mean squared error (RMSE) for each pair. The lowest RMSE value is used to determine the number of clusters.

Note that for the ‘RMSE’ function that we use in this section, loading the packages requires XQuartz to be installed if you’re using a Mac.

require(rgl)
require(qpcR)

#Best-fit line equation:
mod1 <- lm(distance ~ k, data = scenarios)


#Equation from Salvador & Chan (2004):
for(i in 3:(max(scenarios$k)-2)){  #lowest value the 'knee' can be at is 3
  rmse <-  ((i-1)/(max(scenarios$k)-1)*(RMSE(mod1, which = 2:i))) + (((max(scenarios$k)-i)/(max(scenarios$k)-1))*RMSE(mod1, which = (i+1):max(scenarios$k)))
scenarios[i-1,4] <- rmse
}

colnames(scenarios) <- c("k","distance", "marginalDelta", "Lrmse")

Now we can plot RMSE across an increasing number of clusters (k). In our example, the lowest RMSE value where k = 11. This is the ‘knee’ of the graph. Note that there are a roughly balanced number of points on either side of this value.

Total RMSE of possible pairs of best-fit lines

Now it’s time to classify our photos and visualise the categories for good! Go back to ‘A. Distance matrix and clustering’ and re-run the script for the full dataset if a subset of data was used to determine the number of clusters. If not, continue on to the next section…