class: center, middle, inverse, title-slide # R Club - Day 05 ##
https://opueco.github.io/rclub-slides/2019f-05/slides.html
### Kenji Sato ### 2019/6/5 --- ## 復習 - 仮説検定の考え方 - 二項検定 - t検定 --- ## 今日の範囲 - 奥村§3.5, 3.7, 3.8 - ただし 3.7, 3.8 は関連問題を紹介しているので各自で読む - スキップの項 - §3.6 はさんざんやったのでスキップ - 最尤法の話は必要に迫られたらやる --- **t検定の例** ある機械で生産された10個の製品の重量を測定したところ 101.1, 103.2, 102.1, 99.2, 100.5, 101.3, 99.7, 100.5, 98.9, 101.4 のような結果を得た。母平均は100gと考えてよいか。有意水準5%で検定せよ。 (東大統計教室『統計学入門』p.252)) ```r w <- c(101.1, 103.2, 102.1, 99.2, 100.5, 101.3, 99.7, 100.5, 98.9, 101.4) t.test(w, mu = 100) ``` ``` ## ## One Sample t-test ## ## data: w ## t = 1.8909, df = 9, p-value = 0.0912 ## alternative hypothesis: true mean is not equal to 100 ## 95 percent confidence interval: ## 99.8449 101.7351 ## sample estimates: ## mean of x ## 100.79 ``` --- ## 仮説を変えてみる (1/) 帰無仮説 `\(\mu = 100\)` のもとでは,p値は 0.09。5%水準では有意ではない。 `\(\mu = 101\)` にするとどうなる? ```r t.test(w, mu = 101) ``` ``` ## ## One Sample t-test ## ## data: w ## t = -0.50265, df = 9, p-value = 0.6273 ## alternative hypothesis: true mean is not equal to 101 ## 95 percent confidence interval: ## 99.8449 101.7351 ## sample estimates: ## mean of x ## 100.79 ``` --- ## 仮説を変えてみる (2/) `\(\mu = 102\)` にすると? ```r t.test(w, mu = 102) ``` ``` ## ## One Sample t-test ## ## data: w ## t = -2.8962, df = 9, p-value = 0.0177 ## alternative hypothesis: true mean is not equal to 102 ## 95 percent confidence interval: ## 99.8449 101.7351 ## sample estimates: ## mean of x ## 100.79 ``` p値は `\(0.02\)`。データは `\(\mu = 102\)` をサポートしないようだ。 --- ## 信頼区間 confidence interval この数字が変化していないことにお気づきだろうか? ``` ## 95 percent confidence interval: ## 99.8449 101.7351 ``` 95% 信頼区間(confidence interval)というのは,5%水準で帰無仮説が棄却されないような帰無仮説の範囲である。したがって, .box2[ 信頼区間の中におさまるように仮説を選べば p値は有意水準を下回らない。逆に,信頼区間の外側に仮説を選べば p値は有意水準を上回る。 ] --- ## 信頼区間に関する注意 「真の母平均が 95% の確率で信頼区間の中に入る」ということでは**ない**。 通常,95%信頼区間の意味は次のように説明される。 1. 信頼区間は確率変数であり,データを観測するたびに信頼区間が変わる。 2. データを100回観測して 100通りの信頼区間を作ったら,95通りは真の未知母数を含む。 動くのが信頼区間で,母数は動かないというのをイメージできればOK。図を使って見てみよう。 --- **真の母数(未知)** `\(\mu\)` を持つ未知の分布にしたがって生成されたデータを持っている。その **標本平均(固定値)** は `\(\bar{x}\)` <img src="slides_files/figure-html/ci-image-2-1.png" height="650" style="display: block; margin: auto;" /> --- データから `\(\mu\)` について何がわかるか? <img src="slides_files/figure-html/ci-image-1-1.png" height="650" style="display: block; margin: auto;" /> --- 仮説検定の考え方を使う。色付き部分は 5%。 `\(\mu = \mu_0\)` は棄却 <img src="slides_files/figure-html/ci-image-3-1.png" height="650" style="display: block; margin: auto;" /> --- `\(\mu_0\)` を増やしてみる。これも棄却。 <img src="slides_files/figure-html/ci-image-4-1.png" height="650" style="display: block; margin: auto;" /> --- さらに `\(\mu_0\)` を増やしてみる。これはぎりぎり棄却されない。 <img src="slides_files/figure-html/ci-image-5-1.png" height="650" style="display: block; margin: auto;" /> --- さらに `\(\mu_0\)` を増やしてみる。棄却されない。 <img src="slides_files/figure-html/ci-image-6-1.png" height="650" style="display: block; margin: auto;" /> --- さらに `\(\mu_0\)` を増やしてみる。やっぱり棄却されない。 <img src="slides_files/figure-html/ci-image-7-1.png" height="650" style="display: block; margin: auto;" /> --- さらに `\(\mu_0\)` を増やしてみる。やっぱり棄却されない。 <img src="slides_files/figure-html/ci-image-8-1.png" height="650" style="display: block; margin: auto;" /> --- さらに `\(\mu_0\)` を増やしてみる。やっぱり棄却されない。 <img src="slides_files/figure-html/ci-image-9-1.png" height="650" style="display: block; margin: auto;" /> --- さらに `\(\mu_0\)` を増やしてみる。やっぱり棄却されない。 <img src="slides_files/figure-html/ci-image-10-1.png" height="650" style="display: block; margin: auto;" /> --- さらに `\(\mu_0\)` を増やしてみる。やっぱり棄却されない。 <img src="slides_files/figure-html/ci-image-11-1.png" height="650" style="display: block; margin: auto;" /> --- さらに `\(\mu_0\)` を増やしてみる。やっぱり棄却されない。 <img src="slides_files/figure-html/ci-image-12-1.png" height="650" style="display: block; margin: auto;" /> --- 帰無仮説の `\(\mu_0\)` が大きくなってここまで来ると棄却される。 <img src="slides_files/figure-html/ci-image-13-1.png" height="650" style="display: block; margin: auto;" /> --- 赤点の入る区間が信頼区間 <img src="slides_files/figure-html/ci-image-14-1.png" height="650" style="display: block; margin: auto;" /> --- ## 信頼区間の計算方法 さきほどの図からわかったこと - 信頼区間の左端は,その点を仮説に選べば,データがちょうど下側 `\(100(1-\alpha/2)\)`% 点にくるような点。あるいは上側 p値が `\(\alpha/2\)` になる点 - 信頼区間の右端は,その点を仮説に選べば,データがちょうど下側 `\(100(\alpha/2)\)`% 点にくるような点。あるいは下側p値が `\(\alpha/2\)` になる点 --- ## 信頼区間の計算 (1/) 教科書的な計算方法を確認しておく。 $$ t = \frac{\bar{X} - \mu_0}{\sqrt{s^2 / n}} $$ は自由度 `\(n-1\)` の t-分布に従う。 この t-分布の下側 `\(100\alpha\)`%点を `\(t_{n-1}^\alpha\)` と書くとすれば `$$\mathrm{Prob}\left\{ t_{n-1}^{\alpha / 2} \le \frac{\bar{X} - \mu_0}{\sqrt{s^2 / n}} \le t_{n-1}^{1 - \alpha/2} \right\} = 1-\alpha$$` `\(\mu_0\)` がこの範囲に入っていれば帰無仮説は棄却されない。 --- ## 信頼区間の計算 (2/) `$$\mathrm{Prob}\left\{ t_{n-1}^{\alpha / 2} \sqrt{\frac{s^2}{n}} \le \bar{X} - \mu_0 \le t_{n-1}^{1 - \alpha/2} \sqrt{\frac{s^2}{n}} \right\} = 1-\alpha$$` ここから `\(\mu_0\)` の範囲を調べる `$$\mathrm{Prob}\left\{ \bar{X} - t_{n-1}^{1 - \alpha / 2} \sqrt{\frac{s^2}{n}} \le \mu_0 \le \bar {X} - t_{n-1}^{\alpha/2} \sqrt{\frac{s^2}{n}} \right\} = 1-\alpha$$` したがって,信頼区間は `$$\left[ \bar{X} - t_{n-1}^{1 - \alpha / 2}\sqrt{\frac{s^2}{n}}, \bar {X} - t_{n-1}^{\alpha/2} \sqrt{\frac{s^2}{n}} \right]$$` --- ## 計算練習 `t.test()` を使わずに,95% 信頼区間を計算しなさい。 ```r w <- c(101.1, 103.2, 102.1, 99.2, 100.5, 101.3, 99.7, 100.5, 98.9, 101.4) ``` -- ```r > mean(w) - sqrt(var(w) / length(w)) * qt(0.975, length(w) - 1) ``` ``` ## [1] 99.8449 ``` ```r > mean(w) - sqrt(var(w) / length(w)) * qt(0.025, length(w) - 1) ``` ``` ## [1] 101.7351 ``` --- ## 比較 ```r > t.test(w) ``` ``` ## ## One Sample t-test ## ## data: w ## t = 241.25, df = 9, p-value < 2.2e-16 ## alternative hypothesis: true mean is not equal to 0 ## 95 percent confidence interval: ## 99.8449 101.7351 ## sample estimates: ## mean of x ## 100.79 ``` --- ## ここまでのまとめ 帰無仮説を動かして,得られたデータのもとで帰無仮説が棄却されない領域を探す。これが信頼区間。 信頼区間の下限と上限では p値がちょうど, `\(1 -\)` 信頼水準 になっている。 つまり, p値がちょうど 5% になるような点を(2つ)見つけられれば,それが信頼区間の両端になっている。 -- 離散確率変数の場合はp値がジャンプするのでうまくいかない。 --- ## 2項検定を使った実験 (1/) 「硬貨を10回投げて表が4回出た」という観測を使って,「表の出る確率は 0.5」という仮説を棄却できるか? ```r > binom.test(4, 10, 0.5) ``` ``` ## ## Exact binomial test ## ## data: 4 and 10 ## number of successes = 4, number of trials = 10, p-value ## = 0.7539 ## alternative hypothesis: true probability of success is not equal to 0.5 ## 95 percent confidence interval: ## 0.1215523 0.7376219 ## sample estimates: ## probability of success ## 0.4 ``` --- ## 2項検定を使った実験 (2/) 「硬貨を10回投げて表が4回出た」という観測を使って,「表の出る確率は <span style="color:red;">0.6</span>」という仮説を棄却できるか? ```r > binom.test(4, 10, 0.6) ``` ``` ## ## Exact binomial test ## ## data: 4 and 10 ## number of successes = 4, number of trials = 10, p-value ## = 0.2126 ## alternative hypothesis: true probability of success is not equal to 0.6 ## 95 percent confidence interval: ## 0.1215523 0.7376219 ## sample estimates: ## probability of success ## 0.4 ``` --- ## 2項検定を使った実験 (3/) 「硬貨を10回投げて表が4回出た」という観測を使って,「表の出る確率は <span style="color:red;">0.7</span>」という仮説を棄却できるか? ```r > binom.test(4, 10, 0.7) ``` ``` ## ## Exact binomial test ## ## data: 4 and 10 ## number of successes = 4, number of trials = 10, p-value ## = 0.0756 ## alternative hypothesis: true probability of success is not equal to 0.7 ## 95 percent confidence interval: ## 0.1215523 0.7376219 ## sample estimates: ## probability of success ## 0.4 ``` --- ## 2項検定を使った実験 (4/) 2項検定の場合でも信頼区間は変化していないことを確認 ``` ## 95 percent confidence interval: ## 0.1215523 0.7376219 ``` 「仮説が 0.737... を上回れば p値が 0.5 を下回る」はず ```r > binom.test(4, 10, 0.73)$p.value ``` ``` ## [1] 0.02872244 ``` あれ?信頼区間の範囲内で帰無仮説を選んだのに, p値が 5% を下回っている。 --- ## 離散確率変数に関する注意 離散確率変数の場合,帰無仮説の少しの変化に対して(両側)p値が大きくジャンプすることがある。 - 帰無仮説を動かすと「まれな現象」の範囲が変わる。 先週は `\(\theta = 0.5\)` という対象な分布をみていたからあまり意識しなかった。 - まれな現象として新しい事象が加わると有意確率が大きく増える --- ## 問題の概要 10回硬貨を投げて2回表が観測されたとしよう。 `\(\theta = 0.45\)` であれば,平均的に10回中4.5回表が出る。だから,2回表がでることよりも8回表が出ることの方がめずらしい。 `\(\theta = 0.5\)` になると,2回表が出ることと8回表が出ることがちょうど同じくらい珍しい。だからどちらの場合も8回表が出る確率を含める。 しかし, `\(\theta=0.51\)` になると,8回表が出ることは 2回表がでることほど珍しくなくなってしまう。だから,8回表が出る確率の分だけドーンと p値が減ってしまう。 --- ## 両側p値の変化 0.45 → 0.5 の変化はバーの高さが変わるだけなので,p値はなめらかに変化する。0.5 → 0.51 の変化はバーが1つすっぽり抜けてしまうので不連続だ。 <img src="slides_files/figure-html/binom-graph4-1.png" width="700" style="display: block; margin: auto;" /> --- ## 解決策:片側p値を使う 平均(帰無仮説)を中心に, 観測されたデータを裏返して確率を計算するからこんなことが起こる。 - データと同じかより多くの表が出る確率(上側p値) - データと同じかより少なく表が出る確率(下側p値) だけを使えばこんなことは起こらない。 `binom.test()` の信頼区間も片側p値に基づいて計算されている。 --- **上側p値** <img src="slides_files/figure-html/binom-graph5-1.png" width="450" style="display: block; margin: auto;" /> **下側p値** <img src="slides_files/figure-html/binom-graph6-1.png" width="450" style="display: block; margin: auto;" /> --- ## `binom.test()` の信頼区間 (1/) `binom.test()` で出力される信頼区間を自力で計算してみよう。 ```r > binom.test(4, 10) ``` ``` ## ## Exact binomial test ## ## data: 4 and 10 ## number of successes = 4, number of trials = 10, p-value ## = 0.7539 ## alternative hypothesis: true probability of success is not equal to 0.5 ## 95 percent confidence interval: ## 0.1215523 0.7376219 ## sample estimates: ## probability of success ## 0.4 ``` --- ## `binom.test()` の信頼区間 (2/) 0 から 4 の値を取る確率 `sum(dbinom(0:4, 10, x))` が 2.5% になる `x` が95%信頼区間の上端 4 から 10 の値を取る確率 `sum(dbinom(4:10, 10, x))` が 2.5% になる `x` が95%信頼区間の下端 --- ## `binom.test()` の信頼区間 (3/) ```r x <- 0:100 / 100 p_up <- sapply(x, function(t) sum(dbinom(4:10, 10, t))) p_lo <- sapply(x, function(t) sum(dbinom(0:4, 10, t))) plot(x, p_up, type = "l", xlab = "θ", ylab = "p") lines(x, p_lo) abline(h = 0.025, col = "red") polygon(c(0.1215523, 0.7376219, 0.7376219, 0.1215523), c(-0.1, -0.1, 0.025, 0.025), col=rgb(1,0,0,0.1), lty=0) ``` --- ## `binom.test()` の信頼区間 (4/) <img src="slides_files/figure-html/unnamed-chunk-12-1.png" height="500" style="display: block; margin: auto;" /> --- ## uniroot 方程式 `\(f(x) = 0\)` の解を範囲 `\((a, b)\)` の中で見つける関数。 `\(f(a)\)` と `\(f(b)\)` の符号が逆でないとだめ。 ```r > quadratic <- function(x) x ^ 2 - 1 > uniroot(quadratic, c(0, 2), tol = 1e-6) ``` ``` *## $root *## [1] 1 ## ## $f.root ## [1] 4.643496e-07 ## ## $iter ## [1] 7 ## ## $init.it ## [1] NA ## ## $estim.prec ## [1] 5e-07 ``` --- ## `$$\sum_{r=4}^{10} {}_{10}C_r \theta^r (1-\theta)^{10-r} = 0.025$$` ```r uniroot(function(t) sum(dbinom(4:10, 10, t)) - 0.025, interval=c(0.1,0.8), tol=1e-8) ``` ``` ## $root ## [1] 0.1215523 ## ## $f.root ## [1] -1.580923e-12 ## ## $iter ## [1] 13 ## ## $init.it ## [1] NA ## ## $estim.prec ## [1] 5e-09 ``` --- ## `$$\sum_{r=0}^{4} {}_{10}C_r \theta^r (1-\theta)^{10-r} = 0.025$$` ```r uniroot(function(t) sum(dbinom(0:4, 10, t)) - 0.025, interval=c(0.1,0.8), tol=1e-8) ``` ``` ## $root ## [1] 0.7376219 ## ## $f.root ## [1] 4.170969e-14 ## ## $iter ## [1] 12 ## ## $init.it ## [1] NA ## ## $estim.prec ## [1] 5e-09 ``` --- **2項検定の例題: 統計検定2級('18/11)** .pull-left[ 次の表は, 2017年に実施された, JR北海道の利用状況に関する調査結果である。調査は北海道の18歳以上の男女2067人に行われ, 回答数は 1338人であった。この調査結果は, 母集団を北海道の 18歳以上の男女とし, 標本サイズ 1338の単純無作為抽出に基づくものとみなす。 ほぼ毎日利用した人の割合の95% 信頼区間を求めよ。 ] .pull-right[ | 利用頻度 | 割合 | |--------------------|-------| | ほぼ毎日 | 2.0% | | 週に数回程度 | 2.5% | | 月に数回程度 | 9.0% | | 年に数回程度 | 26.8% | | ほとんど利用しない | 58.9% | | わからない,無回答 | 0.8% | ] --- ```r > binom.test(floor(1338 * 0.02), 1338, 0.02)$conf.int ``` ``` ## [1] 0.01273191 0.02834346 ## attr(,"conf.level") ## [1] 0.95 ``` ```r > binom.test(ceiling(1338 * 0.02), 1338, 0.02)$conf.int ``` ``` ## [1] 0.01333942 0.02922492 ## attr(,"conf.level") ## [1] 0.95 ``` 正規近似を使って次のように近似的に計算することもできる。 ```r > 0.02 + qnorm(0.025) * sqrt(0.02 * 0.98 / 1338) ``` ``` ## [1] 0.0124985 ``` ```r > 0.02 + qnorm(0.975) * sqrt(0.02 * 0.98 / 1338) ``` ``` ## [1] 0.0275015 ``` --- ## 内閣支持率(NHK世論調査) NHKは今月10日から3日間、全国の18歳以上の男女を対象にコンピューターで無作為に発生させた固定電話と携帯電話の番号に電話をかける「RDD」という方法で世論調査を行いました。 調査の対象となったのは2333人で、54%に当たる1260人から回答を得ました。 それによりますと、安倍内閣を「支持する」と答えた人は、先月の調査より1ポイント上がって48%だったのに対し、「支持しない」と答えた人は先月より3ポイント下がって32%でした。 <http://www.nhk.or.jp/senkyo/shijiritsu/> ('19/6/1 閲覧) --- ## 支持率は上昇したか? 標本サイズ1260のランダムサンプルとして信頼区間を求める。 -- ```r > binom.test(605, 1260, p = 0.47) ``` ``` ## ## Exact binomial test ## ## data: 605 and 1260 ## number of successes = 605, number of trials = 1260, ## p-value = 0.4805 ## alternative hypothesis: true probability of success is not equal to 0.47 ## 95 percent confidence interval: ## 0.4522495 0.5081607 ## sample estimates: ## probability of success ## 0.4801587 ``` --- ## 不支持率は低下したか? 標本サイズ1260のランダムサンプルとして信頼区間を求める。 -- ```r > binom.test(403, 1260, p = 0.35) ``` ``` ## ## Exact binomial test ## ## data: 403 and 1260 ## number of successes = 403, number of trials = 1260, ## p-value = 0.0248 ## alternative hypothesis: true probability of success is not equal to 0.35 ## 95 percent confidence interval: ## 0.2941364 0.3463896 ## sample estimates: ## probability of success ## 0.3198413 ``` --- ## 練習してみよう ```r devtools::install_github("opueco/R4FunDrill") R4FunDrill::start("Day05") ``` チュートリアルがはじまるよ。