R | Data manipulation (2): arrange, mutate, summarise

這篇介紹另外幾個處理資料的功能，上篇請看這裡：R | Data manipulation (1): select, filter, slice

詳細功能介紹請參考這兩頁：Working with Data - Part 1 & Intro to R - Part 3

練習的部分可以參考 slides: Working with Data - Part 1 (pdf)

這幾個功能在 dplyr 這個 package 裡面，所以要先跑。同時要用內建檔案 mtcars 來練習，所以也要叫出來。下面還會用檔案資料 flights，因為它在 nycflights13 這個 packages 裡面，所以也需要安裝。

library(dplyr) library(nycflights13) mtcars flights

下面先介紹一個觀看檔案的功能。

Glimpse function: columns run down the page, and data runs across, making it possible to see every column in a data frame

glimpse() 和平常看檔案的方式不一樣。我們平常看檔案的時候，variables 是在第一列(rows)，每個欄位(column)是一個 variables。但用 glimpse() 的話，它會橫向顯示每個 variable 的數據值。

下面先用資料檔案 flights 練習。

dim(flights)

dim: 顯示檔案資料的大小 dimension (row x column)

[1] 336776     19

上面顯示的結果告訴我們，共有 336776 個班機(observations)，且用 19 個 variables 呈現各個班機的資訊。

glimpse(flights)

Observations: 336,776 Variables: 19 $ year           <int> 2013, 2013, 2013, 2013, 2013, 2013, ... $ month          <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ... $ day            <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ... $ dep_time       <int> 517, 533, 542, 544, 554, 554, 555, ...

上面用 glimpse() 功能可以看到每班飛機的資訊(variables)包括有 year, month, day, dep_time 等等。

5. Arrange function: order observations (rows)

沒特別指令的話會由小到大排列。

Default setting: order mpg values from small to large. (依 mpg 的大小排列)

arrange(mtcars, mpg)

mpg cyl  disp  hp drat    wt  qsec vs am gear carb 1  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4 2  10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4 3  13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4 4  14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4 5  14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4

如果要由大到小排列的話，則用：desc()

Use desc() to sort in descending order.
Order mpg from large to small.

arrange(mtcars, desc(mpg))

mpg cyl  disp  hp drat    wt  qsec vs am gear carb 1  33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1 2  32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1 3  30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2 4  30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2 5  27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1

如果有兩個 variables 的話，會先排第一個，然後再在第一個裡面排第二個的大小順序。

mtcars %>%    arrange(cyl, gear)

Arrange with cyl first; then within cyl, arrange gear.

上面的例子中，會先依照 cyl 的大小排列，再依 gear 的大小排列。

mpg cyl  disp  hp drat    wt  qsec vs am gear carb 1  21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1 2  22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1 3  24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2 4  22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2 5  32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1

如果想讓兩個 variables 都由大到小排列，可以分開寫，也可以合在一起。

Reorder the mtcars in descending order of mpg and displacement (disp)

mtcars %>%   arrange(desc(carb), desc(mpg))

或寫成這樣：

mtcars %>%   arrange(desc(carb, mpg))

mpg cyl  disp  hp drat    wt  qsec vs am gear carb 1  15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8 2  19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6 3  21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4 4  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4 5  19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4 6  17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4

從上面可以看到是先依 carb 的大小排列，之後再依 mpg 的大小排列（最後四排當 carb 都是 4 的時候，mpg 是由大到小排列），但兩個都是由大到小。

Exercise

這邊我們用檔案資料 flights 來練習。

Q1. Which flights are the most delayed?

哪個班機起飛時間延遲最多，可以用 desc() 由大到小排列延遲的長度，第一個就是延遲最多的。下面先用 glimpse 呈現結果。

flights %>%   arrange(desc(dep_delay)) %>%   glimpse

Observations: 336,776 Variables: 19 $ year           <int> 2013, 2013, 2013, 2013, 2013, ... $ month          <int> 1, 6, 1, 9, 7, 4, 3, 6, 7, 12, 5, 1, 2, ... $ day            <int> 9, 15, 10, 20, 22, 10, 17, 27, 22, ... $ dep_time       <int> 641, 1432, 1121, 1139, 845, 1100, 2321, ... $ sched_dep_time <int> 900, 1935, 1635, 1845, 1600, 1900, 810, ... $ dep_delay      <dbl> 1301, 1137, 1126, 1014, 1005, 960, ...

從上面顯示的最後一行 dep_delay 可以看到是由大到小排列。

如果不確定哪個 variable 是指延遲時間的話，可以用 ?  或 help() 的功能來查看：?flights 或是 help(flights)。

會在 R Studio 右下角的視窗中出現解釋框(如下圖，點圖可以放大)，裡面寫：

dep_delay, arr_delay: Departure and arrival delays, in minutes. Negative times represent early departures/arrivals.

下面是沒用 glimpse() 功能呈現，應該會比較清楚。

arrange(flights, desc(dep_delay))

# A tibble: 336,776 x 19     year month   day dep_time sched_dep_time dep_delay arr_time    <int> <int> <int>    <int>          <int>     <dbl>    <int>  1  2013     1     9      641            900      1301     1242  2  2013     6    15     1432           1935      1137     1607  3  2013     1    10     1121           1635      1126     1239  4  2013     9    20     1139           1845      1014     1457  5  2013     7    22      845           1600      1005     1044  6  2013     4    10     1100           1900       960     1342  7  2013     3    17     2321            810       911      135  8  2013     6    27      959           1900       899     1236  9  2013     7    22     2257            759       898      121 10  2013    12     5      756           1700       896     1058 # ... with 336,766 more rows, and 11 more variables: #   carrier &arr_delay <dbl>,lt;chr>, flight <int>, tailnum <chr>, #   origin <chr>, air_time <dbl>, distance <dbl>, hour <dbl>, #   minute <dbl>, time_hour <dttm>

Q2: Which flights caught up the most time during the flight?

哪班飛機趕上的時間最多，就是說即便是延遲起飛，但是仍然準時抵達的，也就是延遲時間最短的(比預計的飛行時間短)，這邊用延遲起飛的時間減掉延遲抵達的時間，再由大到小排列。

flights %>%   arrange(desc(dep_delay - arr_delay)) %>%   glimpse

Observations: 336,776 Variables: 19 $ year           <int> 2013, 2013, 2013, 2013, 2013, 2013, ... $ month          <int> 6, 2, 2, 5, 2, 7, 7, 12, 5, 11, 5, 5, ... $ day            <int> 13, 26, 23, 13, 27, 14, 17, 27, 2, 13, ... $ dep_time       <int> 1907, 1000, 1226, 1917, 924, 1917, ... $ sched_dep_time <int> 1512, 900, 900, 1900, 900, 1829, 1930, ...

不用 glimpse 顯示：

arrange(flights, desc(dep_delay - arr_delay))

# A tibble: 336,776 x 19     year month   day dep_time sched_dep_time dep_delay arr_time    <int> <int> <int>    <int>          <int>     <dbl>    <int>  1  2013     6    13     1907           1512       235     2134  2  2013     2    26     1000            900        60     1513  3  2013     2    23     1226            900       206     1746  4  2013     5    13     1917           1900        17     2149  5  2013     2    27      924            900        24     1448  6  2013     7    14     1917           1829        48     2109  7  2013     7    17     2004           1930        34     2224  8  2013    12    27     1719           1648        31     1956  9  2013     5     2     1947           1949        -2     2209 10  2013    11    13     2024           2015         9     2251 # ... with 336,766 more rows, and 11 more variables: #   arr_delay <dbl>, carrier <chr>, flight <int>, tailnum <chr>, #   origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>, #   hour <dbl>, minute <dbl>, time_hour <dttm>

可以指定 catchup 為 flights，再帶入 arrange，或是用 head 只顯示前十個 carrier。

flights %>%   arrange(desc(dep_delay - arr_delay)) head(catchup$carrier, 10)

[1] "EV" "HA" "HA" "DL" "HA" "UA" "UA" "UA" "UA" "DL"

如果只看資料的前十個的話，則是下面這樣。

head(flights$carrier, 10)

[1] "UA" "UA" "AA" "B6" "DL" "UA" "B6" "EV" "B6" "AA"

6. Mutate function: 設定新的 variable，新設定的 variable 會顯示在最後一欄。

Make new variables: disp_l and wt_kg

1L = 61.0237 cu.in. (cubic inch)
1kg = 2.2 lbs

檔案資料裡面的 disp 單位是 cu.in.，我們把它轉換成 L，設定其為 disp_l。同時也把裡面原本單位為磅(lbs)的 wt 換算成 kg，指定其為 wt_kg。新訂的 disp_l 和 wt_kg 會顯示在最後兩的 column。

mtcars %>%  mutate(disp_l = disp/61.0237,         wt_kg = wt/2.2)

mpg cyl  disp  hp drat    wt  qsec vs am gear carb   disp_l     wt_kg 1  21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4 2.621932 1.1909091 2  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4 2.621932 1.3068182 3  22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1 1.769804 1.0545455 4  21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1 4.227866 1.4613636 5  18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2 5.899347 1.5636364

把重量換算成噸，指定其為：wt_tonnes
Include a variable for weight in tonnes (1t = 2,204.6 lbs)

mtcars %>%   mutate(wt_tones = wt / 2.2046)

mpg cyl  disp  hp drat    wt  qsec vs am gear carb  wt_tones 1  21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4 1.1884242 2  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4 1.3040914 3  22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1 1.0523451 4  21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1 1.4583144 5  18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2 1.5603738

transmute: leave only the new variables

如果用 transmute() 的話，就只會顯示新訂的 disp_1 和 wt_kg。

mtcars %>%   transmute(disp_l = disp/61.0237, wt_kg = wt/2.2)

disp_l     wt_kg 1  2.621932 1.1909091 2  2.621932 1.3068182 3  1.769804 1.0545455 4  4.227866 1.4613636 5  5.899347 1.5636364 6  3.687092 1.5727273

Exercise

Q1. Compute speed in mph from time (in minutes) and distance (in miles)

同樣可以用 ?flights 查詢：

distance: Distance between airports, in miles air_time: Amount of time spent in the air, in minutes

計算飛機的速度：mph = miles per hour

也就是飛行的距離 distance (in miles) 除以飛行的時間 air_time (in minutes)。

因為 air_time 是以分鐘為單位，所以我們必須先把它轉換成小時，也就是：air_time / 60

下面把速度指定為 speed。

mutate(flights, speed = distance / (air_time / 60))

# A tibble: 336,776 x 20     year month   day dep_time sched_dep_time dep_delay arr_time    <int> <int> <int>    <int>          <int>     <dbl>    <int>  1  2013     1     1      517            515         2      830  2  2013     1     1      533            529         4      850  3  2013     1     1      542            540         2      923  4  2013     1     1      544            545        -1     1004  5  2013     1     1      554            600        -6      812  6  2013     1     1      554            558        -4      740  7  2013     1     1      555            600        -5      913  8  2013     1     1      557            600        -3      709  9  2013     1     1      557            600        -3      838 10  2013     1     1      558            600        -2      753 # ... with 336,766 more rows, and 12 more variables: #   arr_delay <dbl>, carrier <chr> flight <int>, tailnum <chr>, #   origin <chr>, dest <chr>, air_time <dbl>,  distance <dbl>, #   hour <dbl>, minute <dbl>, time_hour <dttm>, speed <dbl>

顯示整個檔案太長，看不到後面計算出來的 speed，所以用 select() 的功能挑出我們想看的幾個 variables。如果忘記 select() 是什麼、怎麼用的話，可以看前一篇。

mutate(flights, speed = distance / (air_time / 60)) %>%   select(carrier, distance, air_time, speed)

# A tibble: 336,776 x 4    carrier distance air_time    speed      <chr>    <dbl>    <dbl>    <dbl>  1      UA     1400      227 370.0441  2      UA     1416      227 374.2731  3      AA     1089      160 408.3750  4      B6     1576      183 516.7213  5      DL      762      116 394.1379  6      UA      719      150 287.6000  7      B6     1065      158 404.4304  8      EV      229       53 259.2453  9      B6      944      140 404.5714 10      AA      733      138 318.6957 # ... with 336,766 more rows

Q2. Which flight flew the fastest?

哪個班機的速度最快？可以用 arrange() 功能由大到小排列出來，跟上面一樣只挑出我們想看的幾個 variables 來看。

flights %>%   mutate(speed = distance / (air_time / 60)) %>%   arrange(desc(speed)) %>%   select(carrier, distance, air_time, speed)

# A tibble: 336,776 x 4    carrier distance air_time    speed      <chr>    <dbl>    <dbl>    <dbl>  1      DL      762       65 703.3846  2      EV     1008       93 650.3226  3      EV      594       55 648.0000  4      EV      748       70 641.1429  5      DL     1035      105 591.4286  6      DL     1598      170 564.0000  7      B6     1598      172 557.4419  8      AA     1623      175 556.4571  9      DL     1598      173 554.2197 10      B6     1598      173 554.2197 # ... with 336,766 more rows

Q3. 把上面班機延遲練習的 Q2，用 select() 的功能挑出我們想看的幾項出來。

flights %>%   mutate(delay_time = dep_delay - arr_delay) %>%   arrange(desc(delay_time)) %>%   select(carrier, dep_time, arr_time, air_time,          dep_delay, arr_delay, delay_time)

# A tibble: 336,776 x 7    carrier dep_time arr_time air_time dep_delay arr_delay delay_time      <chr>    <int>    <int>    <dbl>     <dbl>     <dbl>      <dbl>  1      EV     1907     2134      126       235       126        109  2      HA     1000     1513      584        60       -27         87  3      HA     1226     1746      599       206       126         80  4      DL     1917     2149      313        17       -62         79  5      HA      924     1448      589        24       -52         76  6      UA     1917     2109      274        48       -26         74  7      UA     2004     2224      295        34       -40         74  8      UA     1719     1956      324        31       -42         73  9      UA     1947     2209      300        -2       -75         73 10      DL     2024     2251      311         9       -63         72 # ... with 336,766 more rows

可以和上面排列前十的 carrier 的結果比對，看排列是不是一樣的。


7. Summarise / Summarize

功能和 mutate() 有點像，不一樣的是它會產生一個新的 data frame，並且是那欄（也就是那個 variable）所有觀察資料的總結。

Summarise works in an analogous way to mutate, except instead of adding columns to an existing data frame, it creates a new data frame. This is particularly useful in conjunction with ddply as it makes it easy to perform group-wise summaries.

例如算出 mpg 的平均值，可以用 mean() 的功能這樣算：

mean(mtcars$mpg)

[1] 20.09062

如果用 mutate() 的話會是這樣：

mutate(mtcars, mean(mpg))

mpg cyl  disp  hp drat    wt  qsec vs am gear carb mean(mpg) 1  21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4  20.09062 2  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4  20.09062 3  22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1  20.09062 4  21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1  20.09062 5  18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2  20.09062

同樣是算出 mean，但它會放在最後一欄，所以所有數值都一樣，因為它是全部觀察算出來的平均值。

如果用 summarise() 的話則是這樣，以一個 data frame 的樣式呈現出來：

summarise(mtcars, mean(mpg))

mean(mpg) 1  20.09062

也可以幫平均值指定一個名稱 mpg_mean，直接寫在裡面即可。

summarise(mtcars, mpg_mean = mean(mpg)

mpg_mean 1 20.09062

也可以設定三個，例如把 mpg 的平均值設為 mean_mpg，把重量的中間值設為 median_wt，然後把這兩者的比例設為 ratio。

找出中間值用的功能是：median()

mtcars %>%   summarise(mean_mpg = mean(mpg),             median_wt = median(wt),             ratio = mean_mpg/median_wt)

mean_mpg median_wt    ratio 1 20.09062     3.325 6.042293

也可以全部寫在一行，下面同時把 mpg 的平均值直接設為 mpg，wt 的中間值設為 wt，把兩者的比例設為 ratio。

summarise(mtcars, mpg = mean(mpg),                   wt = median(wt),                   ratio = mpg/wt)

mpg    wt    ratio 1 20.09062 3.325 6.042293

除了算全部的平均值外，也可以分組算，例如算出 cyl 裡各種觀察值的平均值。下面先用 levels() 看 cyl 裡面有哪幾種，因為 cyl 的觀察值是數字，需要先把它換成 factor。

levels(factor(mtcars$cyl))

[1] "4" "6" "8"

由上面結果可以看出 cyl 的觀察值有三種：4, 6, 8。如果要看各個的 mpg 和重量平均值，可以用 group_by() 的功能。

mtcars %>%   group_by(cyl) %>%　   summarise(mean_mpg = mean(mpg), median_wt = median(wt))

# A tibble: 3 x 3     cyl mean_mpg median_wt   <dbl>    <dbl>     <dbl> 1     4 26.66364     2.200 2     6 19.74286     3.215 3     8 15.10000     3.755

也可以在一組裡面再分組算平均值，例如先用 am 分組後，在再裡面用 cyl 分組，先分的那個放前面，所以語法會是：group_by(am, cyl)

mtcars %>%   group_by(am, cyl) %>%   summarise(mpg = mean(mpg), wt = median(wt))

# A tibble: 6 x 4 # Groups:   am [?]      am   cyl      mpg     wt   <dbl> <dbl>    <dbl>  <dbl> 1     0     4 22.90000 3.1500 2     0     6 19.12500 3.4400 3     0     8 15.05000 3.8100 4     1     4 28.07500 2.0375 5     1     6 20.56667 2.7700 6     1     8 15.40000 3.3700

因為用 summarise() 算出來的結果本身就會轉換成一個 data frame，所以可以把這個 data frame 指定一個名稱，例如為 cars_am_cyl，這樣之後想要看的時候便不需要打全部的語法，只要叫出 cars_am_cyl 即可看。

cars_am_cyl <- mtcars %>%   group_by(am, cyl) %>%   summarise(mpg = mean(mpg), wt = median(wt)) cars_am_cyl

因為 cars_am_cyl 已經是一個 data frame，也可以直接用它來做運算，因為上面設定時已經先用 am 分好組了，所以算平均值的時候會依這個分。

cars_am_cyl %>%   summarise(mpg = mean(mpg), wt = median(wt))

# A tibble: 2 x      am      mpg    wt   <dbl>    <dbl> <dbl> 1     0 19.02500  3.44 2     1 21.34722  2.77

也可以相反設試試看，先用 cyl 分組後再依 am 分組運算，也就是 group_by(cyl, am)，然後和上面比較有何不同。

cars_cyl_am <- mtcars %>%   group_by(cyl, am) %>%   summarise(mpg = mean(mpg), wt = median(wt)) cars_cyl_am

# A tibble: 6 x 4 # Groups:   cyl [?]     cyl    am      mpg     wt   <dbl> <dbl>    <dbl>  <dbl> 1     4     0 22.90000 3.1500 2     4     1 28.07500 2.0375 3     6     0 19.12500 3.4400 4     6     1 20.56667 2.7700 5     8     0 15.05000 3.8100 6     8     1 15.40000 3.3700

因為上面是先用 cyl 分組，所以接下來用 cars_cyl_am 來做運算的話，會算出 cyl 裡各組的平均值。

cars_cyl_am %>%   summarise(mpg = mean(mpg), wt = median(wt))

# A tibble: 3 x 3     cyl      mpg      wt   <dbl>    <dbl>   <dbl> 1     4 25.48750 2.59375 2     6 19.84583 3.10500 3     8 15.22500 3.59000

Exercise

Q1. Compute the minimum and maximum displacement for each engine type (vs) by transmission type (am)

在各個引擎的種類中依 am 分組後，找出其中的最大值和最小值，因此要先依引擎分類，再在其中依 am 分組，也就是：group_by(vs, am)

最大值和最小值的功能為：max() 和 min()

mtcars %>%   group_by(vs, am) %>%   summarise(min = min(disp), max = max(disp))

# A tibble: 4 x 4 # Groups:   vs [?]      vs    am   min   max   <dbl> <dbl> <dbl> <dbl> 1     0     0 275.8   472 2     0     1 120.3   351 3     1     0 120.1   258 4     1     1  71.1   121

Q2. Which destinations have the highest average delays?

flights %>%   group_by(dest) %>%   summarise(avg_delay = mean(arr_delay, na.rm = TRUE)) %>%   arrange(desc(avg_delay))

# A tibble: 105 x 2     dest avg_delay    <chr>     <dbl>  1   CAE  41.76415  2   TUL  33.65986  3   OKC  30.61905  4   JAC  28.09524  5   TYS  24.06920  6   MSN  20.19604  7   RIC  20.11125  8   CAK  19.69834  9   DSM  19.00574 10   GRR  18.18956 # ... with 95 more rows

8. Boxplot

接下來試著畫箱圖，如果想畫出在 cyl < 8 的資料中，mpg 對 cyl 的反應，也就是：

x-axis: cyl < 8 y-axis: mpg, response to cyl

boxplot 的語法是：boxplot(formula, data = )

formula 是指兩軸(也就是兩個 variables)的關係，在這裡是： y ~ x 

表示 Y-axis 對 X-axis 的反應，所以也就是：boxplot(y ~ x, data = )

最後用 subset =  的功能挑出你要的，在這個例子裡就是 cyl < 8。

boxplot(mpg ~ cyl, data = mtcars, subset= cyl < 8)

也可以用 %>% 分開寫。

You can use . as a placeholder when the “data” argument is in the second position.

mtcars %>%   filter(cyl < 8) %>%   boxplot(mpg ~ cyl, data = . )

在上面的語法裡，當前面已經表示過用的資料是 mtcars 的時候，後面在 boxplot() 裡面 data 的部分就可以用 "." 替代，也就是：data = . 

也可以用 ggplot2 來畫圖，用 ggplot 裡的 geom_boxplot() 畫的話，語法就是這樣：

ggplot(mtcars %>%          filter(cyl < 8), aes(factor(cyl), mpg)) +   geom_boxplot()

因為 cyl 是數字，所以要需要先把它變成 factor，想更了解的話可以參考這篇：R | ggplot: Point plot & Box plot

也可以先把畫圖要用的資料用 subset() 的功能挑出來，也就是 cyl < 8 的部分。我們可以把挑出來的部分指定成一個新的 data frame，下面我們把指定其為 cyl_sub，然後再用它來畫圖。

cyl_sub <- subset(mtcars, cyl < 8) ggplot(cyl_sub, aes(factor(cyl), mpg)) + geom_boxplot()

好了，這篇就先到這吧。

CBC / Myth or Science: 關於大便的一些事

今天去了一個 Nerd Nite 的活動，因為這次是今年的最後一次，所以除了有科學的講題外，還有個小遊戲，小遊戲很有趣但不是重點，先講主題。

這次的講者是 BC CDC 的一位科學家：Dr. Jennifer Gardy

連結是她參與的 CBC 科學性節目，叫 Myth or Science，然後她今天講的是關於便便的：The Nature of Shit，主要是放不久前錄好的影片，但是聲音還沒後製上去，所以是邊看無聲的影片她邊講，預計在明年二月 on air 的樣子。

她這次關於便便的影片分三個部分，第一個很有趣，就是她吃了一個膠囊攝影機，攝影機全程錄下它從食道進入胃，再到大腸、小腸的過程。藉由膠囊攝影機，你可以清楚看到胃裡面長什麼樣子，它的皺摺等等。對了對了，你可以看到賁門和幽門，我看到的時候還以為是肛門 😅。然後還可以看到大腸裡面的皺摺和微生物等等，超酷的。不過當我清楚看到講者的胃裡面和大腸裡面長什麼樣子的時候，有種「哎呀，秘密都讓人看光惹，好害羞呀！」的感覺，不知道各位如果在公開場合讓大家看你胃裡面的皺摺和大腸裡面的微生物時，會不會有這種尷尬害羞的感覺。

第二個部分是玉米真的是「明天還會再見到」的食物嗎？😆 她為了做這個一連兩天每餐都吃玉米，之後再從自己的便便裡挑出玉米寄給她在多倫多 U of Guelph 的好友分析（她說：反正他一定是叫他的學生分析 XD），結果分析出來後發現便便裡拿出來的玉米上面帶有好多微生物，他們還用 EM 看玉米，好高級。之後她們跑到(美國的) Cambridge 抽取 sewage，看看裡面有哪些微生物和 chemicals。檢驗微生物可以得到一些關於疾病傳染的資料，檢驗 chemicals 可以知道這區的人都在幹麻，例如他們發現禮拜六和禮拜一早上的 sewage 裡面含大量的 Advil 和 Tylenol (就是台灣的普拿疼之類的)。（在場的大家都發出一種會心一笑的笑聲，我猜是因為週五和週日晚上喝太多酒，結果隔天頭痛所以吃止痛藥。）然後也可以靠這樣檢驗出 sewage 裡面是否有 opioid，這區的人使用鴉片的人口多少之類的。

ps. 她有提到 Guelph 的便便重製機 RePoopulator 重製便便(repoopulate)，話說研究重製便便的教授因為和我母校 Queen's 的教授合作，所以還有到本系的 seminar 演講，我是那時知道 repoopulate 的，有興趣的可以看這兩篇：

Nature / Faeces-filled pill stops gut infection (2013)

U of Guelph / Synthetic 'Poop' Can Cure C. difficile Infection, Study Finds

（標題感覺很聳動，好像是真的吃便便，但其實是和 pro-biotics 有關。）

最後一個影片是關於殺人鯨 Orca，因為溫哥華沿海的殺人鯨數量近幾年少很多，保育學專家想知道為什麼，所以要收集 Orca 的便便檢驗。他們收集的方法很有趣，用的是救生犬，他們訓練救生犬辨別 Orca 便便的味道，然後帶牠出海到海中充滿 Orca 的地方等鯨魚便便，鯨魚便便了以後大便會浮在海面上，救生犬聞到以後頭就會轉向便便的位置，船再開過去，當救生犬確定便便的位子後就會叫，然後研究員就可以把鯨魚的便便撈起來。好了，結果是什麼呢？他們分析發現是因為鮭魚變少了，不夠鯨魚吃，於是牠們出現營養不良的現象。

我覺得真的很有趣，很值得一看，明年上線以後一定要找來看啊~~

最後說一下小遊戲是什麼，遊戲叫 PowerPoint karaoke，就是主辦團體會準備好幾份 slides，每個主題都不同，重點是講者們答應了要上台報告了以後，主辦團當天才會給你 slides，也就是等講者上台了，主持人把 ppt 打開了以後，講者才知道自己要講的主題是什麼。但是因為 slides 的內容都是自己不懂的，所以就看你要在台上怎麼掰、掰功如何，笑果非常好，如果有人想辦活動又不知道要玩什麼遊戲，可以玩這個，很搞笑。XD


ps. 其實這是我第一次去 Nerd Night 的活動，因為最近才發現的。

2017 年大家在關注什麼研究？

2017 年最受人關注的前一百篇研究，大家都有發摟到嗎？XD



❖ 1. Associations of fats and carbohydrate intake with cardiovascular disease and mortality in 18 countries from five continents (PURE): a prospective cohort study. The Lancet.

第一名這篇我沒發摟到。（原文看這）


❖ 2. Work organization and mental health problems in PhD students. Research Policy

這篇應該被(菸酒生)轉爆惹。XD
我自己是沒看這篇原文，而是看 Science News 的這篇(這篇應該也被轉到爆惹 XD)

Science / Ph.D. students face significant mental health challenges
Nature / Mental health: Degree and depression
Nature / Graduate survey: A love–hurt relationship


❖ 4. Correction of a pathogenic gene mutation in human embryos. Nature

參考：把 CRISPR 用在人類胚胎上：這次換美國了~


❖ 5. Gender stereotypes about intellectual ability emerge early and influence childrenÍs interests. Science

參考：女孩從六歲開始就覺得男生比較聰明


❖ 10. An extra-uterine system to physiologically support the extreme premature lamb. Nature Comm

參考：人造子宮：Biobag


❖ 23. Pregnancy leads to long-lasting changes in human brain structure. Nature Neuroscience

這篇我有看到，但是沒看原文。是說女人第一次懷孕後大腦的 gray matter 某個部份和海馬迴有縮水了一點。用 fMRI 觀察發現 gray matter 有縮水的部分會在見到自己的小 baby 時特別活躍，顯示這種改變會讓媽媽對周遭環境比較敏感，能夠較快速的接收到嬰兒的訊息，或是感應到環境中的危險，而且這種改變會持續兩年。至於掌管記憶的海馬迴，如果沒有再次懷孕，則在兩年之後恢復。

Science / Pregnancy resculpts women’s brains for at least 2 years


❖ 27. CRISPR-Cas encoding of a digital movie into the genomes of a population of living bacteria. Nature

這篇泛科學有：電影片源哪裡找？科學家利用 CRISPR 把影像儲存在大腸桿菌的DNA！


❖ 43. A female Viking warrior confirmed by genomics

這篇泛科學也有：維京的高階戰士不可能是女生？DNA分析能推翻十九世紀以來的偏見嗎？


❖ 63. A bioprosthetic ovary created using 3D printed microporous scaffolds restores ovarian function in sterilized mice. Nature Comm.

參考：3D 列印卵巢


不知各位發摟到了幾篇呢？



Article:

Nature News / The Great Pyramid’s void, deadly heat and more: the most popular science stories of 2017

THE / These are the 10 most popular academic papers of 2017

兩個人造的 nucleic acids 可以被轉譯成新的氨基酸

除了 A, G, C, T 外的另外兩個人造 nucleic acids：X for NaM 和 Y for TPT3。

The Scripps Research Institute 的 Dr. Floyd Romesberg 實驗室把 sfGFP (superfolder GFP) 裡轉成 Ser151 的 codon 從 TAC 改為 AXC，其 anti-codon 為 GYT。跟之前不同的是這次細菌能夠把這兩個人造 nucleic acids 轉成合成氨基酸 PrK 或是 pAzF 並且活得好好的。


Figure / Y Zhang et al, Nature 2017

NaMTP 和 TPT3 和天然的 nucleic acids 不一樣，不是 H-bonding 相接合，也不會和天然的 nucleic acids 接合成對。他們把 X 放入 GFP 的基因序列裡面後，需要在培養液中加入人造的 tRNA 才能讓細菌執行正常的轉譯(translation)，成功製造出 sfGFP。如果在把 S151 改成 X，那就需要加入人造 Methanosarcina mazei tRNA^Pyl(GYT) 和合成氨基酸 PrK，才能產生 sfGFP。他們也試了其他種的 tRNA：Methanococcus jannaschii TyrRS-tRNA^Tyr (pAzFRS-tRNA^pAzF(GYT))。同樣的把 tRNA^pAzF(GYT) 和合成氨基酸 pAzF 加入培養液後，也可以讓細菌成功製造出 sfGFP。

除此之外，他們也試著換 codon，把 TXC 換成 GXC。同樣的，當他們把對應的 tRNA^Pyl(GYC) 和合成氨基酸 PrK 放進培養液中後，細菌也可以成功的把 PrK 轉進 sfGFP 的氨基酸序列裡，製造出完整的 sfGFP。

原本的 A, T, C, G 只能轉成二十個氨基酸，多了這兩個後，可以組合的方式變多了，理論上可以拓展到共 152 個氨基酸。是說多出來的那一百多個要變成怎樣的氨基酸，可能需要再想想。XD



Nature / ‘Alien’ DNA makes proteins in living cells for the first time (2017)

Science / Scientists just added two functional letters to the genetic code (2017)

The Scientist / Six-Letter DNA Alphabet Produces Proteins in Cells


Paper:

Y Zhang et al, A semi-synthetic organism that stores and retrieves increased genetic information. Nature (2017)

是否有可能取代鴉片類止痛藥的 NFEPP

目前常用的止痛藥很多是鴉片類(opioids)，像是海洛因、嗎啡等等，它們作用在 opioid receptors (ORs)。Opioid receptors 是一種 GPCRs，分佈在中樞神經系統(central nervous system, CNS)和周邊神經系統(peripheral nervious system, PNS)裡的痛覺神經網絡(nociceptive neural circuitry)，opioid receptors 被啟動後，會抑制神經細胞，中止痛覺的傳送，所以有麻醉止痛的效果。鴉片類藥物雖然止痛(analgesic)效果很好，但其副作用也不小，包括暈眩、噁心、呼吸窘迫( respiratory depression)和上癮等症狀。近來北美最為氾濫的此類藥物當屬 fentanyl，其效用比嗎啡高上八十倍，因為可以自己合成，比海洛因便宜，且只要一點點就可以就可以達到效果，常被混在海洛因或夜店藥裡，因此容易讓人不小心過量致死，加拿大 BC 省光今年上半年(一到七月)因 fentanyl 過量致死的就有九百多人。

ORs 分為四種：mu (μ), kappa (κ), and delta (δ) opioid receptors 和 opioid receptor like-1 (ORL1)，最常用來止痛的鴉片類藥物多是作用在 μ opioid receptors (MORs)，但因為 MORs 也參與了腦內回饋系統(reward system)，所以也會產生上癮的現象。

這篇研究的作者們認為，鴉片類藥物之所以會有其他副作用，是因為 ORs 分佈在中樞神經和周邊神經系統，當藥物作用在中樞神經系統的 ORs，就會出現嘔吐或呼吸困窘等症狀，而痛覺接收和傳送是發生在周邊神經系統，如果能讓藥物只作用在疼痛部位的神經細胞，也許便可以避免掉副作用。在慢性疾病造成的發炎和疼痛情況下，會有細胞組織酸化的現象。跟正常的細胞環境相比，受傷的組織其酸鹼值比較低，因此他們比較了在正常情況下和發炎情況下，藥物和 MORs 的反應，希望能找到一個只有在酸性環境下才會啟動 MORs 的藥物，以避免在正常情況下也會啟動 MORs 而產生包括上癮等等的副作用。以 fentanyl 為例，它的 pKa (disscoiation constant) 大於 8，因此在正常生理情況下(pH7.4)和發炎組織(pH5-7)的情況下都會被質子化(protonated)，進而啟動 MORs。作者者們認為在正常情況下啟動 MORs 會產生不必要的副作用，所以想找一個 pKa 小於七，只有在組織發炎、需要止痛的情況下才會啟動 MORs 的物質。

通常要讓一個化學物質的 pKa 降低，可以把它的 H (hydrogen) 改成 F (fluorine)，少了一個可以做 H-bond 的連結，就會降低兩者(ligand & MOR)之間的 interaction。他們改造了 fentanyl 的幾個地方後，選中了這個 pKa = 6.8 的 NFEPP ((±)-N-(3-fluoro-1-phenethylpiperidin-4-yl)-N-phenyl-propionamide)，因為他們認為當 pKa 在 6~7 之間的話，組織發炎時的酸鹼值就足以質子化 NFEPP。


Figure: V Spahn et al, Science (2017)

跟 fentanyl 相比，NFEPP 在正常的生理情況下不會啟動 MORs，靜脈注射(i.v.)只會使發炎的部位依劑量的不同出現不同程度的麻醉、止痛(analgesia)效果，並且高劑量(75g/kg)的注射也沒產生呼吸困窘的副作用。它不像 fentanyl 在低劑量(8-12ug/kg)時即可在健康部位和發炎部位皆出現明顯的止痛效果，並且在其三分之一的量(32ug/kg)下就會造成致死的呼吸窘迫症狀。除此之外，fentanyl 在皮下注射(s.c.)為 30-50ug/kg 的劑量時會產生的依賴藥物的現象，然而 NFEPP 在高劑量(150ug/kg)下並沒有讓老鼠出現上癮的現象。Fentanyl 在 20-60ug/kg 會依劑量高低產生不同程度的四肢麻痺和便秘情形，NFEPP 則在高劑量(150ug/kg)下也未出現這些症狀。

看來降低 pKa 是能解決 fentanyl 產生的問題，不過他們可能需要多加一個情況做比較。在這個研究裡，試驗的老鼠都是有發炎症狀的，只是一隻腳有一隻腳沒有，然後兩隻腳做比較。保險起見，也許要加一個完全健康沒任何發炎症狀的老鼠做比較，看看是否在健康的老鼠身上也不會出現上癮等等的副作用。



Article:

TM Villanueva, Analgesia: Designing out opioid side effects. Nature Review Drug Discov (2017)

PRF / A New Opioid Targets Active Sites of Inflammation to Relieve Pain


Paper:

V Spahn & G Del Vecchio al, A nontoxic pain killer designed by modeling of pathological receptor conformations. Science (2017)

戒酒藥可以殺死癌細胞

有意思哩！六十年來被用來當戒酒藥的 disulfiram（商品名 Antabuse / Antabus / 戒酒硫）也可以抗癌。在 1971 年的時候有篇臨床報告說一位 38 歲的乳癌病患因為癌細胞擴散到脊髓無法醫治而開始酗酒，醫生為了讓她戒酒便給她吃戒酒藥。十年後病患過世了，沒想到解剖後發現，她的腫瘤消失了，脊髓裡只剩一點癌細胞而已。在那之後，有幾個臨床研究和老鼠實驗也顯示 disulfiram 可以阻止癌細胞生長，雖然如此，disulfiram 並沒有因此受到很大個關注。

Drug Bank: disulfiram (tetraethylthiuram disulfide, DSF)

Disulfiram (DSF) 作用在 aldehyde dehydrogenase (ALDH)，ALDH 是酒精代謝過程中需要的酵素(如圖)，負責把 acetaldehyde 轉換成 acetate，disulfiram 抑制了 ALDH 後會造成 acetaldehyde 堆積在體內，會讓人感到極度的不適。



在這期的 Nature 中有篇由美國、捷克和丹麥合作的研究，檢視了在 2000 年到 2013 年間 240,000 位癌症患者和他們的用藥，發現在三千位有用 Antabuse 的患者中，有持續用藥的患者死亡率比後來停用的患者低惹 34%。這研究中檢視的不只是乳癌，還有攝護腺癌、直腸癌和其他種癌症，disulfiram 對它們的效果都是一樣的。此研究也顯示 disulfiram 可以減緩老鼠體內乳癌腫瘤的生長，並且和銅(copper, Cu)一起使用的話效果更好，喂有 disulfiram 的老鼠其腫瘤比沒有吃的小 57%，而和銅(DSF/Cu)一起吃的老鼠其腫瘤比沒吃的小 77%。他們也比較了這兩種飲食在老鼠體內代謝情形，disulfiram 在體內的主要代謝物 ditiocarb (diethyldithiocarbamate, DTC)會和銅結合成 CuET。吃了 DSF/Cu 的老鼠，CuET 在肝臟和血清(plasma)的含量比只吃 DSF 的高(感覺是廢話 XD)，但有趣的是吃 CuET 的老鼠其腫瘤細胞裡的 CuET 量比其他器官高(只吃 DSF 的老鼠也有同樣的情形，但沒有高那麼多)，而在細胞株上的檢測也顯示 CuET 的毒性(toxicity)比代謝物 DTC 本身還高。

他們發現 CuET 會和 NPL4 結合，NPL4 是 ATPase p97 (VCP) 的 adaptor protein，它通常會和另一個 cofactor, Ufd1, 一起和 p97 結合變成一個 complex (Npl4-Ufd1-p97)。p97 有很多不同的 cofactors，和不同的 cofactor 結合會執行不同的功能，其中主要的是把被 ubiquitin (Ub) 標記的蛋白質(Ub-proteins)從細胞內的細胞膜(i.e., ER membrane)中或是 protein complex 中抓出來，讓其能夠被降解(degradation)，它參與的蛋白質降解管道有 proteasomal degradation 和 ERAD (ER-associated degradation)，而和 NPL4 合作後的功用是負責 Ub-proteasome-mediated degradation，但是 CuET 和 NPL4 結合的話便會抑制它的功能。p97 和 NPL4 都是 cytosolic proteins，但當把 CuET 加到細胞上後，NPL4 無法執行功能，Ub-proteins 無法被降解而堆積在細胞中，引起後續的死亡反應。他們發現 CuET 會引發細胞的 UPR (unfolded protein response) 和 HSR (heat-shock response)。

所以目前可能的情況是 disulfiram 的代謝物 DTC 和銅結合後會抑制 p97-NPL4-Ufd1 參與的蛋白質降解功能，大量無法被降解的 Ub-proteins 在細胞裡堆積，引發癌細胞的死亡反應。因為 CuET 只會堆積在癌細胞腫瘤，所以這也可能是為什麼正常人吃 disulfiram 除了不適外沒有其他問題，只是 ..... 為什麼 CuET 只會堆積在癌細胞腫瘤中，原因還不清楚。

另外一個我覺得有趣的點是根據 Drug Bank 的資料，disulfiram 還有一個功能是：inhibits the peripheral benzodiazepine receptor.



Article:

Science / An old drug for alcoholism finds new life as cancer treatment


Paper:

Z Skrott et al, Alcohol-abuse drug disulfiram targets cancer via p97 segregase adaptor NPL4. Nature (2017)

關於偏頭痛(migraine)

（舊文轉過來）

你有偏頭痛嗎？這兩天 AAAS 在聖地牙哥有個年會，今天便在旗下的《Science》發了篇文章，告訴大家他們和 University of Miami Health System 的神經學家 Teshamae Monteith 討論了什麼關於偏頭痛的最新進展。偏頭痛不是普通的頭痛，它是一陣一陣難以忍受的痛，外加通常會噁心反胃，而且對聲音、光和味道很敏感。（嗯，這樣我好像沒偏頭痛過耶。）全世界大概有十分之一的人有偏頭痛，雖然目前無藥可治，但科學家們仍努力解開偏頭痛之謎，看能不能找出有效的治療方法。

目前我們對偏頭痛的了解是什麼？

這是個複雜的疾病，之前研究學者們都以為是跟血管相關的疾病，因為有的患者會在偏頭痛發作時感覺到血管規律的脈動。（阿，這我有過。）現在呢，科學家則認為應該是跟感官有關的疾病，因為很多的感官，包括光線、聲音、味覺、聽覺等等，都有受到影響。當偏頭痛發作時，患者常常無法專心、胃口和心情都會改變，也很難入眠。最引起 Teshamae Monteith 注意的是偏頭痛患者常被各種不同的徵狀困擾，像是在每波的陣痛之間對光的敏感度會增加，表示他們的神經或傳導系統各有不同和改變。大約三分之二的急性偏頭痛患者有異常性頭痛(allodynia)，會對某種刺激特別敏感，例如沖澡時產生的水蒸氣都會令他們疼痛難當。不同的偏頭痛徵狀可能是因為每個人感官對刺激的敏感度不同。

目前偏頭痛的最新療法是什麼？

現在最主要的治療是一種叫 triptans 的藥物，它作用在血清素接收器(serotonin receptors)，血清素被認為是參與偏頭痛的神經傳遞物質(neurotransmitter)，血清素(5-HT)在平常的時候是低的，但偏頭痛發作的時候會升高。憂鬱症和偏頭痛也有強烈的關係，憂鬱症被指出和異常的血清素量有關，而憂鬱症患者較有可能有偏頭痛，偏頭痛患者也比較容易憂鬱。目前還不知道 triptans 是如何作用的，但它的確對某些偏頭痛患者有用，能夠止痛。它是個好藥，但可惜不是對每個人都有用，所以目前在研發各種新藥，希望是不像 triptans 這種會造成血管收縮的藥。

什麼是目前最有可能的新藥？

有個新藥是針對 CGRP (calcitonin gene-related peptide) ，它是一個被認為在急性偏頭痛發作時釋放出來的其中一種 peptide。在一個隨機、雙盲的預防性治療研究中，給患者不同劑量的 CGRP 藥物 telcagepant，雖然數據顯示效果不錯，但因為在一小部分有服用藥物的患者中，發現他們的某些肝臟酵素出現異常，以至於此藥物無法通過安全測試。不過現在有種抗體藥已進入第二階段(Phase II)的測試，因此被認為是目前最熱門的新東西，雖然第一階段的測試數量太小，無法確定它的效用有多大，但仍然是令人興奮的，因為它是第一個顯示安全的 CGRP 藥物。

偏頭痛遺傳的機率有多大？

關於致病基因的部分還需要更多的研究，因為偏頭痛這個病症太複雜，但目前大多數的患者都有家族史。Monteith 在 2008 年當住院醫師的時候，發現三個基因跟偏癱性偏頭痛(hemiplegic migraine)有關，這個病會使身體有一邊無力。而現在，有一堆基因被發現和偏頭痛有關，其中有些和 glutamate 有關。glutamate 是主要興奮性(excitatory)神經傳導物質，其功能是激化神經細胞，但因為腦部系統很多是用 glutamate 來傳遞訊息，所以要抑制它的功能有困難度。

偏頭痛的致病因子是什麼？

偏頭痛通常在青春期或青春期前發作，女性通常在生理期前後發作，男性的頭痛症狀也和他們的賀爾蒙變化有關。較低的社經地位、肥胖和不好的睡眠品質都會增加偏頭痛發作的頻率，腦部組織損傷(infarctions)是偏頭痛的併發症，其他像是抽菸和口服避孕藥則會增加偏頭痛的機率。

研究偏頭痛最困難的地方在哪？

我們(Teshamae Monteith)目前遇到的瓶頸是動物實驗的部分，因為牠們沒辦法和人腦完全一樣，也無法有極類似的偏頭痛症狀。人類的腦部影像技術有比較大的發展空間，在偏頭痛研究的應用上也比較讓人興奮。透過腦部影像，我們能夠更進一步了解跟偏頭痛有關的腦部構造、神經連結和化學物質。另外，研究曾經有過偏頭痛的已逝患者大腦，也可能對了解和偏頭痛有關的腦部損傷有所幫助，不過這種研究有其難解釋的部分，因為造成死亡的原因，例如中風、癌症等等，都會改變大腦狀態。（編按：如果大腦狀態因為真正的死因而改變，就無法知道哪些改變是真的跟偏頭痛相關，或是因為偏頭痛而造成的。）



原文：E Underwood. The Science of Migraines. Science News (Feb 2015)

R | ggplot: summary & error bar

這篇要講的應該是很常會用到的，就是加 error bar。

課堂講義請看這：R Graphics with Ggplot2 - Day 3

這邊我們要用到檔案資料 diamonds。

library(ggplot2) library(dplyr) data(diamonds)

我們先用 xtabs() 來看檔案中 color 的分佈情形。（xtabs() 的功能請看這篇。）

xtabs(~ color, diamonds)

Table 1

color     D     E     F     G     H     I     J  6775  9797  9542 11292  8304  5422  2808

下面用 bar graph 看分佈，然後和上面做比對。

ggplot(diamonds, aes(color)) + geom_bar()

也可以分開設：（這個部分可以參考這篇會比較清楚）

d <- ggplot(diamonds, aes(color)) d + geom_bar()

#1 (color only)


接著，來看看每種 color 的價格，這裡可以用 by() function 來看每種鑽石顏色的價格總整理(summary)。

語法：by(response var, predictor var, summary)

這裡的 predictor variable (又稱 explanatory or independent variable)通常是個 factor (不是數字)，會放在圖的 X-axis，response (or dependent) variable 指的是對應 factor 出現的反應，因此通常是放在圖的 Y-axis。

by(diamonds$price, diamonds$color, summary)

Table 2

diamonds$color: D    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     357     911    1838    3170    4214   18690 ---------------------------------------------------------- diamonds$color: E    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     326     882    1739    3077    4003   18730 ---------------------------------------------------------- diamonds$color: F    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     342     982    2344    3725    4868   18790 ---------------------------------------------------------- diamonds$color: G    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     354     931    2242    3999    6048   18820 ---------------------------------------------------------- diamonds$color: H    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     337     984    3460    4487    5980   18800 ---------------------------------------------------------- diamonds$color: I    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     334    1120    3730    5092    7202   18820 ---------------------------------------------------------- diamonds$color: J    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     335    1860    4234    5324    7695   18710

也可以用 tapply() 功能來看。

tapply(diamonds$price, diamonds$color, summary)

$D    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     357     911    1838    3170    4214   18690 $E    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     326     882    1739    3077    4003   18730 $F    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     342     982    2344    3725    4868   18790 $G    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     354     931    2242    3999    6048   18820 $H    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     337     984    3460    4487    5980   18800 $I    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     334    1120    3730    5092    7202   18820 $J    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     335    1860    4234    5324    7695   18710

如果只想看平均值，就把 summary 改成 mean。

tapply(diamonds$price, diamonds$color, mean)

D        E        F        G        H        I        J 3169.954 3076.752 3724.886 3999.136 4486.669 5091.875 5323.818

在圖裡，可以用 summary_bin 這個功能來看，語法是：stat = "summary_bin"

或是：stat = "summary" （這兩個是一樣的）

summary_bin 裡面的功能有：fun.data, fun.y, fun.ymax (設最大值), fun.ymin (設最小值)

ggplot(diamonds, aes(color, price)) +        geom_bar(stat = "summary_bin", fun.y = "mean")

在上面的語法裡，Y-axis 雖然已經設好是 price，但你需要在 geom_bar() 裡面的設定 Y-axis 是 price 的平均值(mean)或是中間值(median)等等，設定的語法是：fun.y = 

也可以用上面的已設定好的 d <- ggplot(diamonds, aes(color)) 把上面的語法改成下面這樣：（上面的語法是把 Y-axis 設在 ggplot() 裡面，但是 d 裡面只設了 X-axis，所以在下面的語法中，需要把 Y-axis 設在 geom_bar() 裡面。）

d + geom_bar(aes(y = price),              stat = "summary_bin", fun.y = "mean")

#2 (color vs. price mean)


可以把圖和上面的 summary 表格做比較，看畫出來的是不是 mean。（和上面不同的是表一和圖一的指的是每種顏色的鑽石有幾個，所以圖 #1 的 Y-axis 是 count，而表二和圖二不只是每種顏色有幾個，而是每種鑽石顏色的價格，所以它的 Y-axis 指的不只是價格，還要指定是每種鑽石顏色的平均價格，還是每鑽石顏色其價格的中間值等等。）

下面來畫看看 median 的圖。

d + geom_bar(aes(y = price),              stat = "summary_bin", fun.y = "median")

#3 (color vs. price median)


除了用 geom_bar()，我們也可以直接用 stat_summary_bin() 或 stat_summary() 的功能。

stat_summary: summarise y values at distinct x values.

ggplot(diamonds, aes(color, price)) +        stat_summary_bin(fun.y = "mean", geom = "bar")

或是用上面設好的 d。stat_summary_bin() 和 stat_summary() 是一樣的東西，上面和下面的這兩個語法畫出來的圖和上面用 geom_bar() 畫出來的 mean 那張圖是一樣的(圖 #2)。

d + stat_summary(aes(y = price),                  fun.y = "mean", geom = "bar")

我們也可以畫點圖，用 geom_point() 就是下面那樣：

ggplot(diamonds, aes(color, price)) +        geom_point(stat = "summary_bin", fun.y = "mean")

#4


上圖的 Y-axis 不是從 0 開始，而是從三千，我們可以把它設成從 0 開始。下面用上面設好的 d 和 stat_summary來畫圖，但基本上跟上面是一樣的，只是多了 ylim() 的語法。

ylim(min, max): Y-axis 的最小值和最大值

在下面的語法中，Y-axis 的最小值我們設為 0，最大值沒設，用 NA。

d + stat_summary(aes(y = price),                  fun.y = mean, geom = "point") +     ylim(0, NA)

#5


可以看到上圖的 Y-axis 是從 0 開始，而不是像圖 #4 那樣從三千開始。

接下來，試試在圖中加個 error bar，有四種方法。

注意，下面四種的預設都是 standard error (SE)，除非你有指定另外的算法在裡面。

1. geom_pointrange(): vertical line with centre (中間點和其 error bar / standard error)

ggplot(diamonds, aes(color, price)) +        geom_pointrange(stat = "summary_bin")

我們可以改變點和其 error bar 的大小和顏色。

ggplot(diamonds, aes(color, price)) +        geom_pointrange(stat = "summary_bin",                        size = 0.2, color = 'blue')

## No summary function supplied, defaulting to `mean_se()

summary_bin 裡面的 error bar 預設是 standard error (SE)，你也可以把他設成別的，例如把它設成最大值 fun.ymax = max 和最小值 fun.ymin = min。

ggplot(diamonds, aes(color, price)) +        geom_pointrange(stat = "summary",                        fun.y = mean,                        fun.ymin = min,                        fun.ymax = max)

可以把 pointrange() 和其他圖合在一起，這裡要注意一下，放在後面的會覆蓋住前面的，所以需要把 geom_point() 放在 geom_pointrange() 後面，如果相反的話會在圖裡看不出 geom_point() 裡面的設定。下面把兩個用不同顏色和不同大小的點區分。

ggplot(diamonds, aes(color, price)) +        geom_pointrange(stat = "summary_bin",                        size = 0.5, color = 'red') +        geom_point(stat = "summary_bin", fun.y = "mean",                   size = 1.5, color = 'blue')

Error bar 除了以上兩種外，也可以設其他的，例如讓它的範圍為中間 50%，最小值為第 25%，最大值為第 75%。

下面先用 group_by() 的功能把 color 分類，然後再用 summarize() 的功能在每種顏色裡算出其的價格。先把每種顏色的平均價格設為 y = mean(price)，把誤差值(SE)的最小值指定為 ymin，設其為每種鑽石顏色的第 25%；把最大值指定為 ymax，設其為每種顏色的第 75%。

ps. quantile 分為：0 (0%), 0.25 (25%), 0.75 (75%), 1 (100%)

在下面先把算出來的指定為一個 data frame：df_color_price

df_color_price <- group_by(diamonds, color) %>%   summarize(ymin = quantile(price, 0.25),             ymax = quantile(price, 0.75),             y = mean(price)),             n = n())

# A tibble: 7 x 5   color   ymin    ymax        y     n   <ord>  <dbl>   <dbl>    <dbl> <int> 1     D  911.0 4213.50 3169.954  6775 2     E  882.0 4003.00 3076.752  9797 3     F  982.0 4868.25 3724.886  9542 4     G  931.0 6048.00 3999.136 11292 5     H  984.0 5980.25 4486.669  8304 6     I 1120.5 7201.75 5091.875  5422 7     J 1860.5 7695.00 5323.818  2808

上面語法中的 n() 是指在那個分類裡有幾個，在 R 裡面的解釋是："The number of observations in the current group."（也就是說，上面得出的 n 的數字應該是要表 #1 裡的是一樣的。）

從上面的表可以看出 summarize() 算出每種鑽石顏色的平均價格(y)和前後 25% 的價格(ymin, ymax)，還有每種顏色的鑽石有多少個(n)，接著我們用這個 data frame 畫圖。

ggplot(df_color_price, aes(color, y,                        ymin = ymin, ymax = ymax)) +        geom_pointrange() +        geom_point(aes(size = n))

在前面 df_color_price 的語法裡加入 n = n() 是要讓點圖裡的點會依據每種顏色的數量呈現出不同的大小，例如在前面 pointrange() 的例子裡加入 size =  的話可以指定點的大小。之前都是設數字（例如 size = 2），但是在這裡設 size = n 的話，因為在前面已經設了 n = n()，表示 n 是每種顏色的數量的話，那點的大小就會因數量而異。

也可以把上面兩個語法用 %>% 合在一起。

group_by(diamonds, color) %>%          summarize(ymin = quantile(price, 0.25),                    ymax = quantile(price, 0.75),                    y = mean(price),                    n = n()) %>% ggplot(aes(color, y, ymin = ymin, ymax = ymax)) +        geom_pointrange() +        geom_point(aes(size = n))

也可以和 geom_bar() 合用。

ggplot(diamonds, aes(color, price)) +        geom_bar(stat = "summary", fun.y = "mean") +        geom_pointrange(stat = "summary",                        size = 0.1, color = 'blue')

在上面的圖中，因為每種鑽石顏色的 standard error (SE) 很小，所以看不出來。

可以用上面已經設好的 df_color_price 來畫圖：

ggplot(df_color_price, aes(color, y,                        ymin = ymin, ymax = ymax)) +         geom_bar(stat = "summary") +         geom_pointrange(size = 0.2, color = 'blue')

因為在 geom_pointrange() 裡面並沒有設 stat = "summary"，所以圖呈現出來的會是設在第一層 ggplot() 裡的 error bar (ymin, ymax)。

也可以用 %>% 把兩個語法合在一起，變成下面這樣：（這裡沒用 n = n() 是因為並非畫點圖，所以不需要。）

group_by(diamonds, color) %>%   summarize(ymin = quantile(price, 0.25),             ymax = quantile(price, 0.75),             y = mean(price)) %>% ggplot(aes(color, y, ymin = ymin, ymax = ymax)) +        geom_bar(stat = "summary") +        geom_pointrange(size = 0.2, color = 'blue')

下面把 geom_pointrange() 和 geom_line() 合在一起用。

ggplot(diamonds, aes(color, price)) +      geom_pointrange(stat = "summary_bin",                      size = 0.2, color = 'steelblue') +      geom_line(stat = "summary_bin", fun.y = "mean",                group = 1, alpha = 1/2)

geom_line() 裡面的語法可以參考這篇：ggplot | Regression Line & Grouping

2. geom_linerange(): vertical line (只有一條直線 error bar / standard error)

如果單用的話就會像下圖那樣只有一條線，所以通常會和其他圖和用。

ggplot(diamonds, aes(color, price)) +        geom_linerange(stat = "summary_bin")

可以和點圖和 bar graph 合用。和點圖一起用的話，其實畫出來的圖會跟只用 pointrange 差不多。你也可以幫 geom_linerange() 指定顏色，就是那條代表 SE 的直線，可以設成和點不一樣的顏色，如果直接用 pointrange() 的話就無法用不同的顏色。

ggplot(diamonds, aes(color, price)) +        geom_point(stat = "summary_bin", fun.y = "mean",                   size = 2, color = 'skyblue') +        geom_linerange(stat = "summary_bin", color = 'blue')

同樣可以和 bar graph 合用。

ggplot(diamonds, aes(color, price)) +        geom_bar(stat = "summary_bin", fun.y = "mean") +        geom_linerange(stat = "summary_bin", color = 'red')

預設的 SE 跟之前 pointrange 的圖一樣太小看不出來，下面再用之前設的 df_color_price 再畫一次。

ggplot(df_color_price, aes(color, y,                        ymin = ymin, ymax = ymax)) +        geom_bar(stat = "summary") +        geom_linerange(size = 0.5, color = 'steelblue')

geom_linerange() 裡面的 size =  是條線條的粗細。

3. geom_errorbar(): error bars.

如果只有 error bar 的話是這樣：

ggplot(diamonds, aes(color, price)) +        geom_errorbar(stat = "summary_bin")

通常都是和 bar graph 合在一起用。

ggplot(diamonds, aes(color, price)) +        geom_bar(stat = "summary", fun.y = "mean") +        geom_errorbar(stat = "summary",                      width = 0.1, color = 'grey')

geom_errorbar() 裡面的 width =  是用來設 error bar 的寬度。

下面用上面設好的 df_color_price 再畫一次。

ggplot(df_color_price, aes(color, y,                            ymin = ymin, ymax = ymax)) +        geom_bar(stat = "summary") +        geom_errorbar(width = 0.1, color = 'steelblue')

4. geom_crossbar(): vertical bar with center.

終於來到最後一個。只用 geom_crossbar() 的話是這樣：

ggplot(diamonds, aes(color, price)) +        geom_crossbar(stat = "summary_bin")

用上面的 df_color_price 畫圖：

ggplot(df_color_price, aes(color, y,                            ymin = ymin, ymax = ymax)) +        geom_crossbar(width = 0.2, color = 'steelblue')

接下來我想用檔案資料 mpg，用裡面的 class (X-axis) 和 cty (Y-axis) 畫圖。下面跟前面一樣先做總整理。

by(mpg$cty, mpg$class, summary)

mpg$class: 2seater    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    15.0    15.0    15.0    15.4    16.0    16.0 ---------------------------------------------------------- mpg$class: compact    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.   15.00   18.00   20.00   20.13   21.00   33.00 ---------------------------------------------------------- mpg$class: midsize    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.   15.00   18.00   18.00   18.76   21.00   23.00 ---------------------------------------------------------- mpg$class: minivan    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.   11.00   15.50   16.00   15.82   17.00   18.00 --------------------------------------------------------- mpg$class: pickup    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.       9      11      13      13      14      17 ---------------------------------------------------------- mpg$class: subcompact    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.   14.00   17.00   19.00   20.37   23.50   35.00 ---------------------------------------------------------- mpg$class: suv    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    9.00   12.00   13.00   13.50   14.75   20.00

把點圖和 crossbar 合用：

ggplot(mpg, aes(class, cty)) +        geom_point(alpha = 0.5) +        geom_crossbar(stat = "summary",                      width = 0.2, color = 'red')

也可以像上面一樣，把 error bar 設為中間的 50%。

df_mpg <- group_by(mpg, class) %>%   summarize(cty_ymin = quantile(cty, 0.25),             cty_ymax = quantile(cty, 0.75),             y = mean(cty)) ggplot(df_mpg, aes(class, y,                    ymin = cty_ymin,                    ymax = cty_ymax)) +        geom_crossbar(width = 0.2, color = 'steelblue') +        ylim(10, 35)

我把上面的 Y-axis 設了最大值和最小值 ylim(10,35)，讓它和上面點和 crossbar 合圖的 scale 相同，比較好做比較。上面總整理的表格有 1st Qu. 和 3rd Qu.，Qu. 即是 quantile；1st Qu. 就是第一個 quantile，即第 25% ；3rd Qu. 是第三個 quantile，即第 75%。也就是說 1st Qu. 和 3rd Qu. 即是我們在 df_mpg 裡設的 cty_ymin 和 cty_ymax。可以把 crossbar 的最大值和最小值和表格裡的 1st. Qu 和 3rd Qu. 做比對，看看是否一樣。

下面同樣是把點圖和 crossbar 合用，但是用 stat_summary() 來畫。在 error bar 的部分另外設成 sd，用的是 package "Hmisc" 裡面的 mean_sdl。

mean_sdl: smean.sdl computes the mean plus or minus a constant times the standard deviation

library(Hmisc) ggplot(mpg, aes(class, cty)) +        geom_point(alpha = 0.5) +        stat_summary(aes(y = cty), fun.y = mean,                     fun.data = mean_sdl,                     width = 0.2, color = 'pink',                     geom = "crossbar")

mean_cl_normal: computes 3 summary variables: the sample mean and lower and upper Gaussian confidence limits based on the t-distribution.

ggplot(mpg, aes(class, cty)) +        geom_point(alpha = 0.5) +        stat_summary(aes(y = cty), fun.y = mean,                     fun.data = mean_cl_normal,                     width = 0.2, color = 'pink',                     geom = "crossbar")

最後，來練習一下。

比較 cut 和 depth 的關係，用點畫出中間值(median)和中間 50% 的 error bar。

group_by(diamonds, cut) %>%   summarize(ymin = quantile(depth, 0.25),             ymax = quantile(depth, 0.75),             y = median(depth),             n = n()) %>%   ggplot(aes(cut, y, ymin = ymin, ymax = ymax)) +   geom_pointrange(color = 'steelblue') +   geom_point(aes(size = n), color = 'skyblue')

【 其他補充 】

有的地方可以看到 stat = "identity"，這和 stat = "summary_bin" 有什麼不同呢？

stat = "identity" 通常是用在每個項目只有一個，例如在 diamonds 的檔案裡面，每個 color 的觀察只有一個，也就是說當你畫 price vs. color 的圖的時候，每個 color 裡不需要算價格的平均值或找它的中間值，因為就只有一個數(一個觀察)可以畫在 Y-axis。但是當資料變多的時候，也就是每種 color 裡面有好幾個，每個都是不同的價格的時候，那你畫 price vs. color 時就需要算出每種 color 的平均值，才能畫出 Y-axis，這時候就需要用 stat = "summary_bin" 來算出平均值或找出中間值。

另外，當你沒設 stat = " "  的時候，geom_bar() 裡的預設是 stat = "bin"。

在 ggplot2: geom_bar() 網頁裡的解釋是這樣：

The heights of the bars commonly represent one of two things: either a count of cases in each group, or the values in a column of the data frame. By default, geom_bar uses stat="bin". This makes the height of each bar equal to the number of cases in each group, and it is incompatible with mapping values to the y aesthetic. If you want the heights of the bars to represent values in the data, use stat="identity" and map a value to the y aesthetic.



其他網頁參考：

My ggplot2 cheat sheet: Search by task

R document / ggplot2: geom_bar() function

ggplot 2 / Bar charts

ggplot2 / Vertical intervals: lines, crossbars & errorbars

ggplot2 / Summarise y values at unique/binned x