```{r global_options, include=FALSE}
knitr::opts_chunk$set(eval = TRUE, results = FALSE, fig.show = "hide", message = FALSE)
set.seed(1234)
```

このラボでは、データのランダムサンプルからの統計量が、母集団パラメータの点推定値としてどのように機能するかを調査します。推定値の分布などの特性について学習するために、推定値の*標本分布*を定式化することに関心があります。

::: {#boxedtext}
**シードの設定:** このラボではいくつかのランダムサンプルを取り、標本分布を構築します。これは、ラボの開始時にシードを設定する必要があることを意味します。この概念が初めての場合は、確率に関するラボをレビューしてください。
:::

## はじめに

### パッケージの読み込み

このラボでは、**tidyverse** スイートのパッケージを使用してデータを探索および視覚化します。また、リサンプリングには **infer** パッケージを使用します。

パッケージを読み込みましょう。

```{r load-packages, message=FALSE}
library(tidyverse)
library(openintro)
library(infer)
```

### 再現可能なラボレポートの作成

新しいラボレポートを作成するには、RStudioで「New File」->「R Markdown...」に移動し、「From Template」を選択してから、テンプレートのリストから「Lab Report for OpenIntro Statistics Labs」を選択します。

### データ

2019年のギャラップ報告書は次のように述べています。

> 科学の進歩が人々に恩恵をもたらすという前提は、ワクチン開発から過去数十年の技術の爆発的進歩に至るまで、時代を超えた発見に具現化されてきました。その結果、何十億ものスーパーコンピューターが現在、世界中の人々の手やポケットに収まっています。しかし、世界中の誰もが科学が個人的に恩恵をもたらすと感じているわけではありません。
>
> **出典:** [World Science Day: Is Knowledge Power?](https://news.gallup.com/opinion/gallup/268121/world-science-day-knowledge-power.aspx)

ウェルカム・グローバル・モニターの調査では、世界中の人々の20%が、科学者の仕事が自分たちのような人々に恩恵をもたらすとは信じていないことが分かりました。このラボでは、この20%を真の母集団の割合と仮定し、母集団から小さなサンプルを取ることで、サンプル比率がサンプルごとにどのように異なるかを学びます。まず、母集団サイズを100,000と仮定して、母集団を作成します。これは、母集団の20,000人（20%）が科学者の仕事が個人的に恩恵をもたらさないと考えており、残りの80,000人がそう考えていることを意味します。

```{r}
global_monitor <- tibble(
  scientist_work = c(rep("Benefits", 80000), rep("Doesn't benefit", 20000))
)
```

データフレームの名前は`global_monitor`で、「科学者の仕事はあなたのような人々に恩恵をもたらすと思いますか？」という質問への回答を含む変数の名前は`scientist_work`です。

棒グラフを使用して、これらの回答の分布を素早く視覚化できます。

```{r bar-plot-pop, fig.height=2.5, fig.width=10}
ggplot(global_monitor, aes(x = scientist_work)) +
  geom_bar() +
  labs(
    x = "", y = "",
    title = "Do you believe that the work scientists do benefit people like you?"
  ) +
  coord_flip()
```

データフレームが正しく構築されたことを確認するために、要約統計量も取得できます。

```{r summ-stat-pop, results = TRUE}
global_monitor |>
  count(scientist_work) |>
  mutate(p = n /sum(n))
```

## 未知の標本分布

このラボでは、母集団全体にアクセスできますが、これは実生活ではめったにありません。母集団全体に関する情報を収集することは、非常に費用がかかるか、不可能であることがよくあります。このため、私たちは母集団のサンプルを採取し、それを使用して母集団の特性を理解することがよくあります。

科学者の仕事が自分たちに恩恵をもたらさないと考える人々の割合を推定することに関心がある場合、`sample_n`コマンドを使用して母集団を調査できます。

```{r samp1}
samp1 <- global_monitor |>
  sample_n(50)
```

このコマンドは、`global_monitor`データセットからサイズ50の単純ランダムサンプルを収集し、その結果を`samp1`に割り当てます。これは、母集団全員の名前が書かれた帽子からランダムに名前を引くのと似ています。これらの50の名前で作業することは、母集団の100,000人全員で作業するよりもはるかに簡単です。

1.  このサンプルにおける回答の分布を記述してください。母集団における回答の分布とどのように比較されますか？ **ヒント:** `sample_n`関数はデータセットからランダムな観測値（つまり行）をサンプリングしますが、データセット内の変数を同じ名前で参照することはできます。母集団データを視覚化および要約するために以前提示したコードはサンプルにも役立ちますが、今では母集団パラメータではなく標本統計量を計算しているため、割合を`p`とラベル付けしないように注意してください。統計量のラベルをカスタマイズして、それがサンプルから来ていることを示すことができます。

科学者の仕事が自分たちに恩恵をもたらさないと考えるすべての人々の割合を推定することに関心があるが、母集団データにアクセスできない場合、最良の単一の推測は標本比率です。

```{r phat-samp1}
samp1 |>
  count(scientist_work) |>
  mutate(p_hat = n /sum(n))
```

```{r inline-calc, include=FALSE}
# For use inline below
samp1_p_hat <- samp1 |>
  count(scientist_work) |>
  mutate(p_hat = n /sum(n)) |>
  filter(scientist_work == "Doesn't benefit") |>
  pull(p_hat) |>
  round(2)
```

あなたが選んだ50人によって、あなたの推定値は真の母集団比率`r samp1_p_hat`よりわずかに上または下になる可能性があります。しかし、一般的に、標本比率は真の母集団比率のかなり良い推定値であり、母集団の1%未満をサンプリングするだけでそれを得ることができました。

2.  標本比率が他の学生のサンプルの標本比率と一致すると予想しますか？その理由は何ですか？もし答えが「いいえ」の場合、比率は多少異なりますか、それとも大きく異なりますか？学生のチームにあなたの答えを確認させてください。

3.  サイズ50の2番目のサンプルを取り、`samp2`と名付けてください。`samp2`の標本比率は`samp1`のそれとどのように比較されますか？サイズ100とサイズ1000の2つのサンプルをさらに取ったとします。どちらが母集団比率のより正確な推定値を提供すると考えますか？

当然のことながら、別のランダムサンプルを取るたびに、異なる標本比率が得られる可能性があります。この方法で母集団平均を推定する際に、どれくらいの変動性を期待すべきかを知ることは有用です。標本比率の分布は、*標本分布（比率の）*と呼ばれ、この変動性を理解するのに役立ちます。このラボでは、母集団にアクセスできるため、上記の手順を何度も繰り返すことで、標本比率の標本分布を構築できます。ここでは、Rを使用して母集団からサイズ50の異なる15,000個のサンプルを採取し、各サンプル内の応答の割合を計算し、「恩恵を受けない」応答のみをフィルタリングし、各結果を`sample_props50`というベクトルに格納します。標本分布は置換サンプリングによって構築されるため、`replace = TRUE`を指定していることに注意してください。

```{r iterate}
sample_props50 <- global_monitor |>
  rep_sample_n(size = 50, reps = 15000, replace = TRUE) |>
  count(scientist_work) |>
  mutate(p_hat = n / sum(n)) |>
  filter(scientist_work == "Doesn't benefit")
```

そして、ヒストグラムでこれらの割合の分布を視覚化できます。

```{r fig.show="hide"}
ggplot(data = sample_props50, aes(x = p_hat)) +
  geom_histogram(binwidth = 0.02) +
  labs(
    x = "p_hat (Doesn't benefit)",
    title = "Sampling distribution of p_hat",
    subtitle = "Sample size = 50, Number of samples = 15000"
  )
```

次に、このコードのセットがどのように機能するかを確認します。

4.  `sample_props50`にはいくつの要素がありますか？標本分布を記述し、特にその中心に注目してください。回答には必ず分布のプロットを含めてください。

## 間奏：標本分布

`rep_sample_n`関数の背後にある考え方は*反復*です。以前、あなたは母集団内のすべての人々からサイズ`n`（50）の単一のサンプルを取りました。この新しい関数を使用すると、このサンプリング手順を`rep`回繰り返して、一連の標本統計量の分布を構築できます。これを**標本分布**と呼びます。

実際には、真の標本分布を構築することはめったにありません。なぜなら、母集団全体からのデータにアクセスできることはめったにないからです。

`rep_sample_n`関数がなければ、これは大変な作業になるでしょう。以下のコードを15,000回手動で実行する必要がありました。

```{r sample-code}
global_monitor |>
  sample_n(size = 50, replace = TRUE) |>
  count(scientist_work) |>
  mutate(p_hat = n /sum(n)) |>
  filter(scientist_work == "Doesn't benefit")
```

さらに、結果として得られる標本比率を毎回別々のベクトルに保存します。

15,000回それぞれ割合を計算したとき、**異なる**サンプルから計算したことに注意してください！

5.  標本分布がどのように構築されるか、そして`rep_sample_n`関数が正確に何をするかを理解するために、コードを変更して**サイズ10のサンプル**から**25個の標本比率**の標本分布を作成し、それらを`sample_props_small`というデータフレームに入れてみてください。出力を表示してください。`sample_props_small`というオブジェクトにはいくつの観測値がありますか？各観測値は何を表していますか？

## 標本サイズと標本分布

メカニズムはさておき、`rep_sample_n`関数を使用した理由に戻りましょう。標本分布、特に50人のサンプルからの割合の標本分布を計算するためです。

```{r hist, fig.show='hide'}
ggplot(data = sample_props50, aes(x = p_hat)) +
  geom_histogram(binwidth = 0.02)
```

あなたが計算した標本分布は、科学者の仕事が自分たちに恩恵をもたらさないと考える人々の真の割合を推定することについて多くのことを教えてくれます。標本比率が不偏推定量であるため、標本分布は真の母集団比率を中心とし、分布の広がりは、母集団から一度に50人しかサンプリングしない場合に生じる変動性の大きさを表します。

このセクションの残りの部分では、標本サイズが標本分布に与える影響を理解することに取り組みます。

6.  以下のアプリを使用して、サイズ10、50、100のサンプルから「恩恵を受けない」割合の標本分布を作成してください。5,000回のシミュレーションを使用してください。標本分布の各観測値は何を表していますか？標本サイズが増加するにつれて、標本分布の平均、標準誤差、形状はどのように変化しますか？シミュレーションの数を増やした場合、これらの値はどのように（もし変化するなら）変化しますか？（回答にプロットを含める必要はありません。）

```{r shiny, echo=FALSE, eval=TRUE, results = TRUE}
shinyApp(
  ui <- fluidPage(
    
    # Sidebar with a slider input for number of bins 
    sidebarLayout(
      sidebarPanel(
        
        selectInput("outcome",
                    "Outcome of interest:",
                    choices = c("Benefits", "Doesn't benefit"),
                    selected = "Doesn't benefit"),
        
        numericInput("n_samp",
                     "Sample size:",
                     min = 1,
                     max = nrow(global_monitor),
                     value = 30),
        
        numericInput("n_rep",
                     "Number of samples:",
                     min = 1,
                     max = 30000,
                     value = 15000),
        
        hr(),
        
        sliderInput("binwidth",
                    "Binwidth:",
                    min = 0, max = 0.5,
                    value = 0.02,
                    step = 0.005)
        
      ),
      
      # Show a plot of the generated distribution
      mainPanel(
        plotOutput("sampling_plot"),
        textOutput("sampling_mean"),
        textOutput("sampling_se")
      )
    )
  ),
  
  server <- function(input, output) {
    
    # create sampling distribution
    sampling_dist <- reactive({
      global_monitor |>
        rep_sample_n(size = input$n_samp, reps = input$n_rep, replace = TRUE) |>
        count(scientist_work) |>
        mutate(p_hat = n /sum(n)) |>
        filter(scientist_work == input$outcome)
    })
    
    # plot sampling distribution
    output$sampling_plot <- renderPlot({
      
      ggplot(sampling_dist(), aes(x = p_hat)) +
        geom_histogram(binwidth = input$binwidth) +
        xlim(0, 1) +
        labs(
          x = paste0("p_hat (", input$outcome, ")"),
          title = "Sampling distribution of p_hat",
          subtitle = paste0("Sample size = ", input$n_samp, " Number of samples = ", input$n_rep)
        ) +
        theme(plot.title = element_text(face = "bold", size = 16))
    })
    
    ggplot(data = sample_props50, aes(x = p_hat)) +
      geom_histogram(binwidth = 0.02) +
      labs(
        x = "p_hat (Doesn't benefit)",
        title = "Sampling distribution of p_hat",
        subtitle = "Sample size = 50, Number of samples = 15000"
      )
    
    # mean of sampling distribution
    output$sampling_mean <- renderText({
      paste0("Mean of sampling distribution = ", round(mean(sampling_dist()$p_hat), 2))
    })
    
    # mean of sampling distribution
    output$sampling_se <- renderText({
      paste0("SE of sampling distribution = ", round(sd(sampling_dist()$p_hat), 2))
    })
  },
  
  options = list(height = 900) 
)
```

------------------------------------------------------------------------

## さらなる演習

これまで、科学者の仕事が自分たちに恩恵をもたらさないと考える人々の割合を推定することのみに焦点を当ててきました。次に、恩恵をもたらすと考える人々の割合を推定してみます。

これらの質問のいくつかはアプリを使用して答えられるかもしれませんが、必要なコードを記述し、必要なプロットと要約統計量を生成することが期待されます。探索のためにアプリを使用することは歓迎します。

7.  母集団からサイズ15のサンプルを取り、このサンプル中の人々が科学者の仕事が彼らの生活を向上させると考える割合を計算してください。このサンプルを使用して、科学者の仕事が彼らの生活を向上させると考える人々の母集団の割合の最良の点推定値は何ですか？

8.  母集団にアクセスできるので、母集団からサイズ15の2000個のサンプルを取り、2000個の標本比率を計算することで、科学者の仕事が彼らの生活を向上させると考える人々の割合のサイズ15のサンプルの標本分布をシミュレートしてください。これらの割合を`sample_props15`として保存してください。データをプロットし、この標本分布の形状を記述してください。この標本分布に基づいて、科学者の仕事が彼らの生活を向上させると考える人々の真の割合はいくつだと推測しますか？最後に、母集団の割合を計算して報告してください。

9.  サンプルサイズを15から150に変更し、上記と同じ方法で標本分布を計算し、これらの割合を`sample_props150`という新しいオブジェクトに保存してください。この標本分布の形状を記述し、サイズ15のサンプルの標本分布と比較してください。この標本分布に基づいて、科学者の仕事が彼らの生活を向上させると考える人々の真の割合はいくつだと推測しますか？

10. 8と9の標本分布のうち、広がりが小さいのはどちらですか？真の値に近い推定値をより頻繁に行うことを気にする場合、広がりが大きい標本分布と小さい標本分布のどちらを好みますか？

------------------------------------------------------------------------

<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">![Creative Commons License](https://i.creativecommons.org/l/by-sa/4.0/88x31.png){style="border-width:0"}</a><br />この作品は、<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">クリエイティブ・コモンズ 表示 - 継承 4.0 国際ライセンス</a>の下に提供されています。
