R 資料管理與繪圖

許清芳，鄭中平

2024-04-27

R 簡介（一）

R 是什麼

為何要用 R

R 簡介（二）

點視窗vs寫指令

使用 R 的基本場景

RStudio 與 R Markdown

整齊資料（tidy）搭配管線運算子（pipe operator）

整齊資料

管線運算子

#整體設定，含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup.R")
dta <- read.csv("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/Data/HEXACO.csv", 
                na.strings='NA', stringsAsFactors = TRUE)

dta2 <- select(dta, where(is.factor)) 
head(dta2)
dta |> 
  dplyr::select(where(is.factor))|> 
  head()

資料整理

基本概念

dplyr 及 tidyr

dta |> 
 dplyr::filter(國高中 == '國中') |> 
 dplyr::mutate(攻擊 = 攻擊行為 > (mean(攻擊行為) + 1.96*sd(攻擊行為))) |>
 dplyr::group_by(攻擊, 性別)

我需要記指令嗎

資料整理範例

#整體設定，含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup.R")

資料

資料來自於許功餘、張玉鈴（2015），討論性格向度與青少年問題間的關聯。性格模式包括六個向度（HEXACO）：誠實／謙遜，情緒性，外向性，和悅性，嚴謹性與開放性。青少年問題則測量焦慮／憂鬱、社會退縮兩項內化問題行為，與違反規定以及攻擊行為兩類外化問題行為。

#讀檔案
dta <- read.csv("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/Data/HEXACO.csv", na.strings='NA', stringsAsFactors = TRUE)
#檢視資料結構
#程式報表1.1
str(dta)
'data.frame':   1630 obs. of  14 variables:
 $ 國高中      : Factor w/ 2 levels "高中","國中": 1 1 1 1 1 1 1 1 1 1 ...
 $ 性別        : Factor w/ 2 levels "女","男": 2 1 2 2 NA 1 1 1 2 2 ...
 $ 父親教育程度: Factor w/ 5 levels "大學或專科","小學或不識字",..: 5 4 5 4 NA 2 4 4 2 1 ...
 $ 母親教育程度: Factor w/ 5 levels "大學或專科","小學或不識字",..: 5 5 4 4 NA 4 4 4 2 4 ...
 $ 誠實.謙遜   : int  48 60 53 48 56 61 53 66 55 58 ...
 $ 情緒性      : int  54 43 50 52 60 65 45 68 52 44 ...
 $ 外向性      : int  44 39 48 46 52 39 51 60 38 55 ...
 $ 和悅性      : int  50 55 47 50 49 47 51 48 55 54 ...
 $ 嚴謹性      : int  43 59 52 46 53 49 48 51 43 43 ...
 $ 開放性      : int  44 57 44 49 44 43 45 51 38 50 ...
 $ 攻擊行為    : int  13 2 1 12 4 NA 21 1 4 2 ...
 $ 焦慮.憂鬱   : int  11 2 0 9 2 6 14 6 5 1 ...
 $ 違反規定    : int  7 1 0 9 1 2 12 0 1 3 ...
 $ 社會退縮    : int  9 4 0 5 3 7 12 2 4 2 ...
# 看前六筆
head(dta)
國高中 性別 父親教育程度 母親教育程度 誠實.謙遜 情緒性 外向性 和悅性 嚴謹性 開放性 攻擊行為 焦慮.憂鬱 違反規定 社會退縮
高中 男 國中 國中 48 54 44 50 43 44 13 11 7 9
高中 女 高中 國中 60 43 39 55 59 57 2 2 1 4
高中 男 國中 高中 53 50 48 47 52 44 1 0 0 0
高中 男 高中 高中 48 52 46 50 46 49 12 9 9 5
高中 NA NA NA 56 60 52 49 53 44 4 2 1 3
高中 女 小學或不識字 高中 61 65 39 47 49 43 NA 6 2 7

選取變項

選取變項方式一：利用資料型態

#本章主要利用 tidyverse 中的 dplyr 和 tidyr 兩個套件做資料處理  
#程式報表1.2
dta |> 
  dplyr::select(where(is.factor)) |> 
  head()
國高中 性別 父親教育程度 母親教育程度
高中 男 國中 國中
高中 女 高中 國中
高中 男 國中 高中
高中 男 高中 高中
高中 NA NA NA
高中 女 小學或不識字 高中

選取變項方式二：利用變項名稱（部份）文字

#高誠實／謙遜個體出現較少違反規定的行為嗎？
dta |> 
  dplyr::select(contains(c('誠', '規'))) |> 
  cor(use='pair') |> 
  round(3)
誠實.謙遜 違反規定
誠實.謙遜 1.00 -0.38
違反規定 -0.38 1.00

選取變項方式三：利用變項位置

#焦慮／憂鬱、社會退縮屬於內化問題行為
#違反規定、攻擊行為屬於外化問題行為
#我們將類別相同變項重排在一起
#程式報表1.3
dta |> 
  dplyr::select(c(11, 13, 12, 14)) |>
  cor(use='pair') |> 
  round(3)
攻擊行為 違反規定 焦慮.憂鬱 社會退縮
攻擊行為 1.000 0.696 0.610 0.417
違反規定 0.696 1.000 0.423 0.351
焦慮.憂鬱 0.610 0.423 1.000 0.623
社會退縮 0.417 0.351 0.623 1.000

選取觀察值

選取觀察值（slice）方式一：利用列的位置

#程式報表1.4
dta |> 
  dplyr::slice(11:12) |> 
  dplyr::select(1:4)
國高中 性別 父親教育程度 母親教育程度
高中 女 研究所以上 大學或專科
高中 男 高中 國中

選取觀察值（slice）方式二：利用列的位置（隨機若干位）

#程式報表1.5
dta |> 
  dplyr::slice_sample(n=3) |> 
  dplyr::select(contains(c('親')))
父親教育程度 母親教育程度
高中 高中
高中 高中
國中 國中

篩選觀察值（filter）方式：指定變項與條件

dta |> 
  dplyr::filter(性別 == '男') |>
  dplyr::select(contains(c('誠', '規'))) |> 
  cor(use='pair') |> 
  round(3)
誠實.謙遜 違反規定
誠實.謙遜 1.000 -0.318
違反規定 -0.318 1.000

篩選觀察值（filter）方式：可以同時利用多個條件篩選

dta |> 
  dplyr::filter(母親教育程度=='研究所以上',
                父親教育程度=='研究所以上') |>
  dplyr::select(性別) |> 
  table()
女 男
6 8

轉換變項,轉換變項值

#類別變項水準預設排列順序對應編碼大小，常常需要重排
#程式報表1.7前
with(dta, levels(母親教育程度))
[1] "大學或專科"   "小學或不識字" "研究所以上"   "高中"         "國中"

變項轉換一：製造新變項

#程式報表1.7後
dta |> 
  dplyr::mutate(母教育 = forcats::fct_relevel(母親教育程度,
                                              c("小學或不識字", 
                                                "國中", 
                                                "高中", 
                                                "大學或專科", 
                                                "研究所以上")),
                父教育 = forcats::fct_relevel(父親教育程度,
                                              c("小學或不識字", 
                                                "國中", 
                                                "高中", 
                                                "大學或專科", 
                                                "研究所以上")),
                .keep='none') |> table()
小學或不識字 國中 高中 大學或專科 研究所以上
小學或不識字 18 33 18 2 0
國中 21 181 106 20 0
高中 17 147 475 109 5
大學或專科 2 15 60 245 31
研究所以上 0 0 2 7 15

變項轉換二：修改舊變項

dta <- dta |> 
  dplyr::mutate(母親教育程度 = forcats::fct_relevel(母親教育程度,
                                              c("小學或不識字", 
                                                "國中", 
                                                "高中", 
                                                "大學或專科", 
                                                "研究所以上")),
                父親教育程度 = forcats::fct_relevel(父親教育程度,
                                              c("小學或不識字", 
                                                "國中", 
                                                "高中", 
                                                "大學或專科", 
                                                "研究所以上")))

分組

分組，對各組做相同動作

tmp <- dta |> 
  dplyr::select(性別, 攻擊行為, 違反規定) |>
  dplyr::group_by(性別) |>
  dplyr::mutate(攻擊 = case_when(攻擊行為 < quantile(攻擊行為, prob=.33,
                                  na.rm = TRUE) ~ '1_低',
                         攻擊行為 > quantile(攻擊行為, prob=.67,
                                  na.rm = TRUE) ~ '3_高',
                         .default = '2_中'), 
                違規 = case_when(違反規定 < quantile(違反規定, prob=.33,
                                  na.rm = TRUE) ~ '1_低',
                          違反規定 > quantile(違反規定, prob=.67,
                                  na.rm = TRUE) ~ '3_高',
                         .default = '2_中'), .keep='unused') |> as.data.frame()
#剛剛資料作列聯表
#程式報表1.10
tmp |> ftable()
          違規 1_低 2_中 3_高
性別 攻擊                    
女   1_低        81  102   14
     2_中        30  237   72
     3_高         7   86  163
男   1_低       165   52   12
     2_中        72  199   81
     3_高        17   75  151

分組後合計

#合計：算相關
#reframe 會對每一組做同樣動作，並且存成資料框，以便後續利用方便
#剛剛為展示緣故，將連續變項轉成類別變項（通常不必要）做列聯表
#這邊直接針對連續變項求相關
#程式報表1.11
dta |> 
  dplyr::select(性別, 攻擊行為, 違反規定) |>
  dplyr::group_by(性別) |>
  dplyr::reframe(相關係數 = cor(攻擊行為, 違反規定, use='pair'),
                  合計 = n())
性別 相關係數 合計
女 0.6711 792
男 0.7097 824
NA 0.6960 14 
#合計：算基本統計量
#程式報表1.12
dta |> 
  dplyr::group_by(性別) |>
  dplyr::reframe(違規平均 = mean(違反規定, na.rm = TRUE), 
                 違規標準差 = sd(違反規定, na.rm = TRUE),
                 違規中位數 = median(違反規定, na.rm=T), 
                 違規四分位距 = IQR(違反規定, na.rm=T),
                 合計 = n())
性別 違規平均 違規標準差 違規中位數 違規四分位距 合計
女 2.863 2.570 2.0 3.00 792
男 4.335 3.211 4.0 4.00 824
NA 5.214 4.475 4.5 4.75 14

資料形式變換：寬形到長形

#為示範方便，先選出2筆、4個變項，總共8個數值
#程式報表1.13前
dta |> 
  dplyr::select(11:14) |> 
  dplyr::slice(1:2) |> 
  as_tibble(rownames="識別碼")
識別碼 攻擊行為 焦慮.憂鬱 違反規定 社會退縮
1 13 11 7 9
2 2 2 1 4 
#將剛剛資料轉成長形
#請上下對照並注意 id，看看剛剛8個數值現在如何排列
#程式報表1.13後
dta |> dplyr::select(攻擊行為:社會退縮) |> 
  dplyr::slice(1:2) |>
  as_tibble(rownames="識別碼") |>
  tidyr::pivot_longer(cols = -識別碼, 
                      names_to = '行為問題', values_to = '分數')
識別碼 行為問題 分數
1 攻擊行為 13
1 焦慮.憂鬱 11
1 違反規定 7
1 社會退縮 9
2 攻擊行為 2
2 焦慮.憂鬱 2
2 違反規定 1
2 社會退縮 4 
#利用 arrange指令，依性別排序
#程式報表1.14
dta |> dplyr::select(2, 11:14) |> 
  as_tibble(rownames="識別碼") |>
  tidyr::pivot_longer(cols = -c(識別碼, 性別),
                      names_to = '行為問題', values_to = '分數') |>
  dplyr::group_by(行為問題, 性別) |>
  dplyr::reframe(平均分數 = mean(分數, na.rm=T)) |> 
  dplyr::arrange(性別) |>
  dplyr::filter(性別 != 'NA')
行為問題 性別 平均分數
攻擊行為 女 7.991
焦慮.憂鬱 女 7.014
社會退縮 女 4.589
違反規定 女 2.863
攻擊行為 男 9.166
焦慮.憂鬱 男 6.128
社會退縮 男 4.545
違反規定 男 4.335 
#綜合應用，最後面將資料由長形變成寬形
#程式報表1.15
dta |> dplyr::select(2, c(11:14)) |> 
  as_tibble(rownames="id") |>
  tidyr::pivot_longer(cols = -c(id, 性別),
                      names_to = '行為問題', values_to = '分數') |>
  dplyr::group_by(行為問題, 性別) |>
  dplyr::reframe(平均分數 = mean(分數, na.rm=T)) |> 
  dplyr::filter(性別 != 'NA') |>
  tidyr::pivot_wider(names_from = '行為問題',
                     values_from = '平均分數')
性別 攻擊行為 焦慮.憂鬱 社會退縮 違反規定
女 7.991 7.014 4.589 2.863
男 9.166 6.128 4.545 4.335 
#綜合應用，計算以違反規定PR90分類學生時，女生所佔各類比率
dta |>
 dplyr::mutate(違規九 = 違反規定 > quantile(違反規定, prob = 0.9, na.rm=T)) |>
 dplyr::group_by(違規九, 性別) |>
 dplyr::reframe(合計 = n()) |>  
 pivot_wider(names_from=性別, 
             values_from=合計) |>
 dplyr::mutate(違規九分位 = 違規九, 
               女性百分率 = 100*(女 / (女+男)), .keep='none')  |>
 dplyr::filter(違規九分位 != 'NA')
違規九分位 女性百分率
FALSE 51.36
TRUE 27.39

資料形式變換：長形到寬形

#綜合應用，不同母親教育程度、性別下，內化問題的平均
#程式報表1.16
dta |> 
  dplyr::select(性別, 母親教育程度, 社會退縮, 焦慮.憂鬱) |>
  tidyr::pivot_wider(names_from = '性別',
                     values_from = c('社會退縮', '焦慮.憂鬱'),
                     values_fn = ~ mean(.x, na.rm = TRUE)) |>
  dplyr::filter(母親教育程度 != 'NA') |>
  dplyr::select(-contains('NA')) |>
  dplyr::arrange(社會退縮_女)
母親教育程度 社會退縮_男 社會退縮_女 焦慮.憂鬱_男 焦慮.憂鬱_女
研究所以上 4.846 3.800 7.000 7.100
大學或專科 4.611 4.390 6.378 6.622
高中 4.486 4.509 5.942 6.931
國中 4.680 4.647 6.189 6.822
小學或不識字 4.938 5.575 7.613 9.053

利用 ggplot2 繪製統計圖

簡單認識 ggplot2

好複雜！

練習認出角色與幾何成分

knitr::include_graphics("GenderMarriage103.png")

繪圖範例

資料與管理

#讀檔案
dta <- read.csv(file = "https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/Data/HEXACO_HS.csv",header = TRUE, stringsAsFactors = TRUE)
#檢視資料結構
#程式報表2.1
str(dta)
'data.frame':   897 obs. of  13 variables:
 $ 性別        : Factor w/ 2 levels "女","男": 2 1 2 2 1 1 2 2 1 1 ...
 $ 父親教育程度: Factor w/ 5 levels "大學或專科","小學或不識字",..: 5 4 5 4 4 4 2 1 3 1 ...
 $ 母親教育程度: Factor w/ 5 levels "大學或專科","小學或不識字",..: 5 5 4 4 4 4 2 4 1 4 ...
 $ 誠實.謙遜   : int  48 60 53 48 53 66 55 58 51 58 ...
 $ 情緒性      : int  54 43 50 52 45 68 52 44 57 62 ...
 $ 外向性      : int  44 39 48 46 51 60 38 55 54 67 ...
 $ 和悅性      : int  50 55 47 50 51 48 55 54 51 61 ...
 $ 嚴謹性      : int  43 59 52 46 48 51 43 43 50 53 ...
 $ 開放性      : int  44 57 44 49 45 51 38 50 50 52 ...
 $ 攻擊行為    : int  13 2 1 12 21 1 4 2 8 4 ...
 $ 焦慮.憂鬱   : int  11 2 0 9 14 6 5 1 6 4 ...
 $ 違反規定    : int  7 1 0 9 12 0 1 3 4 3 ...
 $ 社會退縮    : int  9 4 0 5 12 2 4 2 5 0 ...
#設定類別變項的順序
dta <- dta |> 
  dplyr::mutate(母親教育程度 = forcats::fct_relevel(母親教育程度, 
                                              c("小學或不識字", 
                                                "國中", 
                                                "高中", 
                                                "大學或專科", 
                                                "研究所以上")),
                父親教育程度 = forcats::fct_relevel(父親教育程度, 
                                              c("小學或不識字", 
                                                "國中", 
                                                "高中", 
                                                "大學或專科", 
                                                "研究所以上")))

一步一步填滿圖層

#ggplot 是一個一個圖層（layer）疊上去
#底下以散佈圖為例，呈現每一步驟的結果
#第一步驟設定圖的框架，注意圖的X軸與Y軸
#圖2.1
g0 <- ggplot(data = dta, 
             aes(x=社會退縮, y=攻擊行為)) 
g0

#在前一步驟結果（圖的框架）上加入點
#圖2.2
g1 <- g0 + geom_point(alpha=.2) 
g1

#在前一步驟結果上加入橢圓（資料的95%區間）與局部迴歸線
#圖2.3
g2 <- g1 + stat_smooth(method='lm', 
                       formula = y ~ x,
                       se=FALSE, 
                       linewidth = 0.5) 
g2

#在前一步驟結果上要求以變項（性別）區分顏色
#圖2.4
g3 <- g2 + aes(color=性別) 
g3

#在前一步驟結果上設定 x 軸與 y 軸刻度
#圖2.5
g4 <- g3 + 
  scale_y_continuous(breaks=seq(0, 25, by=5)) +
  scale_x_continuous(breaks=seq(0, 12, by=2))
g4

#在前一步驟結果上要求以變項（父母教育）分面（facet）
#圖2.6
g5 <- g4 + facet_wrap(vars(母親教育程度),nrow=1)
g5

#在前一步驟結果上加上X軸、Y軸的標示，以及整個圖形的標題
#圖2.7
g6 <- g5 + labs(x='社會退縮分數',
                y='攻擊行為分數',
                title='散布圖：攻擊行為與社會退縮')
g6

#在前一步驟結果上改變主題，並要求圖示位置
#圖2.8
g7 <- g6 + theme_minimal() + 
  theme(legend.position='top')
g7

#前面步驟可以一次執行
ggplot(data = dta, aes(x=社會退縮, y=攻擊行為)) +
   geom_point(alpha=.2) +
 stat_smooth(method='lm', formula = y ~ x,
       se=FALSE, linewidth = 0.5)+
 aes(color=性別) +
 scale_y_continuous(breaks=seq(0, 25, by=5)) +
 scale_x_continuous(breaks=seq(0, 12, by=2)) +
 facet_wrap(vars(母親教育程度),nrow=1)+
 labs(x='社會退縮分數', y='攻擊行為分數',
     title='散布圖：攻擊行為與社會退縮') +
 theme_minimal() + 
 theme(legend.position='top')

繪製統計摘要

以 ggplot 直接繪製統計摘要結果（連續資料平均數與標準誤）

#圖2.9
ggplot(data=dta, 
     aes(x=母親教育程度, y=誠實.謙遜, color=性別)) +
  stat_summary(fun.data = "mean_cl_boot", 
               position=position_dodge(.2)) +
  stat_summary(aes(group=性別), fun = mean, 
               geom="line",
               position=position_dodge(.2))+
  scale_color_grey(end=.7)+
  labs(y='誠實.謙遜平均分數',
       x='母親教育程度',
       title='不同性別的誠實.謙遜平均跟母親教育程度的關係',
       caption="來源: 許功餘")+
  theme(legend.position='top')

以 ggplot 直接繪製統計摘要結果（類別資料的百分比）

#圖2.10
ggplot(dta, 
       aes(x=母親教育程度, group=父親教育程度)) + 
  geom_bar(aes(y=after_stat(prop), 
          fill = factor(after_stat(x))),
          width = .2) + 
  scale_y_continuous(labels=scales::percent) +
  scale_fill_grey()+
  coord_flip()+
  labs(y="百分比",
       title="門當戶對:父母親教育程度") +
  facet_wrap(vars(父親教育程度), ncol=1)+
  theme(legend.position="none")

以 ggplot 直接繪製統計摘要結果（連續資料平均數與標準誤）

#圖2.11
ggplot(data=dta, 
     aes(x=母親教育程度, y=誠實.謙遜)) +
  stat_summary(fun.data = "mean_cl_boot") +
  scale_color_grey(end=.7)+
  facet_wrap(vars(父親教育程度), nrow=1)+
  labs(y='誠實.謙遜平均分數',
       x='母親教育程度',
       title='不同父親教育程度的誠實.謙遜平均跟母親教育程度的關係',
       caption="來源: 許功餘")+
  theme(legend.position='top',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
Warning: Removed 1 rows containing missing values (`geom_segment()`).
Removed 1 rows containing missing values (`geom_segment()`).
Removed 1 rows containing missing values (`geom_segment()`).

多變量圖形，利用寬轉長

#多變量圖形 
#圖2.12,dpi=300
dta |> tidyr::pivot_longer(cols=4:9, 
                           names_to = '人格維度',
                           values_to = '分數') |>
  ggplot()+
  aes(x=違反規定, y=分數, color=性別)+
  geom_point(size=rel(.5))+
  stat_smooth(method='lm', 
              formula = y ~ x,
              se=F, linewidth=.5)+
  facet_grid(人格維度 ~ 母親教育程度) +
  scale_color_grey(start=.1, end=.6)+
  labs(x="違反規定分數",
       y="人格維度分數")+
  theme(legend.position='top')

延伸

#同時繪製兩個變項的散布圖，以及各自的邊際分布
#圖2.13
p <- ggplot(dta, 
            aes(x=誠實.謙遜, y=違反規定)) +
  geom_point(shape=21, alpha=.5) +
  stat_ellipse() +
  geom_vline(xintercept=mean(dta$誠實.謙遜), col="gray") +
  geom_hline(yintercept=mean(dta$違反規定), col="gray") +
  stat_smooth(method="lm", 
              formula = y ~ x,
              linewidth=.7,
              linetype="dotted",
              se=FALSE,
              alpha=.5,
              col=1) +
  labs(y="違反規定分數", 
       x="誠實.謙遜分數") 
pacman::p_load(ggExtra, KernSmooth)
ggMarginal(p, 
           type="histogram",
           xparams=list(binwidth=dpih(dta$誠實.謙遜),
                          fill="gray90"),
           yparams=list(binwidth=dpih(dta$違反規定),
                          fill="gray90"))

小結