library(tidyverse)

1. z 検定と t 検定の関係

  • z 検定は正規分布を用いる統計学的検定法
  • 「標本の平均」と「母集団の平均」が統計学的にみて有意に異なるかどうかを検定する方法
  • z 検定を用いるために最も重要な前提条件
  • 母集団の平均と標準偏差が既知であること（ほとんどの場合、あり得ない前提）
    →　母集団の正しい標準偏差 \(σ\) を知るというのは一般には現実的でない
  • 標本は母集団から抽出された単純ランダム標本でなければならないこと
  • 母集団は正規分布に従うこと
  • z 検定は、母集団が完全に既知の場合にのみ使われる
    →　z 検定では前提条件が現実的ではないことが多い
    →　z 検定の代わりにスチューデントの t 検定を用いる
  • ここでは t 検定を学ぶ前提として、標準正規分布、ｚ 値、z 分布について学ぶ

2. 「正規分布」と「標準正規分布」

2.1 データの生成

  • ここでは、10,000 人の架空の学生の試験結果 score を作ってみる
  • 母集団の標準偏差を 10、母集団の平均が 50点の「正規分布」
  • R のrnorm() 関数を使って、この母集団から分析で使うデータを生成してみる
set.seed(1982)              # 同一の結果を得るためにシードを設定
score <- rnorm(10000,       # 10000人分の架空データを作成する  
               mean = 50,   # 母集団の平均値 = 50  
               sd = 10) |>  # 母集団の標準偏差 = 10
  round(digits = 0)         # 点数は小数点以下を切り捨て
st_id = seq(1:10000)　　　  　# 学生の ID の作成
df <- tibble(score, st_id)  # 作成した二つの変数を df に統合
DT::datatable(df)
  • df の記述統計を表示
sd(df$score)
[1] 9.975603
  • 標準偏差はほぼ 10
summary(df$score)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     16      43      50      50      57      87
  • 試験点数の範囲は 16点から87点
  • 平均は 50点

2.2 基準化 (Standardization)

  • この試験点数 score を基準化 (Standardization) する
  • 基準化 (Standardization) とは、試験点数 score を z 値に変換すること

「標準」正規分布とは

・「正規分布」に変数変換を施した（＝標準化した）後の z 分布のこと
→　「標準」正規分布
・標準化した z 分布（=「標準」正規分布）の特徴
1. 平均値は 0
2. 標準偏差は 1
3. 範囲は −4 から 4

「正規分布 (score)」を「標準正規分布 (z 値)」に変換する ・score の平均値は 50 　　　　　 => z 値の平均値は 0
・score の標準偏差は 10 　　　　=> z 値の標準偏差は 1
・score の範囲は16点から87点 　 => z 値の範囲は −4 から 4

基準化 (Standardization) するための式 \[z = \frac{個々のデータ−平均}{標準偏差}\]

  • 例えば、試験点数 score の 30点を z 値に変換してみる
(30-50)/10
[1] -2
  • 試験点数 score の 30点は z 値の -2 に相当することがわかる

  • 試験点数の平均 50 点に相当する z 値は 0 になる
  • 基準化された後の試験点数 score は標準正規分布に従う

2.3 標準正規分布表の読み方

標準正規分布表


2.3.1 「z = 1 標準偏差 (1\(\sigma\))」の場合

  • 「z = 1.00」と「0.00」が交差する値「0.3413」の意味
  • この数値は「z の値が 0 と 1.0 の間の面積は 0.3413」という言う意味
  • ベルカーブの内側の全体の面積は 1
  • そのうちの 0.3413 が z = 0 と z = 1 の間（つまり、1 標準偏差 = \(1\sigma\)）の面積
  • 確率は面積で表せるので、100% のうち 34.13% が z = 0 と z = 1 の間の確率

標準正規分布と標準偏差 (\(\sigma\))の関係: 1 標準偏差の場合 ・確率変数 \(X\) が \(N(μ, σ^2)\) に従う時
→　つまり「平均値が \(μ\)、標準偏差が \(\sigma\) 」に従う時
・平均 \(μ\) からのずれが \(±1σ\)（± 1 標準偏差）の範囲に \(X\) が含まれる確率：68.26%

  • 平均 \(μ\) からのずれが \(0 と 1σ\) の範囲に \(X\) が含まれる確率は 34.13%
    → 平均 \(μ\) からのずれが \(-1σ と 1σ\) の範囲（±1 標準偏差）に \(X\) が含まれる確率は 34.13 × 2 = 68.26%

2.3.2 「z = 2 標準偏差 (2\(\sigma\))」の場合

  • 「z = 2.00」と「0.00」が交差する値「0.4772」の意味
  • この数値は「z の値が 0 と 2.0 の間の面積は 0.4772」という言う意味
  • ベルカーブの内側の全体の面積は 1
  • そのうちの 0.4772 が z = 0 と z = 2 の間（つまり、1 標準偏差 = \(2\sigma\)）の面積
  • 確率は面積で表せるので、100% のうち 47.72% が z = 0 と z = 2 の間の確率

標準正規分布と標準偏差 (\(\sigma\))の関係: 2 標準偏差の場合 ・確率変数 \(X\) が \(N(μ, σ^2)\) に従う時
→　つまり「平均値が \(μ\)、標準偏差が \(\sigma\) 」に従う時
・平均 \(μ\) からのずれが \(±2σ\)（± 2 標準偏差）の範囲に \(X\) が含まれる確率：95.44%

  • 平均 \(μ\) からのずれが \(0 と 2σ\) の範囲に \(X\) が含まれる確率は 47.72%
    → 平均 \(μ\) からのずれが \(-2σ と 2σ\) の範囲（±2 標準偏差）に \(X\) が含まれる確率は 47.74 × 2 = 95.44%

2.3.3 「z = 3 標準偏差 (3\(\sigma\))」の場合

  • 「z = 3.00」と「0.00」が交差する値「0.4987」の意味
  • この数値は「z の値が 0 と 3.0 の間の面積は 0.4987」という言う意味
  • ベルカーブの内側の全体の面積は 1
  • そのうちの 0.4987が z = 0 と z = 3 の間（つまり、3 標準偏差 = \(3\sigma\)）の面積
  • 確率は面積で表せるので、100% のうち 49.87% が z = 0 と z = 3 の間の確率

標準正規分布と標準偏差 (\(\sigma\))の関係: 3 標準偏差の場合 ・確率変数 \(X\) が \(N(μ, σ^2)\) に従う時
→　つまり「平均値が \(μ\)、標準偏差が \(\sigma\) 」に従う時
・平均 \(μ\) からのずれが \(±3σ\)（± 3 標準偏差）の範囲に \(X\) が含まれる確率：99.74%

  • 平均 \(μ\) からのずれが \(0 と 3σ\) の範囲に \(X\) が含まれる確率は 49.87%
    → 平均 \(μ\) からのずれが \(-3σ と 3σ\) の範囲（±3 標準偏差）に \(X\) が含まれる確率は 49.87 × 2 = 99.74%

2.3.4 標準正規分布と標準偏差の関係（まとめ）

標準正規分布と標準偏差 (\(\sigma\))の関係 ・確率変数 \(X\) が \(N(μ, σ^2)\) に従う時
→　つまり「平均値が \(μ\)、標準偏差が \(\sigma\) 」に従う時
・平均 \(μ\) からのずれが \(±1σ\)（± 1 標準偏差）の範囲に \(X\) が含まれる確率：68.26%
・平均 \(μ\) からのずれが \(±2σ\)（± 2 標準偏差）の範囲に \(X\) が含まれる確率：95.44%
・平均 \(μ\) からのずれが \(±3σ\)（± 3 標準偏差）の範囲に \(X\) が含まれる確率：99.74%

2.4 偏差値

  • 日本の入試制度で使われている「偏差値」は z 分布 の平均値と標準偏差を次のように変更したもの
平均値 標準偏差
z分布 0 1
↓ ↓
偏差値 50 10

基準化 (Standardization) するための式 \[z = \frac{個々のデータ−平均}{標準偏差}\]

偏差値を求める式 \[偏差値 = 10*\frac{個々のデータ−平均}{標準偏差}+ 50\]

偏差値が 60 なら
  • 「平均値（50点）± 1 標準偏差」の間に、全体の 68% が含まれる
    →　トップ (100 - 68.26) / 2 = 32/2 = 約16% に含まれることになる
偏差値が 70 なら
  • 「平均値（50点）± 2 標準偏差」の間に、全体の 95% が含まれる
    →　トップ (100 - 95.44) / 2 = 4.56/2 = 約2.28% に含まれることになる
偏差値が 80 なら
  • 「平均値（50点）± 3 標準偏差」の間に、全体の 99.7% が含まれる
    →　トップ (100 - 99.74) / 2 = 0.26/2 = 約0.13% に含まれることになる

3. t 検定

  • 前セッションで紹介した中心極限定理によれば、母集団が正規分布していなくても、平均 \(\mu\)，分散 \(σ^2\) の母集団からの無作為標本における標準誤差の分布は「標本サイズが十分に大きければ」近似的に正規分布に従う 
  • しかし、標本サイズが十分に大きいケースばかりとは限らず、少ない標本サイズしか入手できないことが多々ある 
  • この問題を解決したのが t 検定 を考え出した William Sealy Gosset という元ギネスビール社員 
  • student's t という名前で論文を投稿したため、t 検定と呼ばれる 
  • Gosset は「標本サイズが小さくても」母集団からの無作為標本における標準誤差の分布は \((n-1)\) の自由度で \(t\) 分布することを発見した 
  • この発見のおかげで、小さい標本サイズの場合でも \(t\) 分布の特徴を使うことで統計的推計が可能となった 
  • 標本数が 100 を超える場合、標準正規分布の特徴を使った \(z\) 検定が使われ、標本数が少ないときに限り\(t\) 検定が使われていた 
  • その理由は、\(t\) 検定の統計量である\(t\) を計算するためには、紙上で大規模な\(t\) 分布表を作成する必要があり、その作業が大変だったから 
  • しかし、R や Stata など統計ソフトの開発により、の\(t\) 分布表から容易に\(t\) 値を計算できるようになったため、\(z\) 検定は使われなくなった 
  • 下の \(t\) 分布図が示すように、自由度（degree of freedom: 標本サイズから 1 を引いた値）が大きくなるにつれ \(t\) 検定の結果と \(z\) 検定の結果は近似していき、標本サイズが無限大だと、両分布は同一になる 

3.1 「1 グループの t 検定」と有意水準

t 検定の手順

  • 母平均 μ = 5.5` の正規母集団から、10 個のデータ1, 2, 3, 4, 5, 6, 7, 8, 9, 10 を標本として抽出したとする 
  • ここで、標本サイズ 10 から得られた標本平均は 5.5 であるが、母平均は 5 なのかどうか確かめたいとする 
    これを確かめるために \(t\) 値が必要となる 
t 検定の手順
  • 帰無仮説の明示 (null hypothesis): \(H_0\)
  • 対抗仮説の明示 (alternative hypothesis): \(H_1\)
  • t 値を計算する
  • 帰無仮説を棄却する臨界値 (critical values) を特定する
  • t 値が棄却域内かどうかを確かめる
  • 結論

    • 具体的な検定プロセス

    • まず帰無仮説を設定する

    • 帰無仮説とは「棄却されるために設定する仮説」

    • 知りたいのは「サンプル平均は 5.5 だが、母集団の平均は 5 なのか」ということ

    • 帰無仮説 \(H_0\)は「母平均 = 5」

    • 次に対抗仮説を設定する

    • 対抗仮説は「母平均は 5 ではない」

    • 帰無仮説と対抗仮説は相互に排他的 (mutually exclusive)

    • 次に \(t\) 値を計算する

    • \(t\) 値は、次の式で求めることができる 

    \[T = \frac{\bar{x} - μ_0}{SE} = \frac{\bar{x} - μ_0}{u_x / \sqrt{n}}\]

    • \(\bar{x}\) : 標本平均（ここでは 5.5）

    • \(μ_0\) : 母集団で推定したい値（ここでは 5）

    • \(n\) : 標本サイズ（ここでは 10)

    • \(u_x\): 不偏標準偏差（= 標本標準偏差）

    • \(SE\) : 標準誤差 (standard Error: SE)

    • 不偏標準偏差\(u_x\) は不偏分散 (unbiased variance)の平方根であるから、まず、不偏分散を求める 

    \(x\) の不偏分散 \(u_x^2\) は次の式で計算できる:

    \[u_x^2 = \frac{\sum_{i=1}^n (x_i - \bar{x})^2}{n-1}\]

    • これより、\(x\) の不偏標準偏差 (unbiased standard diviation) \(u_x^2\) 3.03 が得られる 

    ・母集団の平均は 5 なのかということを確かめたい 
    → 母集団で推定したい値 \(μ\)に 5 を代入して、次の \(t\) 値を得る 

    \[T = \frac{\bar{x} - μ_0}{u_x / \sqrt{n}}\]

    \[ = \frac{{5.5} - 5}{3.03 / \sqrt{10}}\]

    \[ = 0.522\]

    • 点推定（区間推定に関しては「母平均の推定と信頼区間」を参照）
    • ここで得られた 0.522 という\(t\) 値を使って 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 という標本サイズ 10 のサンプルから、母集団の母平均が 5 であるという仮説が妥当かどうかを推定する 

    \[ t 分布表\]

    • \(t\) 分布表内の数値は、検定を行う際の有意水準 (significance leve) に対応した棄却限界値を示している 

    • 標本から得られた \(t\) 値の絶対値が棄却限界値より大きい場合、帰無仮説を棄却する（両側検定の場合） 

    • 分布表の横軸は片側検定における有意水準を示している 

    • 例えば、二列目の「確率95%」は両側検定における有意水準が 5% \((α = 0.05)\) であることを示す 

    • \(t\) 検定では両側検定で 5% \((α = 0.05)\) の有意水準を使うのが一般的 

    • 従って、 \(t\) 分布表二列目の「確率95%」を使う 

    • 分布表の一列目は自由度 (df: degree of freedom) を示している 

    • ここでは「標本サイズ 10 から 1 を引いた 9」 が自由度 

    • 縦軸の 9 と横軸の「確率95%」交差する数字「2.262」が 5% の両側検定における棄却限界値 

    • これを図示すると次のようになる 
      次の図は、自由度 9 の \(t\) 分布における 5% 有意水準\(（α = 0.05）\)の 2 つの棄却限界値 (critical value) -2.26 と 2.26 を示している 

    • サンプルから得られた\(t\) 値が -2.26 と 2.26 の区間内にあれば、帰無仮説 \(H_0\) は「母平均 = 5」は棄却できない 
    • t 値がこの区間外（赤の斜線部分）にあれば、帰無仮説を棄却する 
    • \(t\) = 0.522 なので、帰無仮説は棄却できない 
    • 従って、得られたサンプルから判断する限り「母平均 = 5」である可能性を排除できない 
    • つまり、得られたサンプルから判断する限り「母平均は 5 でありうる」 

    3.2 有意水準と有意確率

    有意水準 (significance level)

    • 「有意水準 5 %」が意味するのは、帰無仮説が正しい時、誤って帰無仮説をたまたま棄却してしまう危険性が 5%（100回のうち 5回）存在すること 
    • 従って、有意水準は「危険率」（\(α\): アルファ）とも呼ばれる 
    • この過ちは「第一種の過誤」(type I error) と呼ばれる 
    • 正しい帰無仮説をたまたま誤って棄却してしまう危険が 5% を下回った場合、帰無仮説を棄却（つまり仮説は無に帰した）し、対抗仮説を受け入れる 
    • 有意水準 = 有意確率と記述している統計学の本やサイトもあるので注意が必要 

    有意確率 (significance probability) = \(p\) 値

    • 「帰無仮説が正しいとき、検定統計量が標本から得られた\(t\) 値以上に分布の中心からかけ離れた値をとる確率」
    • \(p\) 値は帰無仮説の下での\(t\) 値の異常性を表しており、 \(p\) 値が小さいほど標本から得られた\(t\) 値が異常ということ 
    • 従って\(t\) 値の異常性が十分大きい場合（つまり \(p\) 値が十分小さい時に）帰無仮説を棄却する 
    • \(p\) 値は「帰無仮説が正しい確率」でもないし「対抗仮説が間違っている確率」でもない 
    • 「帰無仮説が正しい確率」は 0 もしくは 1 のいずれかであり、 \(p\) 値のように 0 と 1 の間の値をとることはありえない 
    • 真実は帰無仮説が正しいか間違っているかのいずれかであり、帰無仮説が正しいなら「帰無仮説が正しい確率」は 1 であり、そうでないなら「帰無仮説が正しい確率」は 0 である 
    • 我々は真実（つまり母数）を知らないので、この確率が 0 か 1 かを知ることは不可能 

    3.3 R による t 検定

    # 標本サイズ 10 のサンプルに score という名前をつける
score <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)

    母平均 = 5 を検定する (信頼区間 = 95% → 5%で統計的に有意)

    # 帰無仮説（母平均 = 5) を `t 検定`する
t.test(score, mu = 5)
    
    One Sample t-test

data:  score
t = 0.52223, df = 9, p-value = 0.6141
alternative hypothesis: true mean is not equal to 5
95 percent confidence interval:
 3.334149 7.665851
sample estimates:
mean of x 
      5.5
    • 3 行目に \(t\) = 0.52223、p-value が 0.6141 という結果が得られた 
    • p-value = 0.6141 は、両側検定の p 値を表している 
    • 4 行目に alternative hypothesis: true mean is not equal to 5 と対抗仮説が明示されている 
    • p-value (0.6141) が 0.05 より大きいので、帰無仮説（「母平均 = 5」）は棄却できない 
    • もし、p-valueが 0.05 よりも小さければ 5% の有意水準（α = 0.05）で帰無仮説を棄却できる 
    • 従って、得られたサンプルから判断する限り「母平均 = 5」である可能性を排除できない 
    • つまり、得られたサンプルから判断する限り「母平均は 5 でありうる」 

    母平均 = 3 を検定する (信頼区間 = 95% → 5%で統計的に有意)

    # 帰無仮説（母平均 = 3) を `t 検定`する
t.test(score, mu = 3)
    
    One Sample t-test

data:  score
t = 2.6112, df = 9, p-value = 0.02822
alternative hypothesis: true mean is not equal to 3
95 percent confidence interval:
 3.334149 7.665851
sample estimates:
mean of x 
      5.5
    • 3 行目に \(t\) = 2.6112、p-value が 0.02822 という結果が得られた 
    • p-value = 0.02822 は、両側検定の p 値を表している 
    • 4 行目に alternative hypothesis: true mean is not equal to 3 と対抗仮説が明示されている 
    • p-value (0.02822) が 0.05 より小さいので、帰無仮説（「母平均 = 3」）は棄却できる 
    • 従って、得られたサンプルから判断する限り「母平均 = 3」ではない 
    • つまり、得られたサンプルから判断する限り「母平均は 3 でありえず、それよりも大きい）」 

    母平均 = 3 を検定する (信頼区間 = 99% → 1%で統計的に有意)

    • 信頼区間を変更したい場合は conf.level を指定する
    # 帰無仮説（母平均 = 3) を `t 検定`する
t.test(score, mu = 3, conf.level = 0.99)
    
    One Sample t-test

data:  score
t = 2.6112, df = 9, p-value = 0.02822
alternative hypothesis: true mean is not equal to 3
99 percent confidence interval:
 2.388519 8.611481
sample estimates:
mean of x 
      5.5
    • 3 行目に \(t\) = 2.6112、p-value が 0.02822 という結果が得られた 
    • p-value = 0.02822 は、両側検定の p 値を表している 
    • 4 行目に alternative hypothesis: true mean is not equal to 3 と対抗仮説が明示されている 
    • p-value (0.02822) が 0.01 より小さいので、帰無仮説（「母平均 = 3」）は棄却できる 
    • 従って、得られたサンプルから判断する限り「母平均 = 3」ではない 
    • つまり、得られたサンプルから判断する限り「母平均は 3 でありえず、それよりも大きい）」 

    4. 割合に関する推定

    内閣支持29%、不支持52%（NHK世論調査）記事
    NHKの世論調査によりますと、菅内閣を「支持する」と答えた人は、先月より4ポイント下がって29%で、去年9月の内閣発足以降最低を更新しました。・・・調査の対象となったのは、2115人で、57%にあたる1214人から回答を得ました。2021年8月（8月10日更新）

    Q:この世論調査から「菅内閣の支持率は30%を切っている」といえるのか？

    • 割合に関する推定を実行してみる
    • 「支持する」と答えたのは、回答した人数である1214人の 29%なので 352 人
    prop.test(c(352), c(1214))
    
    1-sample proportions test with continuity correction

data:  c(352) out of c(1214), null probability 0.5
X-squared = 213.41, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.2647212 0.3165265
sample estimates:
        p 
0.2899506

    • 内閣支持率の 95%信頼区間は 26.47% ~ 31.66%

    • この調査で内閣支持率 29%という推定値を得たからといって「有権者全体における内閣支持率は 30%を切っている」とは 95% の確率ではいえない

    • 仮に、菅内閣を「支持する」と答えたのが 352人ではなく300人だったとすると

    prop.test(c(300), c(1214))
    
    1-sample proportions test with continuity correction

data:  c(300) out of c(1214), null probability 0.5
X-squared = 309.53, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.2232793 0.2725771
sample estimates:
       p 
0.247117
    • 内閣支持率の 95%信頼区間は 22.33% ~ 27.26%
    • 「有権者全体における内閣支持率は 30%を切っている」と 95% の確率でいえる

    5. Excercise

    Q5.1: 「2. z 分布」を参考にして、次の問題を解きなさい　　

    • B 県の高校 1 年生全員がある予備校の数学のテストを受けました
    • 採点したところ 「数学のテスト結果」は平均が 45 点で標準偏差が 10 の正規分布に従うとみなせることがわかりました
      ・63 点以上得点した受験生の割合(⇒確率)をこたえなさい

    Q5.2: 「2. z 分布」を参考にして、次の問題を解きなさい　　

    • 高校 1 年生全員がある予備校の数学のテストを受けました　　
    • 採点したところ「数学のテス ト結果」は平均が 45 点で標準偏差が 10 の正規分布に従うとみなせることがわかりました
      ・Q5.2.1: 70 点以上得点した受験生の割合は?
      ・Q5.2-2: 50 点以上得点した受験生の割合は?
      ・Q5.2-3: 40 点以下得点した受験生の割合は?
      ・Q5.2-4: トップ 5% に入るためには何点必要か?

    Q5.3: 「3. t 検定」を参考にして、次の問題を解きなさい

    • 早稲田大学の A 君は、大隈商店街で喫茶店「ルオー」を開いた
    • A 君は「ルオー」の売上額を正規母集団から観測されるデータをみなし、その母平均 \(\mu\) を代表的な売り上げとして推定しようとした
    • 毎日の売り上げ総額を示す伝票の中からランダムに8 枚を抜き出してみると次のような数字が出てきた

    \[45, 39, 42, 57, 28, 33, 40, 52（単位は万円）\]
    ・母集団では一日あたりの「ルオー」の売り上げが 50 万円である、という仮説を検定しなさい

    Q5.4: 「3. t 検定」を参考にして、次の問題を解きなさい

    • 母平均 \(μ\) = 5.5 の正規母集団から、10 個のデータ1, 2, 3, 4, 5, 6, 7, 8, 9, 10 を標本として抽出したとする
      ・ここで、標本サイズ 10 から得られた標本平均は 5.5 であるが、母平均は 3 なのかどうか R を使って確かめなさい

    Q5.5: 「2. z 分布」を参考にして、次の問題を解きなさい

    • 早稲田大学政治経済学部を 20,000 人が受験したところ、その成績は、平均値 65 点、標準偏差 10 点の正規分布に従った

    ・Q5.5.1: ある受験生が，75 点以上 85 点以下である確率を求めなさい
    ・Q5.5.2: この入学試験において，上位 10 %に入るためには，何点以上あればよいか?
    ・Q5.5.3: この入学試験において，上位 1000 人が合格する。合格するためには何点以上必要か?

    参考文献
  • 宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」
  • 土井翔平（北海道大学公共政策大学院）「Rで計量政治学入門」
  • 矢内勇生（高知工科大学）授業一覧
  • 浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年
  • 浅野正彦, 中村公亮.『初めてのRStudio』オーム社、2018年
  • Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.
  • Kieran Healy, DATA VISUALIZATION, Princeton, 2019
  • Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017