![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter9/lm.png?raw=true\")
\n",
"\n",
"この$\\epsilon_i$**が正規分布に従っている**ということを仮定においたモデルを**一般線形モデル**(**線形モデル**)と呼びます。\n",
"\n",
"分散分析や$t$検定なども一般線形モデルの枠組みで扱う事が出来ます。"
],
"metadata": {
"id": "yvQIppd1QxRm"
}
},
{
"cell_type": "markdown",
"source": [
"ではこの仮定($\\epsilon_i$**が正規分布に従っている**)に合わない場合はどうなるでしょうか。\n",
"\n",
"例えば下図のような場合、左の様に正規分布を仮定した場合よりも、右の様にポアソン分布を仮定した方が当てはまりがよさそうです。\n",
"\n",
"しかし、正規分布を仮定しない場合は、線形モデルでは扱う事が出来ません。\n",
"\n",
"![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter9/lm_poisson.png?raw=true\")
\n",
"\n",
"こうした問題に対応するために、一般線形モデルを拡張する必要があります。\n",
"\n",
"一般線形モデルは残差の分布が正規分布に限られていましたが、\n",
"\n",
"これを「正規分布以外」の二項分布やポアソン分布等も扱える様に拡張したモデルが**一般\"化\"線形モデル**(Genelarized Linear Model, GLM)になります。\n",
"\n",
"この時、最小二乗法ではパラメータ$\\beta$の推定が出来ないことが多いので、\n",
"\n",
"確率分布に基づいた**最尤法**という手法を用いて回帰モデルを推定する必要があります。"
],
"metadata": {
"id": "0GtmA3FZEqWW"
}
},
{
"cell_type": "markdown",
"source": [
"### 二項分布の場合\n",
"\n",
"まずは二項分布の例を使って一般化線形モデルについて触れていきます。\n",
"\n",
"例えば下の様なデータがあったとします。\n",
"\n",
"ある土壌栄養の条件(`nutrition`)における種子の発芽率(`rate`)を示しています。\n",
"\n",
"調査は全て10個の種を使用し(`size`)、発芽した種の数は`germination`で表されています。"
],
"metadata": {
"id": "HYB1LsxpGtYt"
}
},
{
"cell_type": "code",
"source": [
"# 各肥料条件における発芽率のデータ\n",
"data <- read.csv(\"https://raw.githubusercontent.com/slt666666/biostatistics_text_wed/refs/heads/main/source/_static/data/chapter9_germination1.csv\")\n",
"# データの一部を表示\n",
"head(data)"
],
"metadata": {
"id": "mz7pGOOkMS9F",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 286
},
"outputId": "84a8c471-facf-4288-8f2c-d0911cf6cc85"
},
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/html": [
"![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter9/over100.png?raw=true\")
\n",
"\n",
"では右図の様な回帰をどのように行っていくかについて学んでいきます。"
],
"metadata": {
"id": "1FWaJ8lhJ0pT"
}
},
{
"cell_type": "markdown",
"source": [
"#### ロジット変換\n",
"\n",
"この時用いられるのが**ロジット変換**と呼ばれる方法になります。\n",
"\n",
"今回の発芽率の様な確率$p$は$0 \\sim 1$の範囲の値しかとらないですが、\n",
"\n",
"この値をロジット$logit(p) = log\\dfrac{p}{1-p}$に変換することで、$(-\\infty \\sim \\infty)$の値をとるように変換することが出来ます。\n",
"\n",
"そのため、このロジットには線形モデルを仮定しても先ほどの様な不合理性が生じないことになります。\n",
"\n",
"$logit(p) = log(\\dfrac{p}{1-p}) = \\beta_0+\\beta_1x$\n",
"\n",
"新たな目的変数に変換することで線形モデルで説明することが可能になるので、\n",
"\n",
"あとはこれまでと同じようにパラメータ$\\beta$を推定すれば良いというわけです。\n",
"\n",
"このように、一般の回帰曲線の式を線形の式に変形するような関数、\n",
"\n",
"今回であればロジット関数$logit(p) = log(\\dfrac{p}{1-p})$を**リンク関数**と呼びます。\n",
"\n",
"また、$\\beta_0+\\beta_1x$の部分を**線形予測子**と呼びます。\n",
"\n",
"![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter9/link_function.png?raw=true\")
\n",
"\n",
"ロジット関数の逆関数を求めると、\n",
"\n",
"$p = \\dfrac{e^{\\beta_0+\\beta_1x}}{1 + e^{\\beta_0+\\beta_1x}} = \\dfrac{1}{1 + e^{-(\\beta_0+\\beta_1x)}}$ となるので、\n",
"\n",
"$g(x) = \\dfrac{1}{1 + e^{-(\\beta_0+\\beta_1x)}}$と表される関数を**ロジスティック関数**と呼びます。\n",
"\n",
"この関数の形を見てみましょう。\n",
"\n",
"Rの`stat_function`という関数を使うと、任意の数式をグラフ化する事が出来ます。"
],
"metadata": {
"id": "W9RzTswn6z-s"
}
},
{
"cell_type": "code",
"source": [
"# ロジスティック関数を可視化する\n",
"library(ggplot2)\n",
"\n",
"beta0 <- -1\n",
"beta1 <- 1\n",
"\n",
"g <- ggplot(data=data.frame(X=c(-5,5)), aes(x=X))\n",
"g <- g + stat_function(fun=function(x) 1/(1+exp(-(beta0+beta1*x))))\n",
"g"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 437
},
"id": "2PgnGgQN2kLR",
"outputId": "30baa2bf-39ee-4ad1-f22c-a3cd354b9e20"
},
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
"plot without title"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAADAFBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4\nuLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnK\nysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc\n3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u\n7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////i\nsF19AAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3deXxU5bnA8SchiSwKigu4UHdc\nb2uLaF0Re71uJNgCisgi1hUV0UqxAgoi4AUBhWorvYh1vaLXDdwQ92pBEJHViGFNyIIoOyQh\nOXeyMiGTnDNnngN53/P7/ZE5c86bJ+cz5PvJZJjMiENECSf7+gSIbAhIRAoBiUghIBEpBCQi\nhYBEpBCQiBQCEpFCepC2/+RWSaHrEh9t2xrE1OLiIKZucb+RfFRY8nMAUzftDGDoTztKNgUw\n9eeiAIb+tK1ki9uSjQFA2rreLafYdYmPtm0JYuqukiCmbt4exNQi58cApv5cGMDQ9TucnwOY\n+mMg31lbnc1uS34CkktAAhKQFAISkICkEJCABCSFgAQkICkEJCABSSEgAQlICgEJSEBSCEhA\nApJCQAISkBQCEpCApBCQgAQkhYAEJCApBCQgAUkhIAEJSAoBCUhAUghIQAoEUvafOldtbhnX\np8fw/N2XQPIckMIO6bPeE6ohjRi0Imfs7SXVl0DyHJDCDunDgn9XQVqfkRX5aXTVgqpLIHkP\nSGGH5DjVkL7sUhr5eMfLVZeRDzs/iPTdZrecXa5LfLRzRxBTS0qDmLq9MIipxc6WAKZuLQ5g\n6OYiZ2sQY4P5znK2uy3ZkgCk964v+zh4ctVl5EN+u0hPe51EZEvVv9j4gdS37GMEUuVl5MO2\nZyJ9vdUtp8R1iY+KCoOYWloaxNSdRUFM3eVsC2Dq9l0BDN1a5GwPYOq2QL6zCp2drl84AUiz\nK+7SvVJ1WbWA35Hc43ckfkeqhrQhY7njbOq8uOoSSN4DUtgh/bR+ZueyW2HmW44zesCK7GH3\nlFZfAslzQAo7pD+ml/WmM2ZI5DeiCb2vG/XT7ksgeQ5IYYfkISC5ByQgAUkhIAEJSAoBCUhA\nUghIQAKSQkACEpAUAhKQgKQQkAyClDcvD0gJB6QwQ8r9/K83ndVUHgdSwgEppJAKZj9542/S\nRKTRSde8B6SEA1IIIa16deDFB5UZOrXHI++s5nckjYAUMkhLp95yRkoE0RGdHnxzdcUuICkE\npBBByny674kRQyln3Dx1UdRuICkEpJBAWvXibacnizQ5/95XVu1xCEgKASkMkD578Pw0kdSz\n730zJ8ZRICkEJNshrXym15GR+3On3vHy6jpWAEkhIFkNacm430V+FDXvNGFhPYuApBCQ7IU0\n58H2kd+K2vZ/M7f+dUBSCEiWQpo3+BSR5PYPzHZfCiSFgGQjpIUjfiOSevGjSzytBpJCQLIO\nUtb4c5Il+fxxmV6nAkkhINkFKf/Vbo0l6cyRi9wWRgUkhYBkE6S5fzpKpM3AefFNBZJCQLIG\nUu7UC5OkydWv5cc7FUgKAckSSIsHHS7S/rEVPqYCSSEgWQHpjYxUadr7Y39TgaQQkMyHtO5v\np4u0HZ3ldyqQFAKS6ZBWDD9Skq98rcD/VCApBCSzIS0ZeKCkXf1FQlOBpBCQTIb0dc80OXjg\ndwlOBZJCQDIX0rweqXL02LUJTwWSQkAyFdL83inSZpzLE7s9BSSFgGQmpPndU+S4SRqMgKQS\nkEyElHlbmpzwRJ7SVCApBCTzIK0d2kIOf0yLEZBUApJpkPKntJFmAxN/iGF3QFIISIZBev0U\nSbvV858aeQpICgHJKEgrr5KkLvOVpwJJISAZBGndyAPk9BnqY4GkEJDMgfR/J8lBI/UeY6gO\nSAoByRRIi66Q5N4BMAKSSkAyA1LB+OZy5ize+lIhIIUY0tcdpMnQPN5DViMghRZS/rhm0rHs\nsTogKQSksEL67Ndy0KTyLSApBKRwQip4KE0yllZsA0khIIUS0sILpeXUqitAUghIYYQ0taV0\n2P2WLEBSCEjhg7Syt+w3NOoFH4GkEJBCB+mdo+W0z6N3AEkhIIUN0sjU5H413/MVSAoBKVyQ\nVmZIy2l77AOSQkAKFaQvT5Yzav29BJAUAlKYIP2tqfTOqbUXSAoBKTyQcm6SZpNj7AeSQkAK\nDaQl7eTkL2MdAJJCQAoLpI+PlIzVMY8ASSEghQTSSwck9a/jXfeApBCQwgFpZHLak3UdA5JC\nQAoDpJzu0ur9Oo8CSSEghQDSsvbyy2/rPgwkhYBkP6S5x0qnNfUcB5JCQLIe0keHSb9638IS\nSAoByXZIbzRPGlj/CiApBCTLIT2zX8pjLkuApBCQ7Ib0SHLTl9zWAEkhIFkNaagc9I7rIiAp\nBCSLIRXcKG1iPruuZkBSCEj2QiroK20Xui8DkkZAshZSwQ1y4mIvC4GkEJBshZR/rZy0xNNK\nICkEJEsh5XWX0z2+lyWQFAKSnZDyrpZfen1PWCApBCQrIa3LkHY/eJ0KJIWAZCOk/C7SfoXn\nqUBSCEgWQiroI7/27ghIGgHJQkj95RSvvx+VBSSFgGQfpPvlmEXxTAWSQkCyDtIjcvjXcU0F\nkkJAsg3SE8kt/xXfVCApBCTLIP0z5YBZcU4FkkJAsgvSG2lNZsQ7FUgKAckqSP86MOXluKcC\nSSEg2QTpu2NlTPxTgaQQkCyCtKad3ONjqlWQdmx1yylxXeKjwp1BTC0tDWLqzqIgpu5ytgUw\ndXtxAEO3Fjnb6zm6qZN03eJj6rZgvrMc12+tbQFA2r7FLWeX6xIfFe4MYmpJaRBTdxQGMXWX\nszWAqduLAxi6pcjZVs/RW+S8H/1M3Vri83Tqbaezw/ULBwCJu3bucdeu3rt2Q+Wk5b6mWnXX\nDkjuAak+SFOSW9V6d1hvAUkhIFkC6f39mn3ocyqQFAKSHZAWt05+3u9UICkEJCsgrfut3O97\nKpAUApIVkHpKp3rfcKLegKQQkGyANEpOjf0+y54CkkJAsgDS9LSD5iYwFUgKAcl8SPMPTnk9\nkalAUghIxkNa+ysZndBUICkEJNMhFWTItYlNBZJCQDId0ihpl5PYVCApBCTDIc1MO9DnM4Oq\nA5JCQDIb0g+/SHo20alAUghIRkMquFL6JzwVSAoByWhIw+XMdQlPBZJCQDIZ0vtpByX6C9J6\nIKkEJIMhLW+THP9rBtUOSAoByVxIBZf7eq2TWgFJISCZC+kBOS9PYyqQFAKSsZBmpR7q6U3L\nXQOSQkAyFdKaE5Ne0JkKJIWAZCqk6+WPSlOBpBCQDIX0clLbNUpTgaQQkMyE9N1hqTO1pgJJ\nISCZCelSGaI2FUgKAclISGPkbJVHvssDkkJAMhHSV82aKzw1qCogKQQkAyHl/kb+pjgVSAoB\nyUBId8sfNKcCSSEgmQdpVsqRP2hOBZJCQDIOUu5/yP+qTgWSQkAyDtK90l13KpAUApJpkD5N\na+Xv/cTqDEgKAckwSLm/kmeUpwJJISAZBul+6ao9FUgKAcksSF/t13KZ9lQgKQQkoyCVnCVP\nq08FkkJAMgrSf0uG/lQgKQQkkyAtbNJyqf5UICkEJIMgFZwv/whgLJAUApJBkCbKpXW9q3ki\nAUkhIJkD6fuDm/wAJCAlXNgh9ZBhdb2reUIBSSEgGQNpRtLxG4EEpMQLN6Tc05Jeq+NdzRMM\nSAoByRRIQ6R7He9qnmhAUghIhkD6pulBy4AEJI1CDelSeSz2mzEnHpAUApIZkKbK2QVAApJK\nIYa08oiUT9YDCUgqhRjS7XJn2QWQgKRQeCH9O+2I1WWXQAKSQuGF9DuZXH4JJCApFFpIL0j7\ngvINIAFJobBCyjk++YOKLSABSaGwQnpAelduAQlICoUU0tLmLb6r3AQSkBQKKaTuMrJqE0hA\nUiickGYlt11XtQ0kICkUSkgF7eTl6itAApJCoYQ0Sa7YfQVIQFIojJBWtkqbs/sakICkUBgh\n3S0Doq4BCUgKhRDSwiaHrIi6CiQgKRRCSNfK2OirQAKSQuGD9FmjE9ZFXwcSkBQKH6SO8lyN\n60ACkkKhg/S6nFNzB5CApFDYIOX/MumdmnuABCSFwgZpkvx+jz1AApJCIYOUfVTaV3vsAhKQ\nFAoZpCFy6567gAQkhcIFKbN5i8w99wEJSAqFC9KNMqzWPiABSaFQQZqb1ia71k4gAUmhUEHq\nKk/U3gkkICkUJkifJ5+cV3svkICkUJggXSFTY+wFEpAUChGkD5J+VRBjN5CApFCIIF0k02Lt\nBhKQFAoPpOlydsz9QAKSQuGBdI68GXM/kICkUGggvSi/i30g7JC2jOvTY3h++ebC9PJmOHeW\nXXQDUhyFBVLBfyS9H/tI2CGNGLQiZ+ztJWWbRWUDlnRb4/SdHtnYAKQ4Cguk/5H0Oo6EHNL6\njKzIT6WrFlTvGPKi43SdW2MNkNwLCaS8to3+VcehkEP6sktp5OMdL1dd/+yPxU5R+sS7bhiV\nDaQ4CgmkidK9rkMhh/Te9WUfB0+uvFpy6weOs7HX+MzMYb22lhmaGOnf291ySlyX+KioMIip\npaVBTC0sDmJqiRPE1J27/H7mxqPTltZ1rNjZ6XdsfQXzneW4f2vFC6lvDUifXb+rcmt7t5mR\nj/ntIj3taRLZ31Ny274+hb1VSfWWN0izK+7avVJ5dfjk6iP9Ir8sOYVzIq3Y6Jazy3WJj3Zu\nD2JqSWkQU7ftDGJqsbMpgKlbin1+4vqjUxfVebDQ2eJzbH1tCuQ7a4ezzW3J5jghbchY7jib\nOi+uuLa1/FGHVZOKHWdHt4+q1vA7knuh+B1pglxf98GQ/47kjB6wInvYPaXOzLciVxakl/2P\n0uYeE3KzR/XdCSTvhQHSuqNT59d9NOyQtk3ofd2oyGeNGRK58nFGcdm+rCHX9ByRV70ESO6F\nAdLj1W+8HKuwQ/IQkNwLAaS841K/rucwkICkUAggTZKe9R0GEpAUsh9S3vGp8+o7DiQgKWQ/\npL/KtfUeBxKQFLIeUt4JjWbXuwBIQFLIekh/q/tZdhUBCUgK2Q4p78RG/65/BZCApJDtkJ6S\nq11WAAlIClkOKf8ktx9IQAKSRpZDmip/cFsCJCApZDmkM5M+clsCJCApZDekV+VS1zVAApJC\ndkO6QN52XQMkIClkNaT35Vz3RUACkkJWQ7o89qt91wxIQFLIZkhfJJ8e6+0n9ghIQFLIZkhX\ny9MeVgEJSApZDOmb1GNivEFfrYAEJIUshnSDPO5lGZCApJC9kJY1PiLHyzogAUkheyHdJQ97\nWgckIClkLaQVLVqu8rQQSEBSyFpIg+U+bwuBBCSFbIWUfej+33tbCSQgKWQrpHHSz+NKIAFJ\nIUshFbRN/cbjUiABSSFLIT0r3bwuBRKQFLIU0m/F9Q/6qgISkBSyE9JM6eB5KpCApJCdkNLl\nZc9TgQQkhayE9HWjUzz8/URlQAKSQlZCulEmeZ8KJCApZCOk75u19vR01YqABCSFbIR0nwyN\nYyqQgKSQhZByDtv/hzimAglIClkIabzcFs9UIAFJIfsgFZyUUs97mNcOSEBSyD5Iz0uXuKYC\nCUgK2QfpXJkV11QgAUkh6yB95OXVVaMDEpAUsg7SNfLP+KYCCUgK2QZp2X5tvLyYXVRAApJC\ntkEa6PG1g3YHJCApZBmknMP2z4pzKpCApJBlkCbKzfFOBRKQFLIM0unJc+OdCiQgKWQXpDfk\n8rinAglICtkF6Qp5Pe6pQAKSQlZBmh/PX8ZWBSQgKWQVpFtkYvxTgQQkhWyCtLL5wdnxTwUS\nkBSyCdLDcq+PqUACkkIWQco/Jm2xj6lAApJCFkF6Vq7xMxVIQFLIIkgXxvmHSJUBCUgK2QPp\ni6SzfE0FEpAUsgdSX3nK11QgAUkhayCtPOCwOF4VMiogAUkhayCNlIH+pgIJSArZAqnghNSF\n/qYCCUgK2QJpmvze51QgAUkhWyBdKjN8TgUSkBSyBNL8Rqf5nQokIClkCaTb5TG/U4EEJIXs\ngLS25YFr/E4FEpAUsgPSBLnD91QgAUkhOyD9Knme76lAApJCVkCaIZf6nwokIClkBaSrZJr/\nqUACkkI2QFqcely+/6lAApJCNkD6s4xMYCqQgKSQBZByD28a7+t9RwckIClkAaSp0juRqUAC\nkkIWQOogHyYyFUhAUsh8SF8lnZnQVCABSSHzId0qf01oKpCApJDxkNYe5P9pduUBCUgKGQ9p\nkvRLbCqQgKSQ8ZB+kzQ7salAApJCpkP6WDomOBVIQFLIdEg95Z8JTgUSkBQyHFJW09brEpwK\nJCApZDik0TIo0alAApJChkM6JcXnq9ntDkhAUshsSG9Ip4SnAsm1na45Je5r4q+4OIippaVB\nTC3aFcTUEg+3ffwV1vrX6irvJjx1l1OY8IzaFQbyr1XsFLmuCQDStp/dcopdl/hoh/sX9lFJ\naRBTt+0MYmqxszGAqZuL9tiRmXbcTwlP3elsTnhG7TbuCmDoz9udrW5LNgUAibt27hl9126w\nPJT4VO7aAUkhkyEVHLvf94lPBRKQFDIZ0jS5WmEqkICkkMmQOvl+4fzogAQkhQyGtCT1xAKF\nqUACkkIGQxosD2tMBRKQFDIXUsExGg81AAlIKpkLSeehBiABSSVzIek81AAkIKlkLCSlhxqA\nBCSVjIU0OKHXKY4KSEBSyFRIWg81AAlIKpkKSeuhBiABSSVTIWk91AAkIKlkKKQlqW11HmoA\nEpBUMhTS/VoPNQAJSCqZCangmMZKDzUACUgqmQlJ76EGIAFJJTMhpct0talAApJCRkLKTDtB\n66EGIAFJJSMhPSjD9aYCCUgKGQmpbdoyvalAApJCJkJ6SzIUpwIJSAqZCKmbTFOcCiQgKWQg\npKwmR+UpTgUSkBQyENJouU9zKpCApJCBkE5v9I3mVCABSSHzIM2U/1KdCiQgKWQepN7yrOpU\nIAFJIeMgrT7gsETf7LJmQAKSQsZBekwG6E4FEpAUMg7SmUlf6U4FEpAUMg3SZ3Kh8lQgAUkh\n0yDdLJOVpwIJSAoZBimnZcts5alAApJChkGaLLdoTwUSkBQyDFIH+VR7KpCApJBZkOYnt1Of\nCiQgKWQWpHtlvPpUIAFJIaMg7WjTdIX6VCABSSGjIM2Q7vpTgQQkhYyC1EXtBb+jAhKQFDIJ\nUtZ+x+u9Cld1QAKSQiZBGikPBDAVSEBSyCRIp6UsDmAqkICkkEGQ3pN0/aFAApJKBkHqJa/p\nDwUSkFQyB9LqAw5z/wf1EZCApJA5kCbJ3YXuq+IPSEBSyBxIZ8tXQAJS4oUc0uykc38GEpAS\nL+SQ7pQngAQkhcINad1hzdcACUgKhRvSs9JnPZCApFC4IV0m7wMJSBqFGtLS1JPXAwlIGoUa\n0gPyEJCApFKoIbVNWQIkIKkUZkjvSKf1QAKSSmGG1FNeWA8kIKkUYkirDzgsdz2QgKRSiCFN\nkrvKLoAEJIVCDOkc+bLsAkhAUii8kOYlnV1+CSQgKRReSAPk8fJLIAFJodBCyj+y6cryDSAB\nSaHQQnpJelRsAAlICoUWUnrV66sCCUgKhRVSZlrV66sCCUgKhRXSwzKkcgtIQFIorJBOTVlU\nuQUkICkUUkiz5JKqTSABSaGQQvqjTK3aBBKQFAonpOwDW+ZUbQMJSAqFE9JkubV6G0hAUiic\nkDrKJ9XbQAKSQqGE9G2jM3ZfARKQFAolpPvkv3dfARKQFAojpIJj9lu++xqQgKRQGCG9IX+I\nugYkICkURkjXyCtR14AEJIVCCGlVsyPyoq4CCUgKhRDSePlT9FUgAUmhEEJqnzQ3+iqQgoC0\nZVyfHsPzK7bvTI/UreY+IHmqIUOanXR+jetACgLSiEErcsbeXlK+3Xd6ZMKGmvuA5KmGDOlO\neaLGdSAFAGl9RlbkJ9BVC8qvdJ1bex+QPNWAIeW2OmBNjR1ACgDSl11KIx/veLlsuyh94l03\njMqusQ9I3mrAkJ6XXjV3ACkASO9dX/Zx8OSyjxt7jc/MHNZra9S+n2+L9HaRW06p6xIf7doV\nxNRSJ4ipu0qCmKpyrp3ls5o7igM51xKnOIixwXxnOa7fWoWxIRXVCanvbkjlbe82M2pffrtI\nT3siSQ2z/NST9/UpmFj1IwQ1IR3S/+vY62dX3I17Zfeefi9G7SvZFOmnH91yil2X+Gj71iCm\n7ioJYuqWHUFMLXI2JDzjIRm2x56NhQkPjdEOZ2MAUzcE8p21zdnstuTn2JAuSpbTx6yLAWlD\nxnLH2dR5cdn2qknFjrOj20fR+8ridyT3Gu7vSKemLN5jD78jJfA7Uu6kC5IaXfbS9lqSRg9Y\nkT3snlJn5lvO5h4TcrNH9d1ZvQ9InmuwkGbKpXvuAlJiDzbkPNZemt84Z4+92yb0vm5U5LPG\nDHGcrCHX9ByRt3sfkDzXYCH1lWf23AWkRB+1W9xDRM6dG+MeXr0Byb2GCin6RU+qAlJCkPLG\n/VIaXfl/b53Z6D0gqddQIf1dbqu1D0j+IRW+mp4iJ40ue7ih6LLjgaReQ4XUQT6ttQ9I/iG1\nlP37fl65/UYSkNRroJC+bfTr2juB5B/S+VO2Vm+vmQIk9RoopD/LmNo7gcTfIykUJkg1X/Sk\nKiABSaEwQXpNusTYCyQgKRQmSN3k1Rh7gQQkhUIEaUWTGi96UhWQgKRQiCA9KgNj7QYSkBQK\nEaR2NV/0pCogAUmh8ED6t1wQcz+QgKRQeCDdIU/G3A8kICkUGkjrDmu+JuYBIAFJodBAelb6\nxD4AJCApFBpIl8l7sQ8ACUgKhQXS0tST6zgCJCApFBZID8hDdRwBEpAUCguktilL6jgCJCAp\nFBJI70inug4BCUgKhQRST3mhrkNAApJC4YC0+oDDcus6BiQgKRQOSJPkrjqPAQlICoUD0jny\nZZ3HgAQkhUIBaV7S2XUfBBKQFAoFpAHyeN0HgQQkhcIAKf/IpivrPgokICkUBkgvSY96jgIJ\nSAqFAVK6zKjnKJCApFAIIGWmHV9Qz2EgAUmhEEB6WIbUdxhIQFIoBJBOTVlU32EgAUkh+yHN\nlMvqPQ4kIClkP6Q+8ly9x4EEJIWsh7S2xaHr6l0AJCApZD2kv0r/+hcACUgKWQ+pvuerlgck\nIClkO6R5Sb91WQEkIClkO6S7ZJLLCiABSSHLIeW2brbKZQmQgKSQ5ZCel95uS4AEJIUsh3S5\nvO+2BEhAUshuSEtTT3FdAyQgKWQ3pAfkYdc1QAKSQnZDapv2nesaIAFJIashzZAM90VAApJC\nVkPqLtPcFwEJSArZDGlVsyPy3FcBCUgK2QxpvPzJwyogAUkhmyH9JvlrD6uABCSFLIb0mVzk\nZRmQgKSQxZBulClelgEJSArZCyn7wJY5XtYBCUgK2QvpSennaR2QgKSQvZDOlX95WgckIClk\nLaTZ9b2VS3RAApJC1kK60/VPYysDEpAUshVSbqvma7ytBBKQFLIV0lS5weNKIAFJIVshXSwf\nelwJJCApZCmkBY1+5XUpkICkkKWQBspYr0uBBCSF7ISU36bxD16nAglICtkJ6SXp7nkqkICk\nkJ2QrpR3PE8FEpAUshLS0tQT63vX2JoBCUgKWQlpqDzkfSqQgKSQjZAKjvPwKlzVAQlICtkI\n6RXpGsdUIAFJIRshpctbcUwFEpAUshDS0tQTvD/UACQgqWQhpCEyIp6pQAKSQvZBKjg2noca\ngAQkleyDNC2uhxqABCSV7IPUSabHNRVIQFLIOkhxPauhLCABSSHrIA328N5iNQISkBSyDVLB\nsWmZ8U0FEpAUsg3SNOkW51QgAUkh2yDF+1ADkICkkmWQlsT7UAOQgKSSZZDifqgBSEBSyS5I\nBcfu9328U4EEJIXsghT/Qw1AApJKdkG6XGbEPRVI+whS4S63nFLXJT4qKQliquMEMbUkkFug\n1O1cVzY6Nf6pQZ1rIP9cwXxnuZ9rcQCQ+Ink3j76iXSPjIl/Kj+RuGunkE2Qcls3WxH/VCAB\nSSGbIE2Rvj6mAglICtkE6Xz5yMdUIAFJIYsgfZl0lp+pQAKSQhZBuln+5mcqkICkkD2Q1h7U\nMtvPVCABSSF7ID0ud/qaCiQgKWQPpDOS5viaCiQgKWQNpI/kYn9TgQQkhayB1EOe8zcVSEBS\nyBZIPzQ9Ms/fVCABSSFbII2Uv/icCiQgKWQLpJNTFvmcCiQgKWQJpNck3e9UIAFJIUsgXSlv\n+J0KJCApZAekb1JOivfFg6oDEpAUsgPSXfKo76lAApJCVkDKOaTFat9TgQQkhayANElu9T8V\nSEBSyApIv06a7X8qkICkkA2Q3pXfJTAVSEBSyAZIXeXFBKYCCUgKWQDpu7Rj8hOYCiQgKWQB\npPvkoUSmAglICpkPKfeIJnG/cH50QAKSQuZDmiK9E5oKJCApZD6kc+TThKYCCUgKGQ/p86Rz\nEpsKJCApZDyk3jIlsalAApJCpkNa3rT1usSmAglICpkO6UHff2JeFZCApJDhkPJ+kbY0walA\nApJChkP6H+mZ6FQgAUkhwyG1T/Cx7/VAApJKZkOaJR0SngokIClkNqSu8kLCU4EEJIWMhrQo\n7bhEnvddEZCApJDRkO6W0YlPBRKQFDIZUs4hB/h4F/M9AxKQFDIZ0njppzAVSEBSyGRIpzaa\npzAVSEBSyGBIr0gnjalAApJCBkO6RKZrTAUSkBQyF9Kc5F+qTAUSkBQyF9KN8oTKVCABSSFj\nIf3Q7LAclalAApJCxkIaLIN1pgIJSAqZCimndZNMnalAApJCpkKaKDcpTQUSkBQyFFLBaY3m\nKk0FEpAUMhTSS9JZayqQgKSQoZAukHe1pgIJSAqZCemTpHPVpgIJSAqZCamLPK82FUhAUshI\nSN+mnpD4X8ZWBSQgKWQkpFtlvN5UIAFJIRMhrWh+yFq9qUACkkImQhqa8MsURwckIClkIKR1\nR2o9O6g8IAFJIQMhTZIbNacCCUgKmQep4BS1ZweVByQgKWQepBckQ3UqkICkkHmQzpJZqlOB\nBCSFjIP0pnTUnQokIClkHKSO8qbuVCABSSHTIH2c1E55KpCApJBpkNIV3silZkACkkKGQZrb\n6CS9p6tWBCQgKWQYpB7yd+2pQAKSQmZBWph2dK72VCABSSGzIN0ij6pPBRKQFDIK0o/NDlX8\n+4nKgAQkhYyCNFQe1J8KJAG0RZcAABJSSURBVCApZBKknw9unhXAVCABKfFMgvSI/CmAqUAC\nkkIGQVrbupnmH/RVBSQgKWQQpIfl3h/dV8UdkICkkDmQsg9vvA5IAQwFkkbmQBolAxwgBTAU\nSBoZAynn8P1WAwlIiRdySKPlliIgASnxwg0p54j9FgIJSAqFG9J/y03rgQQkhUINaV2btIVA\nMgjSlnF9egzPr9jeMLbn1fdlOs6d6ZG6ASmO1CGNLXtRSCCZA2nEoBU5Y28vKd++e1DWukev\n2+H0nR6ZtAFIcaQNKfIDaT6QDIK0PiMr8lPpqgVl25tHrXGcgvTvna5za6wBknvakMbJDeuB\nZBCkL7uURj7e8XL1jmWdfypKn3jXDaOygRRHypDW/SI18gMJSOZAeu/6so+DJ1dd39xvqrOx\n1/jMzGG9tpb9wOoY6flStxzHdUmDyYhznSy3lF0EdK6BTDXpXD1M3RUvpL41IK29+cnSiq3t\n3WZGPm7IiPTyLrecUtclPiotCWKq4wQxtUT1Fth2dFpW2WWpAedaVakTyD9XIOda4n6uxXFC\nml1x1+6VimsLekyvPtLvxaot7tq5p3vXbnTl+7hw186Yu3YbMpY7zqbOi8uvLLl2XtnFqkkR\njju6fQQk76lCWnt444XlG0AyBpIzesCK7GH3lDoz33IKb3qpbMSOzT0m5GaP6rsTSN5ThTRE\n+ldsAMkcSNsm9L5uVOSzxgxxFqSXN8PJGnJNzxF51UuA5J4mpKyDmlf+YSyQzIHkISC5pwnp\nHhlUuQUkICkUUkiZ+7dcUbkJJCApFFJIt8rwqk0gAUmhcEJa1LjVmqptIAFJoXBC6iNjq7eB\nBCSFQglpflqbnOorQAKSQqGE1E2e2H0FSEBSKIyQPk8+KW/3NSABSaEwQrpMpkZdAxKQFAoh\npLekXUHUVSABSaHwQSpoLzOirwMJSAqFD9JT0qnGdSABSaHQQco5JuWLGjuABCSFQgfpocq/\n56sOSEBSKGyQfmi5/9Kae4AEJIXCBqmfDN5jD5CApFDIIM1PO3zNHruABCSFQgapi0zccxeQ\ngKRQuCB9nHxK3p77gAQkhcIF6UKZVmsfkICkUKggvSAdau8EEpAUChOknOMafVx7L5CApFCY\nIA2VPjH2AglICoUI0tIDDvwuxm4gAUmhEEG6RkbF2g0kICkUHkizktuui7UfSEBSKDSQCtrF\neOi7LCABSaHQQJooV8Y+ACQgKRQWSCtbpc2JfQRIQFIoLJBul7vrOAIkICkUEkhfpR2+qo5D\nQAKSQiGB9J/y97oOAQlICoUD0gvSvqCuY0ACkkKhgLSmTaOP6jwIJCApFApI/eWWug8CCUgK\nhQHS56mtsuo+CiQgKRQCSAXn1Xit7z0DEpAUCgGkCXJxfYeBBCSF7IeU2bLJvPqOAwlICtkP\nqbsMrfc4kICkkPWQ3ko6OeZfT1QHJCApZDuknLZJb9S/AkhAUsh2SH+Rni4rgAQkhSyHNKfx\nwZkuS4AEJIXshpR/jjzptgZIQFLIbkij5RLXNUACkkJWQ5rfrPm3rouABCSFbIZU0FEed18F\nJCApZDOkR+XCOv8KaXdAApJCFkNa2KJpvc8NqgxIQFLIYki/k7FelgEJSArZC2minOfhjh2Q\ngKSStZAWH9jkK08LgQQkhayFdJmM9LYQSEBSyFZIj8tZ+d5WAglIClkKad7+Tet4heJaAQlI\nCtkJKbe9TPS6FkhAUshOSAPreuuJGAEJSApZCemD1NZufzyxOyABSSEbIa0+Pul/vU8FEpAU\nshHSdfW9sGqtgAQkhSyE9E85aW0cU4EEJIXsg7TskLRP4pkKJCApZB2k/I4yPK6pQAKSQtZB\nGiwdPD6loTIgAUkh2yC9mXLo4vimAglIClkGaVnr5FfjnAokIClkF6TIL0j3xzsVSEBSyC5I\nA+WCvHinAglIClkF6fVGhy2JeyqQgKSQTZCWtkqZHv9UIAFJIYsg5XdweSek2AEJSApZBGmA\nXOLp1U72CEhAUsgeSFOT2nj/24mogAQkhayB9MUB+33gayqQgKSQLZC+PybpKX9TgQQkhSyB\nlHex9Pc5FUhAUsgSSDdLx7j/J7YyIAFJITsgPSHH/+B3KpCApJAVkD5svP/nvqcCCUgK2QBp\nyRHJz/ufCiQgKWQBpDW/kb8kMBVIQFLIfEh5l8vv/TyjoSogAUkh8yHdJL/NSWQqkOyC9KNb\nTrHrEh9tc//CPtpVEsTULdtj7HxITlye0NQiZ0NCnx+7jYUBDP1xh7MxgKkbgvnOcja7Lfk5\nAEiFxW45pa5LfFRSEsTUUieIqbtinOv/Jh/yXWJTAzrXYP61nF1BjN1X51oUACTu2rkX467d\n+00av5vgVO7a2XXXzvXEgFQb0rxDkv+Z6FQgAUkhoyF9d6w8kvBUIAFJIZMhrThDbk98KpCA\npJDBkNacI13je1HVmAEJSAqZC2ndJXLZOoWpQAKSQsZCyrtKLsjWmAokIClkKqSCXtJulcpU\nIAFJIVMh3S6nfK8zFUhAUshQSPfKsfG/pmrsgAQkhcyENFSO+kZrKpCApJCRkAZL6zlqU4EE\nJIVMhDRUjtRzBCQgaWQgpHvlqLmKU4EEJIXMg3SPtJmnORVIQFLIOEgDpM3XqlOBBCSFDINU\ncLMc/63uVCABSSGzIOV2l7Zxvmm5a0ACkkJGQfrxEjl9qfZUIAFJIZMgrTlXzstSnwokIClk\nEKRFp8sVa/XHAglICpkD6YujpHeu/lggAUkjYyC931LudnlXc38BCUgKmQLp6SaNHnV5M2af\nAQlICpkBqeAvSU2ecXtXc58BCUgKGQEp52ppPdPtXc39BiQgKWQCpMXt5PSyPz8CEpCA5L9P\n20j6mrINIAEJSL6b2jRpYMX7HwEJSEDyWd6ApMb/qNwGEpCA5K+l50ubD6quAAlIQPLV9NZy\nwbLqa0ACEpB8VDAyNal/1Kt7AwlIQIq/lRnSclr0DiABCUhx994x0r7mH8MCCUhAirPcQSlJ\nt+3xZhNAAhKQ4mv+b+XQF/fcCSQgASmuprSQKzNr7QUSkIAUR1ldpPHIGPuBBCQgee+ZVtLu\nq1gHgAQkIHltabqk/jn2n5QDCUhA8tiUlnL6h3UcAxKQgOSp+R2l8dC8uo4CCUhA8lDuiGZy\nwdy6jwMJSEBy7/VTpMWEgnoWAAlIQHJr0dVJSRkL6l0CJCABqf7WjdxfTp/hsghIQAJSvT1/\nnBw0ps4HGaoCEpCAVE/vnSPJfWo/I6hWQAISkOpsTnqSdPzYy0ogAQlIdZTZP01OfsHbWiAB\nCUgxy7y7mfzi7/nuC8sDEpCAFKPMAftLyxE5nqcCCUhAqlXmwObScmA878AHJCABaY+W3NFM\nDnlwVVxTgQQkINXoX9elyaHDV8c5FUhAAlJUb1ySJL8YvSbuqUACEpCqWjf5DJFf/4/r0xhi\nBCQgAamihfe2luRL3/Q3FUhAAlKkgtc6pcj+13/hdyqQgASk9ZkjTxQ5+ZEV/qcCCUhhh5T7\nQqc0Se3s8z5dZUACUrgh/euOViInDF6c4FQgASnEkBaO+LVI817vJD4VSEAKK6TMMecmS/KF\nT67VmAokIIUSUubjHVMk6cyRi5SmAglI4YP09YhzG4mcNmS+3lQgASlckArmPnCaSNIZg79U\nHQskIIUI0vf/6H6oSGqHR751XxtfQAJSSCDlzbz/rMgdupbXPLdcc2xlQAJSCCAVfD7q8hYi\nyWf86Z08nXc13zMgAclySPmfjflD5P6cHNn9yWVl14EEJCDFWc7bQ//rwAiigzs/OqdqH5CA\nBKQ4+nryre33iyA6/PejP41++XsgAQlI3lr28l8uPSRiKPmk3k98vedBIAEJSK598+y9lx4e\nMSSHXnb//8X8swggAQlI9ZT19ti+v21RZuiQjgOm1P2kBSABCUgxy5w+/taObcoIJR19xZ+f\nd/nPViABCUg1y/5s6tDrzjm4jJC0PO/G8e+u9DAVSEACUmXL3n3q/h7nH5VUTujIi25+9M3v\nPE8FEpDCDqlg0TuTH7jhv05pVg5IDjn72sFPfxLva9ABCUghhbRmzvR/jLil81lHpVYAanzi\nf/7xoX9+4uV+XIyABKQwQcpZ+OFrfx96a9fz2u5fwUeSWp1x+U0PTZ25NLHBQAKS5ZDyln35\n9nOThvfvcVn741tIVfufeG7XW0f84+353t9jpd6ABCTbIGUt+GzGi38fM+TOXukX/vIXu+2I\ntDzhrCtuuH/8M2/Pif9luF0CEpCCgLRlXJ8ew/NrbkfvU4KUu3zBnFlvTpsycfTQ/n/sntGh\nXdsjmkt0aYee2P7S7rfdP/bpNz9fklv+Ofv6rS/jCUhhhzRi0IqcsbeX1NiO3ucdUs7y5V/P\n+3TW269Om/KPcWOHDunfr3f3jEs6tPuPo1sdmCx71uKok8+8KKPX7X8Z/cRzb332bay3UQES\nkIyBtD4jK/IT6KoF0dvR+7xBWn3QgSm1qFTV7KCjT/rV+RdnXHtD/78MGzf5+ddmzVnk6cWB\ngQQkYyB92aU08vGOl6O3o/btyo6U/5Nb+ce2PeOMDh0u69y1z/V33TXkwXETpkx99fUPP56/\ncMU610+us+3b/H9u3ZWUBjF1644gphY5PwcwdVNRAEN/2ulsCmDqz7sCGPrTdmer25KNcUJ6\n7/qyj4MnR29H7ctvF+lpT5OILKr6FxuvkPqWfayEVLkdtW/ToEizdrrllLgu8dGu4iCmlpYG\nMbV4VxBTS5wgphYF86/lFAYwtTCYfy2nyHVNnJBmV9yNeyV6O3pfWfv6mQ268TsSvyMF8DvS\nhozlkZ87nRdHb0fvA5K3gBRySM7oASuyh91T6sx8a/d21SWQPAeksEPaNqH3daMinzVmyO7t\nqksgeQ5IYYfkISC5ByQgAUkhIAEJSAoBCUhAUghIQAKSQkACEpAUAhKQgKQQkIAEJIWABCQg\nKQQkIAFJISABCUgKAQlIQFIISEACkkJAAhKQFAISkICkEJCABCSFgAQkICkEJCABSSEgAQlI\nCgEJSEBSCEhAApJCQAISkBQCEpCApBCQgAQkhYAEJCApBCQgAUkhIAEJSAoBCUhAUghIQEq8\n0pHP7L0vlmhPTNjXZ+C9F0cW7utT8Ny7I/P29Sl4bvbIpd4X701I7fruvS+WaBmX7Osz8N7t\n7bbu61Pw3Oh23+3rU/Dcc+0+8L4YSLEDUjABKfGAFFBACiYgKQSkYAISEdUdkIgUAhKRQkAi\nUmjvQLozPVK3iu0t4/r0GJ6/V76svzaM7Xn1fZkV29En3uCKviUb+q1qzI3q+Pxu3TuQ+k5f\nv379hortEYNW5Iy9vWSvfF1f3T0oa92j1+0o344+8QZX9C3Z0G9VY25Ux+d3696B1HVu9eb6\njKyI86sW7JWv66fNo9Y4TkH69+VXok68wRV9Szb0W9WYG7UsX9+tewVSUfrEu24YlV2+/WWX\n0sjHO17eG1/Xf8s6lz8dMfrEG1zRt6QRt6oJN6rj97t1r0Da2Gt8ZuawXuX/bfje9WUfB0/e\nG1/Xd5v7TS2/jD7xBlf0LWnCrWrEjer4/W4NGtLnnTt3rngO7fZuM8su3it/ekPD/CevOtm1\nNz9Zuntv5Yk3uKJvyYZ8q1Zmxo1aVdzfrUFD2rZq1aqdFZv9Xiz7OLvih+UrAX9dX1We7IIe\n02vsrjjxBlf0LdmQb9WKDLlRq4v3u3Wv3LVbNanYcXZ0+6hse0PGcsfZ1Hnx3vi6/lpy7byq\nzegTb3BF35IN/lY15UZ1/H637hVIm3tMyM0e1XenM/Mtxxk9YEX2sHtK3T9rH1V400tlf/u4\no+xkq0+8QVZ1Sxpwq5pzo/r9bt07D39nDbmm54g8xxkzJHL/aULv60b95P45+6oF6eXNKD/Z\nqhNvkFXdkgbcqubcqI7P71aeIkSkEJCIFAISkUJAIlIISEQKAYlIISARKQQkIoWAZHQ7Tj2k\n/H82Sy9qsmxfn0u4A5LZLUhLL7uYJJP29ZmEPCAZ3qMyxXGyml3aYJ9mF5KAZHilv2u+qvSi\nljn7+jzCHpBML7vlRY/LtH19FqEPSMb3iiT32tfnQEAyvh+bS0P+E/OQBCTju+rAi1qs3tcn\nEfqAZHpPyXM5B17Mg3b7OCAZXmaz3zvOVHlsX59H2AOS2RW1O7QgcnEFT2zYxwHJ7AbJa2UX\n2S3aF+/rUwl3QDK6T5J7Vmw8LcP27ZmEPSARKQQkIoWARKQQkIgUAhKRQkAiUghIRAoBiUgh\nIBEpBCQihYBEpBCQiBT6f/7FasUso+yxAAAAAElFTkSuQmCC"
},
"metadata": {
"image/png": {
"width": 420,
"height": 420
}
}
}
]
},
{
"cell_type": "markdown",
"source": [
"ロジスティック関数によって、求めていた$0 \\sim 1$の範囲に留まる曲線が描けているのが分かるかと思います。\n",
"\n",
"パラメータ$\\beta_0, \\beta_1$の値によって形が変わってくるので、あとは得られたデータに合うように$\\beta_0, \\beta_1$を推定していきます。"
],
"metadata": {
"id": "ygzcwuYuNFNZ"
}
},
{
"cell_type": "markdown",
"source": [
"#### 最尤法 (most likelihood method)\n",
"\n",
"線形モデルでは、予測からのデータの残差が最小になるようにパラメータを推定する最小二乗法を用いていました。\n",
"\n",
"しかしこの方法は、残差が正規分布に従うことを仮定した方法だったので、今回は使用できません。\n",
"\n",
"そこで**最尤法**という手法でパラメータを推定します。\n",
"\n",
"最尤法は、手元の観測データを最も高い確率で再現できるモデルを探索してパラメータを決める方法になります。\n",
"\n",
"得られた観測データに対し最も尤もらしい(もっともらしい)モデルを求める、という意味になります。\n",
"\n",
"具体的に見ていくと、例えば今回下記の様なデータが得られていたとします。\n",
"\n",
"\\begin{array}{c|ccc} \\hline\n",
" 観測データ & 発芽種子数 & 未発芽種子数 & 発芽率 & 栄養成分 \\\\ \\hline\n",
" A & 1 & 9 & 0.1 & 1.411 \\\\\n",
" B & 6 & 4 & 0.6 & 5.229 \\\\\n",
" C & 9 & 1 & 0.9 & 9.45 \\\\ \\hline\n",
"\\end{array}\n",
"\n",
"栄養成分$1.411$の時に発芽する確率を$p_A$とおくと、観測データAが生じる確率は、二項分布より\n",
"\n",
"* $_1C_{10}\\times p_A^1\\times(1-p_A)^9$\n",
"\n",
"と計算できます。同様に、観測データBやCが生じる確率はそれぞれ\n",
"\n",
"* $_6C_{10}\\times p_B^6\\times(1-p_B)^4$\n",
"\n",
"* $_9C_{10}\\times p_C^9\\times(1-p_C)^1$\n",
"\n",
"と計算できます。\n",
"\n",
"先ほどのロジット関数の逆関数より、$p = \\dfrac{1}{1 + e^{-(\\beta_0+\\beta_1x)}}$となるので、各観測値における施肥量から\n",
"\n",
"$p_A = \\dfrac{1}{1 + e^{-(\\beta_0+\\beta_1 \\times 1.14)}}$、$p_B = \\dfrac{1}{1 + e^{-(\\beta_0+\\beta_1 \\times 5.229)}}$、$p_C = \\dfrac{1}{1 + e^{-(\\beta_0+\\beta_1 \\times 9.45)}}$\n",
"\n",
"となり、観測データA,B,Cが生じる確率は$\\beta_0, \\beta_1$の関数となることが分かります。\n",
"\n",
"観測データAが生じる確率$\\times$観測データBが生じる確率$\\times$観測データCが生じる確率が最大になる、つまり\n",
"\n",
"$p_A^1\\times(1-p_A)^9 \\times p_B^6\\times(1-p_B)^4 \\times p_C^9\\times(1-p_C)^1$\n",
"\n",
"が最大になる様な$\\beta_0, \\beta_1$を求めるのが最尤法となります。\n",
"\n",
"\n", "\n", "もう少し一般化すると、観測データが$n$得られている場合、\n", "\n", "各観測データにおける発芽した種子数を$y_i$, 発芽しなかった種子数を$10 - y_i$, 発芽する確率を$p_i$とすると、\n", "\n", "各観測データが生じる確率の積は\n", "\n", "$L = \\prod_{i=1}^{n}p_i^{y_i}(1-p_i)^{10 - y_i}$\n", "\n", "と表すことが出来ます。各観測値が起きうる確率の積$L$を**尤度**と呼びます。\n", "\n", "
\n", "\n", "この尤度が最大になるパラメータを求めることになりますが、\n", "\n", "積の数がとても多いので対数をとり、和の形に変換します。\n", "\n", "$log(L) = log(\\prod_{i=1}^{n}p_i^{y_i}(1-p_i)^{10 - y_i}) $\n", "\n", "$= \\sum_{i=1}^{n}\\lbrace y_i log(p_i)+(10-y_i)log(1-p_i)\\rbrace$\n", "\n", "この対数を取った尤度$log(L)$を**対数尤度**と呼びます。\n", "\n", "
\n", "\n", "$p = \\dfrac{1}{1 + e^{-(\\beta_0+\\beta_1x)}}$を対数尤度の関数に代入し、尤度(対数尤度)を最大にするパラメータ(今回は$\\beta_0, \\beta_1$)の値を計算します。\n", "\n", "計算によって算出されたパラメータ($\\hat{\\beta_0}, \\hat{\\beta_1}$)を**最尤推定量**と呼びます。\n", "\n", "最尤推定量は、標本サイズが大きくなるにつれて、真のパラメータに確率収束していく一致性という性質や、\n", "\n", "線形不偏推定量の中で分散が最も小さい推定量という性質を持ち、観測されたデータから元となる母集団を推定しようとする際に広く用いられる推定量になります。\n", "\n", "
\n", "\n", "最尤推定量は偏微分$\\dfrac{\\partial}{\\partial \\beta_0}log(L) = 0$や、$\\dfrac{\\partial}{\\partial \\beta_1}log(L) = 0$を解くことで得る事が出来ますが、\n", "\n", "(※厳密にはロジスティック関数の場合は最尤推定量が解析的に求められないので、ニュートンラフソン法等を使い数値的に求める必要があります。)\n", "\n", "これもRの関数を用いて簡単に算出することが可能です。\n", "\n", "\n", "\n" ], "metadata": { "id": "ppOfXtgQQAqg" } }, { "cell_type": "code", "source": [ "# データの表示\n", "head(data)" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 286 }, "id": "sxPR5cOJq9rp", "outputId": "ac9b69d8-2349-4897-f39d-b981cdaa9b45" }, "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/html": [ "
| germination | size | germination_rate | nutrition | |
|---|---|---|---|---|
| <int> | <int> | <dbl> | <dbl> | |
| 1 | 0 | 10 | 0.0 | 1.113 |
| 2 | 1 | 10 | 0.1 | 1.411 |
| 3 | 0 | 10 | 0.0 | 1.162 |
| 4 | 0 | 10 | 0.0 | 1.181 |
| 5 | 0 | 10 | 0.0 | 1.422 |
| 6 | 0 | 10 | 0.0 | 2.237 |
| germination | size | germination_rate | solar | nutrition | |
|---|---|---|---|---|---|
| <int> | <int> | <dbl> | <chr> | <dbl> | |
| 1 | 0 | 10 | 0 | shade | 1.213 |
| 2 | 0 | 10 | 0 | shade | 1.190 |
| 3 | 0 | 10 | 0 | shade | 1.199 |
| 4 | 0 | 10 | 0 | shade | 1.070 |
| 5 | 0 | 10 | 0 | shade | 1.425 |
| 6 | 0 | 10 | 0 | shade | 2.117 |
| germination | size | germination_rate | solar | nutrition | rain | temperature | wind | |
|---|---|---|---|---|---|---|---|---|
| <int> | <int> | <dbl> | <chr> | <dbl> | <int> | <int> | <int> | |
| 1 | 0 | 10 | 0 | shade | 1.213 | 5 | 20 | 5 |
| 2 | 0 | 10 | 0 | shade | 1.190 | 51 | 25 | 5 |
| 3 | 0 | 10 | 0 | shade | 1.199 | 61 | 15 | 1 |
| 4 | 0 | 10 | 0 | shade | 1.070 | 12 | 15 | 2 |
| 5 | 0 | 10 | 0 | shade | 1.425 | 99 | 20 | 4 |
| 6 | 0 | 10 | 0 | shade | 2.117 | 78 | 25 | 0 |
\n", "\n", "実際の解析では、この**良くないモデル**というものを定量的に評価する必要があります。\n", "\n", "例えば、今回得られたデータからは、\n", "\n", "* 1つ目のモデル: $logit(p) = log(\\dfrac{p}{1-p}) = \\beta_0+\\beta_1 \\times 栄養成分$\n", "\n", "* 2つ目のモデル: $logit(p) = log(\\dfrac{p}{1-p}) = \\beta_0+\\beta_1 \\times 栄養成分+\\beta_2 \\times 日照条件$\n", "\n", "* 3つ目のモデル: $logit(p) = log(\\dfrac{p}{1-p}) = \\beta_0+\\beta_1 \\times 栄養成分+\\beta_2 \\times 日照条件 + \\beta_3 \\times 降水量 + \\beta_4 \\times 気温 + \\beta_5 \\times 風速$\n", "\n", "* ...\n", "\n", "といった形で様々なモデルを考える事が出来ました。\n", "\n", "このうちどのモデルが良いモデルなのか選択する基準として**AIC** (**Akaike's Information Criterion**)と呼ばれる指標があります。\n", "\n" ], "metadata": { "id": "e66p1QQ9Vy1e" } }, { "cell_type": "markdown", "source": [ "まず、モデルの良さを考えるときに、統計モデルのパラメータを最尤推定した時の**対数尤度**$log(L)$というものがありました。\n", "\n", "対数尤度は「各観測データが起きうる確率の積」によって計算されたので、尤度が高ければ高いほど「\"観測データ\"へのあてはまりの良さ」が高いと考えられます。\n", "\n", "最尤法では、対数尤度が最大になる様なパラメータ$\\beta$が得られますが、この最大対数尤度を$log(L)^*$と表記しておきます。\n", "\n", "この時、\n", "\n", "$D = -2log(L)^*$\n", "\n", "で計算される値が**逸脱度**と定義されます。\n", "\n", "これは、値が大きいほどあてはまりが悪いことを示す値になります。\n", "\n", "
\n", "\n", "この逸脱度が小さいほど良いモデルになるのでは？という気がしてきますが、\n", "\n", "逸脱度はどんな変数であっても変数を加えれば加えるほど小さくなっていくという特徴があります。\n", "\n", "例えば右図の変数を大量に組み込んだ式の方がたしかに\"観測データ\"に対するあてはまりは良くなっていますが、その事象を説明する統計モデルとして正しいか？というと良くない可能性が高そうですね。\n", "\n", "
![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter9/overfit.png?raw=true\")
\n",
"\n",
"よって「ある程度変数の数を抑えつつ、逸脱度は小さい」様なモデルを選択する必要があります。\n",
"\n",
"この選択基準として、変数の数を$k$としたときに、\n",
"\n",
"$AIC = -2(log(L)^* - k) = D + 2k$\n",
"\n",
"と定義される**AIC**が良く用いられます。\n",
"\n",
"このAICが一番小さいモデルが良いモデルとして選択されます。\n",
"\n",
"AICは観測データへの当てはまりの良さだけでなく、モデルのパラメータ数でペナルティをつけてモデル選択することで、観測データに過剰に適合しないモデルを選択する事が出来る指標です。\n",
"\n",
"Rの`glm`関数では、最尤推定を実施する際に、このAICも同時に計算してくれています。"
],
"metadata": {
"id": "-A6gnUwSaP2-"
}
},
{
"cell_type": "code",
"source": [
"# GLMを実施し、AICを確認する\n",
"result <- glm(cbind(germination, size - germination) ~ nutrition, family = binomial, data = data)\n",
"summary(result)"
],
"metadata": {
"id": "nQnN7B8hlvl0",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 372
},
"outputId": "f51f1f4a-55c4-4c6d-a58d-a1ad3e97352c"
},
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
"\n",
"Call:\n",
"glm(formula = cbind(germination, size - germination) ~ nutrition, \n",
" family = binomial, data = data)\n",
"\n",
"Coefficients:\n",
" Estimate Std. Error z value Pr(>|z|) \n",
"(Intercept) -3.82415 0.25237 -15.15 <2e-16 ***\n",
"nutrition 0.44661 0.03369 13.26 <2e-16 ***\n",
"---\n",
"Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
"\n",
"(Dispersion parameter for binomial family taken to be 1)\n",
"\n",
" Null deviance: 656.62 on 99 degrees of freedom\n",
"Residual deviance: 413.64 on 98 degrees of freedom\n",
"AIC: 544\n",
"\n",
"Number of Fisher Scoring iterations: 5\n"
]
},
"metadata": {}
}
]
},
{
"cell_type": "markdown",
"source": [
"では先程提示した今回のデータから考えられる3つのモデル\n",
"\n",
"* 1つ目のモデル: $logit(p) = log(\\dfrac{p}{1-p}) = \\beta_0+\\beta_1 \\times 栄養成分$\n",
"\n",
"* 2つ目のモデル: $logit(p) = log(\\dfrac{p}{1-p}) = \\beta_0+\\beta_1 \\times 栄養成分+\\beta_2 \\times 日照条件$\n",
"\n",
"* 3つ目のモデル: $logit(p) = log(\\dfrac{p}{1-p}) = \\beta_0+\\beta_1 \\times 栄養成分+\\beta_2 \\times 日照条件 + \\beta_3 \\times 降水量 + \\beta_4 \\times 気温 + \\beta_5 \\times 風速$\n",
"\n",
"これらのAICをそれぞれ確認してみます。\n",
"\n",
"`glm`関数の結果を表示する際に、`summary(result)$aic`という形でAICのみを抽出する事が出来ます。"
],
"metadata": {
"id": "dd-5fwEVlwgn"
}
},
{
"cell_type": "code",
"source": [
"# 3つのモデルそれぞれのAICを計算する\n",
"# 1つ目のモデル\n",
"result <- glm(cbind(germination, size - germination) ~ nutrition, family = binomial, data = data)\n",
"summary(result)$aic\n",
"\n",
"# 2つ目のモデル\n",
"result <- glm(cbind(germination, size - germination) ~ nutrition + solar, family = binomial, data = data)\n",
"summary(result)$aic\n",
"\n",
"# 3つ目のモデル\n",
"result <- glm(cbind(germination, size - germination) ~ nutrition + solar + rain + temperature + wind, family = binomial, data = data)\n",
"summary(result)$aic"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 69
},
"id": "UXNJL9jFcjVB",
"outputId": "7b3db652-e180-49a0-90e4-9fea11943ce7"
},
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/html": [
"543.997350753646"
],
"text/markdown": "543.997350753646",
"text/latex": "543.997350753646",
"text/plain": [
"[1] 543.9974"
]
},
"metadata": {}
},
{
"output_type": "display_data",
"data": {
"text/html": [
"203.490122016901"
],
"text/markdown": "203.490122016901",
"text/latex": "203.490122016901",
"text/plain": [
"[1] 203.4901"
]
},
"metadata": {}
},
{
"output_type": "display_data",
"data": {
"text/html": [
"207.499829547914"
],
"text/markdown": "207.499829547914",
"text/latex": "207.499829547914",
"text/plain": [
"[1] 207.4998"
]
},
"metadata": {}
}
]
},
{
"cell_type": "markdown",
"source": [
"2つ目の栄養成分`nutrition`と日照条件`solar`のみ組み込んだモデルのAICが最も小さいことが分かったので、\n",
"\n",
"今回は2つ目のモデルが最も良いモデルとして選択できます。"
],
"metadata": {
"id": "d7TE1fAamG3L"
}
},
{
"cell_type": "markdown",
"source": [
"#### 他のモデル選択の指標\n",
"\n",
"他によく使用されるモデル選択のための情報量基準には、**BIC**(**Bayesian Informaiton Criterion**)、ベイズ情報量基準というものがあります。\n",
"\n",
"BICは最大対数尤度を$log(L)^*$、パラメータ数$k$,サンプルサイズ$n$とすると、\n",
"\n",
"$BIC = -2log(L)^* + klog(n)$\n",
"\n",
"で定義されます。\n",
"\n",
"AICと同じように、この値を最小化するモデルを良いモデルとして選択します。\n",
"\n",
"AICが新しく得られるデータの予測精度を高くするモデルを選ぶという指標であるのに対し、BICは候補となる複数の統計モデルに対し、真のモデルである確率が高いモデルが良いという指標のもと算出されます。\n",
"\n",
"\n", "\n", "AICやBICは他の生物学的な解析でも登場することがあります。\n", "\n", "例えば系統解析の際に、最適な塩基置換モデルをAICやBICで選択することがあります。" ], "metadata": { "id": "Jf1Wvz3hmdiD" } }, { "cell_type": "markdown", "source": [ "### 尤度比検定\n", "\n", "統計モデルを比較する方法として、AICを基準とした比較以外に、**尤度比検定**と呼ばれるものもあります。\n", "\n", "これは2つのモデルを比較し、データへの当てはまりがどれだけ改善されたかに基づいて検定を行う手法になります。\n", "\n", "例として、先ほどの観測データに対し下の2つのモデルを比較してみます。\n", "\n", "* 1つ目のモデル: $logit(p) = log(\\dfrac{p}{1-p}) = \\beta_0+\\beta_1 \\times 栄養成分$\n", "\n", "* 2つ目のモデル: $logit(p) = log(\\dfrac{p}{1-p}) = \\beta_0+\\beta_1 \\times 栄養成分+\\beta_2 \\times 日照条件$\n", "\n", "このモデルを帰無仮説の枠組みで考えると、\n", "\n", "帰無仮説: $H_0: \\beta_2 = 0$\n", "\n", "対立仮説: $H_1: \\beta_2 \\ne 0$\n", "\n", "と考える事が出来ます。\n", "\n", "この検定における統計量では、**尤度比**というものを扱います。\n", "\n", "各モデルの最大尤度を$L_1^*, L_2^*$とすると、尤度比は$\\dfrac{L_1^*}{L_2^*}$です。\n", "\n", "この尤度比の対数を取り$-2$をかけたもの、\n", "\n", "$𝛥D_{1,2}=-2\\times(log(L_1^*)-log(L_2^*))$\n", "\n", "を検定統計量として使用します。(**逸脱度の差**を計算していることになります。)\n", "\n", "この$𝛥D_{1,2}$は十分なサンプルサイズがある時、自由度$1$の$\\chi^2$分布に従うと近似できることが知られています。\n", "\n", "Rでは分散分析を実施するのに用いられる`anova`関数を使用することで検定を行う事が出来ます。\n" ], "metadata": { "id": "mGiHs_BFwKxq" } }, { "cell_type": "code", "source": [ "# anova関数を用いてモデルの尤度比検定を実施する\n", "# 1つ目のモデル\n", "model1 <- glm(cbind(germination, size - germination) ~ nutrition, family = binomial, data = data)\n", "\n", "# 2つ目のモデル\n", "model2 <- glm(cbind(germination, size - germination) ~ nutrition + solar, family = binomial, data = data)\n", "\n", "# 尤度比検定\n", "anova(model1, model2, test=\"LRT\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 160 }, "id": "35cXWaiL3ibY", "outputId": "be45eafc-fcca-41a7-c315-44dfbb0c414a" }, "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/html": [ "
| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) | |
|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 98 | 413.63965 | NA | NA | NA |
| 2 | 97 | 71.13242 | 1 | 342.5072 | 1.814875e-76 |
\n",
"\n",
"誤差項の従う分布に応じて、リンク関数も様々な関数に変わっていきます。\n",
"\n",
"\\begin{array}{cc} \\hline\n",
" 分布 & リンク関数 \\\\ \\hline\n",
" 正規分布 & identity: g(\\mu) = \\mu \\\\\n",
" 二項分布 & logit: g(\\mu) = log\\frac{\\mu}{1-\\mu} \\\\\n",
" ポアソン分布 & log: g(\\mu) = log\\mu \\\\ \\hline\n",
"\\end{array}\n",
"\n"
],
"metadata": {
"id": "I9Y7W_c4Bnst"
}
},
{
"cell_type": "markdown",
"source": [
"### ポアソン分布の場合\n",
"\n",
"では、別の一般化線形モデルの例として、ポアソン分布の場合を見ておきます。\n",
"\n",
"今回はある地域の一定区画$1m^2$にある植物の数とその土壌水分量のデータを使用します。"
],
"metadata": {
"id": "g0CZJ2knKmy4"
}
},
{
"cell_type": "code",
"source": [
"# 植物数のデータ読み込み\n",
"data <- read.csv(\"https://raw.githubusercontent.com/slt666666/biostatistics_text_wed/refs/heads/main/source/_static/data/chapter9_plant_number.csv\")\n",
"# データの一部を表示\n",
"head(data)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 286
},
"id": "Im7h73O3Li1B",
"outputId": "b556862f-869d-41bd-cdd0-10f6fc8a9dbf"
},
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/html": [
"| num | moisture | |
|---|---|---|
| <int> | <dbl> | |
| 1 | 1 | 13.0 |
| 2 | 0 | 10.7 |
| 3 | 0 | 8.1 |
| 4 | 0 | 5.2 |
| 5 | 5 | 22.9 |
| 6 | 1 | 5.7 |
![\"title\"](\"https://github.com/slt666666/biostatistics_text_wed/blob/main/source/_static/images/chapter9/fix_random.png?raw=true\")
\n",
"\n",
"例えば、ある作物を圃場Aと圃場Bで栽培した際の種子の発芽率の違いを調査するとします。\n",
"\n",
"この時、単純に考えると、最初にやった様に種子の発芽率$p$をロジット変換して、\n",
"\n",
"$logit(p) = \\beta_0 + 圃場の効果 \\times \\beta_1 + \\epsilon$\n",
"\n",
"というモデルで圃場の効果を考えられそうです。\n",
"\n",
"圃場Aの時は$\\beta_1 = 20$で圃場Bの時は$\\beta_1 = 50$といった形です。\n",
"\n",
"しかし実際のデータにおいては、圃場の真ん中あたりと圃場の端っこの方では環境が異なり、圃場の効果にもばらつきがあると考えられます。\n",
"\n",
"この様に個体ごとに変動すると考えられる効果のことを**変量効果**と呼びます。\n",
"\n",
"同じ個体を同じ条件・環境で栽培してもデータからは読み取れない何らかの原因でばらつきは生じるものです(**個体差**)。\n",
"\n",
"この個体差に関しても変量効果として考える事が出来ます。\n",
"\n",
"\n", "\n", "一方、ある作物に肥料を与えない時と、肥料を100g与えた時の種子の発芽数の違いを調査する場合、\n", "\n", "「肥料を与えない」「肥料を100g与える」という処理は全個体共通で行う事が出来ると考えられるので、\n", "\n", "肥料を100g与えた効果は$\\beta_1 = 50$といった形で表せそうです。\n", "\n", "この様な、効果を固定して考えてよいものを**固定効果**と呼びます。\n", "\n", "
\n", "\n", "この固定効果と変量効果、両者が混合しているモデルを、一般化線形**混合**モデル（Generalized Linear Mixed Model: GLMM）と呼びます。\n", "\n", "例えば、圃場の違い($S_1$)や肥料の有無($\\beta_1$)、個体差($S_0$)を組み込んだモデルを考えると\n", "\n", "$logit(p) = \\beta_0 + S_0 + X\\beta_{1} + ZS_{1} + \\epsilon$\n", "\n", "といったモデルになります。\n", "\n", "個体差$S_0$や圃場効果$S_1$は変量効果で、各個体毎に独立にばらついています。\n", "\n", "一方で、肥料の効果$\\beta_1$はどの個体でも共通の固定効果と考えています。\n", "\n", "
\n", "\n", "種子の発芽数は種子数$n$、発芽率$p$の二項分布$Bi(n, p)$に従います。\n", "\n", "しかし、$logit(p) = \\beta_0 + S_0 + X\\beta_{1} + ZS_{1} + \\epsilon$\n", "\n", "には変動効果が混ざっているので、$logit(p)$もばらつくことになります。\n", "\n", "簡単にするために、$logit(p)$が正規分布に従ってばらついていると仮定すると、$p$はロジット正規分布と呼ばれる分布に従います。(詳細は省く)\n", "\n", "つまり、種子の発芽数は二項分布とロジット正規分布という２つの分布が混ざった分布に従うことになります。\n", "\n", "GLMMではこの様な確率分布のパラメータをデータから最尤推定することになります。\n" ], "metadata": { "id": "wOtRtCRzkFAp" } }, { "cell_type": "markdown", "source": [ "### (参考) RでのGLMM\n", "\n", "細かい説明は省きますが、Rでは`lme4`や`glmmML`といったパッケージでGLMMのパラメータ推定が出来ます。\n", "\n", "今回は`lme4`を使用して見ます。Google Colabにはデフォルトでは`lme4`が無いので、インストールする必要があります。" ], "metadata": { "id": "13CnIQtdZR6K" } }, { "cell_type": "code", "source": [ "# lme4パッケージのインストール\n", "install.packages(\"lme4\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "4hcBSOfp7hny", "outputId": "4a19da18-919a-4b31-e9e8-c88837c975a1" }, "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stderr", "text": [ "Installing package into ‘/usr/local/lib/R/site-library’\n", "(as ‘lib’ is unspecified)\n", "\n", "also installing the dependencies ‘rbibutils’, ‘minqa’, ‘nloptr’, ‘reformulas’, ‘Rdpack’, ‘RcppEigen’\n", "\n", "\n" ] } ] }, { "cell_type": "markdown", "source": [ "サンプルデータとして種子の発芽データを使用します。\n", "\n", "* 個体`plant`:1~10\n", "* 施肥量`fertilizer`:100~500\n", "\n", "のそれぞれの条件で種子が発芽したかどうかの反復を4回とっています。\n" ], "metadata": { "id": "6mzAbY0-Zilr" } }, { "cell_type": "code", "source": [ "# GLMM用のデータ読み込み\n", "data <- read.csv(\"https://raw.githubusercontent.com/slt666666/biostatistics_text_wed/refs/heads/main/source/_static/data/chapter9_GLMM.csv\")\n", "# データの一部を表示\n", "head(data)" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 286 }, "id": "wG4lJReCZg8W", "outputId": "7bb4e901-3363-4346-d195-be59af5438fb" }, "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/html": [ "
| germinate | replicate | plant | fertilizer | |
|---|---|---|---|---|
| <int> | <int> | <int> | <int> | |
| 1 | 0 | 1 | 1 | 100 |
| 2 | 0 | 2 | 1 | 100 |
| 3 | 0 | 3 | 1 | 100 |
| 4 | 0 | 4 | 1 | 100 |
| 5 | 0 | 1 | 2 | 100 |
| 6 | 1 | 2 | 2 | 100 |