16 以 moderndive 進行單變數迴歸
本章為 Ismay and Kim (2019) 第 5 章內容。
需要的套件
library(tidyverse)
library(moderndive)
library(skimr)
library(gapminder)其中，moderndive 使用 tidyverse style 編寫，可以從事基本的線性迴歸；skimr 則可以幫助我們快速地計算敘述統計。
16.1 解釋變數為連續變數
以下，我們想了解的是老師的外貌與其教學評鑑的分數有關係嗎？是什麼關係？因此，解釋變數（explanatory variable） x 為老師的外貌分數，而結果變數（outcome variable）為老師的教學評鑑分數。
16.1.1 探索性資料分析
在 moderdive 中，有一 data frame 為 eval，其包含了 UT Austin 的 463 門課程。我們使用 select() 選取我們接下來將要使用的變數，包括 ID、score、bty_avg 與 age，並且存在名為 evals_ch5 裡頭。
evals_ch5 <- evals %>%
select(ID, score, bty_avg, age)回憶第 7 章所提及的，當我們拿到資料要分析或建模之前，第一件事情就是進行探索性資料分析，即 EDA。EDA 讓我們了解個別變數的分配、變數之間是否有潛在的關係與是否有 outliers 或缺漏值。EDA 常見的三個步驟：5
觀察原始資料的值。
計算概括的統計量，如平均數、中位數或四分位距等。
資料視覺化。
16.1.1.1 EDA：觀察原始資料的值
我們可以使用 tibble 中的 glimpse() 查看資料，如：
glimpse(evals_ch5)## Rows: 463
## Columns: 4
## $ ID <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,…
## $ score <dbl> 4.7, 4.1, 3.9, 4.8, 4.6, 4.3, 2.8, 4.1, 3.4, 4.5, 3.8, 4.5, 4.…
## $ bty_avg <dbl> 5.000, 5.000, 5.000, 5.000, 3.000, 3.000, 3.000, 3.333, 3.333,…
## $ age <int> 36, 36, 36, 36, 59, 59, 59, 51, 51, 40, 40, 40, 40, 40, 40, 40…也能使用 dplyr 中的 sample_n() 函數，隨機挑選給定數量的觀察值，如：
evals_ch5 %>%
sample_n(size = 5)## # A tibble: 5 × 4
## ID score bty_avg age
## <int> <dbl> <dbl> <int>
## 1 275 4.8 6.5 38
## 2 302 4.5 3.33 43
## 3 462 4.4 5.33 42
## 4 388 3.5 3 57
## 5 249 3.6 3.17 50由上，我們可知有 463 個觀察值。其中，各個變數分別代表：
ID：在這份資料中用以編號觀察值的變數。score：連續變數；老師的平均教學評鑑分數；最高 5 分，最低 1 分；為我們所關切的結果變數 y。bty_avg：連續變數；老師的平均外貌分數，最高 10 分，最低 1 分；為我們所關切的解釋變數 x。age：連續變數；老師的年齡；為另一個解釋變數 x。
16.1.1.2 EDA：簡單的敘述統計
既然 score 與 bty_avg 是連續的數值變數，那我們可以計算他們的平均值與中位數，如：
evals_ch5 %>%
summarise(mean_score = mean(score),
mean_bty_avg = mean(bty_avg),
median_score = median(score),
median_bty_avg = median(bty_avg))## # A tibble: 1 × 4
## mean_score mean_bty_avg median_score median_bty_avg
## <dbl> <dbl> <dbl> <dbl>
## 1 4.17 4.42 4.3 4.33不過，如果要使用 summarize 來計算標準差、極值與分量的話，並非不可，只是比較麻煩。我們可以使用套件 skimr 中的 skim() 函數。把 data frame 丟入這個函數，將會回傳常用的敘述統計量，如：
evals_ch5 %>% select(score, bty_avg) %>% skim()── Data Summary ────────────────────────
Values
Name Piped data
Number of rows 463
Number of columns 2
_______________________
Column type frequency:
numeric 2
________________________
Group variables None
── Variable type: numeric ───────────────────────────────────────────────────────────────────
skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
1 score 0 1 4.17 0.544 2.3 3.8 4.3 4.6 5 ▁▁▅▇▇
2 bty_avg 0 1 4.42 1.53 1.67 3.17 4.33 5.5 8.17 ▃▇▇▃▂其中，回傳的有：
missing：缺漏值的數量。complete：非缺漏值的數量。n：所有值變數值的數量。mean：平均值。sd：標準差（standard deviation）。p0：0th 分量，即最小值。p25：25th 分量。p50：50th 分量。p75：75th 分量。p100：100th 分量，即最大值。
不過，skim() 回傳的這些都是所謂的 univariate summary statistics。但我們還會想知道 binary summary statistics。例如，如果兩個變數都是數值變數，那我們可以計算相關係數（correlation coefficients） 。
想要計算相關係數的話，我們有兩種方式。
- 使用
moderndive中的get_correlation()函數，使用 formula~，把解釋變數放在~右邊，結果變數放在~左邊，如：
evals_ch5 %>%
get_correlation(formula = score ~ bty_avg)## # A tibble: 1 × 1
## cor
## <dbl>
## 1 0.187- 使用
summarize()中的cor函數，如：
evals_ch5 %>%
summarise(correlation = cor(score, bty_avg))## # A tibble: 1 × 1
## correlation
## <dbl>
## 1 0.18716.1.1.3 EDA：視覺化
加上 position = "jitter"（或者不用 geom_point()，改用 hgeom_jitter()）如同第 5 章所說的，是為了避免 overplotting 的情況。
evals_ch5 %>%
ggplot(aes(x = bty_avg, y = score)) +
geom_point(position = "jitter") +
labs(x = "Beauty Score",
y = "Teaching Score",
title = "Scatterplot of realtionship of traching and beauty scores")
我們也可以使用 geom_smooth(method = "lm")，就能以最小平方法來適配模型，並畫出迴歸線（regression line），如：
evals_ch5 %>%
ggplot(aes(x = bty_avg, y = score)) +
geom_jitter() +
geom_smooth(method = "lm", se = FALSE, formula = 'y ~ x') +
labs(x = "Beauty Score",
y = "Teaching Score",
title = "Relationship between teaching and beauty scores")
事實上，藍色線段的斜率恰等於上節所算出的相關係數 0.0187。不過，這其實是一個特例。
參考文獻
可以發現 Ismay and Kim (2019) 的說法與 Wickham and Grolemund (2016) 不太相同。↩︎