עד כה בקורס עסקנו בקשר שבין משתנים מסבירים (\(X_1,\ldots,X_p\)) לבין משתנה תלוי \(y\), ואיך מייצרים מודלים שונים שמתארים את הקשר שביניהם (בין אם מודלי רגרסיה, או מודלים מתקדמים יותר).

לעיתים הבעיה הניצבת בפנינו כחוקרים היא לא לתאר קשר בין משתנה תלוי למשתנים בלתי תלויים, אלא ללמוד על מבני הנתונים מתוך משתנים בלתי תלויים. לדוגמה:

בגנטיקה של אוכלוסיות לדוגמה, עוסקים בדיוק בבעיות כאלו: כיצד לסווג מינים שונים או כיצד להבין ממחקר גנטי על הקשר ההיסטורי שבין אוכלוסיות שונות. במקרה הבא, הגדרה של שני משתנים בלבד (PC1, PC2) מצליחה להפריד חלק מהאוכלוסיות באופן טוב, שדי מזכיר את מפת אירופה. PCA used in genetic studies https://stats.stackexchange.com/questions/8777/in-genome-wide-association-studies-what-are-principal-components

המפה מתארת איך ממגוון משתנים המתארים שינויים גנטיים, ניתן לחלץ שני משתנים בלבד המסבירים את מירב השונות בין אנשים שונים, ואת הדימיון הרב המתקבל בין משתנים אלו לבין מפת הפיזור של האנשים באירופה.

המינים הבאים סווגו לפי מידת הקרבה שלהן אחד לשני, בתרשים שנקרא Dendrogram. התרשים מסודר כעץ שבו ענפים כאשר ענפים סמוכים הם תצפיות שיש ביניהן דמיון. Clustering of species https://www.researchgate.net/figure/Dendrogram-showing-the-genetic-diversity-of-the-genomic-selection-training-population_fig2_317632929

סוג הניתוחים הללו רלוונטי למקרים בהם אנחנו מנסים לפענח כיצד פרטים שונים מתקבצים יחדיו באוכלוסיה (Clustering) או כיצד משתנים שונים “מתנהגים” בהשוואה אחד לשני (Principle Component Analysis).

שימו לב שיש הבדל מהותי בין יצירת אשכולות - המתייחס לפרטים שונים באוכלוסיה (תצפיות) לבין ניתוח גורמים (המתייחס למשתנים עצמם).

הרבה פעמים נשתמש בכלים אלו של unsupervised learning, כשלב ביניים לפני שנעבור למידול supervised.

נדון ראשית בניתוח גורמים (Principle Component Analysis - PCA).

ניתוח גורמים / Principle Component Analysis

בניתוח מסוג PCA, אנחנו משתמשים בכלים של אלגברה לינארית כדי לסובב ולהזיז את מערכת הצירים של הנתונים. הזזה זו מתבצעת באופן כזה, שבו מערכת הצירים החדשה היא מערכת צירים שבה הציר הראשון הוא בעל השונות הגבוהה ביותר, הציר השני בעל שונות פחות גבוהה וכן הלאה (הכוונה לפיזור הנתונים בכל ציר חדש).

אם אנחנו בוחרים תת-קבוצה של צירים אלו, המשמעות היא שאנחנו בוחרים תת-קבוצה שמסבירה “X% מהשונות שיש בנתונים”.

library(tidyverse)
## -- Attaching packages --------------------------------------------------------------------- tidyverse 1.2.1 --
## <U+221A> ggplot2 3.0.0     <U+221A> purrr   0.2.5
## <U+221A> tibble  1.4.2     <U+221A> dplyr   0.7.6
## <U+221A> tidyr   0.8.1     <U+221A> stringr 1.2.0
## <U+221A> readr   1.1.1     <U+221A> forcats 0.2.0
## -- Conflicts ------------------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
# lets try to reduce the dimension of the iris dataset
ggplot(iris, aes(Sepal.Length, Sepal.Width, color = Species)) + geom_point()

ggplot(iris, aes(Sepal.Length, Petal.Length, color = Species)) + geom_point()

ggplot(iris, aes(Sepal.Length, Petal.Width, color = Species)) + geom_point()

ggplot(iris, aes(Sepal.Width, Petal.Width, color = Species)) + geom_point()

ggplot(iris, aes(Sepal.Width, Petal.Length, color = Species)) + geom_point()

ggplot(iris, aes(Petal.Length, Petal.Width, color = Species)) + geom_point()

# It looks like every two axis have their interaction (or not)
# Some are very correlated, for example the Petal.Width and Petal.Length seem to have a very strong linear relationship.
# To a lesser extent, the same can be said for Petal.Length<->Sepal.Length, and Petal.Width<->Sepal.Length

# Now, lets run the PCA
iris_pca <- prcomp(iris %>% select(-Species))
iris_pca
## Standard deviations (1, .., p=4):
## [1] 2.0562689 0.4926162 0.2796596 0.1543862
## 
## Rotation (n x k) = (4 x 4):
##                      PC1         PC2         PC3        PC4
## Sepal.Length  0.36138659 -0.65658877  0.58202985  0.3154872
## Sepal.Width  -0.08452251 -0.73016143 -0.59791083 -0.3197231
## Petal.Length  0.85667061  0.17337266 -0.07623608 -0.4798390
## Petal.Width   0.35828920  0.07548102 -0.54583143  0.7536574

הפקודה מדווחת לנו על ארבעה רכיבים (תמיד מספר הרכיבים יהיה כמספר המשתנים הכללי שהכנסנו, במקרה זה ארבע משתנים).

כל משתנה מקבל “מקדם” שהוא המקדם שבו צריך להכפיל את המשתנה בשביל לקבל את ערך הרכיב החדש, הסכום שלהם הוא הרכיב החדש. לפעמים ניתן לתת פרשנות לרכיב החדש (כמעין ממוצע משוקלל), ולפעמים לא. הרבה פעמים התוצר של PCA הוא פשוט “קופסה שחורה”.

כאשר משתמשים בפקודת summary על אובייקט PCA, מקבלים את השונות של כל משתנה חדש, ואת שיעור השונות המוסברת (שלו והמצטברת).

summary(iris_pca)
## Importance of components%s:
##                           PC1     PC2    PC3     PC4
## Standard deviation     2.0563 0.49262 0.2797 0.15439
## Proportion of Variance 0.9246 0.05307 0.0171 0.00521
## Cumulative Proportion  0.9246 0.97769 0.9948 1.00000
# the following matrix contains the variables after they were rotated by the PCA
head(iris_pca$x)
##            PC1        PC2         PC3          PC4
## [1,] -2.684126 -0.3193972  0.02791483  0.002262437
## [2,] -2.714142  0.1770012  0.21046427  0.099026550
## [3,] -2.888991  0.1449494 -0.01790026  0.019968390
## [4,] -2.745343  0.3182990 -0.03155937 -0.075575817
## [5,] -2.728717 -0.3267545 -0.09007924 -0.061258593
## [6,] -2.280860 -0.7413304 -0.16867766 -0.024200858
# to rotate new data, just use the standard predict form
predict(object = iris_pca, newdata = iris %>% slice(1:5))
## Warning: package 'bindrcpp' was built under R version 3.4.4
##            PC1        PC2         PC3          PC4
## [1,] -2.684126 -0.3193972  0.02791483  0.002262437
## [2,] -2.714142  0.1770012  0.21046427  0.099026550
## [3,] -2.888991  0.1449494 -0.01790026  0.019968390
## [4,] -2.745343  0.3182990 -0.03155937 -0.075575817
## [5,] -2.728717 -0.3267545 -0.09007924 -0.061258593
# `slice` selects the rows 1:5
ggplot(as_tibble(iris_pca$x) %>%
         mutate(Species = iris$Species), aes(x = PC1, y = PC2, color = Species)) + 
  geom_point() + 
  coord_equal() + 
  ggtitle("The first two components of PCA on the iris dataset")

שימו לב לטווחים של ציר ה-PC1 לעומת הטווחים של ציר ה-PC2.

לעיתים, שלב ה-PCA מהווה שלב מקדים בשביל לבנות רגרסיה, והמשתנים המסבירים ברגרסיה הופכים להיות חלק מה-PCs.

תרגיל PCA

  1. קראו את נתוני המטופלים (קובץ נטישת מטופלים), ובצעו PCA על כל המשתנים שהם משתני 0-1 (לא כולל משתנה הנטישה).
  2. כמה משתנים צריך בשביל להסביר 80% מהשונות שבנתונים? (כמה רכיבי PC).
  3. בצעו PCA שוב, והפעם כללו גם את משתנה הגיל. מה ההשפעה של הוספת המשתנה לתוצאת ה-PCA?
  4. לפקודת PCA ישנם שני פרמטרים: scale. ו-center. פרמטרים אלו “מכינים” את הנתונים לפני ביצוע ה-PCA על ידי מרכוז ונרמול. הריצו שוב את פקודת ה-PCA כולל שימוש במרכוז ונרמול של המשתנים. האם קיבלתם תוצאה שונה? כיצד אתם מסבירים זאת?
  5. בנו מודל רגרסיה לוגיסטית המבוסס על הרכיבים החדשים, וחוזה את אי-ההגעה לפגישה. השוו בין מודל זה, לבין המודל המקורי (ללא PCA). מהם ביצוע המודל המבוסס על PCA לעומת המודל ללא PCA?

באפשרותכם להיעזר בקוד הבא:

appointments <- read_csv("data-files/Medical_Appointments_No_Shows_KaggleV2-May-2016.csv") %>%
  mutate(no_show = `No-show` == "Yes")

appointments_pca_all <- prcomp(appointments %>% 
                                 select(XXX),
                               XXX = XXX, XXX = XXX)
summary(appointments_pca_all)

# Add the results of PCA as new variables into a dataset
appointments_dataset_pca <- appointments %>% 
  bind_cols(as_tibble(appointments_pca_all$x))

# Try various options of # of PCs. When can you say that a smaller number of variables achieves similar results to using the full 7 variables?
appointments_pcaX_glm <- glm(formula = 
                            no_show ~ PC1 + PC2 + ... + PCX,
                          family = binomial,
                          data = appointments_dataset_pca %>% filter(is_train))

# To examine performance you can either look at the residuals or visualize using an ROC...

עד כה עסקנו בטיפול במשתנים, כעת נעסוק בטיפול בתצפיות באמצעות חלוקה לאשכולות (clustering).

חלוקה לאשכולות (KMeans clustering, Hierarchical clustering)

KMeans clustering

פקודת kmeans מחלקת את הנקודות במרחב לפי המרחק אחת מהשניה. במובן זה היא מזכירה את knn (שגם משתמשת במרחק בין נקודות), אך היא אינה נותנת תחזית או מתייחסת למשתנה תלוי כלשהו.

איך עובדת KMeans?

השיטה מנסה למצוא קבוצות אשר השונות בתוכן היא מינימלית (התצפיות בכל קבוצה דומות זו לזו), והשונות בין הקבוצות גבוהה (הקבוצות נמצאות רחוק אחת מהשניה).

נוסחה מקובלת לKMeans עושה שימוש במרחק אוקלידי בין הנקודות בכל קבוצה:

\[ W(C_k)=\frac{1}{|C_k|}\sum_{i,i'\in C_k}\sum_{j=1}^p(x_{ij}-x_{i'j})^2 \]

כאשר \(C_k\) הוא מספר התצפיות בקבוצה \(k\). נוסחה זו מייצגת את השונות בתוך אשכול \(k\).

פונקציית המטרה שלנו תהיה למצוא חלוקה לאשכולות \(C_1,\ldots,C_K\) כך שתקטין את סך השונות:

\[ \min_{C_1,\ldots,C_k}\sum_{k=1}^K{W(C_k)} \]

מי קובע את \(K\)? אתם.

איך? שאלה טובה.

לפעמים זה ניסוי וטעיה, לפעמים זה נגזר מדרישות או מטרות העבודה.

איך עובד האלגוריתם למציאת החלוקה \(C_1,\ldots, C_k\)?

  1. הגרל סיווג בין \(1,\ldots,K\) לכל אחת מהתצפיות.
  2. חזור על הפעולות הבאות עד שאין שינוי השמה של אשכול:
    1. לכל אשכול, חשב את הcentroid (נקודת מרכז המסה של האשכול, ממוצע על פני כלל המשתנים)
    2. שנה את הקצאות התצפיות לפי הcentroid אליו הם הכי קרובים.
# The original split of iris:
ggplot(iris, aes(Sepal.Width, Petal.Length, color = Species)) + geom_point()

# Now, illustrating KMeans on the two variables
iris_kmeans2 <- kmeans(iris %>% select(Sepal.Width, Petal.Length),
                       centers = 2)
iris_kmeans2
## K-means clustering with 2 clusters of sizes 51, 99
## 
## Cluster means:
##   Sepal.Width Petal.Length
## 1    3.409804     1.492157
## 2    2.875758     4.925253
## 
## Clustering vector:
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [71] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2
## [106] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [141] 2 2 2 2 2 2 2 2 2 2
## 
## Within cluster sum of squares by cluster:
## [1] 11.68196 74.62869
##  (between_SS / total_SS =  82.5 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"
iris_kmeans2$centers
##   Sepal.Width Petal.Length
## 1    3.409804     1.492157
## 2    2.875758     4.925253
glimpse(iris_kmeans2)
## List of 9
##  $ cluster     : int [1:150] 1 1 1 1 1 1 1 1 1 1 ...
##  $ centers     : num [1:2, 1:2] 3.41 2.88 1.49 4.93
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "1" "2"
##   .. ..$ : chr [1:2] "Sepal.Width" "Petal.Length"
##  $ totss       : num 493
##  $ withinss    : num [1:2] 11.7 74.6
##  $ tot.withinss: num 86.3
##  $ betweenss   : num 406
##  $ size        : int [1:2] 51 99
##  $ iter        : int 1
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"
iris_kmeans3 <- kmeans(iris %>% select(Sepal.Width, Petal.Length),
                       centers = 3)
iris_kmeans5 <- kmeans(iris %>% select(Sepal.Width, Petal.Length),
                       centers = 5)

iris_kmeans <- iris %>%
  mutate(kmeans2 = iris_kmeans2$cluster,
         kmeans3 = iris_kmeans3$cluster,
         kmeans5 = iris_kmeans5$cluster) %>%
  select(starts_with("kmeans"), Sepal.Width, Petal.Length) %>%
  gather(clusters, cluster_num, -Sepal.Width, -Petal.Length)

iris_kmeans_centers <- as_tibble(iris_kmeans2$centers) %>%
  rownames_to_column() %>%
  mutate(clusters = "kmeans2") %>%
  bind_rows(
    as_tibble(iris_kmeans3$centers) %>%
  rownames_to_column() %>%
  mutate(clusters = "kmeans3")
  ) %>%
  bind_rows(
    as_tibble(iris_kmeans5$centers) %>%
  rownames_to_column() %>%
  mutate(clusters = "kmeans5")
  ) %>%
  rename(cluster_num = rowname)

ggplot(iris_kmeans, aes(x = Sepal.Width, y = Petal.Length, color = factor(cluster_num))) +
  geom_point(size = 2) +
  facet_wrap(~ clusters) + 
  guides(color = guide_legend(title = "Cluster classification")) + 
  geom_point(data = iris_kmeans_centers, size = 5, shape = 21, 
             aes(fill = factor(cluster_num)), color = "black", alpha = 0.5) + 
  scale_fill_discrete(guide = "none")

# and here is the original, just for comparison
ggplot(iris, aes(Sepal.Length, Sepal.Width, color = Species)) + geom_point()

תרגיל

  1. השתמש ב-kmeans כדי לבנות קבוצות שונות של מטופלים בקובץ המטופלים. לשם כך השתמש במשתני ה-0-1 לא כולל ה-no-show. בחנו מספר אפשרויות שונות לערכי k.
  2. האם אתם מזהים מאפיינים כלשהם לקבוצות? לדוגמה האם הקבוצות מאובחנות ביניהן בפרמטרים מסוימים (יכול להיות פרמטרים שנכנסו למודל ה-KMeans או פרמטרים אחרים כגון גיל/מגדר שלא היו חלק מהחלוקה ל-KMeans).
  3. האם חלוקה זו לקבוצות יכולה לסייע כדי להבחין בין מטופלים שמגיעים לבין מטופלים שהולכים להבריז מפגישה?

כדי לבחון את התפלגות הפרמטרים, באפשרותכם להיעזר בקוד הבא:

appointments_kmeansXXX <- kmeans(appointments %>% select(Scholarship:SMS_received), centers = XXX)

appointments_cluster <- appointments %>%
  mutate(kmeansXXX_cluster = XXX) %>%
  mutate(male = XXX == XXX)

# The following code will help you show the averages of the different variables,
# by using the mechanism we talked about in previous units (summarize_at):

appointments_cluster %>%
  select(Scholarship:SMS_received, kmeansXXX_cluster, Age, male, no_show) %>%
  add_count(kmeansXXX_cluster) %>%
  group_by(kmeansXXX_cluster) %>%
  summarize_at(.vars = vars(1:6, Age, male, n, no_show, -kmeansXXX_cluster), 
               .funs = funs(mean(.)))

חלוקה לאשכולות בשיטת Hierarchical Clustering

שיטה זו משתמשת בעקרונות שונים משיטת KMeans clustering. במקום לחפש “מרכזי מסה”, ואז לסווג אילו תצפיות שייכות לכל מרכז מסה, השיטה ההירארכית מחפשת בכל נקודה אילו תצפיות “למזג” אחת לשניה לקבוצה אחת.

  1. התחל עם \(n\) תצפיות ופונקצית מרחק בין כל שתי תצפיות. סה“כ ישנם \(n*(n-1)/2\) מרחקים אפשריים. כל תצפית נחשבת לאשכול (האלגוריתם מתחיל עם \(n\) אשכולות).
  2. כל עוד מספר האשכולות גדול מ-1 בצע את הצעדים הבאים:
    1. בדוק את כל המרחקים בין האשכולות הקיימים. בחר את שני האשכולות שהכי קרובים זה לזה. מזג את שני האשכולות הללו. המרחק ביניהם ישמש אותנו כמדד ל“אי-הדימיון” ביניהם.
    2. חשב מחדש את המרחקים בין כל האשכולות, לפי פונקציית מרחק בין-אשכולות.

האלגוריתם המתואר פועל בשיטת bottom-up (מתחיל מהתצפית הבודדת). ישנם אלגוריתמים הפועלים בשיטה הפוכה (פיצול).

בהתייחס לשלב 2b באלגוריתם, ישנן פונקציות מרחק שונות. בדומה ל-KMeans אפשר להשתמש במרחק בין נקודות מרכז מסה (centroid), אך אפשר גם להשתמש במרחק בין שתי הנקודות הקרובות ביותר (single linkage), במרחק הממוצע (average linkage), ובמרחק מירבי (complete linkage).

# prepare the dataset with observation numbering
iris_prep <- iris %>% 
  mutate(obs = paste0(seq_along(Species), "-", Species)) %>%
  column_to_rownames("obs") %>%
  select(-Species)

# compute the distance matrix
iris_dist <- dist(iris_prep, method = "euclidean")

# generate the h-clustering
iris_clust <- hclust(iris_dist, method = "complete")

# We can generate the default plot which shows the tree
plot(iris_clust)

# The following will show us the order in which observations were merged:
iris_clust$order
##   [1] 108 131 103 126 130 119 106 123 118 132 110 136 141 145 125 121 144
##  [18] 101 137 149 116 111 148 113 140 142 146 109 104 117 138 105 129 133
##  [35] 150  71 128 139 115 122 114 102 143 135 112 147 124 127  73  84 134
##  [52] 120  69  88  66  76  77  55  59  78  87  51  53  86  52  57  75  98
##  [69]  74  79  64  92  61  99  58  94 107  67  85  56  91  62  72  68  83
##  [86]  93  95 100  89  96  97  63  65  80  60  54  90  70  81  82  42  30
## [103]  31  26  10  35  13   2  46  36   5  38  28  29  41   1  18  50   8
## [120]  40  23   7  43   3   4  48  14   9  39  17  33  34  15  16   6  19
## [137]  21  32  37  11  49  45  47  20  22  44  24  27  12  25
# The height object shows us a measure for the dissimilarity between the merged clusters
iris_clust$height
##   [1] 0.0000000 0.1000000 0.1000000 0.1000000 0.1000000 0.1000000 0.1414214
##   [8] 0.1414214 0.1414214 0.1414214 0.1414214 0.1414214 0.1414214 0.1414214
##  [15] 0.1414214 0.1414214 0.1414214 0.1414214 0.1414214 0.1414214 0.1414214
##  [22] 0.1732051 0.1732051 0.1732051 0.1732051 0.1732051 0.1732051 0.2000000
##  [29] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2236068 0.2236068
##  [36] 0.2449490 0.2449490 0.2449490 0.2449490 0.2449490 0.2449490 0.2449490
##  [43] 0.2449490 0.2645751 0.2645751 0.2645751 0.2645751 0.2645751 0.2645751
##  [50] 0.2645751 0.2645751 0.2645751 0.2828427 0.2828427 0.3000000 0.3000000
##  [57] 0.3000000 0.3000000 0.3000000 0.3162278 0.3162278 0.3162278 0.3316625
##  [64] 0.3316625 0.3316625 0.3316625 0.3316625 0.3316625 0.3464102 0.3464102
##  [71] 0.3464102 0.3464102 0.3605551 0.3605551 0.3741657 0.3741657 0.3741657
##  [78] 0.3872983 0.3872983 0.3872983 0.4000000 0.4000000 0.4123106 0.4123106
##  [85] 0.4242641 0.4242641 0.4358899 0.4472136 0.4472136 0.4582576 0.4582576
##  [92] 0.4690416 0.4690416 0.4690416 0.5099020 0.5196152 0.5196152 0.5196152
##  [99] 0.5196152 0.5196152 0.5477226 0.5477226 0.5477226 0.5567764 0.5567764
## [106] 0.5567764 0.6000000 0.6082763 0.6082763 0.6164414 0.6164414 0.6480741
## [113] 0.6480741 0.6557439 0.6782330 0.6855655 0.7071068 0.7211103 0.7348469
## [120] 0.7549834 0.7615773 0.7874008 0.8062258 0.8062258 0.8306624 0.9055385
## [127] 0.9219544 1.0049876 1.0099505 1.0295630 1.0677078 1.0677078 1.1090537
## [134] 1.1747340 1.2124356 1.2247449 1.3892444 1.4071247 1.4177447 1.4491377
## [141] 1.4525839 1.4628739 1.6613248 1.7058722 2.2360680 2.4289916 3.2109189
## [148] 4.0249224 7.0851958
# We can even cut the tree at a specific height or a specific number of clusters
cutree(iris_clust, k = 3)
##       1-setosa       2-setosa       3-setosa       4-setosa       5-setosa 
##              1              1              1              1              1 
##       6-setosa       7-setosa       8-setosa       9-setosa      10-setosa 
##              1              1              1              1              1 
##      11-setosa      12-setosa      13-setosa      14-setosa      15-setosa 
##              1              1              1              1              1 
##      16-setosa      17-setosa      18-setosa      19-setosa      20-setosa 
##              1              1              1              1              1 
##      21-setosa      22-setosa      23-setosa      24-setosa      25-setosa 
##              1              1              1              1              1 
##      26-setosa      27-setosa      28-setosa      29-setosa      30-setosa 
##              1              1              1              1              1 
##      31-setosa      32-setosa      33-setosa      34-setosa      35-setosa 
##              1              1              1              1              1 
##      36-setosa      37-setosa      38-setosa      39-setosa      40-setosa 
##              1              1              1              1              1 
##      41-setosa      42-setosa      43-setosa      44-setosa      45-setosa 
##              1              1              1              1              1 
##      46-setosa      47-setosa      48-setosa      49-setosa      50-setosa 
##              1              1              1              1              1 
##  51-versicolor  52-versicolor  53-versicolor  54-versicolor  55-versicolor 
##              2              2              2              3              2 
##  56-versicolor  57-versicolor  58-versicolor  59-versicolor  60-versicolor 
##              3              2              3              2              3 
##  61-versicolor  62-versicolor  63-versicolor  64-versicolor  65-versicolor 
##              3              3              3              2              3 
##  66-versicolor  67-versicolor  68-versicolor  69-versicolor  70-versicolor 
##              2              3              3              2              3 
##  71-versicolor  72-versicolor  73-versicolor  74-versicolor  75-versicolor 
##              2              3              2              2              2 
##  76-versicolor  77-versicolor  78-versicolor  79-versicolor  80-versicolor 
##              2              2              2              2              3 
##  81-versicolor  82-versicolor  83-versicolor  84-versicolor  85-versicolor 
##              3              3              3              2              3 
##  86-versicolor  87-versicolor  88-versicolor  89-versicolor  90-versicolor 
##              2              2              2              3              3 
##  91-versicolor  92-versicolor  93-versicolor  94-versicolor  95-versicolor 
##              3              2              3              3              3 
##  96-versicolor  97-versicolor  98-versicolor  99-versicolor 100-versicolor 
##              3              3              2              3              3 
##  101-virginica  102-virginica  103-virginica  104-virginica  105-virginica 
##              2              2              2              2              2 
##  106-virginica  107-virginica  108-virginica  109-virginica  110-virginica 
##              2              3              2              2              2 
##  111-virginica  112-virginica  113-virginica  114-virginica  115-virginica 
##              2              2              2              2              2 
##  116-virginica  117-virginica  118-virginica  119-virginica  120-virginica 
##              2              2              2              2              2 
##  121-virginica  122-virginica  123-virginica  124-virginica  125-virginica 
##              2              2              2              2              2 
##  126-virginica  127-virginica  128-virginica  129-virginica  130-virginica 
##              2              2              2              2              2 
##  131-virginica  132-virginica  133-virginica  134-virginica  135-virginica 
##              2              2              2              2              2 
##  136-virginica  137-virginica  138-virginica  139-virginica  140-virginica 
##              2              2              2              2              2 
##  141-virginica  142-virginica  143-virginica  144-virginica  145-virginica 
##              2              2              2              2              2 
##  146-virginica  147-virginica  148-virginica  149-virginica  150-virginica 
##              2              2              2              2              2
cutree(iris_clust, h = 1)
##       1-setosa       2-setosa       3-setosa       4-setosa       5-setosa 
##              1              1              2              2              1 
##       6-setosa       7-setosa       8-setosa       9-setosa      10-setosa 
##              3              2              1              2              1 
##      11-setosa      12-setosa      13-setosa      14-setosa      15-setosa 
##              3              4              1              2              5 
##      16-setosa      17-setosa      18-setosa      19-setosa      20-setosa 
##              5              5              1              3              4 
##      21-setosa      22-setosa      23-setosa      24-setosa      25-setosa 
##              3              4              2              4              4 
##      26-setosa      27-setosa      28-setosa      29-setosa      30-setosa 
##              1              4              1              1              1 
##      31-setosa      32-setosa      33-setosa      34-setosa      35-setosa 
##              1              3              5              5              1 
##      36-setosa      37-setosa      38-setosa      39-setosa      40-setosa 
##              1              3              1              2              1 
##      41-setosa      42-setosa      43-setosa      44-setosa      45-setosa 
##              1              6              2              4              4 
##      46-setosa      47-setosa      48-setosa      49-setosa      50-setosa 
##              1              4              2              3              1 
##  51-versicolor  52-versicolor  53-versicolor  54-versicolor  55-versicolor 
##              7              8              7              9              7 
##  56-versicolor  57-versicolor  58-versicolor  59-versicolor  60-versicolor 
##             10              8             11              7              9 
##  61-versicolor  62-versicolor  63-versicolor  64-versicolor  65-versicolor 
##             11             10             12              8              9 
##  66-versicolor  67-versicolor  68-versicolor  69-versicolor  70-versicolor 
##              7             10             10             13              9 
##  71-versicolor  72-versicolor  73-versicolor  74-versicolor  75-versicolor 
##             14             10             15              8              8 
##  76-versicolor  77-versicolor  78-versicolor  79-versicolor  80-versicolor 
##              7              7              7              8              9 
##  81-versicolor  82-versicolor  83-versicolor  84-versicolor  85-versicolor 
##              9              9             10             15             10 
##  86-versicolor  87-versicolor  88-versicolor  89-versicolor  90-versicolor 
##              8              7             13             10              9 
##  91-versicolor  92-versicolor  93-versicolor  94-versicolor  95-versicolor 
##             10              8             10             11             10 
##  96-versicolor  97-versicolor  98-versicolor  99-versicolor 100-versicolor 
##             10             10              8             11             10 
##  101-virginica  102-virginica  103-virginica  104-virginica  105-virginica 
##             16             14             17             18             18 
##  106-virginica  107-virginica  108-virginica  109-virginica  110-virginica 
##             19             20             17             18             21 
##  111-virginica  112-virginica  113-virginica  114-virginica  115-virginica 
##             22             15             22             14             14 
##  116-virginica  117-virginica  118-virginica  119-virginica  120-virginica 
##             22             18             23             19             13 
##  121-virginica  122-virginica  123-virginica  124-virginica  125-virginica 
##             16             14             19             15             16 
##  126-virginica  127-virginica  128-virginica  129-virginica  130-virginica 
##             17             15             14             18             17 
##  131-virginica  132-virginica  133-virginica  134-virginica  135-virginica 
##             17             23             18             15             15 
##  136-virginica  137-virginica  138-virginica  139-virginica  140-virginica 
##             21             16             18             14             22 
##  141-virginica  142-virginica  143-virginica  144-virginica  145-virginica 
##             16             22             14             16             16 
##  146-virginica  147-virginica  148-virginica  149-virginica  150-virginica 
##             22             15             22             16             14
# For better visualizations, here is a ggplot2 wrapper, also visualizing single and average linkage
#install.packages("ggdendro")
library(ggdendro)
## Warning: package 'ggdendro' was built under R version 3.4.4
ggdendrogram(hclust(iris_dist, method = "complete"))

ggdendrogram(hclust(iris_dist, method = "average"))

ggdendrogram(hclust(iris_dist, method = "single"))

שיטה זו דורשת כוח חישוב חזק, משום שנדרש לחשב בה פונקציית מרחק בין נקודות רבות, לעומת KMeans (שם המרחק מחושב רק מ“מרכזי מסה”).

תרגיל

  1. באמצעות פונקציית sample_n דגמו 5,000 תצפיות מקובץ הנתונים של המטופלים (למה מדגם ולא כולם?)
  2. בצעו קיבוץ ל-3 אשכולות באמצעות אלגוריתם hierarchical clustering, השוו את תוצאות קיבוץ זה למול הקיבוץ שנתקבל בתרגיל הקודם באמצעות אלגוריתם KMeans. האם הקיבוצים נותנים תוצאות דומות או שונות?
    1. זכרו! אם אתם משתמשים במשתנה Age, מומלץ לנרמל אותו, לדוגמה באמצעות הפונקציה scale.
    2. אם אתם משתמשים בפונקציה scale השתמשו גם בפונקציה as.numeric כדי להפוך את התוצאה לוקטורית ולא לרשימה.
  3. לצורך ההשוואה בסעיף הקודם, מומלץ לבדוק מספר פונקציות Linkage (נניח complete, single, average).
  4. איזה סוג קיבוץ מצליח ליצור קבוצות המובחנות זו מזו בהיבט המשתנים הדמוגרפיים מגדר וגיל?
# If you decide to use Age, don't forget to scale it otherwise it might pull all the results in unwanted direction - it will have a lot of weight on the Euclidean distance function
# You can either scale using the scale function (it uses mean and standard deviation) or scale to 0-1 range.

appointments_sample <- appointments %>%
   sample_n(5000) %>%
   select(Age, Scholarship:SMS_received) %>%
   mutate(age_scaled = as.numeric(scale(Age))) %>%
   select(-Age)

# another option is to use a 0-1 scaling as follows
appointments_sample <- appointments %>%
   sample_n(5000) %>%
   select(Age, Scholarship:SMS_received) %>%
   mutate(age_scaled = (Age - min(Age))/(max(Age)-min(Age))) %>%
   select(-Age)

# Now continue to use dist() and hclust() as in the example above.