進化学実習 2026 牧野研 東北大学

1. 導入: データ解析の全体像。Rの基本。
2. データの可視化。
3. データ構造の処理1: 抽出、集約など。
4. データ構造の処理2: 結合、変形など。
5. データ内容の処理: 数値、文字列など。
6. データ入力、レポート作成
7. 統計モデリング1: 確率分布、尤度
8. 統計モデリング2: 一般化線形モデル
9. 発表会

作業再開、出欠確認

1. r-training-2026.Rproj をダブルクリック。
   正しい working directory でRStudioが起動される。
2. 初日と同じ手順で出欠確認。
3. 今日のぶんのスクリプトを新規作成し、好きな名前で保存。
4. tidyverseの読み込みやパレット設定など、 おまじないをまず実行。

この実習の目標

✅ 生物学研究にはデータとモデルが必須だと認識

✅ 再現可能な解析を楽にやりたい気持ちになる

✅ 必要な方法を調べ、実践する力をつける

⬜ データ解析の基本に触れる

個々の方法は覚えなくても大丈夫！
忘れては調べ、を何度も繰り返しながら染み込ませていこう。

データ解析のおおまかな流れ

1. コンピュータ環境の整備 ✅
2. データの取得、読み込み ✅ readr簡単
3. 探索的データ解析
   • 前処理、加工 ✅ dplyr, tidyr, stringr手放せない
   • 可視化、仮説生成 ✅ ggplot楽しい
   • 統計解析、仮説検証 ⬜ 👈今回
4. 報告、発表 ✅ Quarto楽しい

データを使ってやりたいこと

• 現象を理解したい
• 将来を予測したい
• ものを分類・判別したい
• 挙動を制御したい
• 新しい何かを生成したい

そのために解析は必要？ 未加工の生データこそ宝？

初日を振り返る

データ科学における数理モデル

データ生成をうまく真似できそうな仮定の数式表現。
e.g., 大きいほど高く売れる: $\text{price} = A \times \text{carat} + B + \epsilon$

新しく採れたダイヤモンドの価格予想とかにも使える。

このように「YをXの関数として表す」ようなモデルを回帰と呼ぶ。

よりよい回帰をめざして

久保先生の緑本に沿ってちょっとずつ線形モデルを発展させていく。

「データ解析のための統計モデリング入門」久保拓弥 2012 より改変

回帰モデルの2段階

1. Define a family of models: だいたいどんな形か、式をたてる

   • 直線: $y = a_1 + a_2 x$
   • 対数: $\log(y) = a_1 + a_2 x$
   • 二次曲線: $y = a_1 + a_2 x^2$

2. Generate a fitted model: データに合うようにパラメータを調整

   • $y = 3x + 7$
   • $y = 9x^2$

https://r4ds.had.co.nz/model-basics.html

たぶん身長が高いほど体重も重い

なんとなく $y = a x + b$ でいい線が引けそう
 

たぶん身長が高いほど体重も重い

なんとなく $y = a x + b$ でいい線が引けそう
じゃあ傾き a と切片 b、どう決める？

最小二乗法 (Ordinary Least Square: OLS)

回帰直線からの残差平方和(RSS)を最小化する。

残差平方和(RSS)が最小となるパラメータを探せ

ランダムに試してみて、上位のものを採用。
この程度の試行回数では足りなそう。

残差平方和(RSS)が最小となるパラメータを探せ

グリッドサーチ: パラメータ空間の一定範囲内を均等に試す。
さっきのランダムよりはちょっとマシか。

こうした最適化の手法はいろいろあるけど、ここでは扱わない。

これくらいなら一瞬で計算してもらえる

par_init = c(intercept = 0, slope = 0)
result = optim(par_init, fn = rss_weight, data = df_weight)
result$par

intercept     slope 
-69.68394  78.53490 

上記コードは最適化一般の書き方。覚えなくていい。
回帰が目的なら次ページのようにするのが楽 →

lm() で直線あてはめしてみる

fit = lm(data = mpg, formula = hwy ~ displ)
broom::tidy(fit)

         term  estimate std.error statistic       p.value
1 (Intercept) 35.697651 0.7203676  49.55477 2.123519e-125
2       displ -3.530589 0.1945137 -18.15085  2.038974e-46

mpg_aug = broom::augment(fit, type.predict = "response")
ggplot(mpg_aug) + aes(displ, hwy) + geom_point() +
  geom_line(aes(y = .fitted), linewidth = 1, color = "#3366ff")

🔰 diamonds と iris でも lm() を試してみよう。

何でもかんでも直線あてはめではよろしくない

• 観察データは常に正の値なのに予測が負に突入してない？
• 縦軸は整数。しかものばらつきが横軸に応じて変化？

何でもかんでも直線あてはめではよろしくない

• 観察データは常に正の値なのに予測が負に突入してない？
• 縦軸は整数。しかものばらつきが横軸に応じて変化？
• データに合わせた統計モデルを使うとマシ

よりよい回帰をめざして

久保先生の緑本に沿ってちょっとずつ線形モデルを発展させていく。

確率分布

発生する事象(値)と頻度の関係。

手元のデータを数えて作るのが経験分布
e.g., サイコロを12回投げた結果、学生1000人の身長

一方、少数のパラメータと数式で作るのが理論分布。
(こちらを単に「確率分布」と呼ぶことが多い印象）

分布を特徴づける代表値 central tendency

平均値 mean

和を観察数で割る

中央値 median

順に並べて真ん中

最頻値 mode

最も頻度が高い値

外れ値に対する応答

もし総資産額20兆円の大富豪が鳥取県に引っ越してきたら
→ 県民の平均資産は4000万円上昇。中央値・最頻値はほぼそのまま。

目的や状況に応じて使い分けよう。

ばらつきを捉える記述統計量

分散 variance

平均値からの差の自乗の平均。 $\frac 1 n \sum _i ^n (X_i - \bar X)^2$

これの平方根が標準偏差 (standard deviation)。

Percentile, Quantile (四分位)

小さい順にならべて上位何%にあるか。

中央値 = 50th percentile = 第二四分位(Q2)

記述統計量に頼りすぎず分布を可視化する

同じデータでも見せ方で印象・情報量が変わる。

確率変数$X$はパラメータ$\theta$の確率分布$f$に従う…?

$X \sim f(\theta)$

e.g.,
コインを3枚投げたうち表の出る枚数 $X$ は二項分布に従う。
$X \sim \text{Binomial}(n = 3, p = 0.5)$

一緒に実験してみよう。（日本の硬貨は年号ありが裏）

試行を繰り返して記録してみる

コインを3枚投げたうち表の出た枚数 $X$

試行1: 表 裏 表 → $X = 2$
試行2: 裏 裏 裏 → $X = 0$
試行3: 表 裏 裏 → $X = 1$ 続けて $2, 1, 3, 0, 2, \ldots$

コイントスしなくても $X$ らしきものを生成できる

• コインを3枚投げたうち表の出る枚数 $X$
• $n = 3, p = 0.5$ の二項分布からサンプルする乱数 $X$

$X \sim \text{Binomial}(n = 3, p = 0.5)$

   ↓ サンプル

{2, 0, 1, 2, 1, 3, 0, 2, …}

これらはとてもよく似ているので
「コインをn枚投げたうち表の出る枚数は二項分布に従う」
みたいな言い方をする。逆に言うと
「二項分布とはn回試行のうちの成功回数を確率変数とする分布」
のように理解できる。

統計モデリングの一環とも捉えられる

コイン3枚投げを繰り返して得たデータ {2, 0, 1, 2, 1, 3, 0, 2, …}

   ↓ たった2つのパラメータで記述。情報を圧縮。

$n = 3, p = 0.5$ の二項分布で説明・再現できるぞ

こういうふうに現象と対応した確率分布、ほかにもある？

有名な確率分布、それに「従う」もの

離散一様分布

コインの表裏、サイコロの出目1–6

負の二項分布 (幾何分布 if n = 1)

成功率pの試行がn回成功するまでの失敗回数

二項分布

成功率p、試行回数nのうちの成功回数

ポアソン分布

単位時間あたり平均$\lambda$回起こる事象の発生回数

ガンマ分布 (指数分布 if k = 1)

ポアソン過程でk回起こるまでの待ち時間

正規分布

確率変数の和、平均値など。

離散一様分布

同じ確率で起こるn通りの事象のうちXが起こる確率

e.g., コインの表裏、サイコロの出目1–6

🔰 一様分布になりそうな例を考えてみよう

幾何分布 $~\text{Geom}(p)$

成功率pの試行が初めて成功するまでの失敗回数

e.g., コイントスで表が出るまでに何回裏が出るか

\[ \Pr(X = k \mid p) = p (1 - p)^k \]

「初めて成功するまでの試行回数」とする定義もある。

🔰 幾何分布になりそうな例を考えてみよう

負の二項分布 $~\text{NB}(n, p)$

成功率pの試行がn回成功するまでの失敗回数X。 n = 1 のとき幾何分布と一致。

\[ \Pr(X = k \mid n,~p) = \binom {n + k - 1} k p^n (1 - p)^k \]

失敗回数ではなく試行回数を変数とする定義もある。連続である必要はない。

🔰 負の二項分布になりそうな例を考えてみよう

二項分布 $~\text{Binomial}(n,~p)$

確率$p$で当たるクジを$n$回引いてX回当たる確率。平均は$np$。

\[ \Pr(X = k \mid n,~p) = \binom n k p^k (1 - p)^{n - k} \]

🔰 二項分布になりそうな例を考えてみよう

ポアソン分布 $~\text{Poisson}(\lambda)$

平均$\lambda$で単位時間(空間)あたりに発生する事象の回数。

e.g., 1時間あたりのメッセージ受信件数、メッシュ区画内の生物個体数

\[ \Pr(X = k \mid \lambda) = \frac {\lambda^k e^{-\lambda}} {k!} \]

二項分布の極限 $(\lambda = np;~n \to \infty;~p \to 0)$。
めったに起きないことを何回も試行するような感じ。

指数分布 $~\text{Exp}(\lambda)$

ポアソン過程の事象の発生間隔。平均は $1 / \lambda$ 。

e.g., メッセージの受信間隔、道路沿いに落ちてる手袋の間隔

\[ \Pr(x \mid \lambda) = \lambda e^{-\lambda x} \]

幾何分布の連続値版。

🔰 ポアソン分布・指数分布になりそうな例を考えてみよう

ガンマ分布 $~\text{Gamma}(k,~\lambda)$

ポアソン過程の事象k回発生までの待ち時間

e.g., メッセージを2つ受信するまでの待ち時間

\[ \Pr(x \mid k,~\lambda) = \frac {\lambda^k x^{k - 1} e^{-\lambda x}} {\Gamma(k)} \]

指数分布をkのぶん右に膨らませた感じ。
shapeパラメータ $k = 1$ のとき指数分布と一致。

正規分布 $~\mathcal{N}(\mu,~\sigma)$

平均 $\mu$、標準偏差 $\sigma$ の美しい分布。よく登場する。
e.g., $\mu = 50, ~\sigma = 10$ (濃い灰色にデータの95%, 99%が含まれる):

\[ \Pr(x \mid \mu,~\sigma) = \frac 1 {\sqrt{2 \pi \sigma^2}} \exp \left(\frac {-(x - \mu)^2} {2\sigma^2} \right) \]

正規分布に近づくものがいろいろある

標本平均の反復(中心極限定理); e.g., 一様分布 [0, 100) から40サンプル

大きい$n$の二項分布

正規分布に近づくものがいろいろある

大きい$\lambda$のポアソン分布

平均値固定なら$k$が大きくなるほど左右対称に尖るガンマ分布

有名な確率分布対応関係ふりかえり

離散一様分布

コインの表裏、サイコロの出目1–6

負の二項分布 (幾何分布 if n = 1)

成功率pの試行がn回成功するまでの失敗回数

二項分布

成功率p、試行回数nのうちの成功回数

ポアソン分布

単位時間あたり平均$\lambda$回起こる事象の発生回数

ガンマ分布 (指数分布 if k = 1)

ポアソン過程でk回起こるまでの待ち時間

正規分布

確率変数の和、平均値。使い勝手が良く、よく登場する。

現実には、確率分布に「従わない」ことが多い

植物100個体から8個ずつ種子を取って植えたら全体で半分ちょい発芽。
親1個体あたりの生存数はn=8の二項分布になるはずだけど、
極端な値(全部死亡、全部生存)が多かった。

「それはなぜ？」と考えて要因を探るのも統計モデリングの仕事。
「普通はこれに従うはず」を理解してこそできる思考。

疑似乱数生成器 Pseudo Random Number Generator

コンピューター上でランダムっぽい数値を出力する装置。
実際には決定論的に計算されているので、
シード(出発点)と呼び出し回数が同じなら出る数も同じになる。

set.seed(42)
runif(3L)
# 0.9148060 0.9370754 0.2861395
runif(3L)
# 0.8304476 0.6417455 0.5190959
set.seed(42)
runif(6L)
# 0.9148060 0.9370754 0.2861395 0.8304476 0.6417455 0.5190959

シードに適当な固定値を与えておくことで再現性を保てる。
ただし「このシードじゃないと良い結果が出ない」はダメ。

さまざまな「分布に従う」乱数を生成することもできる。

いろんな乱数を生成・可視化して感覚を掴もう

n = 100
x = sample.int(6, n, replace = TRUE)  # 一様分布(整数)
x = runif(n, min = 0, max = 1)        # 一様分布
x = rgeom(n, prob = 0.5)              # 幾何分布
x = rbinom(n, size = 5, prob = 0.5)   # 二項分布
x = rpois(n, lambda = 2.1)            # ポアソン分布
x = rnorm(n, mean = 50, sd = 10)      # 正規分布
print(x)

p1 = ggplot(data.frame(x)) + aes(x)
p1 + geom_histogram() # for continuous values
p1 + geom_bar()       # for discrete values

🔰 正規分布の n, mean, sd を変えて作図し、それぞれの影響を確認しよう。

🔰 ポアソン分布の n, lambda を変えて作図し、それぞれの影響を確認しよう。

🔰 5%の当たりを狙って20連ガチャを回す人が100万人いたら、
    全部はずれ、1つ当たり、2つ当たり… の人はどれくらいいるか？

(Quartoでどうまとめるか、腕の見せ所)

データに分布をあてはめたい

ある植物を50個体調べて、それぞれの種子数Xを数えた。
個体Aは種2個、個体Bは種4個、、、サンプルサイズ n = 50 のデータ。

カウントデータだし形もポアソン分布っぽい。
分布のパラメータ $\lambda$ はどれくらいがいいだろう？

データに分布をあてはめたい

ある植物を50個体調べて、それぞれの種子数Xを数えた。
個体Aは種2個、個体Bは種4個、、、サンプルサイズ n = 50 のデータ。

カウントデータだし形もポアソン分布っぽい。
分布のパラメータ $\lambda$ はどれくらいがいいだろう？

黒が観察データ。青がポアソン分布。 よく重なるのは $\lambda \approx 3$ くらいか。

尤ゆう度 (likelihood)

尤もっともらしさ。 モデルのあてはまりの良さの尺度のひとつ。

あるモデル$M$の下でそのデータ$D$が観察される確率。
定義通り素直に書くと
$\Pr(D \mid M)$

データ$D$を固定し、モデル$M$の関数とみなしたものが尤度関数:
$L(M \mid D)$

モデルの構造も固定してパラメータ$\theta$だけ動かす場合はこう書く:
$L(\theta \mid D)$ とか $L(\theta)$ とか

尤度を手計算できる例

コインを5枚投げた結果 $D$: 表 4, 裏 1

表が出る確率 $p = 0.5$ と仮定:

表が出る確率 $p = 0.8$ と仮定:

$L(0.8 \mid D) > L(0.5 \mid D)$

$p = 0.8$ のほうがより尤もらしい。

種子数ポアソン分布の例でも尤度を計算してみる

$n = 50$個体ぶん、且つ、且つ、且つ、と確率を掛けていく:

この中では $\lambda = 3$ がいいけど、より尤もらしい値を求めたい。

最尤推定 Maximum Likelihood Estimation

扱いやすい 対数尤度 (log likelihood) にしてから計算する。
一階微分が0になる $\lambda$ を求めると…標本平均と一致。

最尤推定を使っても“真のλ”は得られない

今回のデータは真の生成ルール“$X \sim \text{Poisson}(\lambda = 3.0)$”で作った。
「50個体サンプル→最尤推定」を1,000回繰り返してみると:

サンプルの取れ方によってはかなりズレた推定をしてしまう。
(標本データへのあてはまりはかなり良く見えるのに！)

サンプルサイズを増やすほどマシにはなる

“$X \sim \text{Poisson}(\lambda = 3.0)$”からnサンプル→最尤推定を1,000回繰り返す:

Q. じゃあどれくらいのサンプル数nを確保すればいいのか？
A. 推定したい統計量とか、許容できる誤差とかによる。

すべてのモデルは間違っている

確率分布がいい感じに最尤推定できたとしても、
それはあくまでモデル。仮定。近似。

All models are wrong, but some are useful. — George E. P. Box

統計モデリングの道具 — まとめ

• 何はともあれ作図して俯瞰
• 確率変数 $X$
• 確率分布 $X \sim f(\theta)$
  • 少ないパラメータ $\theta$ でばらつきの様子を表現
  • この現象はこの分布を作りがち(〜に従う) という知見がある
• 尤度
  • あるモデルでこのデータになる確率 $\Pr(D \mid M)$
  • データ固定でモデル探索 → 尤度関数 $L(M \mid D),~L(\theta \mid D)$
  • 対数を取ったほうが扱いやすい → 対数尤度 $\log L(M \mid D)$
  • これを最大化するようなパラメータ $\hat \theta$ 探し ＝ 最尤法

🔰 4日目の課題1: 尤度

サイコロを10回振ったら6の目が3回出た。

1. 6の目の出る確率が1/6だとした場合の尤度は?
2. 6の目の出る確率が0.2だとした場合の尤度は?
3. 横軸を6の目の出る確率、縦軸を対数尤度とするグラフを描こう。
4. このサイコロで6の目が出る確率を最尤推定しよう。
   数学で解ければ優。Rで見つければ良。目分量・勘で可。

ヒント

確率pで当たるクジをn回引いてk回当たる確率、と同じ計算を使う。

数学の $\binom 5 2 = {}_5 \mathrm{C} _2 = 10$ はRでは choose(5, 2) 。

🔰 4日目の課題2: 分布を当てはめる

植物25個体から8個ずつ種をとって植え、生き残ったものを数えた。

1. データの分布を描いてみて、当てはまりそうな確率分布を検討する
2. 理論分布を適当なパラメータで描いてみる
3. パラメータや分布を変えてみて、データの分布にすり寄せる
4. 対数尤度の変化を可視化し、パラメータを最尤推定する

trials,survived
8,6
8,6
8,4
8,7
8,8
8,7
8,7
8,3
8,6
8,5
8,8
8,7
8,4
8,8
8,7
8,2
8,8
8,8
8,6
8,7
8,5
8,4
8,7
8,8
8,7

参考文献

• データ解析のための統計モデリング入門 久保拓弥 2012
• StanとRでベイズ統計モデリング 松浦健太郎 2016
• RとStanではじめる ベイズ統計モデリングによるデータ分析入門 馬場真哉 2019
• データ分析のための数理モデル入門 江崎貴裕 2020
• 分析者のためのデータ解釈学入門 江崎貴裕 2020
• 統計学を哲学する 大塚淳 2020
• 科学とモデル—シミュレーションの哲学 入門 Michael Weisberg 2017
  (原著: Simulation and Similarity 2013)