第2講 R言語速習

練習問題

準備

データの操作と視覚化のために tidyverse パッケージを読み込む．

```{r}
#| warning: true
library(tidyverse)
```

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Quartoファイルの中で使われているパッケージがインストールされていない場合は， エディタ上部にインストールを促す警告が出るので，その指示に従ってインストールすればよい． 手動でインストールする場合は以下のいずれかを実行する．

• Package タブからインストール
• コンソール上で次のコマンドを実行 install.packages("パッケージ名")

tidyverse は第1講のサンプルで使っているので通常はインストールされているが， R言語のバージョンがアップグレードされた場合などに対応が必要となる．

library(tidyverse) を実行すると読み込まれたパッケージが表示される． 同名の関数が存在する場合には Conflicts として衝突する関数名が表示される． 衝突する場合はパッケージ名を明示的に付ける必要がある． 例えば filter() には以下の2つがある．

• dplyer::filter() データフレームの抽出のための関数
• stats::filter() 時系列処理のための線形フィルタ関数

名前の衝突による不具合を避けたい場合は conflicted パッケージの利用を推奨する．

library(conflicted)

データフレームの作成

データフレームは 同じ長さのベクトル(関数 base::c() で作成)を関数 tibble::tibble() に渡して作成する．

```{r}
#| eval: false
#' (... <- ...) は代入した結果を表示
(foo <- tibble(one = c(1,2,3),two = c("AB","CD","EF"))) 
(bar <- tibble(three = c("x","y","z"),four = c(0.9,0.5,-0.3)))
```

データフレームを結合するには関数 dplyr::bind_cols() を用いる．

```{r}
#| eval: false
(baz <- bind_cols(foo,bar))
```

問題

次の表に対応するデータフレームを作成しなさい．

name	math	phys	chem	bio
Alice	90	25	65	70
Bob	80	50	100	50
Carol	70	75	70	30
Dave	60	100	40	80
Eve	50	80	75	100

解答欄

```{r}
#' ここにコードを書く
```

解答例

各項目が同じ長さのベクトルなので，これらを並べればよい．

```{r}
#| code-fold: true
(grade_data <- tibble( # 変数名は自由に決めてよい
   name = c("Alice", "Bob", "Carol", "Dave", "Eve"),
   math = c(90, 80, 70, 60, 50),
   phys = c(25, 50, 75,100, 80),
   chem = c(65,100, 70, 40, 75),
   bio  = c(70, 50, 30, 80,100)))
```

# A tibble: 5 × 5
  name   math  phys  chem   bio
  <chr> <dbl> <dbl> <dbl> <dbl>
1 Alice    90    25    65    70
2 Bob      80    50   100    50
3 Carol    70    75    70    30
4 Dave     60   100    40    80
5 Eve      50    80    75   100

行や列の名前を操作するには以下のようにする．

```{r}
#| code-fold: true
names(grade_data)   # 列の名前を表示する
print(grade_data)   # データフレームの内容を表示(printは省略できる)
glimpse(grade_data) # データフレームの内容を別の形式で表示
```

[1] "name" "math" "phys" "chem" "bio" 
# A tibble: 5 × 5
  name   math  phys  chem   bio
  <chr> <dbl> <dbl> <dbl> <dbl>
1 Alice    90    25    65    70
2 Bob      80    50   100    50
3 Carol    70    75    70    30
4 Dave     60   100    40    80
5 Eve      50    80    75   100
Rows: 5
Columns: 5
$ name <chr> "Alice", "Bob", "Carol", "Dave", "Eve"
$ math <dbl> 90, 80, 70, 60, 50
$ phys <dbl> 25, 50, 75, 100, 80
$ chem <dbl> 65, 100, 70, 40, 75
$ bio  <dbl> 70, 50, 30, 80, 100

データの取り出し方は後ほど詳しく述べるが，例えば以下のように操作することができる．

```{r}
#| code-fold: true
grade_data[2,3]      # 特定の要素を数値で参照する
grade_data[2,"phys"] # 列を名前で参照する (上記と同じ結果)
grade_data[3,]       # 特定の行を表示 (データフレームになる)
grade_data["bio"]    # 特定の列を表示 (データフレームになる)
grade_data[,"bio"]   # 上記と同じ結果
grade_data[["bio"]]  # ベクトルとして取り出す (リストとしての処理)
grade_data$bio       # 上記と同じ結果
```

# A tibble: 1 × 1
   phys
  <dbl>
1    50
# A tibble: 1 × 1
   phys
  <dbl>
1    50
# A tibble: 1 × 5
  name   math  phys  chem   bio
  <chr> <dbl> <dbl> <dbl> <dbl>
1 Carol    70    75    70    30
# A tibble: 5 × 1
    bio
  <dbl>
1    70
2    50
3    30
4    80
5   100
# A tibble: 5 × 1
    bio
  <dbl>
1    70
2    50
3    30
4    80
5   100
[1]  70  50  30  80 100
[1]  70  50  30  80 100

変数や関数の名称については各自で命名規則を決めておくと良い． 例えば以下のサイトが参考になる．

• https://style.tidyverse.org/index.html
• https://google.github.io/styleguide/Rguide.html

ファイルの読み書き

ファイルの書き出しには関数 readr::write_csv() を用いる．

write_csv(x, file = "ファイル名") # データフレームxをファイルに書き出す

ファイルの読み込みには関数 readr::~read_csv() を用いる．

y <- read_csv(file = "ファイル名") # 変数yにCSVファイルの内容を読み込む

問題

以下の問いに答えなさい．

1. 前の演習で作成したデータフレームを適当なファイルに書き出しなさい．
2. 書き出したファイルから別の変数に読み込みなさい．
3. pcr_case_daily.csv (厚労省からダウンロードしたファイル)を変数 pcr_data に読み込みなさい．

解答欄

```{r}
#' ここにコードを書く
```

解答例

前の練習問題で作ったデータフレームを作業ディレクトリの data 以下に grade_data.csv として保存する．

```{r}
#| code-fold: true
write_csv(grade_data, file = "data/grade_data.csv")
```

保存したファイルは File タブからその中身を含め確認することができる．

grade_copy というオブジェクトにファイルの内容を代入して，その中身を確認する．

```{r}
#| code-fold: true
(grade_copy <- read_csv(file = "data/grade_data.csv"))
```

# A tibble: 5 × 5
  name   math  phys  chem   bio
  <chr> <dbl> <dbl> <dbl> <dbl>
1 Alice    90    25    65    70
2 Bob      80    50   100    50
3 Carol    70    75    70    30
4 Dave     60   100    40    80
5 Eve      50    80    75   100

Files タブの操作で読み込みことも可能である．

ダウンロードしたファイルをファイル名 pcr_case_daily.csv として作業ディレクトリの data 以下に保存しておく．

```{r}
#| code-fold: true
pcr_data <- read_csv("data/pcr_case_daily.csv") 
print(pcr_data) # 中身の一部を表示
```

# A tibble: 1,175 × 10
   日付    国立感染症研究所 検疫所 地方衛生研究所・保健所…¹… 民間検査会社（主に行政検査）…²…
   <chr>              <dbl>  <dbl>                  <dbl>                  <dbl>
 1 2020/2…              472     75                    398                      0
 2 2020/2…               15     68                    609                      0
 3 2020/2…               20     15                    758                      0
 4 2020/2…              261    188                    902                    132
 5 2020/2…              341    127                    677                      2
 6 2020/2…               53     72                    529                      0
 7 2020/2…               22    103                    460                     17
 8 2020/2…              195     38                    695                      0
 9 2020/2…              267     19                    970                    149
10 2020/2…              237     61                   1101                      0
# ℹ 1,165 more rows
# ℹ abbreviated names: ¹​`地方衛生研究所・保健所`,
#   ²​`民間検査会社（主に行政検査）`
# ℹ 5 more variables: 大学等 <dbl>, 医療機関 <dbl>, 小計 <dbl>,
#   `民間検査会社（主に自費検査）` <dbl>, 計 <dbl>

列名を確認して pcr_colnames に保存しておく．

```{r}
#| code-fold: true
(pcr_colnames <- names(pcr_data)) # colnames(pcr_data) でも良い
```

 [1] "日付"                         "国立感染症研究所"            
 [3] "検疫所"                       "地方衛生研究所・保健所"      
 [5] "民間検査会社（主に行政検査）" "大学等"                      
 [7] "医療機関"                     "小計"                        
 [9] "民間検査会社（主に自費検査）" "計"                          

列名を扱い易いように英語略記に変更する．

略称はそれぞれ以下の意味．

• National Institute of Infectious Diseases
• Customs-Immigration-Quarantine
• Health Center
• Administrative Inspection
• University
• Medical Institution
• subtotal
• Self Inspection
• total

```{r}
#| code-fold: true
names(pcr_data) <- # 
  c("date","niid","ciq","hc","ai","univ","mi","sub","si","total")
pcr_data # 中身を確認(10行だけ表示される)
names(pcr_colnames) <- names(pcr_data) # 和英の列名の対応づけができるようにしておく
pcr_colnames
```

# A tibble: 1,175 × 10
   date       niid   ciq    hc    ai  univ    mi   sub    si total
   <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020/2/18   472    75   398     0    79    NA  1024    NA  1024
 2 2020/2/19    15    68   609     0     0    NA   692    NA   692
 3 2020/2/20    20    15   758     0     0    NA   793    NA   793
 4 2020/2/21   261   188   902   132   108    NA  1591    NA  1591
 5 2020/2/22   341   127   677     2    19    NA  1166    NA  1166
 6 2020/2/23    53    72   529     0     0    NA   654    NA   654
 7 2020/2/24    22   103   460    17     0    NA   602    NA   602
 8 2020/2/25   195    38   695     0     0    NA   928    NA   928
 9 2020/2/26   267    19   970   149     0    NA  1405    NA  1405
10 2020/2/27   237    61  1101     0     0    NA  1399    NA  1399
# ℹ 1,165 more rows
                          date                           niid 
                        "日付"             "国立感染症研究所" 
                           ciq                             hc 
                      "検疫所"       "地方衛生研究所・保健所" 
                            ai                           univ 
"民間検査会社（主に行政検査）"                       "大学等" 
                            mi                            sub 
                    "医療機関"                         "小計" 
                            si                          total 
"民間検査会社（主に自費検査）"                           "計" 

以降の処理のために date 列を関数 lubridate::date() で date 型に変換する． 列の変換・追加などには関数 dplyr::mutate() を用いる．

```{r}
#| code-fold: true
(pcr_data <- mutate(pcr_data, date = date(date)))
```

# A tibble: 1,175 × 10
   date        niid   ciq    hc    ai  univ    mi   sub    si total
   <date>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020-02-18   472    75   398     0    79    NA  1024    NA  1024
 2 2020-02-19    15    68   609     0     0    NA   692    NA   692
 3 2020-02-20    20    15   758     0     0    NA   793    NA   793
 4 2020-02-21   261   188   902   132   108    NA  1591    NA  1591
 5 2020-02-22   341   127   677     2    19    NA  1166    NA  1166
 6 2020-02-23    53    72   529     0     0    NA   654    NA   654
 7 2020-02-24    22   103   460    17     0    NA   602    NA   602
 8 2020-02-25   195    38   695     0     0    NA   928    NA   928
 9 2020-02-26   267    19   970   149     0    NA  1405    NA  1405
10 2020-02-27   237    61  1101     0     0    NA  1399    NA  1399
# ℹ 1,165 more rows

以下，この項で用いた関数に関する補足．

#' 関数 print() を用いると表示する行数を指定できる
print(pcr_data, n = 5) # 全ては n = Inf
#' 日本語を含むファイルでは文字化けが起こった場合は以下で対応する
#' 関数 readr::guess_encoding() でファイルの文字コードを推測する
guess_encoding("data/pcr_case_daily.csv")
#' "UTF-8" であると 1 の信頼度で認識される
#' 文字コードを指定して読み込む場合は以下のように記述する
pcr_data <- # 文字コードとして UTF-8 を指定
  read_csv(file = "data/pcr_case_daily.csv",
           locale = locale(encoding = "utf-8"))
#' その他の文字コードとしては "sjis", "shift-jis", "shift_jis", 
#' "cp932"(拡張文字を含む)などを大文字小文字は区別せず指定できる
#' URLを指定して読み込むこともできる 
pcr_data <- # 更新される情報を追跡する場合に利用を推奨
  read_csv("https://www.mhlw.go.jp/content/pcr_case_daily.csv")
#' 列名の変更にはいろいろな方法があるので適宜使用する
#' 読み込み時に行う方法
pcr_data <- read_csv("data/pcr_case_daily.csv",
                     skip = 1, # 列名の行を読み飛ばす
                     col_names = c("date","niid","ciq","hc","ai",
                                   "univ","mi","sub","si","total"))
pcr_data <- mutate(pcr_data, date = date(date)) # date型に変更
#' 関数 dplyr::rename() を使う方法
pcr_data <- read_csv("data/pcr_case_daily.csv") # そのまま読み込む
(pcr_colnames <- set_names(names(pcr_data), # 新旧の列名に対応するベクトルを作成
                           c("date","niid","ciq","hc","ai",
                             "univ","mi","sub","si","total")))
pcr_data <- rename(pcr_data, all_of(pcr_colnames)) # 列名を変更
pcr_data <- mutate(pcr_data, date = date(date)) # date型に変更

データフレームの操作

行の選択には関数 dplyer::filter() を， 列の選択には関数 dplyer::select() を用いる．

```{r}
#| eval: false
#' 前に作成したデータフレーム z を用いた例
(foo <- filter(z, three >= 7))      # 列 three の値が7以上の行を選択
(bar <- select(foo, c(one, three))) # 列 one,three を選択
#' パイプを用いると以下のように簡潔に書ける
z |> filter(three >= 7) |> select(one, three)
```

問題

pcr_case_daily.csv から以下の条件を満たすデータを取り出しなさい．

1. 医療機関 (mi) での検査件数が2000を越えたときの国立感染症研究所 (niid) と医療機関 (mi) のデータ．
2. 大学等 (univ) と医療機関 (mi) でともに検査件数が2000を越えたデータ．
3. 2020年3月の各機関(sub,total は集計なので除く)の検査件数データ．

解答欄

```{r}
#' ここにコードを書く
```

解答例

それぞれ以下のようにすれば条件に合致したデータを抽出することができる．

```{r}
#| code-fold: true
pcr_data |>            # データフレーム
  filter(mi > 2000) |> # 行の条件による絞り込み
  select(c(niid,mi))   # 列の選択 select(niid,mi) としても良い
pcr_data |>
  filter(univ > 2000 & mi > 2000) # 複合的な条件の指定
pcr_data |>
  filter(date >= "2020-03-01" & date < "2020-04-01") |> # 日付の範囲での指定
  select(!c(sub,total)) # 列の削除は select(-c(sub,total)), select(-sub,-total) も可
```

# A tibble: 977 × 2
    niid    mi
   <dbl> <dbl>
 1     0  3510
 2     0  2636
 3     0  3459
 4     0  3694
 5     0  3437
 6    83  2579
 7     0  4247
 8    36  3469
 9     0  3725
10    20  3643
# ℹ 967 more rows
# A tibble: 746 × 10
   date        niid   ciq    hc    ai  univ    mi   sub    si total
   <date>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020-08-12     0   107  5805 15264  2152  3694 27022    NA 27022
 2 2020-08-17     0    11  3193 14606  2303  4247 24360    NA 24360
 3 2020-08-18    36     4  5631 15884  2322  3469 27346    NA 27346
 4 2020-08-19     0    77  5505 16463  2412  3725 28182    NA 28182
 5 2020-08-20    20    37  4940 17632  2397  3643 28669    NA 28669
 6 2020-08-24     0     3  3123  6599  2482  3454 15661    NA 15661
 7 2020-08-25     0     0  4740 14602  2291  3512 25145    NA 25145
 8 2020-08-26     0   111  4611 16481  2496  3488 27187    NA 27187
 9 2020-08-27    27     0  5183 16751  2071  3412 27444    NA 27444
10 2020-08-31     0     0  2948 16133  2714  4578 26373    NA 26373
# ℹ 736 more rows
# A tibble: 31 × 8
   date        niid   ciq    hc    ai  univ    mi    si
   <date>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020-03-01     0    22   544   116     0    NA    NA
 2 2020-03-02   125    24  1038    15     0    NA    NA
 3 2020-03-03    20    27  1718    11     0    NA    NA
 4 2020-03-04    71    12  1398     7     0    NA    NA
 5 2020-03-05    83    11  1579     5     0    NA    NA
 6 2020-03-06   107    35  1786     5    62    68    NA
 7 2020-03-07    18    13  1616    11    20    44    NA
 8 2020-03-08     0    10   579     2     8     9    NA
 9 2020-03-09     0    27  1072     4    97    48    NA
10 2020-03-10     8    57  1763    97    65    62    NA
# ℹ 21 more rows

データフレームの集約

列の集計には関数 dplyr::summarise() を用いる．

```{r}
#| eval: false
#' 練習問題のデータフレームを用いた例
grade_data |>
  summarise(math_mean = mean(math), nums = n()) # 数学の平均を求める
#' 複数の列にまたがる操作には関数 dplyr::across() を利用
grade_data |>
  summarise(across(!name, mean))                # 名前の列以外の平均を求める
```

集計に必要な関数が用意されていない場合は， 無名関数として定義して利用すれば良い．

問題

pcr_case_daily.csv について以下の集計を行いなさい．

1. 各機関でのPCR検査件数の最大値．
2. 2021年の各機関でのPCR検査件数の月ごとの最大値．

ヘルプを用いて datasets::mtcars の内容を調べた上で以下の集計を行いなさい．

1. 気筒数 (cyl) ごとに排気量 (disp) の最大値，最小値．
2. 気筒数 (cyl) とギア数 (gear) ごとの燃費 (mpg) の平均値．

解答欄

```{r}
#' ここにコードを書く
```

解答例

pcr_case_daily.csv において各機関でのPCR検査件数の最大値を求めるには基本的には以下のようにすればよい．

```{r}
#| code-fold: true
pcr_data |> summarise(across(!date, max))
```

# A tibble: 1 × 9
   niid   ciq    hc     ai  univ    mi    sub    si total
  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
1   517  1733 16973 161077 12731    NA 272018    NA    NA

NAが含まれる列については最大値求めることができないので，最大値の計算で NA (欠損データ) を除く必要がある．

```{r}
#| code-fold: true
pcr_data |> # max の無名関数を利用する方法
  summarise(across(!date, \(x) max(x, na.rm = TRUE)))
pcr_data |> # package::purrr の lambda 式を利用する方法
  summarise(across(!date, ~ max(.x, na.rm = TRUE)))
```

# A tibble: 1 × 9
   niid   ciq    hc     ai  univ     mi    sub     si  total
  <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1   517  1733 16973 161077 12731 136570 272018 131822 371225
# A tibble: 1 × 9
   niid   ciq    hc     ai  univ     mi    sub     si  total
  <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1   517  1733 16973 161077 12731 136570 272018 131822 371225

2021年の月ごとの各機関でのPCR検査件数の最大値は以下のようになる．

```{r}
#| code-fold: true
pcr_data |>
  filter(year(date) == 2021) |>
  group_by(month(date)) |>
  summarise(across(!date, max))
```

# A tibble: 12 × 10
   `month(date)`  niid   ciq    hc     ai  univ    mi    sub     si  total
           <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>  <dbl>  <dbl>  <dbl>
 1             1     0    23 12042  58398  6556 23417  91261     NA  91261
 2             2     0    38  7032  39353  5533 22878  68847     NA  68847
 3             3     0     3  5117  66274  5005 23983  95939     NA  95939
 4             4     1     4  8866  75440  5121 30378 115724     NA 115724
 5             5     0    20 11181  89381  6392 35379 130410     NA 130410
 6             6     0    28  5433  68911  5439 28268  95330     NA  95330
 7             7    10    21  7739  63796  6500 36674 104740     NA 174233
 8             8     1    23 13200 109734  7397 54257 168857 110172 272564
 9             9     0    20 11602  94414  7472 43944 150671  90013 239984
10            10     0    21  2294  56724  5521 29978  86137  76467 159101
11            11     0    21  1022  31346  5706 27344  55361  71597 123910
12            12     0     6  2452  24376  5733 26647  51718  66888 112438

datasets::mtcars の集計は以下のとおりである．

```{r}
#| code-fold: true
mtcars |>
  group_by(cyl) |>
  summarise(max_disp = max(disp))
mtcars |>
  group_by(cyl) |>
  summarise(min_disp = min(disp))
mtcars |> # まとめて計算することも可能
  group_by(cyl) |> # 列名の作られ方に注意
  summarise(across(disp, list(max = max, min = min)))
mtcars |>
  group_by(cyl, gear) |>
  summarise(mpg = mean(mpg))
```

# A tibble: 3 × 2
    cyl max_disp
  <dbl>    <dbl>
1     4     147.
2     6     258 
3     8     472 
# A tibble: 3 × 2
    cyl min_disp
  <dbl>    <dbl>
1     4     71.1
2     6    145  
3     8    276. 
# A tibble: 3 × 3
    cyl disp_max disp_min
  <dbl>    <dbl>    <dbl>
1     4     147.     71.1
2     6     258     145  
3     8     472     276. 
# A tibble: 8 × 3
# Groups:   cyl [3]
    cyl  gear   mpg
  <dbl> <dbl> <dbl>
1     4     3  21.5
2     4     4  26.9
3     4     5  28.2
4     6     3  19.8
5     6     4  19.8
6     6     5  19.7
7     8     3  15.0
8     8     5  15.4

グループは既定値では順次解除されるので以下のような集計も可能である．

```{r}
#| eval: false
mtcars |> 
  group_by(cyl, gear) |>
  summarise(mpg = mean(mpg)) |> # cyl のグループは残っている
  summarise(mpg = max(mpg)) # cyl ごとに平均値の最大値を求める
```

基本的なグラフの描画

グラフの描画には 関数 ggplot2::ggplot() で描画に用いるデータフレームを指定し， 描きたいグラフの種類に応じて関数 ggplot2::geom_XXX() を加える．

```{r}
#| eval: false
pcr_data |> # パイプ演算子でデータフレームを関数 ggplot2::ggplot() に渡す
  ggplot(aes(x = date)) + # date をx軸に指定
  geom_line(aes(y = ai), colour = "blue") + # 行政検査を青
  geom_line(aes(y = mi), colour = "red") +  # 医療機関を赤
  labs(y = "number of tests") # y軸のラベルを変更
```

散布図行列を描くには関数 GGally::ggpairs() を用いると良い．

```{r}
#| eval: false
library(GGally)
pcr_data |>
    select(!c(date,sub,total)) |> # 日付と集計値を除いて必要なデータフレームに整形
    ggpairs() # 標準の散布図行列
```

問題

pcr_case_daily.csv を用いて以下の描画を行いなさい．

1. 検疫所 (b)，地方衛生研究所.保健所 (c)，民間検査会社 (d) における検査件数の推移．
2. 民間検査会社 (d)，大学等 (e)，医療機関 (f) での検査件数の関係 (散布図)．

解答欄

```{r}
#' ここにコードを書く
```

解答例

書き方はいろいろあるので，以下はあくまで一例である．

```{r}
#| code-fold: true
if(Sys.info()["sysname"] == "Darwin") { # MacOSの場合には日本語フォントを指定する
  theme_update(text = element_text(family = "HiraginoSans-W4"))}
pcr_data |>
  select(c(date,ciq,hc,ai)) |> # 描画対象の列を抽出
  pivot_longer(!date, names_to = "organ", values_to = "nums") |> # 
  ggplot(aes(x = date, y = nums, colour = organ)) +
  geom_line() +
  labs(x = "日付", y = "検査件数")
```

y軸を対数表示にする場合は以下のようにすればよい．

```{r}
#| code-fold: true
pcr_data |>
  select(c(date,ciq,hc,ai)) |> # 描画対象の列を抽出
  pivot_longer(!date, names_to = "organ", values_to = "nums") |> # 
  ggplot(aes(x = date, y = nums, colour = organ)) +
  geom_line() +
  scale_y_log10() + # y軸を対数表示 (log10(0)=-Inf の警告が出る場合がある)
  labs(x = "日付", y = "検査件数")
```

散布図を並べるために GGally パッケージを利用する．

```{r}
#| code-fold: true
#| fig-width: 7
#| fig-height: 7
library(GGally)
pcr_data |>
  select(c(ai,univ,mi)) |> # 描画対象の列を抽出
  ggpairs(columnLabels = pcr_colnames[c("ai","univ","mi")]) # ラベルを渡す
```

さまざまなグラフ

行政検査での検査件数の分布

行政検査(ai)での検査件数の分布を描画する．

```{r}
pcr_data |>
  ggplot(aes(x = ai)) + # 分布を描画する列を指定
  geom_histogram(bins = 30, fill = "lightblue", colour = "blue") +
  labs(x = pcr_colnames["ai"], y = "頻度", title = "検査件数のヒストグラム")
```

各ビンの頻度を表示するためには例えば以下のようにすればよい．

```{r}
pcr_data |>
  ggplot(aes(x = ai)) + # 分布を描画する列を指定
  geom_histogram(bins = 30, fill = "lightblue", colour = "blue") +
  geom_text(stat="bin", bins = 30, colour = "darkblue", size = 3, 
            aes(label = after_stat(count), y = after_stat(count) + 2)) +
  labs(x = pcr_colnames["ai"], y = "頻度", title = "検査件数のヒストグラム")
```

月ごとの大学等での検査件数の分布(分位点)

大学等(univ)での検査件数の分布(2021年分)を描画する．

```{r}
pcr_data |>
  filter(year(date) == 2021) |> # 2021年を抽出
  mutate(date = as_factor(month(date))) |> # 月を因子化する
  ggplot(aes(x = date, y = univ)) + # 月毎に集計する
  geom_boxplot(fill = "orange") + # 塗り潰しの色を指定
  labs(title = "月ごとの検査件数 (2021年)", x = "月", y = pcr_colnames["univ"])
```

機関ごとの月の検査件数の推移

機関ごとの月の検査件数の推移(2021年分)を描画する．

```{r}
pcr_data |>
  filter(year(date) == 2021) |>
  mutate(month = as_factor(month(date))) |> # 月を作成
  select(!c(date,sub,total)) |> # 機関に限定
  group_by(month) |> # 月でグループ化
  summarize(across(everything(), sum)) |> # 全て(月以外)を集計
  pivot_longer(!month, names_to = "organ", values_to = "nums",
               names_transform = list(organ = as_factor)) |>
  ## 最後のオプションは organ 列のラベルを出てきた順で因子化して元の列の並びにしている
  ggplot(aes(x = organ, y = nums, fill = month)) +
  geom_bar(stat = "identity", position = "dodge", na.rm = TRUE) +
  theme(legend.position = "top") + guides(fill = guide_legend(nrow = 1))
```

擬似乱数

問題

ヘルプを用いて以下の関数を調べよ．

1. 関数 base::sample()
2. 関数 stats::rbinom()
3. 関数 stats::runif()
4. 関数 stats::rnorm()
5. 関数 base::set.seed()

以下の試行を実装してみよ．

1. サイコロを10回振る．
2. 4枚のコインを投げたときの表の枚数．

解答欄

```{r}
#' ここにコードを書く
```

解答例

それぞれの関数の基本的な使い方は以下の例のようになる(問題の試行も含まれる)．

```{r}
#| code-fold: true
#' 関数sampleの使い方
(x <- 1:10)   # サンプリング対象の集合を定義
set.seed(123) # 乱数のシード値(任意に決めてよい)を指定
sample(x, 5)                 # xから5つの要素を重複なしでランダムに抽出
sample(x, 5, replace = TRUE) # xから5つの要素を重複ありでランダムに抽出
sample(x, length(x))         # xの要素のランダムな並べ替え
sample(1:6, 10, replace = TRUE)             # サイコロを10回振る実験の再現
sample(1:6, 10, prob = 6:1, replace = TRUE) # 出る目の確率に偏りがある場合

#' 関数rbinomの使い方
rbinom(10, size = 4, prob = 0.5) # 表(1)の出る確率が0.5にコインを4枚投げる試行を10回
rbinom(20, size = 4, prob = 0.2) # 個数を20, 確率を0.2に変更

#' 関数runifの使い方
runif(5, min = -1, max = 2) # 区間(-1,2)上の一様乱数を5個発生
runif(5)                    # 指定しない場合は区間(0,1)が既定値

#' 関数rnormの使い方
rnorm(10, mean = 5, sd = 3) # 平均5，分散3^2の正規乱数を10個発生
rnorm(10)                   # 指定しない場合は mu=0, sd=1 が既定値

#' 関数set.seedについて
set.seed(1) # 乱数の初期値をseed=1で指定
runif(5) 
set.seed(2) # 乱数の初期値をseed=2で指定
runif(5)    # seed=1の場合と異なる結果
set.seed(1) # 乱数の初期値をseed=1で指定
runif(5)    # 初めのseed=1の場合と同じ結果
```

 [1]  1  2  3  4  5  6  7  8  9 10
[1]  3 10  2  8  6
[1]  5  4  6  9 10
 [1]  5  3  9  1  4  7 10  6  8  2
 [1] 2 2 1 6 3 4 6 1 3 5
 [1] 1 1 2 2 2 1 1 1 2 1
 [1] 3 0 2 3 1 2 1 1 3 3
 [1] 0 1 0 0 0 1 1 1 1 1 1 1 1 1 0 1 0 0 1 0
[1] -0.6665937 -0.2691416  1.0041668  0.2529403  1.3645875
[1] 0.1028646 0.4348927 0.9849570 0.8930511 0.8864691
 [1]  2.196845  6.181126  6.210894  2.340690  1.043187  5.086532  3.703611
 [8] 10.069618  8.685178  5.828070
 [1] -1.0489755 -0.5208693  1.6232025 -1.0700682  1.6858872 -0.2416898
 [7] -0.4682005 -0.7729782  2.1499193 -1.3343536
[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
[1] 0.1848823 0.7023740 0.5733263 0.1680519 0.9438393
[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819

モンテカルロ法

例 : 中心極限定理

“多数の独立な確率変数の和はほぼ正規分布に従う” ことを主張している．これをモンテカルロ法(確率シミュレーション)で確認する．

確率変数として区間[-1,1]の一様乱数を設定する．

```{r}
mc_rand <- function(n) { # n個の乱数を生成
  return(runif(n, min = -1, max = 1))
}
```

標本平均の計算を行う関数を定義する．

```{r}
mc_mean <- function(n) { # n個のデータで計算
  return(mean(mc_rand(n)))
}
```

実験を行う．

```{r}
set.seed(123) # 実験を再現したい場合はシード値を指定する
mu <- 0; sigma <- sqrt(1/3) # 理論平均と標準偏差
mc_num <- 5000 # 実験の繰り返し回数
for(n in c(1,2,4,8,16)){ # nを変えて実験
  p <- tibble(x = replicate(mc_num, mc_mean(n))) |> # 繰り返し実験し標本平均を記録
    ggplot(aes(x)) + 
    geom_histogram(aes(y = after_stat(density)), # 密度表示
                   fill = "orchid", alpha = 0.5, # 塗り潰しの色
                   colour = "purple") + # 境界線の色
    geom_function(fun = \(x) dnorm(x, mean = mu, sd = sigma/sqrt(n)),
                  colour = "orange", linewidth = 1.5) + # 理論曲線を重ねる
    labs(x = expression(bar(X)), # x軸の表示
         title = paste0("n=", n)) # タイトルにnを記載
  print(p) # for 文の中では明示的に print する必要がある
}
```

例 : コイン投げの賭け

AとBの二人で交互にコインを投げ，最初に表が出た方を勝ちとするとき，AとBそれぞれの勝率はいくつとなるかをモンテカルロ法で実験する．

コイン投げは関数 sample(), rbinom() などを用いて模擬できる

```{r}
#| eval: false
sample(0:1, 1) # 0と1が入った壺からから1つ選ぶ
rbinom(1, size = 1, prob = 0.5) # 表裏が等確率で出る1枚のコインを1回投げる
```

コイン投げの試行を定義する．いろいろな書き方があるので以下は一例である．

```{r}
mc_trial <- function(){
  while(TRUE){ # 永久に回るループ
    if(rbinom(1, size = 1, prob = 0.5)==1){return("A")} # Aが表で終了
    if(rbinom(1, size = 1, prob = 0.5)==1){return("B")} # Bが表で終了
    #' どちらも裏ならもう一度ループ
  }
}
```

実験を行う．

```{r}
set.seed(8888) # 実験を再現したい場合はシード値を指定する
mc_num <- 10000 # 実験回数を設定 
mc_data <- replicate(mc_num, mc_trial()) 
```

実験結果の簡単な集計は以下のようになる．

```{r}
table(mc_data)        # 頻度
table(mc_data)/mc_num # 確率(推定値)
```

mc_data
   A    B 
6810 3190 
mc_data
    A     B 
0.681 0.319 

問題

以下の簡単な双六ゲームを考える．

• ゴールまでのます目は100とする．
• さいころを振り出た目の数だけ進む．
• ゴールに辿り着くまで繰り返す．

さいころを振る回数の分布がどうなるか実験を行いなさい．

解答欄

```{r}
#' ここにコードを書く
```

解答例

双六の試行を行う関数を作成する．

```{r}
#| code-fold: true
mc_trial <- function(){
  step <- 0 # 最初の位置
  num <- 0  # さいころを振る回数
  while(TRUE){ # 永久に回るループ
    step <- step + sample(1:6, 1) # さいころを振る
    num <- num + 1 # 回数を記録
    if(step >= 100) { # ゴールしたか?
      return(num) # 回数を出力して関数を終了
    }
  }
}
```

同じ試行でも関数の作り方はいろいろある． 例えば100回サイコロを振れば必ずどこかで100を越えるので，計算は無駄があるが条件分岐のない関数を考えることもできる．

mc_trial <- function() {
  which.max(cumsum(sample(1:6, 100, replace = TRUE)) >= 100)
} # 関数 which.max() は初めて TRUE(1) になった場所を返す

試しに10回の試行を実行してみる．

```{r}
#| code-fold: true
for(i in 1:10) print(mc_trial())
```

[1] 32
[1] 33
[1] 31
[1] 36
[1] 34
[1] 24
[1] 25
[1] 30
[1] 26
[1] 27

確率シミュレーション(Monte-Carlo実験)を行う．

```{r}
#| code-fold: true
set.seed(12345) # 必要に応じて乱数のシードを設定する
mc_num <- 10000 # 実験回数を設定 
mc_data <- replicate(mc_num, mc_trial()) 
summary(mc_data) # 簡単な集計
```

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  20.00   27.00   29.00   29.09   31.00   41.00 

同じ操作を繰り返すには関数 purrr::map() を用いることもできる．

```{r}
#| eval: false
mc_data <-
  1:mc_num |>                        # ベクトルによる実験回数の指定
  map(\(x)tibble(x = mc_trial())) |> # 単一試行のデータフレームのリスト
  list_rbind()                       # リストを結合 
```

試行が複数の出力を持つ場合は， こちらを利用した整理が簡単な場合もある．

ヒストグラムを作成して結果を視覚化する．

```{r}
#| code-fold: true
#| fig-width: 7
#| fig-height: 7
tibble(x = mc_data) |> # ヒストグラムを出力
  ggplot(aes(x)) + 
  geom_histogram(binwidth = 1,
                 fill = "slateblue", alpha = 0.5, # 塗り潰しの色
                 colour = "slateblue")  # 縁の色
```