---
title: "第2講 R言語速習"
subtitle: "統計データ解析"
date: "`r Sys.Date()`"
format:
    html: 
      toc: true
      html-math-method: katex
      self-contained: true
      grid: 
        margin-width: 350px
execute: 
  echo: true
  warning: false
reference-location: margin
citation-location: margin
tbl-cap-location: margin
fig-cap-location: margin
editor: visual
---

## 準備

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

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

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

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

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

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

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

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

``` r
library(conflicted)
```
:::

## データフレームの作成

::: callout-note
### 参考

データフレームは 同じ長さのベクトル(関数 `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}
(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)))
```

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

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

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

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

```

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

-   <https://style.tidyverse.org/index.html>
-   <https://google.github.io/styleguide/Rguide.html>
:::

## ファイルの読み書き

::: callout-note
### 参考

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

```{r}
#| eval: false
write_csv(x, file = "ファイル名") # データフレームxをファイルに書き出す
```

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

```{r}
#| eval: false
y <- read_csv(file = "ファイル名") # 変数yにCSVファイルの内容を読み込む
```
:::

### 問題

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

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

### 解答例

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

```{r}
write_csv(grade_data, file = "data/grade_data.csv")
```

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

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

```{r}
(grade_copy <- read_csv(file = "data/grade_data.csv"))
```

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

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

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

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

```{r}
(pcr_colnames <- names(pcr_data)) # colnames(pcr_data) でも良い
```

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

::: column-margin
略称はそれぞれ以下の意味．

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

```{r}
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
```

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

```{r}
(pcr_data <- mutate(pcr_data, date = date(date)))
```

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

``` r
#' 関数 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型に変更
```
:::

## データフレームの操作

::: callout-note
### 参考

行の選択には関数 `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}
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) も可
```

## データフレームの集約

::: callout-note
### 参考

列の集計には関数 `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) の平均値．

### 解答例

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

```{r}
pcr_data |> summarise(across(!date, max))
```

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

```{r}
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)))
```

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

```{r}
pcr_data |>
    filter(year(date) == 2021) |>
    group_by(month(date)) |>
    summarise(across(!date, max))
```

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

```{r}
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))
```

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

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

## 基本的なグラフの描画

::: callout-note
### 参考

グラフの描画には 関数 `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}
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}
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}
#| fig-width: 7
#| fig-height: 7
pcr_data |>
    select(c(ai,univ,mi)) |> # 描画対象の列を抽出
    GGally::ggpairs(columnLabels = pcr_colnames[c("ai","univ","mi")]) # ラベルを渡す
```

## 擬似乱数

### 問題

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

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

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

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

### 解答例

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

```{r}
#' 関数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の場合と同じ結果
```

## 双六ゲーム

### 問題

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

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

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

### 解答例

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

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

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

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

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

```{r}
for(i in 1:10) print(mc_trial())
```

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

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

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

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

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

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

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