![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter7/baratuki.png?raw=true\")
\n",
"\n",
"因子による効果は**主効果**と呼ばれます。\n",
"\n",
"観測値全体のばらつきを**総変動**と呼び、総変動は、因子の水準の違いに起因するばらつきである**級間変動**と、同じ水準内で偶然生じるばらつき**級内変動**とに分解できます。\n",
"\n",
"$総変動 = 級間変動 + 級内変動$\n",
"\n",
"つまり、分散分析とは級間変動が級内変動に比べて有意に大きいか否かを調べる手法、となります。\n",
"\n",
"具体的な手順を順にみていきます。"
],
"metadata": {
"id": "Qe8XMeQPbPol"
}
},
{
"cell_type": "markdown",
"source": [
"#### 1. 観測値の総平均からのズレを分解する\n",
"\n",
"先ほど示した通り、総変動を水準の効果による級間変動と偶然による級内変動に分解して捉える必要があります。\n",
"\n",
"総変動や級間変動、級内変動を具体的にデータから考えていくと、\n",
"\n",
"まず個々の観測値$y_{ij}$について考えてみます。\n",
"\n",
"観測値$y_{ij}$の、データ全体の平均(総平均$\\bar{y_{\\cdot\\cdot}}$)からのズレ($d_1$)は、\n",
"\n",
"水準の効果に起因する水準平均$\\bar{y_{i\\cdot}}$の総平均$\\bar{y_{\\cdot\\cdot}}$からのズレ($d_2$)と、\n",
"\n",
"偶然誤差による観測値$y_{ij}$の水準平均$\\bar{y_{i\\cdot}}$からのズレ($d_3$)、\n",
"\n",
"に分解して捉えることができます。\n",
"\n",
"つまり、\n",
"\n",
"$(y_{ij}-\\bar{y_{\\cdot\\cdot}})=(\\bar{y_{i\\cdot}}-\\bar{y_{\\cdot\\cdot}})+(y_{ij}-\\bar{y_{i\\cdot}})$\n",
"\n",
"$d_1 = d_2 + d_3$ と表すことができます。\n",
"\n",
"![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter7/zure.png?raw=true\")
"
],
"metadata": {
"id": "fEhLNnYLyWmV"
}
},
{
"cell_type": "markdown",
"source": [
"#### 2. 変動を偏差平方和で表す\n",
"\n",
"先ほどの観測値の総平均からの偏差($d_1$)、水準平均の総平均からの偏差($d_2$)、偶然による観測値の水準平均からの偏差($d_3$)は、各観測値ごと・各水準ごとに値が出るので、\n",
"\n",
"全ての観測値について足し合わせて、データ全体としてまとめて捉える必要があります。\n",
"\n",
"ただし、平均からの差を合計すると0になってしまうので、\n",
"\n",
"例) $\\sum_{i}\\sum_{j}(y_{ij}-\\bar{y_{\\cdot\\cdot}})=0$\n",
"\n",
"分散分析では、個々の偏差を2乗したものの和、**偏差平方和**によってデータの変動を表します。\n",
"\n",
"* 総変動を表す**総平方和**$S_T=\\sum_{i}\\sum_{j}(y_{ij}-\\bar{y_{\\cdot\\cdot}})^2$\n",
"\n",
"* 級間変動を表す**級間平方和**$S_A = \\sum_{i}\\sum_{j}(\\bar{y_{i\\cdot}}-\\bar{y_{\\cdot\\cdot}})^2$\n",
"\n",
"* 級内変動を表す**残差平方和**$S_e = \\sum_{i}\\sum_{j}(\\bar{y_{ij}}-\\bar{y_{i\\cdot}})^2$\n",
"\n",
"\n",
"![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter7/hendou.png?raw=true\")
\n",
"\n",
"実際に下記の様なデータが得られた場合に、\n",
"\n",
"\\begin{array}{cccc} \\hline\n",
" 施肥量 & 10g & 20g & 30g \\\\ \\hline\n",
" & 30 & 31 & 38 \\\\\n",
" 反復 & 25 & 26 & 35 \\\\\n",
" & 29 & 30 & 35 \\\\ \\hline\n",
"\\end{array}\n",
"\n",
"このデータから変動を計算してみると..."
],
"metadata": {
"id": "IJhG0Sc72KRw"
}
},
{
"cell_type": "code",
"source": [
"# 総平方和、級間平方和、残差平方和を求める\n",
"all <- c(30, 25, 29, 31, 26, 30, 38, 35, 35)\n",
"A1 <- c(30, 25, 29)\n",
"A2 <- c(31, 26, 30)\n",
"A3 <- c(38, 35, 35)\n",
"\n",
"all_mean <- mean(all)\n",
"A1_mean <- mean(A1)\n",
"A2_mean <- mean(A2)\n",
"A3_mean <- mean(A3)\n",
"\n",
"ST <- sum((all - all_mean)^2)\n",
"SA <- 3*(A1_mean - all_mean)^2 + 3*(A2_mean - all_mean)^2 + 3*(A3_mean - all_mean)^2\n",
"Se <- sum((A1 - A1_mean)^2) + sum((A2 - A2_mean)^2) + sum((A3 - A3_mean)^2)\n",
"\n",
"print(ST)\n",
"print(SA)\n",
"print(Se)"
],
"metadata": {
"id": "oilIvEtCLJDN",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "65474607-5a71-49bb-8d7e-ede080771bc3"
},
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"[1] 148\n",
"[1] 114\n",
"[1] 34\n"
]
}
]
},
{
"cell_type": "markdown",
"source": [
"ということで、$S_T = 148, S_A = 114, S_e = 34$となり、\n",
"\n",
"\n",
"総平方和$S_T$ = 級間平方和$S_A$ $+$ 残差平方和$S_e$が成り立っています。\n",
"\n",
"これでデータ全体のばらつきを級間変動と級内変動に分割できました。"
],
"metadata": {
"id": "rRd5UNqIS4pp"
}
},
{
"cell_type": "markdown",
"source": [
"#### 3. 級間変動と級内変動の大きさを比較\n",
"\n",
"平方和$S_A$と$S_e$はそれらを構成している変数の数が多くなれば、大きな値をとりやすくなります。\n",
"\n",
"そのため、級間変動が級内変動と比べて意味のある大きさなのかどうかを比較するためには、その部分を考慮する必要があります。\n",
"\n",
"よって、級間平方和と残差平方和をそれぞれの**自由度**で割った**平均平方**を求めて、その比の大きさを調べます。\n",
"\n",
"級間平方和$S_A$は$\\bar{y_{1\\cdot}}-\\bar{y_{\\cdot\\cdot}}, \\bar{y_{2\\cdot}}-\\bar{y_{\\cdot\\cdot}}, ..., \\bar{y_{a\\cdot}}-\\bar{y_{\\cdot\\cdot}}$という$a$個の要素を持ち、そのうち$a-1$個が決まれば残りの1個も決まるため、\n",
"\n",
"級間平方和$S_A$の自由度$\\nu_A$は$a-1$になります。\n",
"\n",
"級内平方和$S_e$は各水準毎に$y_{i1} - \\bar{y_{i\\cdot}}, y_{i2} - \\bar{y_{i\\cdot}}, ... y_{in} - \\bar{y_{i\\cdot}}$という$n$個の要素から成る平方和なので、各水準毎に$n-1$個が決まれば残りの1個も決まるため、\n",
"\n",
"級内平方和$S_e$の自由度$\\nu_{e}$は$a\\times(n-1)$となります。\n",
"\n",
"総平方和$S_T$の自由度$\\nu_{T}$は$N-1 = a\\times n - 1$なので、\n",
"\n",
"自由度についても$\\nu_T = \\nu_A + \\nu_e$が成り立ちます。\n",
"\n",
"\\begin{array}{cccc} \\hline\n",
" 施肥量 & 10g & 20g & 30g \\\\ \\hline\n",
" & 30 & 31 & 38 \\\\\n",
" 反復 & 25 & 26 & 35 \\\\\n",
" & 29 & 30 & 35 \\\\ \\hline\n",
"\\end{array}\n",
"\n",
"今回のデータだと、\n",
"\n",
"$\\nu_A = 3 - 1 = 2$\n",
"\n",
"$\\nu_e = 3 \\times (3 - 1) = 6$\n",
"\n",
"となります。\n",
"\n",
"級間平方和と残差平方和をそれぞれの**自由度**で割った**平均平方**を計算すると、\n",
"\n",
"$V_A = S_A / \\nu_A = 114 / 2 = 57$\n",
"\n",
"$V_e = S_e / \\nu_e = 34 / 6 = 5.666...$\n",
"\n",
"となります。\n",
"\n",
"分散分析では、この様に求めた級間変動と級内変動の2つの平均平方の比$V_A / V_e$(これを**分散比**と呼ぶ)の大きさで、級間変動が級内変動と比較して大きいか評価します。"
],
"metadata": {
"id": "Mhlieeb_TVbp"
}
},
{
"cell_type": "markdown",
"source": [
"#### 4. 帰無仮説と検定\n",
"\n",
"この計算した分散比を評価する方法になります。\n",
"\n",
"分散分析での帰無仮説は、$a$個ある水準の母平均を$\\mu_1, \\mu_2, ..., \\mu_a$とすると、\n",
"\n",
"帰無仮説$H_0: \\mu_1 = \\mu_2 = ... = \\mu_a$ (全ての水準の間に差はない)\n",
"\n",
"対立仮説$H_1:$ 少なくともある2つの水準間には差がある\n",
"\n",
"となります。\n",
"\n",
"このとき、平方和$S_A, S_e$はそれぞれ自由度$\\nu_A, \\nu_e$の$\\chi^2$分布に従うことが知られています。\n",
"\n",
"(導出は「自然科学の統計学」(東京大学出版)等を参照)\n",
"\n",
"よって、母分散の比検定で扱った通り、帰無仮説のもとで、\n",
"\n",
"$F = \\dfrac{S_A/\\nu_A}{S_e/\\nu_e}$\n",
"\n",
"が自由度$\\nu_A, \\nu_e$の$F$分布$F(\\nu_A, \\nu_e)$に従うので、\n",
"\n",
"このことを利用して$F$検定で仮説を検定できます。\n",
"\n",
"これが**分散分析検定**となります。\n",
"\n",
"今回の場合、\n",
"\n",
"$V_A = S_A / \\nu_A = 114 / 2 = 57$\n",
"\n",
"$V_e = S_e / \\nu_e = 34 / 6 = 5.666...$\n",
"\n",
"から分散比を計算して、\n",
"\n",
"$F = V_A/V_e = 57/5.666... = 10.0588$\n",
"\n",
"と$F$値を求められたので、\n",
"\n",
"$F$分布$F(\\nu_A, \\nu_e)$のパーセント点$F_{0.05}(2,6)$を求めると、"
],
"metadata": {
"id": "HMsdi-3qboIV"
}
},
{
"cell_type": "code",
"source": [
"# 自由度2, 6のF分布における95%パーセント点をqf関数で求める\n",
"qf(0.95, 2, 6)"
],
"metadata": {
"id": "CW0k4XFmhNFA",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"outputId": "118c237d-cca4-421d-98eb-692770bff72a"
},
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/html": [
"5.14325284978472"
],
"text/markdown": "5.14325284978472",
"text/latex": "5.14325284978472",
"text/plain": [
"[1] 5.143253"
]
},
"metadata": {}
}
]
},
{
"cell_type": "markdown",
"source": [
"ということで、$F_{0.05}(2,6) < F$となり、有意水準5%のもとで帰無仮説は棄却される。\n",
"\n",
"つまり、施肥量の効果は有意水準5%で有意である、と言えます。"
],
"metadata": {
"id": "fD2Z6-9Wh5t9"
}
},
{
"cell_type": "markdown",
"source": [
"#### 5. 分散分析表\n",
"\n",
"以上の様な分散分析の結果は、最終的に下の様な**分散分析表**としてまとめることが多いです。\n",
"\n",
"\\begin{array}{c|cccc} \\hline\n",
" 要因 & 平方和 & 自由度 & 平均平方 & F値 \\\\ \\hline\n",
" 級間A & S_A & \\nu_A & V_A=S_A/\\nu_A & F=V_A/V_e \\\\\n",
" 誤差E & S_e & \\nu_e & V_e=S_e/\\nu_e & \\\\ \\hline\n",
" 計 & S_T & \\nu_T & & \\\\ \\hline\n",
"\\end{array}\n",
"\n",
"今回の場合だと、下表の様に値が求まりました。\n",
"\n",
"\\begin{array}{c|cccc} \\hline\n",
" 要因 & 平方和 & 自由度 & 平均平方 & F値 \\\\ \\hline\n",
" 施肥量 & 114 & 2 & 57 & 10.0588^* \\\\\n",
" 誤差 & 34 & 6 & 5.6667 & \\\\ \\hline\n",
" 計 & 148 & 8 & & \\\\ \\hline\n",
"\\end{array}\n",
"\n",
"このとき、$F$検定を行った結果、5%水準で有意であれば$*$を、1%水準で有意であれば$**$を分散比$F$の右肩に記す場合が多いです。\n",
"\n"
],
"metadata": {
"id": "MdfdG6Mtig31"
}
},
{
"cell_type": "markdown",
"source": [
"#### 多重比較\n",
"\n",
"さて、今回のデータでは施肥量の効果が5%の有意水準で有意であると判定されました。\n",
"\n",
"ただし、これはあくまでも少なくともある2つの水準間には差があるということが分かっただけで、どの水準間に差があるのかは分かりません。\n",
"\n",
"そこで**多重比較**により平均値間の差の検定を行います。\n",
"\n",
"最初の方でも述べた通り、検定を繰り返してしまうと、第一種の過誤が生じる確率が上がってしまうので、第一種の過誤を抑え込む方法を利用する必要があります。\n",
"\n",
"Tukey-Kramer法やDunnet法、Bonferroni法など様々な手法があります。\n",
"\n",
"目的や扱うデータの特徴に応じて適切なものを選ぶ必要があります。\n",
"\n",
"* Dunnet法 ... ある基準となる群と他の全ての群を比較する検定方法 (A vs B, A vs C)\n",
"* Tukey-Kramer法 ... 全ての群間を比較する検定方法 (A vs B, B vs C, C vs A)\n",
"* Bonferroni法 ... 有意水準を補正する方法、検出力はかなり低くなるので保守的な検定方法\n",
"* 他、Scheffe法，Games–Howell法，Fisher PLSD法など。\n",
"\n",
"一部のみ紹介しましたが、実際に自分の研究データで多重比較をする際には、類似する論文等を読み、適切な多重比較手法を調べて使用する必要がります。\n",
"\n",
"今回は一例としてTukey-Kramer法を使用してみます。\n",
"\n",
"RでTukey-Kramer法を使用するには`TurkeyHSD`関数を使用します。\n"
],
"metadata": {
"id": "URgH3oAHld3C"
}
},
{
"cell_type": "code",
"source": [
"# TukeyHSD関数によるTurkey法を用いた多重比較\n",
"vx <- c(30, 25, 29, 31, 26, 30, 38, 35, 35)\n",
"fx=factor(rep(c(\"A_10g\", \"B_20g\", \"C_30g\"), c(3, 3, 3)))\n",
"TukeyHSD(aov(vx~fx))"
],
"metadata": {
"id": "uhkr-JVbhlZw",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 204
},
"outputId": "684f974b-1280-44bd-aa8b-a06dbc72761c"
},
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
" Tukey multiple comparisons of means\n",
" 95% family-wise confidence level\n",
"\n",
"Fit: aov(formula = vx ~ fx)\n",
"\n",
"$fx\n",
" diff lwr upr p adj\n",
"B_20g-A_10g 1 -4.963654 6.963654 0.8672718\n",
"C_30g-A_10g 8 2.036346 13.963654 0.0147233\n",
"C_30g-B_20g 7 1.036346 12.963654 0.0263670\n"
]
},
"metadata": {}
}
]
},
{
"cell_type": "markdown",
"source": [
"ということで、30gと10g、30gと20gを比較した際に、有意水準5%より低いp値が見られたので、これらの2群間に有意差があると判定できました。\n",
"\n",
"#### 分散分析を多重比較検定の前にやる必要はあるのか？\n",
"\n",
"こうしてみると、分散分析なんてやらずに、いきなり多重検定の手法で複数群間の比較をすれば良いのでは？という気がしてきます。\n",
"\n",
"これはある意味正解で、分散分析をやらないといけない場合とそうでない場合とがあり、これは用いる多重比較法によって異なってきます。\n",
"\n",
"例えばDunnett法やTukey-Kramer法では、分散分析で用いる$F$統計量とは異なる統計量を基に検定を行うため、\n",
"\n",
"「分散分析では有意だが多重比較では有意ではない」や「分散分析では有意ではないが、多重比較では有意」みたいなことが起きる可能性があります。\n",
"\n",
"そのため、これらの手法は分散分析とは独立した検定だと考えておく必要があります。"
],
"metadata": {
"id": "il-47jqpqPtd"
}
},
{
"cell_type": "markdown",
"source": [
"### 二元配置の分散分析\n",
"\n",
"先ほどの実験では施肥量のみを因子として扱いましたが、他にも温度や湿度など、様々な因子が影響を与えます。\n",
"\n",
"複数の因子を扱う実験では、それぞれの各因子単独の効果(**主効果**)だけではなく、それらの因子間で生じる**交互作用**も考慮に入れる必要があります。\n",
"\n",
"![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter7/kougo.png?raw=true\")
\n",
"\n",
"複数の因子を考慮した実験を完全無作為化法で行う場合、それを**多元配置**と呼びます。\n",
"\n",
"今回はその中でも最も簡単な2因子の場合である**二元配置**の場合における分散分析を扱います。"
],
"metadata": {
"id": "Y5b2b8pQqfef"
}
},
{
"cell_type": "markdown",
"source": [
"#### 繰り返しのある二元配置分散分析\n",
"\n",
"まず、下の様な、施肥量と温度の2因子を考慮する場合を考えます。\n",
"\n",
"施肥量3水準、温度2水準の6種類の処理に、それぞれ反復が2回ずつあるデータになっています。\n",
"\n",
"\\begin{array}{cccc} \\hline\n",
" & 施肥量 & 10g & 20g & 30g \\\\ \\hline\n",
" 温度 & 30度 & 30 & 31 & 38 \\\\\n",
" & & 29 & 28 & 36 \\\\ \\hline\n",
" 温度 & 20度 & 23 & 26 & 36 \\\\\n",
" & & 22 & 27 & 35 \\\\ \\hline\n",
"\\end{array}\n",
"\n",
"一般的に、水準$A$(施肥量)と水準$B$(温度)の組み合わせ$A_iB_j$での$k$番目の観測値を$y_{ijk}$とします。\n",
"\n",
"$A$の水準数は$a$、$B$の水準数は$b$、反復は$r$とします。\n",
"\n",
"全観測値の数は$n = a \\times b \\times r$と表せます。\n",
"\n",
"今回のデータの場合は施肥量の水準数$a=3, 温度の水準数b=2, 反復r=2$ということになります。"
],
"metadata": {
"id": "av_T-YUpv7hf"
}
},
{
"cell_type": "markdown",
"source": [
"#### 観測値の総平均からのズレを分解する\n",
"\n",
"一元配置の時と同じように、総変動($S_T$)を級間変動($S_{AB}$)や誤差変動$S_e$に分解していきます。\n",
"\n",
"今回は2因子あるので、級間変動$S_{AB}$は因子$A$による変動$S_A$と因子$B$による変動$S_B$、そして因子$A$と因子$B$の交互作用による変動$S_{A\\times B}$に分解できます。\n",
"\n",
"![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter7/hendou2.png?raw=true\")
\n",
"\n",
"つまり、\n",
"\n",
"$S_T = S_{AB}+S_e$\n",
"\n",
"$S_{AB} = S_A + S_B + S_{A\\times B}$\n",
"\n",
"と表すことができます。\n",
"\n",
"一元配置の時と同じように、データ全体におけるそれぞれの変動を表す偏差平方和は、\n",
"\n",
"総変動を表す総平方和$S_T = \\sum\\limits_{i}\\sum\\limits_{j}\\sum\\limits_{k} (y_{ijk} - \\bar{y_{\\cdot\\cdot\\cdot}})^2$\n",
"\n",
"級間変動を表す級間平方和$S_{AB} = \\sum\\limits_{i}\\sum\\limits_{j}\\sum\\limits_{k} (\\bar{y_{ij\\cdot}} - \\bar{y_{\\cdot\\cdot\\cdot}})^2$\n",
"\n",
"因子Aによる変動を表すA間平方和$S_A = \\sum\\limits_{i}\\sum\\limits_{j}\\sum\\limits_{k} (\\bar{y_{i\\cdot\\cdot}} - \\bar{y_{\\cdot\\cdot\\cdot}})^2$\n",
"\n",
"因子Bによる変動を表すB間平方和$S_B = \\sum\\limits_{i}\\sum\\limits_{j}\\sum\\limits_{k} (\\bar{y_{\\cdot j\\cdot}} - \\bar{y_{\\cdot\\cdot\\cdot}})^2$\n",
"\n",
"交互作用による変動を表す交互作用平方和$S_{A\\times B} = S_{AB} - S_A - S_B$\n",
"\n",
"と求めることができます。\n",
"\n",
"一例として下表のようなデータがあった場合、\n",
"\n",
"\\begin{array}{cccc} \\hline\n",
" & 施肥量 & 10g & 20g & 30g \\\\ \\hline\n",
" 温度 & 30度 & 30 & 31 & 38 \\\\\n",
" & & 29 & 28 & 36 \\\\ \\hline\n",
" 温度 & 20度 & 23 & 26 & 36 \\\\\n",
" & & 22 & 27 & 35 \\\\ \\hline\n",
"\\end{array}\n",
"\n",
"計算過程は省きますが(理由は後程説明)、上のデータから偏差平方和を計算すると、\n",
"\n",
"$S_T = 304.9167$ $S_{AB} = 296.4167$ $S_A = 236.1667$ $S_B = 44.0833$ $S_{A\\times B} = 16.1667$\n",
"\n",
"と求めることができ、誤差変動は$S_e = S_T - S_{AB} = 8.5$ となります。\n",
"\n",
"\n",
"\n"
],
"metadata": {
"id": "mMGd_oo72Sn8"
}
},
{
"cell_type": "markdown",
"source": [
"#### 級間変動と級内変動の大きさを比較\n",
"\n",
"一元配置の時と同じように、平方和$S_A, S_B, S_{A \\times B}$の大きさを$S_e$を基準に判定するためには、自由度を考慮する必要があります。\n",
"\n",
"一元配置の時と同じ考え方で、各平方和の自由度は、\n",
"\n",
"総平方和の自由度$\\nu_T=n-1 = a \\times b \\times r - 1$\n",
"\n",
"級間平方和の自由度$\\nu_{AB}=a\\times b - 1$\n",
"\n",
"A間平方和の自由度$\\nu_A=a-1$\n",
"\n",
"B間平方和の自由度$\\nu_B=b-1$\n",
"\n",
"交互作用平方和の自由度$\\nu_{A\\times B}=\\nu_{AB} - \\nu_A - \\nu_B = (a-1)\\times(b-1)$\n",
"\n",
"残差平方和の自由度$\\nu_e = \\nu_T - \\nu_{AB} = a\\times b \\times (r-1)$\n",
"\n",
"と求める形になります。\n",
"\n",
"先ほど求めた各平方和を自由度で割って平均平方を求めると、\n",
"\n",
"因子$A$の平均平方$V_A = S_A / \\nu_A = 118.0833$\n",
"\n",
"因子$B$の平均平方$V_B = S_B / \\nu_B = 44.0833$\n",
"\n",
"交互作用の平均平方$V_{A\\times B} = S_{A\\times B} / \\nu_{A\\times B} = 8.0833$\n",
"\n",
"残差の平均平方$V_e = S_e/\\nu_e = 1.4167$\n",
"\n",
"と計算できます。"
],
"metadata": {
"id": "5eAOQOjCBYo3"
}
},
{
"cell_type": "markdown",
"source": [
"#### 帰無仮説と検定\n",
"\n",
"今回の検定における帰無仮説は以下の3つになります。\n",
"\n",
"$H_{0(A)}:$ 因子$A$の水準間に差はない(施肥量の効果は無い)\n",
"\n",
"$H_{0(B)}:$ 因子$B$の水準間に差はない(温度の効果は無い)\n",
"\n",
"$H_{0(A\\times B)}:$ 因子$A,B$の間に交互作用はない(施肥量と温度の間に交互作用は無い)\n",
"\n",
"それぞれの検定に用いられる$F$統計量は、分散比を計算して\n",
"\n",
"$H_{0(A)}: F_A = V_A / V_e = 118.0833/1.4167 = 83.3529$\n",
"\n",
"$H_{0(B)}: F_B = V_B / V_e = 44.0833/1.4167 = 31.1176$\n",
"\n",
"$H_{0(A\\times B)}: F_{A\\times B} = V_{A\\times B} / V_e = 8.0833/1.4167 = 5.7095$\n",
"\n",
"と求められます。\n",
"\n",
"それぞれの検定の自由度における$F$分布のパーセント点$F_{0.05}(2, 6), F_{0.05}(1, 6), F_{0.05}(2, 6)$を見ると\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"metadata": {
"id": "vblvQmQEG7IX"
}
},
{
"cell_type": "code",
"source": [
"# それぞれの検定におけるF分布のパーセンタイル点をqf関数で求める\n",
"qf(0.95, 2, 6)\n",
"qf(0.95, 1, 6)\n",
"qf(0.95, 2, 6)"
],
"metadata": {
"id": "w_TeFUSLLufM",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 69
},
"outputId": "0c00e056-6461-40f7-9b46-8fd814c765c1"
},
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/html": [
"5.14325284978472"
],
"text/markdown": "5.14325284978472",
"text/latex": "5.14325284978472",
"text/plain": [
"[1] 5.143253"
]
},
"metadata": {}
},
{
"output_type": "display_data",
"data": {
"text/html": [
"5.9873776072737"
],
"text/markdown": "5.9873776072737",
"text/latex": "5.9873776072737",
"text/plain": [
"[1] 5.987378"
]
},
"metadata": {}
},
{
"output_type": "display_data",
"data": {
"text/html": [
"5.14325284978472"
],
"text/markdown": "5.14325284978472",
"text/latex": "5.14325284978472",
"text/plain": [
"[1] 5.143253"
]
},
"metadata": {}
}
]
},
{
"cell_type": "markdown",
"source": [
"ということで、いずれの検定においても$F > F_{0.05}$となっているので、\n",
"\n",
"因子$A$, 因子$B$, 交互作用いずれの効果も5%有意水準の下で有意であると判定されます。"
],
"metadata": {
"id": "BZNPhNriLz2P"
}
},
{
"cell_type": "markdown",
"source": [
"以上の結果は、一元配置の時と同じように分散分析表として書くことができます。\n",
"\n",
"\\begin{array}{c|cccc} \\hline\n",
" 要因 & 平方和 & 自由度 & 平均平方 & F値 \\\\ \\hline\n",
" 級間A & S_A & \\nu_A & V_A=S_A/\\nu_A & F_A=V_A/V_e \\\\\n",
" 級間B & S_B & \\nu_B & V_B=S_B/\\nu_B & F_B=V_B/V_e \\\\\n",
" 交互作用A\\times B & S_{A\\times B} & \\nu_{A\\times B} & V_{A\\times B}=S_{A\\times B}/\\nu_{A\\times B} & F_{A\\times B}=V_{A\\times B}/V_e \\\\\n",
" 誤差E & S_e & \\nu_e & V_e=S_e/\\nu_e & \\\\ \\hline\n",
" 計 & S_T & \\nu_T & & \\\\ \\hline\n",
"\\end{array}\n",
"\n",
"実際に計算した結果は、\n",
"\n",
"\\begin{array}{c|cccc} \\hline\n",
" 要因 & 平方和 & 自由度 & 平均平方 & F値 \\\\ \\hline\n",
" 施肥量 & 236.1667 & 2 & 118.0833 & 83.3529^{*} \\\\\n",
" 温度 & 44.0833 & 1 & 44.0833 & 31.1176^{*} \\\\\n",
" 施肥量\\times 温度 & 16.1667 & 2 & 8.0833 & 5.7059^{*} \\\\\n",
" 誤差 & 8.5 & 6 & 1.4167 & \\\\ \\hline\n",
"\\end{array}\n",
"\n",
"となります。"
],
"metadata": {
"id": "-fB6WyszK8TO"
}
},
{
"cell_type": "markdown",
"source": [
"#### RによるANOVAの実施\n",
"\n",
"さて、二元配置の例では計算を省略しましたが、理由があります。\n",
"\n",
"当然ですが、他の検定等と同様に、分散分析を実行する`aov`という関数がRにはあります。\n",
"\n",
"この関数を使用することで、煩わしい偏差平方和や平均平方の計算を行ってくれます。\n",
"\n",
"`aov`関数を使用するには、データフレームにデータを整理しておく必要があります。\n"
],
"metadata": {
"id": "w5NQhhmHOE-b"
}
},
{
"cell_type": "code",
"source": [
"# 分散分析を実施するaov関数のためにデータフレームを作成する\n",
"size <- c(30, 29, 23, 22, 31, 28, 26, 27, 38, 36, 36, 35)\n",
"fertilizer <- factor(rep(c(\"A_10g\", \"A_20g\", \"A_30g\"), c(4, 4, 4)))\n",
"temperature <- factor(rep(rep(c(\"B_30\", \"B_20\"), c(2, 2)), 3))\n",
"df <- data.frame(fertilizer=fertilizer, temperature=temperature, size=size)\n",
"df"
],
"metadata": {
"id": "TTgrMH3PM--R",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 474
},
"outputId": "788602ef-3ecf-4b03-ba00-131d17f6721d"
},
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/html": [
"![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter6/false_positive.png?raw=true\")
\n",
"\n",
"この第一種の過誤が有意水準$\\alpha$のことを指していましたが、\n",
"\n",
"有意水準と対するものとして、**検出力**があります。\n",
"\n",
"第二種の過誤の確率を$\\beta$とすると、検出力は$1-\\beta$で表されます。つまり、帰無仮説が正しくないときに、正しく帰無仮説を棄却する確率を表します。\n",
"\n",
"検出力は、実施した検定の能力がどのくらいだったのかを表す指標になります。\n",
"\n",
"検出力を計算する際には、対立仮説をより具体的に考える必要があります。(なのでややこしい)\n",
"\n",
"例えばある島に草丈の標本平均が25㎝の植物がいたとします。\n",
"\n",
"この草丈の母平均が10cmかどうかを検定する場合、\n",
"\n",
"帰無仮説$H_0: \\mu = 10$\n",
"\n",
"対立仮説$H_1: \\mu \\ne 10$\n",
"\n",
"となります。この場合は15㎝の差が有意かどうか調べることになるので、$H_0$が簡単に棄却されそうです。\n",
"\n",
"では草丈の母平均が24cmかどうかを検定する場合、\n",
"\n",
"帰無仮説$H_0: \\mu = 24$\n",
"\n",
"対立仮説$H_1: \\mu \\ne 24$\n",
"\n",
"となり、この場合は1cmの差が有意かどうか調べることになるので、$H_0$はなかなか棄却されなさそうな感じがします。\n",
"\n",
"平均値の差を検定する場合だと、この差によって$H_1$を棄却できる確率は変わってきます。\n",
"\n",
"この検出したい差のことを**効果量**と呼びます。効果量が大きいほど当然ですが検出力は高くなります。\n",
"\n",
"![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter7/power.png?raw=true\")
\n",
"\n",
"また、当然ですがサンプル数も大きく検出力に影響します。\n",
"\n",
"帰無仮説が$H_0: \\mu = 24$だったとしても、サンプル数1,000,000個体位集めれば、標本平均は限りなく25に近付いていくので、$H_0$も棄却できそうです。\n",
"\n",
"この様に、検出力は有意水準$\\alpha$、効果量、サンプルサイズによって決まります。\n",
"\n",
"\n", "\n", "検出力はRの`power.t.test`等の関数によって計算する事が出来ます。\n", "\n", "例えば、$t$検定において、効果量$5$、標本の標準偏差が$1.5$程度、有意水準$\\alpha=0.05$、標本数$10$の時の検出力は…" ], "metadata": { "id": "lpq8OilVYmot" } }, { "cell_type": "code", "source": [ "# ある条件下での検出力をpower.t.test関数で計算する\n", "power.t.test(n = 10, delta = 5, sd = 1.5, sig.level = 0.05, power = NULL, type = \"two.sample\", alternative = \"two.sided\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 223 }, "id": "j79dsraCu5EZ", "outputId": "cff71807-4881-4e57-a53f-60d9dfde14db" }, "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "\n", " Two-sample t test power calculation \n", "\n", " n = 10\n", " delta = 5\n", " sd = 1.5\n", " sig.level = 0.05\n", " power = 0.9999998\n", " alternative = two.sided\n", "\n", "NOTE: n is number in *each* group\n" ] }, "metadata": {} } ] }, { "cell_type": "markdown", "source": [ "power = 0.9999998ということで、とても検出力が高いです。(基準は大体0.8や0.9が使用されることが多いです。)\n", "\n", "効果量(`delta`の部分)をもっと小さくして、効果量$1$にしてみると" ], "metadata": { "id": "qaHXz2h_vsSt" } }, { "cell_type": "code", "source": [ "# ある条件下での検出力をpower.t.test関数で計算する\n", "power.t.test(n = 10, delta = 1, sd = 1.5, sig.level = 0.05, power = NULL, type = \"two.sample\", alternative = \"two.sided\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 223 }, "id": "npfaqR7qwDZM", "outputId": "0d0db27c-ae4f-47a1-ff9b-0132b21c9854" }, "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "\n", " Two-sample t test power calculation \n", "\n", " n = 10\n", " delta = 1\n", " sd = 1.5\n", " sig.level = 0.05\n", " power = 0.2917418\n", " alternative = two.sided\n", "\n", "NOTE: n is number in *each* group\n" ] }, "metadata": {} } ] }, { "cell_type": "markdown", "source": [ "power = 0.2917418ということで、一気に検出力が下がりました。" ], "metadata": { "id": "VBR0fwZmwJcY" } }, { "cell_type": "markdown", "source": [ "### サンプルサイズの設計\n", "\n", "同じ関数`power.t.test`を使用することで、効果量が$1$の時に、検出力0.9を達成するのに必要な標本サイズは？ということも計算できます。\n", "\n", "先ほどの$t$検定において、効果量$1$、標本の標準偏差が$1.5$程度、有意水準$\\alpha=0.05$において、検出力$0.9$を達成するために必要な標本サイズは…" ], "metadata": { "id": "cCPTQWPtwYKu" } }, { "cell_type": "code", "source": [ "# ある検出力を達成するための標本サイズをpower.t.test関数で計算する\n", "power.t.test(n = NULL, delta = 1, sd = 1.5, sig.level = 0.05, power = 0.9, type = \"two.sample\", alternative = \"two.sided\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 223 }, "id": "eRgE0q8TwzNp", "outputId": "8590c74f-f13d-4330-d538-c3e418968116" }, "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "\n", " Two-sample t test power calculation \n", "\n", " n = 48.26431\n", " delta = 1\n", " sd = 1.5\n", " sig.level = 0.05\n", " power = 0.9\n", " alternative = two.sided\n", "\n", "NOTE: n is number in *each* group\n" ] }, "metadata": {} } ] }, { "cell_type": "markdown", "source": [ "$n = 48.26431$ということで、約$50$の標本サイズが必要ということが分かります。" ], "metadata": { "id": "GnjFBI1Vw2iX" } }, { "cell_type": "code", "source": [ "# ある検出力を達成するための標本サイズをpower.t.test関数で計算する\n", "power.t.test(n = NULL, delta = 15, sd = 1.5, sig.level = 0.05, power = 0.9, type = \"two.sample\", alternative = \"two.sided\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 219 }, "id": "z-GeGgxruHxR", "outputId": "6285be27-be51-4f4a-8392-6486359fb6a3" }, "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "\n", " Two-sample t test power calculation \n", "\n", " n = 1.762908\n", " delta = 15\n", " sd = 1.5\n", " sig.level = 0.05\n", " power = 0.9\n", " alternative = two.sided\n", "\n", "NOTE: n is number in *each* group\n" ] }, "metadata": {} } ] }, { "cell_type": "markdown", "source": [ "この様な形で、目的とする効果量をある程度の検出力を持って検出するために必要なサンプルサイズの設計ができます。\n", "\n", "$t$検定以外にも、`power.prop.test`や`power.anova.test`などの関数が準備されており、様々な検定の検出力・サンプルサイズ計算が可能です。\n", "\n", "実験計画を立てる際に、大まかにこの様な計算をしてサンプル数等を設計できれば理想的です。\n", "\n", "サンプルサイズが小さすぎると、本来検出出来る差が検出できず実験自体に意味がなくなってしまう可能性があります。\n", "\n", "また、サンプルサイズが大きすぎると、コストの問題や統計的に有意になりやすくなるという問題があります。\n", "\n", "したがって、実験は適切なサンプルサイズに対して行う必要があります。" ], "metadata": { "id": "593URUd_uKlE" } }, { "cell_type": "markdown", "source": [ "## 研究倫理・注意点(p値ハッキング)\n", "\n", "最後に、検定を行う際に避けるべき行為である**p値ハッキング**について説明します。\n", "\n", "p値ハッキングとは、統計的に有意でないデータを、有意と示してしまう誤ったデータ分析の行為です。\n", "\n", "多重比較で扱った検定を繰り返す行為もこれにあたります。\n", "\n", "(本来有意ではないものが、検定を繰り返すことによって有意と出てしまう。)\n", "\n", "他にも、\n", "\n", "* 検定の結果有意では無かったので、新たなデータを取得・追加した\n", "* 恣意的に外れ値と判断し、いくつかの測定値を除外した\n", "* 都合が良い測定値を選択的に抽出した\n", "* 有意になるように検定手法を変えたり、検定に微調整を加えた\n", "\n", "などがこの行為に当たります。\n", "\n", "p値ハッキングは意図的に行う場合もあれば、意図せず行ってしまっていた場合もあります。\n", "\n", "この様な不正行為を避けるためにも、使用する統計手法について正しく理解し、下記の様な点について常に気を付けておく必要があります。\n", "\n", "* 事前に研究計画を立てる\n", " * 結果に合わせて検定手法を変えてはいけない\n", "* p < 0.05はあくまでも1つの基準なので、実際にデータが意味する所は何かを意識する\n", " * p = 0.051のため有意差が無い=本当に意味がない？\n", "\n", "* データの選択を行う場合、理論的に正しい基準で実施されているかどうか確かめる\n", "* 解析したデータの透明性を確保する" ], "metadata": { "id": "3AVJ3F1NXb9O" } } ] }