1 Cluster Analysis

1.1 Important Characteristics of Our Data - 0.6 points

library(readr)
library(ggplot2)
library(dplyr)
library(rvest)
library(stringr)
library(tidytext)
library(tidyr)
library(ggwordcloud)
library(psych)
library(knitr)
library(kableExtra)
library(factoextra)
library(cluster)

dataforproject2 <- read_csv("/srv/store/students/vvsuschevskiy/datanal/4year/clusters/dataforproject2.csv")

## Warning: Missing column names filled in: 'X1' [1]

head(dataforproject2[,-1:-2]) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

dataforproject2 %>% 
  filter(str_detect(LanguageWorkedWith, "R")) %>%
  mutate(LearnR = ifelse(str_detect(LanguageDesireNextYear, "R"), "WantR", "Nope")) -> Use_R

Use_R[sapply(Use_R, is.character)] <- lapply(Use_R[sapply(Use_R, is.character)], as.factor)

#describeBy(Use_R, group = Use_R$LearnR)

names_text = names(Use_R[,c(3:7, 9:11, 13, 27, 35,  37:40)])

b = 1
for( i in names(Use_R[,c(3:7, 9:11, 13, 27, 35,  37:40) ] )){
  print(
    Use_R %>% 
    ggplot(aes_string(x = i))+
    geom_bar(aes(fill = LearnR),color = "black", stat="count", position = "dodge")+
      ggtitle(names_text[b])+
      geom_text(aes( label =  paste0(round((..count..)/sum(..count..)*100), "%"),
                   y=  (..count..)/sum(..count..)), stat= "count", vjust = -.5)+
      theme_minimal()+
      scale_x_discrete(guide = guide_axis(n.dodge = 3))
  )
  b = b + 1
}

R users are mostly developers by profession, or at least develop something as a part of their work. Also, they usually code for a hobby, and there is no big difference between people who code in R and people who want to learn R. However, that might be a proble of the data, because it was pre-filltered.

1.2 Justification of Variables (you may want to use not all the variables) - 1 point

To justify of variables lets take a look what is happening with R right now. And collect data from UseR 2019. Unfortunately, I founf only keynotes not proceedings, but I am not ready to care too much.

first_page <- read_html("https://user2019.r-project.org/program/")

first_page %>% 
  html_nodes(".speaker-bio , ul") %>% 
  html_text() %>% 
  as_tibble() -> wikipedia_text

## Warning: Calling `as_tibble()` on a vector is discouraged, because the behavior is likely to change in the future. Use `tibble::enframe(name = NULL)` instead.
## This warning is displayed once per session.

head(wikipedia_text)%>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

      Program
      Program overview
        
        Talk schedule
        
        Important dates
        
        Keynotes
        
        Tutorials
        
        Datathon
        
        Information for presenters 
        
        Social Program
        
        Side events
        
        Posters
        
      
    
    
    
    
      Registration
      Registration
        
        Abstract submission
        
        Scholarships
        
      
    
    
    
    
      Venue
      Travel
        
        Toulouse
        
        Gala dinner
        
        Around Toulouse
        
      
    
    
    
    
      About
      Local organization committee
        
        Scientific committee
        
        Past events
        
        FAQ
        
        Carbon footprint
        
        Legal information
        
      
    
    
    
    
      Code of Conduct
    
    
    
    
    
      Contact </td>

        Talk schedule
        
        Important dates
        
        Keynotes
        
        Tutorials
        
        Datathon
        
        Information for presenters 
        
        Social Program
        
        Side events
        
        Posters </td>

        Abstract submission
        
        Scholarships </td>

        Toulouse
        
        Gala dinner
        
        Around Toulouse </td>

        Scientific committee
        
        Past events
        
        FAQ
        
        Carbon footprint
        
        Legal information </td>

stop_words = get_stopwords("en")

stop_words = rbind(stop_words, c("packages"))

draw_wc = function(text){
  text %>%
  unnest_tokens(bigram, value, token = "skip_ngrams", n = 2, k = 5) %>%
  separate(bigram, c("word1", "word2"), sep = " ") %>%
  filter(!word1 %in% stop_words$word) %>%
  filter(!word2 %in% stop_words$word) %>%
  filter(!word1 %in% c(1:20)) %>%
  filter(!word2 %in% c(1:20)) %>% 
  count(word1, word2, sort = TRUE) %>% na.omit() %>% 
  mutate(n = n, word = paste(word1, word2, sep = " ")) %>% 
  select(word, n) %>% 
  filter(n > 1) %>% 
  ggplot(aes(label = word, size = n)) +
  geom_text_wordcloud(rm_outside = T) +
  scale_size_area(max_size = 20) +
  scale_color_manual(values = c("red", "skyblue", "black"))+
  theme_minimal()
}

draw_wc(wikipedia_text)

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 2954 rows [1,
## 8, 15, 22, 29, 36, 43, 50, 57, 64, 71, 78, 85, 92, 99, 106, 113, 120, 127,
## 134, ...].

In 2019 R was mostly about DS, Science, Statistics, creating applications, little bit about AI and education. So lets try to find different R users.

Thus, there are different reasons to code in R, so I will try to capture this reasons with following variables:

1. OpenSourcer – because people who work in idustry usually could not contribute
2. Employment – work is an important for any R user, and part time workers might be more interesting
3. Students – just to control for students
4. YearsCode – because R is not so old
5. JobFactors – beacuse that is an answer for my question
6. WorkPlan – cool programmers do not have shedule

1.3 the distance metric matches variable types - 1 point (if this is incorrect, interpretation will fail)

Use_R %>% 
  select(Respondent, OpenSourcer, Employment, Student, YearsCode, JobFactors, WorkPlan) %>% na.omit() -> df_cluster

1.4 k-means

df_cluster_num <- mutate_all(df_cluster, function(x) as.numeric(x))


fviz_nbclust(df_cluster_num[,-1], kmeans, method = "wss")

fviz_nbclust(df_cluster_num[,-1], kmeans, method = "silhouette")

fviz_nbclust(df_cluster_num[,-1], kmeans, method = "gap")

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 185350)

## Warning: did not converge in 10 iterations

wss and silhouette say that 2 clusters is an optimal amout for this data, but when we are talking about our data, we should take into consideration, that convert factors into numeric should be illigal.

clusters <- kmeans(df_cluster_num[,-1],
                      2, # how many groups to locate
                      nstart = 20 # R will try 20 different random starting assignments 
                      # and then select the one with the lowest within cluster variation.
)


fviz_cluster(clusters, data = df_cluster_num[,-1],
  ellipse.type = "convex",
  palette = "jco",
  repel = TRUE)

As we see our data could not be properly plotted, bacause our dementions do not take enough load.

so we need to create a proper distance matrix

I will use gower, since I have mostly factors

Use_R %>% 
  select(Respondent, OpenSourcer, Employment, Student, YearsCode, JobFactors, WorkPlan) %>% mutate_if(is.factor, addNA) -> df_cluster

df_cluster$YearsCode = as.numeric(as.character(df_cluster$YearsCode))

## Warning: NAs introduced by coercion

# to perform different types of hierarchical clustering
# package functions used: daisy(), diana(), clusplot()
gower.dist <- daisy(df_cluster[ ,-1], metric = c("gower"))
#class(gower.dist) 
## dissimilarity , dist

1.5 DIVISIVE

divisive.clust <- diana(as.matrix(gower.dist),
                  diss = TRUE, keep.diss = F)

?diana
plot(divisive.clust)

works too bad, could not compute((

1.6 AGGLOMERATIVE

aggl.clust.c <- hclust(gower.dist, method = "complete")

plot(aggl.clust.c,
     main = "Agglomerative, complete linkages", hang = -1, cex = 0.6, ylim =c(0.5,1))

1.7 Vizualization

 #install.packages("ape")
library("ape")
# Default plot
plot(as.phylo(aggl.clust.c), cex = 0.6, label.offset = 0.5)

plot(as.phylo(aggl.clust.c), type = "unrooted", cex = 0.6,
     no.margin = TRUE, show.tip.label = F)

It looks that we have around 12 clusters what is not good, but not really bad, lets cut them, and try to look at them

clus4 = cutree(aggl.clust.c, 12)

plot(as.phylo(aggl.clust.c), type = "fan", 
     label.offset = 1, cex = 0.3, show.tip.label = F)

lets put these clusters back to data

Use_R$cluster = clus4 

and draw all graphs

names_text = names(Use_R[,c(3:7, 9:11, 13, 27, 35,  37:40)])

b = 1
for( i in names(Use_R[,c(3:7, 9:11, 13, 27, 35,  37:40) ] )){
  print(
    Use_R %>% 
    ggplot(aes_string(x = i))+
    geom_bar(aes(fill = as.factor(cluster)), color = "black", stat="count", position = "dodge")+
      ggtitle(names_text[b])+
      geom_text(aes( label =  paste0(round((..count..)/sum(..count..)*100), "%"),
                   y=  (..count..)/sum(..count..)), stat= "count", vjust = -.5)+
      theme_minimal()+
      scale_x_discrete(guide = guide_axis(n.dodge = 3))
  )
  b = b + 1
}

Ok, I could not really understand why are they are different, we need a computer to make this decision for me. All credits to https://towardsdatascience.com/hierarchical-clustering-on-categorical-data-in-r-a27e578f2995

# Cluster stats comes out as list while it is more convenient to look at it as a table
# This code below will produce a dataframe with observations in columns and variables in row
# Not quite tidy data, which will require a tweak for plotting, but I prefer this view as an output here as I find it more comprehensive 
library(fpc)
cstats.table <- function(dist, tree, k) {
clust.assess <- c("cluster.number","n","within.cluster.ss","average.within","average.between",
                  "wb.ratio","dunn2","avg.silwidth")
clust.size <- c("cluster.size")
stats.names <- c()
row.clust <- c()
output.stats <- matrix(ncol = k, nrow = length(clust.assess))
cluster.sizes <- matrix(ncol = k, nrow = k)
for(i in c(1:k)){
  row.clust[i] <- paste("Cluster-", i, " size")
}
for(i in c(2:k)){
  stats.names[i] <- paste("Test", i-1)
  
  for(j in seq_along(clust.assess)){
    output.stats[j, i] <- unlist(cluster.stats(d = dist, clustering = cutree(tree, k = i))[clust.assess])[j]
    
  }
  
  for(d in 1:k) {
    cluster.sizes[d, i] <- unlist(cluster.stats(d = dist, clustering = cutree(tree, k = i))[clust.size])[d]
    dim(cluster.sizes[d, i]) <- c(length(cluster.sizes[i]), 1)
    cluster.sizes[d, i]
    
  }
}
output.stats.df <- data.frame(output.stats)
cluster.sizes <- data.frame(cluster.sizes)
cluster.sizes[is.na(cluster.sizes)] <- 0
rows.all <- c(clust.assess, row.clust)
# rownames(output.stats.df) <- clust.assess
output <- rbind(output.stats.df, cluster.sizes)[ ,-1]
colnames(output) <- stats.names[2:k]
rownames(output) <- rows.all
is.num <- sapply(output, is.numeric)
output[is.num] <- lapply(output[is.num], round, 2)
output
}
# I am capping the maximum amout of clusters by 7
# I want to choose a reasonable number, based on which I will be able to see basic differences between customer groups as a result
stats.df.divisive <- cstats.table(gower.dist, divisive.clust, 7)
stats.df.divisive

##                    Test 1  Test 2  Test 3  Test 4  Test 5  Test 6
## cluster.number       2.00    3.00    4.00    5.00    6.00    7.00
## n                 5048.00 5048.00 5048.00 5048.00 5048.00 5048.00
## within.cluster.ss  780.28  732.45  643.05  628.96  626.63  592.96
## average.within       0.53    0.51    0.48    0.48    0.48    0.46
## average.between      0.71    0.70    0.71    0.71    0.71    0.71
## wb.ratio             0.75    0.73    0.68    0.67    0.67    0.66
## dunn2                1.13    0.98    1.08    1.09    1.10    1.16
## avg.silwidth         0.25    0.18    0.21    0.21    0.21    0.22
## Cluster- 1  size  2934.00 2934.00 2934.00 2934.00 2934.00 2934.00
## Cluster- 2  size  2114.00 1679.00 1342.00 1299.00 1293.00 1094.00
## Cluster- 3  size     0.00  435.00  435.00  435.00  435.00  435.00
## Cluster- 4  size     0.00    0.00  337.00  337.00  337.00  337.00
## Cluster- 5  size     0.00    0.00    0.00   43.00   43.00  199.00
## Cluster- 6  size     0.00    0.00    0.00    0.00    6.00   43.00
## Cluster- 7  size     0.00    0.00    0.00    0.00    0.00    6.00

stats.df.aggl <- cstats.table(gower.dist, aggl.clust.c, 7)
stats.df.aggl

##                    Test 1  Test 2  Test 3  Test 4  Test 5  Test 6
## cluster.number       2.00    3.00    4.00    5.00    6.00    7.00
## n                 5048.00 5048.00 5048.00 5048.00 5048.00 5048.00
## within.cluster.ss 1001.66  981.15  924.83  920.76  897.20  864.59
## average.within       0.60    0.59    0.57    0.57    0.56    0.55
## average.between      0.67    0.69    0.64    0.64    0.65    0.65
## wb.ratio             0.90    0.87    0.90    0.89    0.87    0.85
## dunn2                1.11    1.11    1.01    0.98    0.98    0.98
## avg.silwidth         0.10    0.07   -0.01   -0.01   -0.02   -0.02
## Cluster- 1  size  4913.00 4785.00 4259.00 4234.00 4124.00 3826.00
## Cluster- 2  size   135.00  128.00  526.00   25.00  110.00  298.00
## Cluster- 3  size     0.00  135.00  128.00  526.00   25.00  110.00
## Cluster- 4  size     0.00    0.00  135.00  128.00  526.00   25.00
## Cluster- 5  size     0.00    0.00    0.00  135.00  128.00  526.00
## Cluster- 6  size     0.00    0.00    0.00    0.00  135.00  128.00
## Cluster- 7  size     0.00    0.00    0.00    0.00    0.00  135.00

and plot them

stats.df.divisive["method",] = "divi"
stats.df.aggl["method",] = "aggl"

df_elbow = as.data.frame(rbind(t(stats.df.divisive),t(stats.df.aggl)))

ggplot(data = df_elbow, aes(x=cluster.number, y=within.cluster.ss, group = method)) + 
  geom_point()+
  geom_line(aes(color = method))+
  ggtitle("clustering") +
  labs(x = "Num.of clusters", y = "Within clusters sum of squares (SS)") +
  theme(plot.title = element_text(hjust = 0.5))

according to our picture there is no elbow, (probably because it happen later, but I could no dercribe more than 7 clusters)

ggplot(data = df_elbow, aes(x=cluster.number, y=avg.silwidth, group = method)) + 
  geom_point()+
  geom_line(aes(color = method))+
  ggtitle(" clustering") +
  labs(x = "Num.of clusters", y = "Average silhouette width") +
  theme(plot.title = element_text(hjust = 0.5))

We could se it here, but for JK, I could not run that script with test one more time, my PC will explodes. & is fine, just fine. Divi is better in all metrics

clus7 = cutree(divisive.clust, 7)

Use_R$cluster = clus7 

df_cluster$cluster = clus7

names_text = names(df_cluster[,c(-1,-6,-5)])

b = 1
for( i in names(df_cluster[,c(-1,-6,-5)] )){
  print(
    Use_R %>% 
    ggplot(aes_string(x = i))+
    geom_bar(aes(fill = as.factor(cluster)), color = "black", stat="count", position = "dodge")+
      ggtitle(names_text[b])+
      geom_text(aes( label =  paste0(round((..count..)/sum(..count..)*100), "%"),
                   y=  (..count..)/sum(..count..)), stat= "count", vjust = -.5)+
      theme_minimal()+
      scale_x_discrete(guide = guide_axis(n.dodge = 3))
  )
  b = b + 1
}

df_cluster %>% 
ggplot(aes(YearsCode, fill = as.factor(cluster)))+
  geom_histogram()+
  theme_minimal()+
  facet_grid(. ~ cluster)

## Warning: Removed 112 rows containing non-finite values (stat_bin).

1.8 Cluster names

1. cluster 1 is mostly developers, or people who develop, they have enough time for hobbies, definately not a students. They work for sure. mostly from office

2. cluster 2 students, but or just graduated

3. • 5 office workers, but develompent is not their primary job

6- 7. mostly noise. less then 1% of data

df_cluster$JobFactors = as.factor(df_cluster$JobFactors)

Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

df_cluster %>% group_by(cluster) %>% summarise(mode = Mode(JobFactors))

## # A tibble: 7 x 2
##   cluster mode                                                                  
##     <int> <fct>                                                                 
## 1       1 Languages, frameworks, and other technologies I'd be working with;Off…
## 2       2 <NA>                                                                  
## 3       3 Office environment or company culture;Opportunities for professional …
## 4       4 Office environment or company culture;Opportunities for professional …
## 5       5 Office environment or company culture;Opportunities for professional …
## 6       6 Languages, frameworks, and other technologies I'd be working with;Off…
## 7       7 Industry that I'd be working in;Financial performance or funding stat…