kur0cky

こんにちは

「大規模計算時代の統計推論」 を全部Rでやってみる ~第1章~

かの有名な「カステラ本」の姉妹編?「大規模計算時代の統計推論―原理と発展―」の和訳が発売されました

著者はブラッドリー・エフロン,トレヴァー・ヘイスティという超レジェンド研究者達ですが,訳者にもそうそうたる名前が並んでいます.

せっかく新たなバイブルに出会えたので,この本で行われている解析を全てRでやってみる,ということをやっていきます. 解析には随時適当なライブラリを使用し,作図・作表には基本的にggplot2, gtパッケージを使います.tidyverseも可能な限り活用していきます.

今回の記事では第1章「アルゴリズムと推論」の再現をします.

準備

まずは,みんな大好きtidyverseを読み込み,ggplot2で日本語フォントが使えるように設定しておきます.

library(tidyverse)
theme_set(theme_classic(base_family = "HiraKakuPro-W3"))

1.1 回帰の例

本節では,157人の被験者について調査された腎機能のデータを用いた分析が示されています. データは2変量であり,ひとつは目的変数である腎機能の健全度tot,もうひとつは説明変数である年齢ageです.

データは公式サイトに上がっているので,tidyverseにも含まれているreadrパッケージで読み込みます. urlに飛べばわかるのですが,元ファイルがスペース区切りなのでread_delim()delim = " "を指定します.

kidney <- read_delim(
  file = "https://web.stanford.edu/~hastie/CASI_files/DATA/kidney.txt",
  delim = " ",
  col_types = "dd")
head(kidney)
# A tibble: 6 x 2
    age   tot
  <dbl> <dbl>
1    18  2.44
2    19  3.86
3    19 -1.22
4    20  2.3 
5    21  0.98
6    21 -0.5 

fig 1.1

それでは分析を行っていきます.ここでは2つの作業が行われています.

  • totageで線形回帰し,回帰直線を求める.
  • 20歳から80歳まで,10歳おきに予測値の標準誤差を求め,エラーバーを追加する.

とりあえず線形回帰し,切片と傾きを得ましょう.

fit_lm_kidney <- lm(tot ~ age, data = kidney)
beta0 <- fit_lm_kidney$coefficients[1]
beta1 <- fit_lm_kidney$coefficients[2]

続いて予測値の標準誤差を求めます.
これはあまり知られていないかもしれませんが,lm()で出力されるモデルでは,predict()のなかでse.fit = TRUEとしてやることで,説明変数の各値での標準誤差を計算することができます

予めage = 20, 30, ..., 80のデータフレームを作っておき,それぞれに対する予測値,標準誤差を求めましょう.

ages <- tibble(age = seq(20, 80, 10))
df = predict(fit_lm_kidney, ages, se.fit = TRUE) %>% 
  data.frame() %>% 
  bind_cols(ages) %>% 
  select(age, 
         fit_lm = fit,
         se_lm = se.fit)
print(df)
  age     fit_lm     se_lm
1  20  1.2882585 0.2066481
2  30  0.5023743 0.1549648
3  40 -0.2835098 0.1473983
4  50 -1.0693940 0.1893143
5  60 -1.8552781 0.2575947
6  70 -2.6411623 0.3365586
7  80 -3.4270465 0.4202260

結果が出たので,図を作成します.

agetotの散布図の上に,先ほど求めた回帰係数で直線を引き,さらにエラーバーを描きます. エラーバーはgeom_segment()で線分を描きますが,本書では標準誤差の±2倍の長さで引いているらしいので,そのようにします.

できるだけ元の図を再現するため,いろいろなオプションを追加しています.

p <- kidney %>% 
  ggplot(aes(age, tot))+
  geom_point(colour = "steelblue", shape = 8)+
  geom_abline(intercept = beta0, slope = beta1, colour = "red")+
  geom_segment(data = df,
               aes(x = age, xend = age,
                   y = fit_lm - 2*se_lm, yend = fit_lm + 2*se_lm))+
  labs(x = "年齢", y = "複合尺度tot")

f:id:kur0cky:20200806000055p:plain

fig 1.2

次は,lowess (locally weighted scatterplot smoother)という手法を用いた平滑化曲線の推定が行われています. 線形回帰と異なり,lowessでは標準誤差を陽に表すことができません.そこで,反復計算で推論を行うブートストラップ法で計算することになります.

lowessのアルゴリズムの詳細については割愛しますが,詳しくはwikipediaなどを見てください. en.wikipedia.org

Rにはstats::lowess()があるので,早速実行していきましょう.本書にあるとおりlowess(x, y, 1/3)を実行してみます.

lowess(x = kidney$age,
       y = kidney$tot,
       f = 1/3)
$x
[1] 18 19 19 20 21 21 ...

$y
[1] 1.338870 1.239242 1.239242 1.145344 1.057087 1.057087 ...

これはまずいですね. 入力に対する回帰結果だけが返ってきています.10歳ごとのエラーバーが欲しいのに,データに存在しない40歳や50歳の標準誤差を得ることができません

方針を変え,loess(locally estimated scatterplot smoother) でできる限り同じように解析したいと思います.loessについてもアルゴリズムの詳細は上記リンクなどをご覧ください.

まずは,データを全て使い曲線を推定します.回帰値もさきほどのdfに追加しておきましょう

fit_loess_kidney <- loess(tot ~ age, kidney, 
                          span = 1/3, 
                          degree = 1,
                          family = "symmetric")
df <- df %>% 
  mutate(fit_loess = predict(fit_loess_kidney, age))

続いてブートストラップ法を用いて標準誤差を推定します.
ブートストラップ法では,重複ありで同じ大きさのサンプルを抽出し,複製されたサンプルに対し推定を行います. この作業を大量に行うことにより,推定値自体の経験分布を得ます.

今回は,ブートストラップの反復回数は1000回で行いました.
ブートストラップ法にはbootパッケージという便利なものもあるのですが,今回は簡単な処理なのでforループで書いてしまいます.

result_boot <- tibble() # 空のdata.frameを用意しておく

for(i in 1:1000){
  # ブートストラップサンプルの抽出
  kidney_boot <- kidney %>% 
    sample_n(size = n(), # 同じ大きさ
             replace = T) # 重複あり

  # 抽出したブートストラップサンプルに対し曲線を推定
  fit_boot <- loess(tot ~ age, 
                    data = kidney_boot, 
                    span = 1/3,
                    degree = 1,
                    family = "symmetric")

  # 10, 20, ..., 80歳で当てはめた値を計算
  result_boot <- ages %>% 
    mutate(fit_loess_boot = predict(fit_boot, ages)) %>% 
    bind_rows(result_boot, .)
}

それでは,ブートストラップ標準誤差を計算し,結果に追加します

df <- result_boot %>% 
  group_by(age) %>% 
  summarise(se_loess = sd(fit_loess_boot, na.rm = T)) %>% 
  left_join(df,., by = "age")

これで念願の図を書くことができます!!

kidney %>% 
  ggplot(aes(age, tot))+
  geom_point(colour = "steelblue", shape = 8)+
  geom_line(data = kidney %>% 
              mutate(fit_loess = predict(fit_loess_kidney)),
            aes(y = fit_loess),
            colour = "red")+
  geom_segment(data = df,
               aes(x = age, xend = age,
                   y = fit_loess - 2*se_loess, yend = fit_loess + 2*se_loess))+
  labs(x = "年齢", y = "複合尺度tot")

f:id:kur0cky:20200809203756p:plain

10歳ごとにエラーバーを引かなければいけなかったので少々面倒くさい作業になってしまいました.
脳死な平滑化で良ければ以下のコードで一発だったのですけれどね.笑

kidney %>% 
  ggplot(aes(age, tot))+
  geom_point()+
  geom_smooth(method = "loess", span = 1/3)

f:id:kur0cky:20200809222520p:plain

tab 1.1

ここまできたら結果はdfにまとまっているので,サクッと表にしてしまいましょう.
そこまで込み入った表ではないのでgt は使わず,knitr::kable()で出力してしまいます.

df %>% 
  rename(線形回帰 = fit_lm,
         標準誤差 = se_lm,
         loess = fit_loess,
         ブートストラップ標準誤差 = se_loess) %>% 
  gather(key, value, -age) %>% 
  spread(age, value) %>% 
  rename(age = key) %>% 
  slice(c(3, 4, 1, 2)) %>% 
  mutate_if(is.numeric, round, 2) %>% 
  knitr::kable(format = "html")
age 20 30 40 50 60 70 80
線形回帰 1.29 0.50 -0.28 -1.07 -1.86 -2.64 -3.43
標準誤差 0.21 0.15 0.15 0.19 0.26 0.34 0.42
loess 1.65 0.65 -0.61 -1.24 -1.92 -2.68 -3.49
ブートストラップ標準誤差 0.67 0.24 0.36 0.35 0.37 0.49 0.67

fig 1.3

上で行ったブートストラップ推定のうち,はじめの25個を図示します.また,全データを用いたのものも合わせて示します.
飽きてきて若干コード汚くなってきてますがご容赦ください

result_boot <- tibble()
for(i in 1:25){
  kidney_boot <- kidney %>% 
    sample_n(size = n(),
             replace = T)
  fit_boot <- loess(tot ~ age,
                    data = kidney_boot,
                    span = 1/3,
                    degree = 1,
                    family = "symmetric")
  result_boot <- kidney_boot %>% 
    mutate(fit_loess = predict(fit_boot),
           iter = i) %>% 
    bind_rows(result_boot,.)
}

result_boot %>% 
  ggplot(aes(age, fit_loess, colour = factor(iter)))+
  geom_line(linetype = "dashed", size = .3)+
  geom_line(data = kidney %>% 
              mutate(fit_loess = predict(fit_loess_kidney)),
            aes(y = fit_loess),
            colour = "black", 
            size = 1)+
  theme(legend.position = "none")+
  labs(x = "年齢", y = "複合尺度tot")

f:id:kur0cky:20200809211721p:plain

1.2 仮説検定

2つ目の例では,72名の白血病患者について,7128個の遺伝子群について調べたデータを扱っています. データには2種類の白血病が混ざっており,47名はALL,25名はAMLの患者です. まずはデータを手に入れましょう.飽きてきたので整形もやってしまいます.

leukemia_big <- read.csv("http://web.stanford.edu/~hastie/CASI_files/DATA/leukemia_big.csv")

gene136AML <- leukemia_big %>% 
  select(starts_with("AML")) %>% 
  slice(136) %>% 
  unlist() %>% 
  unname()
gene136ALL <- leukemia_big %>% 
  select(starts_with("ALL")) %>% 
  slice(136) %>% 
  unlist() %>% 
  unname()

fig 1.4

7128個の遺伝子のうち,136番目のものの値をALL, AMLそれぞれでヒストグラムにします.

tibble(x = c(gene136AML, gene136ALL),
       type = rep(c("AML", "ALL"), 
                  times = c(length(gene136AML), length(gene136ALL)))) %>% 
  ggplot(aes(x))+
  geom_histogram(bins = 11, colour="black", fill = "steelblue")+
  facet_wrap(~type,
             ncol = 1)

f:id:kur0cky:20200809213716p:plain

2標本t検定

遺伝子136の平均は,ALLとAMLでそれぞれ0.75, 0.95です.
この2つの分布の平均は異なっていると言えるかどうか,t検定を行います.

> mean(gene136ALL)
[1] 0.7524794

> mean(gene136AML)
[1] 0.9499731

> t.test(gene136AML, gene136ALL, var.equal = T)

    Two Sample t-test

data:  gene136AML and gene136ALL
t = 3.014, df = 70, p-value = 0.003589
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.06680742 0.32817995
sample estimates:
mean of x mean of y 
0.9499731 0.7524794 

結果はp-value = 0.003589となり,本書と同様の値を出すことができました.

もし本当になにの違いも無ければ,非常にレアな事象が起こっている」ということになります

fig 1.5

最後に図1.5です.

全7128個の遺伝子に対して同様のt検定を行い,自由度70のt分布と比較してみます.
つまり7128回のt検定を行わなくてはならないのですが,ただ単にforループするだけではつまらないので,ここではtidyr::nest()purrr:map()を使い,オシャレに解析を管理する方法を共有したいと思います.

nest()した段階で一度出力を見てみます

> df <- leukemia_big %>% 
+   rownames_to_column("gene_index") %>% # 遺伝子番号を示す行名をカラムへ
+   gather(type, value, -gene_index) %>% # カラム名でALL, AMLを切り分けたいので縦長形式に
+   mutate(type = str_sub(type, 1, 3)) %>% # stringrパッケージを活用し,病種を抽出
+   group_by(gene_index, type) %>% 
+   nest()
> print(df)

# A tibble: 14,256 x 3
# Groups:   gene_index, type [14,256]
   gene_index type  data             
   <chr>      <chr> <list>           
 1 1          ALL   <tibble [47 × 1]>
 2 2          ALL   <tibble [47 × 1]>
 3 3          ALL   <tibble [47 × 1]>
 4 4          ALL   <tibble [47 × 1]>
 5 5          ALL   <tibble [47 × 1]>
 6 6          ALL   <tibble [47 × 1]>
 7 7          ALL   <tibble [47 × 1]>
 8 8          ALL   <tibble [47 × 1]>
 9 9          ALL   <tibble [47 × 1]>
10 10         ALL   <tibble [47 × 1]>
# … with 14,246 more rows

nest()で生まれた新たなdataカラムの各要素には,47 × 1のtibbleが格納されていることがわかります!
このように,tibbleはリストを持つことができるのです.

つまりあとは,spread()で横長形式に戻し,t検定を各レコードに対して実行してやれば良いことになるのです.

レコードごとに処理を行いたいので,purrr::map()を活用します.map2()は,2つの同じ長さのリストを利用した処理を行うための関数です.

> df <- df %>% 
+   spread(type, data) %>% 
+   mutate(t_test = map2(AML, ALL, ~ t.test(unlist(.x), unlist(.y))))
> print(df)

# A tibble: 7,128 x 4
# Groups:   gene_index [7,128]
   gene_index ALL               AML               t_test 
   <chr>      <list>            <list>            <list> 
 1 1          <tibble [47 × 1]> <tibble [25 × 1]> <htest>
 2 10         <tibble [47 × 1]> <tibble [25 × 1]> <htest>
 3 100        <tibble [47 × 1]> <tibble [25 × 1]> <htest>
 4 1000       <tibble [47 × 1]> <tibble [25 × 1]> <htest>
 5 1001       <tibble [47 × 1]> <tibble [25 × 1]> <htest>
 6 1002       <tibble [47 × 1]> <tibble [25 × 1]> <htest>
 7 1003       <tibble [47 × 1]> <tibble [25 × 1]> <htest>
 8 1004       <tibble [47 × 1]> <tibble [25 × 1]> <htest>
 9 1005       <tibble [47 × 1]> <tibble [25 × 1]> <htest>
10 1006       <tibble [47 × 1]> <tibble [25 × 1]> <htest>
# … with 7,118 more rows

このような形にしたことで,それぞれ遺伝子に対応する結果を同じレコードで管理できます.
tidyverse,おそろしいですねぇ

それでは最後に,作図して終えたいと思います.

dens = tibble(x = seq(-10, 10, length.out = 1000)) %>% 
  mutate(dens = dt(x, df = 70))
df %>% 
  mutate(p = map_dbl(t_test, ~ .x$statistic)) %>% 
  ggplot(aes(p))+
  geom_histogram(bins = 70, colour = "black", fill = "green")+
  geom_line(aes(x = x, y = dens*3500), data = dens)+
  scale_y_continuous(limits = c(0, 700))+
  labs(x = "t統計量", y = "度数")

f:id:kur0cky:20200806000109p:plain

おわりに

少ないはずの1章から,結構な量になってしまいました...若干飽き始めてますが,最後までやりきりたいと思います.

本書については私はまだ読み途中ですが,歴史的敬意や,主義的な対立を踏まえた俯瞰など,「お気持ち」が込められておりニヤニヤしながら読んでいます.

特に,ブートストラップ法への愛がにじみ出ている気がします.