V04_Distribution.Rmd

```{r include=FALSE}
knitr::opts_chunk$set(cache = T, warning = F, message = F, 
                      class.output = "output", out.width='100%', 
                      fig.asp = 0.5, fig.align = 'center')
library(tidyverse)
options(scipen = 999)
```

# DISTRIBUTION: Histogram & Density

本章節將介紹與資料分布相關的視覺化方法。資料分布是指數據中每個值出現的頻率或概率。在統計學中，了解資料分布是非常重要的，因為它可以幫助我們判斷數據是否為正態分佈，或者是否存在異常值或極端值。本章節將涵蓋常見的資料分布視覺化方法，包括直方圖、密度圖、箱形圖和金字塔圖等。

以下是R語言ggplot2套件中，用於資料分布視覺化的一些常用函式：

-   `geom_histogram()`：用於創建直方圖。

-   `geom_density()`：用於創建密度圖。

-   `geom_boxplot()`：用於創建箱形圖。

-   `geom_bar()`：用於創建柱狀圖。

-   `geom_freqpoly()`：用於創建頻率多邊形圖。

註：本節的設計概念不少是參考 Claus O. Wilke 所著的「Foundations of Data Visualization」一書的章節，同時也參考臺灣和資料新聞的案例進行了改編。

```{r include=FALSE}
library(printr)
library(knitr)
th <- theme_minimal() + 
  theme(title = element_text(family="Heiti TC Light"),
        text = element_text(family="Heiti TC Light"), 
        axis.text.y = element_text(family="PingFang TC"),
        axis.text.x = element_text(family="Heiti TC Light"),
        legend.text = element_text(family="Heiti TC Light"),
        plot.title = element_text(family="Heiti TC Light")
        )
```

接下來我們將使用Histogram和Density Plot這兩種資料視覺化方法來探索台灣村里長的年齡和性別分布情況。我們所使用的資料來源包括內政部和中選會的投票資料，這些資料能夠提供具有代表性的統計樣本，幫助我們更好地了解村里長的整體特徵。在進行資料視覺化的過程中，我們將會運用R語言中的ggplot2套件，並根據不同的視覺化需求進行相應的設置和調整。 <https://www.moi.gov.tw/LocalOfficial.aspx?n=577&TYP=KND0007>。

```{r Fig.v1.clean-data}





vilmaster <- readr::read_csv("data/tw_vil2018_elccand.csv") %>%
  drop_na(當選註記)
```

## Density plot

密度圖（Density Plot)是一種展示數據集分佈情況的圖表，它可以幫助我們更好地理解數據集中數值出現的概率。圖表的 X 軸代表數據集的數值範圍，Y 軸則代表每個數值的出現概率。與直方圖不同，密度圖的曲線是光滑的，因為它是通過連續的數值範圍估算出的概率密度函數。通過比較不同數據集的密度圖，我們可以更好地了解它們之間的差異。在ggplot2中，可以用`geom_density()`函數來創建密度圖。

```{r Fig.v1.geom_density, fig.asp=0.4, message=FALSE, warning=FALSE, out.width='100%'}
p1 <- vilmaster %>%
  ggplot()  + aes(年齡) + 
  geom_density() + th

p2 <- vilmaster %>%
  ggplot()  + aes(年齡, fill=factor(性別)) + 
  geom_density(alpha=0.5) + th + 
  scale_fill_manual(
    limits=c('1', '2'),           # original chart group
    values=c("gold", "skyblue"),  # map to color
    name="性別",                  # legend title
    breaks=c(1, 2),               # original legend group labels
    labels=c("Male","Female"),    # map to new labels
    na.value = "lightgrey"        # color for other groups
    )

cowplot::plot_grid(
  p1, p2,
  labels = c("(a) Overall", "(b) Group by gender"),
  nrow = 1, rel_widths = c(1, 1)
)

```

### Density with different bandwidth

參數`bw`指的是bnadwidth，為繪製histogram時的bar所涵蓋的資料寬度。以step-plot來說，`bw`越大，則梯距越寬；以density-plot來說，若bw越大則越是平滑。

```{r Fig.v1.geom_density_line, fig.asp=1}
library(ggridges) # for geom_density_line()
p.b05 <- vilmaster %>% ggplot()  + aes(年齡) + 
  geom_density_line(fill='gold', bw=0.5, kernel='gaussian') + th

p.b1 <- vilmaster %>% ggplot()  + aes(年齡) + 
  geom_density_line(fill='gold', bw=1, kernel='gaussian') + th

p.b5 <- vilmaster %>% ggplot()  + aes(年齡) + 
  geom_density_line(fill='gold', bw=5, kernel='gaussian') + th

p.rect <- vilmaster %>% ggplot()  + aes(年齡) + 
  geom_density_line(fill='gold', bw=10, kernel='rectangular') + th

cowplot::plot_grid( p.b05, p.b1, p.b5, p.rect,
  labels = c("(a) bw=.5", "(b) bw=1", "(c) bw=2", "(b) rect"),
  nrow = 2, rel_widths = c(1, 1)
)
```

## Histogram

直方圖（Histogram）是一種用於展示數據集分佈的圖表。它通過將數據範圍分成若干個區間（稱為 "bins" 或 "buckets"），然後計算落在每個區間內的數據的數量（稱為 "frequency"），來展示數據集的分佈情況。直方圖的 X 軸表示數據範圍，Y 軸表示每個區間中的頻數。直方圖可以幫助我們快速了解數據的分佈情況，特別是數據的中心趨勢、數據的離散程度和是否存在異常值等。

### Histogram with different number of bins

```{r Fig.v1.geom_histogram, fig.asp=1}
p10 <- vilmaster %>%
  ggplot()  + aes(年齡) + 
  geom_histogram(bins=10, fill='royalblue') + th

p20 <- vilmaster %>%
  ggplot()  + aes(年齡) + 
  geom_histogram(bins=20, fill='royalblue') + th

p30 <- vilmaster %>%
  ggplot()  + aes(年齡) + 
  geom_histogram(bins=30, fill='royalblue') + th

p40 <- vilmaster %>%
  ggplot()  + aes(年齡) + 
  geom_histogram(bins=40, fill='royalblue') + th

cowplot::plot_grid(
  p10, p20, p30, p40,
  labels = c("(a) bins=10", "(b) bins=20", "(c) bins=30", "(b) bins=40"),
  nrow = 2, rel_widths = c(1, 1)
)
```

### Density vs histogram

Histogram通常用來顯示數據的分佈情況，它會把數據區間分成若干個等寬的區間，然後計算每個區間內數據的頻率，再將這些頻率表示在y軸上。因此，histogram顯示的是數據的頻率，而不是數據的密度。

Density plot則是用來顯示數據的概率密度函數，它會通過核密度估計（Kernel Density Estimation, KDE）方法，將數據點周圍的密度估計出來，然後將這些估計值表示在y軸上。因此，density plot顯示的是數據的密度，而不是數據的頻率。

```{r fig.asp=0.5, warning=F}
pd <- vilmaster %>%
  ggplot()  + aes(年齡, fill=factor(性別)) + 
  geom_density(alpha=0.5) + th + 
  scale_fill_manual(
    values=c("1"='gold', '2'="skyblue"),
    labels=c('1'="Male",'2'="Female"),
    name='Sex'
    )

ph <- vilmaster %>%
  ggplot()  + aes(年齡, fill=factor(性別)) + 
  geom_histogram(bins=20, position="dodge") + th + 
  scale_fill_manual(values=c("1"='gold', '2'="skyblue ")) + 
  theme(legend.position="none")

cowplot::plot_grid(
  pd, ph,
  labels = c("(a) geom_density()", "(b) geom_histogram()"),
  nrow = 1, rel_widths = c(6, 4)
)
```

### Positions of bar chart

```{r Fig.v1.hist.style}
p.hist.dodge <- vilmaster %>%
  ggplot()  + aes(年齡, fill=factor(性別)) + 
  geom_histogram(bins=20, position="dodge") + th + 
  scale_fill_manual(
    values=c("1"='gold', '2'="skyblue "),
    labels=c('1'="Male",'2'="Female"),
    name='Sex'
  )

p.hist.stack <- vilmaster %>%
  ggplot()  + aes(年齡, fill=factor(性別)) + 
  geom_histogram(bins=20, position="stack") + th +
  scale_fill_manual(values=c("1"='gold', '2'="skyblue ")) + 
  theme(legend.position="none")

cowplot::plot_grid(
  p.hist.dodge, p.hist.stack,
  labels = c("(a) position:dodge", "(b) position:stack"),
  nrow = 1, rel_widths = c(6, 4)
)
```

### Display two groups histogram by facet_wrap()

1.  `geom_histogram(bins=20, position="dodge")` 用於繪製直方圖， `bins=20`表示將數據分成20個區間， `position="dodge"`表示將不同性別的數據分開顯示。

2.  `th` 是本範例在最早先所建立的ggplot主題，用於設置圖表的樣式（例如背景顏色、字體等）。

3.  `scale_fill_manual()` 用於手動設置填充顏色， `values=c("1"='gold', '2'="skyblue")` 表示性別為1時填充金色，性別為2時填充天藍色。 `labels=c('1'="Male",'2'="Female")` 表示將性別1標記為Male，性別2標記為Female。 `name='Sex'` 表示設置顏色圖例的標題為`Sex`。

4.  `facet_wrap(.~性別, nrow=1)` 表示將不同性別的數據分開顯示，每直行顯示一個性別。`.~性別` 表示將數據按性別分組。

```{r Fig.v1.hist.facet, fig.asp=0.5}
vilmaster %>%
  ggplot()  + aes(年齡, fill=factor(性別)) + 
  geom_histogram(bins=20, position="dodge") + th + 
  scale_fill_manual(
    values=c("1"='gold', '2'="skyblue "),
    labels=c('1'="Male",'2'="Female"),
    name='Sex'
  ) + 
  facet_wrap(.~性別, nrow=1)
```

## Pyramid Plot {#pyramid}

金字塔圖（Pyramid plot）是一種用於比較兩個群體的統計圖表。它的形狀像一座金字塔，通常用於展示男女或年齡分佈等相關的數據。金字塔圖以垂直線為軸線，其中一側代表一個群體（如男性），另一側代表另一個群體（如女性）。圖表的左右兩側是對稱的，並以一條中心線分開。圖表中的每一行表示一個年齡段，而每一列則表示一個群體的比例或頻數。金字塔圖的高度表示總人數或總比例，並且可以用不同的顏色區分不同的群體。金字塔圖可以直觀地顯示兩個群體之間的差異，特別是在不同年齡段之間。

### Modify geom_col() to pyramid plot

```{r Fig.v1.hist.pyramid}

vilmaster %>%
  group_by(性別) %>%
  mutate(age_group = cut(年齡, 0:20*5+.01)) %>%
  count(age_group) %>%
  ungroup() %>%
  ggplot()  + aes(x=age_group, 
                  y=ifelse(性別=='1', -1, 1)*n,
                  fill=factor(性別)) + 
  geom_col() + 
  scale_y_continuous(name = "Count", breaks = 250*(-6:2), labels = c("1500", "1250", "1000", "750", "500", "250", "0", "250", "500")) + 
  coord_flip() + 
  scale_fill_manual(
    values=c("1"='gold', '2'="skyblue "),
    labels=c('1'="Male",'2'="Female"),
    name='Sex'
  ) + th + labs(y="Count", x="Age Group")
```

## Box plot: Muitiple Distrubution

箱形圖（Box plot）是一種用於展示數據分佈情況的統計圖表。它通常顯示數據的中位數、四分位數、極值和異常值等統計量。箱形圖的中間線表示數據的中位數，箱子的上下邊界則分別表示數據的上四分位數和下四分位數。箱子的高度表示數據的變異程度，而箱子外的線段則表示數據的最大值和最小值。如果數據中存在異常值，則通常使用圓圈或星號等符號來標記。箱形圖可以用來比較不同數據集之間的分佈情況，以及檢查數據是否存在異常值。

### TW-Salary (boxplot) {#twsalary}

Inspired by [Six Myths About Choosing a College Major - The New York Times (nytimes.com)](https://www.nytimes.com/2017/11/03/education/edlife/choosing-a-college-major.html) and [What's Going On in This Graph? \| Jan. 9, 2018 - The New York Times (nytimes.com)](https://www.nytimes.com/2018/01/04/learning/whats-going-on-in-this-graph-jan-9-2018.html)

![](https://lh5.googleusercontent.com/sl38jyzhqytJUX60SZypH6ZcmoOyGCDk0-6lk3oOS-YxyXJpRgioK_iF5Y8ZkjESzkuFXPbl6k5qovTuHaJhX-_KdH1sHGWl9nJFy3ixt5sryeEhXmiRfEhrjZjvZi50RtKKTm7UceUL)

```{r tw-109-salary}
library(readxl)
raw <- read_excel("data/tw_salary109.xlsx", sheet=1, trim_ws = T)
raw
```

```{r plot-salary-dist}
raw %>%
    slice(-(1:12)) %>%
    mutate(Category = reorder(Category, desc(Median))) %>%
    ggplot() + aes(y = Category, 
                   xlower=Q1, xmiddle=Median, xupper=Q3, xmin=0, xmax=150) +
    geom_boxplot(stat = "identity", color="white", fill="skyblue") +
    geom_point(aes(x = Mean)) + 
    th + 
    theme(panel.grid.minor = element_blank(), 
          panel.grid.major = element_blank())

```

### TW-Income (boxplot) {#twincome}

本案例用BoxPlot來呈現某個行政區（鄉鎮市區）各村里的所得中位數、平均數、四分位數的分佈。如果在箱型圖中，平均數高於第三分位數，這代表數據集呈現右偏分佈。也就是說，數據中的大部分觀測值都分佈在第一、二分位數之間，但存在一些較大的極端值，使平均值被往右偏移。從這樣的分佈中可以察覺某些里因為有少數極端高收入住戶，而使得平均高於四分位數。

```{r plot-box-plot, message=FALSE, warning=FALSE}
library(gghighlight)
toplot <- read_csv("data/tw_income_107.csv", ) %>%
    filter(!`村里` %in% c("合計", "其他", "福住里")) %>%
    filter(鄉鎮市區 %in% c("信義區")) %>%
    mutate(村里 = reorder(村里, desc(中位數)))

toplot %>% 
    mutate(group = if_else((平均數>第三分位數), "highlight", "none")) %>%
    ggplot() + aes(y = 村里, 
                   xlower=第一分位數, xmiddle=中位數, xupper=第三分位數, 
                   xmin= min(第一分位數), xmax=max(第三分位數), fill=group) +
    geom_boxplot(stat = "identity", color="white") +
    scale_fill_manual(values = c("highlight"="orangered", "none"="skyblue")) + guides(fill=FALSE) + 
    geom_point(aes(x = 平均數)) + 
    xlab("年所得（單位：千元）") + 
    th + 
    theme(panel.grid.minor = element_blank(), 
          panel.grid.major = element_blank())
```

## Likert plot

```{r}
raw <- read_rds("data/tfc_survey.rds")
dt <- raw %>%
  mutate(QA3_lv = ordered(QA3, 
                          levels=c("20-24", "25-29", "30-34", "35-39", 
                                   "40-44", "45-49", "50-54", "55-59", 
                                   "60-64", "65-69", "70及以上"),
                          labels = c("青年", "青年", "壯年", "壯年", 
                                     "壯年", "中年", "中年", "中年", 
                                     "中年", "老年", "老年"))) %>%
  mutate(Q7 = ordered(Q7, levels=c("一點也不會", "不會", "會", "絕對會"))) %>%
  mutate(Q8 = ordered(Q8, levels=c("一點也不會", "不會", "會", "絕對會")))
```

### Stacked or dodged bar

要比較不同年齡層在某個題項的填答結果時，最常見的是用Stacked或Dodged長條圖。Stacked是便於看到各組的總數大小但難以比較各組之間回應的比例，而Dodged是便於比較各組之間每個項目的比例，而不容易觀察總數大小。但這兩種呈現方法，以上面這個例子來說，從視覺化上均難以閱讀出來，哪個年齡層的填答比較靠近「會或絕對會」，又哪個年齡層比較靠近「絕對不會或不會」。

這時候我們可以用一種繪製方法來表達這類Likert問卷的結果，這種圖表稱為Likert Plot（Graph）。

```{r}
p1 <- dt %>%
    count(QA3_lv, Q7) %>%
    ggplot() + 
    aes(QA3_lv, n, fill=Q7) + 
    geom_col(position = position_stack(reverse = TRUE)) + 
    coord_flip() + th

p2 <- dt %>%
    count(QA3_lv, Q7) %>%
    ggplot() + 
    aes(QA3_lv, n, fill=Q7) + 
    geom_col(position="dodge") + 
    th

cowplot::plot_grid(
  p1, NULL, p2,
  labels = c("(a) Stacked", "", "(b) Dodged"),
  ncol = 1, rel_heights = c(1, 0.1, 1)
)

```

### Likert Graph

Likert Graph繪製重點有幾個：

1.  **要轉用比例來繪製**。例如下圖就是用絕對的數值來繪製，因為年齡層人數的不同，例如壯年人數比較多，而老年人數少非常多，反而難以跨組比較。
2.  `ggstats`的套件有`gglikert()`可以用（請見[Plot Likert-type items with \`gglikert()\` • ggstats (larmarange.github.io)](https://larmarange.github.io/ggstats/articles/gglikert.html)）的說明，但也可以用`geom_segment()`來自己刻。
3.  用`geom_segment()`時在`aes()`多了幾個參數，為該資料在X軸的起始點與終點（`x`, `xend`）與Y軸的起始點與終點（`y`, `yend`）。要自己運算。

```{r}
color <- c("#9393C6", "#A8A8A8","#FFA166", "#FF6200")
dt %>%
    count(QA3_lv, Q7) %>%
    mutate(y_acc = cumsum(n)) %>%
    group_by(QA3_lv) %>%
    mutate(y_end = y_acc - min(y_acc) - n[[2]]) %>%
    mutate(y_start = y_end - n) %>%
    ungroup() %>% 
    ggplot() + 
    aes(x = QA3_lv, xend = QA3_lv, y = y_start, yend = y_end, , color=Q7) + 
    geom_segment(linewidth = 18) +
    coord_flip() + theme_bw() + 
    scale_color_manual("", 
                       labels = c("一點也不會", "不會", "會", "絕對會"), 
                       values = color, guide = "legend") +
    th
```

正確用比例繪製的結果如下。Likert Graph和本節所提到的[Pyramid Graph](#pyramid)在數位敘事上的效果很類似，都是對應到一般的Stacked或Dodged長條圖不易做組間比較。Pyramid Graph適於做兩組間的數值左右對照，Likert Graph則有助於快速看出不同題項或不同組別間的填答差異。

```{r}

library(scales)

dt %>% 
    count(QA3_lv, Q7) %>% 
    group_by(QA3_lv) %>%
    mutate(perc = n/sum(n)) %>%
    mutate(y_acc = cumsum(perc)) %>%
    mutate(y_end = y_acc - y_acc[[2]]) %>%
    # mutate(y_end = y_acc - perc[[1]] - perc[[2]]) %>%
    # mutate(y_end = y_acc - min(y_acc) - perc[[2]]) %>%
    mutate(y_start = y_end - perc) %>%
    ungroup() %>% 
    ggplot() + 
    aes(x = QA3_lv, xend = QA3_lv, y = y_start, yend = y_end, , color=Q7) + 
    geom_segment(linewidth = 18) +
    scale_y_continuous(labels = percent_format()) + 
    coord_flip() +
    scale_color_manual("", 
                       labels = c("一點也不會", "不會", "會", "絕對會"), 
                       values = color, guide = "legend") +
    ylab("Perc(%)") + xlab("Age group") + 
    th

```