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

## Quartoの使い方

### Quarto とは

Quarto は Pandoc を利用した出版システムで，文章とプログラムをまとめて記述した上で，レポートやスライドあるいはWebページのような文書の作成を行うことができる． 同様な働きをするものとしては RStudio では R Markdown という形式が用意されているが，R 以外の言語 (Python や Julia) も利用できるという意味では R Markdown の発展形と捉えることもできる． R のプログラムを主として扱う形式である R Script でもプログラム中にコメントを残すことはできるが，コメントが多いとコードそのものが読み難いといった問題がある．

### Quarto の構造

Quarto は，文書全体の体裁を整えるための YAMLヘッダー，プログラムを記述するチャンク(chunk)，テキスト(markdown文)の3つの領域から構成される．

YAMLヘッダーは文書の先頭に置き，`---` と `---` の間に必要な事項を記述する． この文書のYAMLヘッダーは以下のように書かれている．

``` yaml
---
title: "統計データ解析"
subtitle: "第1講 サンプルコード"
date: "`r Sys.Date()`"
format:
    html:                     # <1>
      toc: true
      html-math-method: katex # <2>
      self-contained: true
      grid: 
        margin-width: 350px   # <3>
execute: 
  echo: fenced                # <4>
reference-location: margin    # <5>
citation-location: margin
tbl-cap-location: margin
fig-cap-location: margin
warning: false                # <6>
editor: visual
---
```

1.  HTML形式の文書を作成することを宣言する．
2.  TeXのレンダリングにKaTeXを利用する．
3.  右側に大き目のマージンを指定する．
4.  チャンク内の情報を全て表示するように指定する．
5.  脚注，文献，表や図の説明をマージンに表示する．
6.  作成した文書中に警告を表示しない．

チャンクは ```` ```{r} ```` (R言語の場合) と ```` ``` ```` の間にコードを記述する． 例えば，後の解析で共通に用いるパッケージ [^1] を読み込むには以下のように記述する．

[^1]: 以下ではこの他に `ggrepel` `broom` `MASS` などパッケージに含まれる関数の一部を利用するが，これらは個別に指定して利用する．

```{r}
library(conflicted) # 名前の衝突の警告を出してもらう
conflicts_prefer(
  dplyr::filter(),
  dplyr::select(),
)
library(tidyverse)  # tidyverseの基本パッケージ群の読み込み
library(ggfortify)  # 分析結果の視覚化にggplotを利用するためのパッケージ
library(GGally)     # ggplotで散布図などを作成するためのパッケージ
library(gt)         # 表を作成するためのパッケージ
#' macOSのための日本語表示の設定
if(Sys.info()["sysname"] == "Darwin") { # macOSか調べる
    theme_update(text = element_text(family = "HiraMaruProN-W4"))
    label_family <- "HiraMaruProN-W4"    
} else {label_family <- NULL}
```

Rのコードはパッケージ `knitr` [^2] を利用して処理される． また，文書全体を通して1つのプロセスが処理し，チャンク間で情報は共有されるため，前のチャンクでの計算結果を後のチャンクでも使うことができる．

[^2]: knitr のオプションと Quarto の関係は [Code Cells: Knitr](https://quarto.org/docs/reference/cells/cells-knitr.html) に詳しく記載されている．

テキストは一般的な markdown 文法で解釈される． 数式はTeXの基本的な記述を理解するので，簡単な数式をそのまま記述することができる． テキスト中の数式は `$ 数式 $` で，独立した数式は `$$ 数式 $$` で記述する．

::: callout-tip
### 数式の記述の例

正規分布は $\mathbb{R}=(-\infty,\infty)$ 上の確率分布で，平均 $\mu$ 分散 $\sigma^{2}$ である分布の密度関数は $$
    \phi(x) = \frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}
$$ となる．
:::

### RStudio での利用

Quarto ファイルを開くと左上のペインに表示される． エディタの左上の `Source/Visual` で表示を切り替えることができる． また `Render` の右の設定で，文書の出力先 `Window/Viewer Pane`，チャンクの実行結果の出力先 `inline/console` を選択することができる．

## 講義で扱う内容

### 単回帰分析

例として `MASS::Animals` を用いる[^3]． 拡張したデータフレーム形式である `tibble` に変換して `bb_data` で参照する．

[^3]: オリジナルのデータは `MASS::Animals` あるいは `data(Animals, package = "MASS")` とすれば `Animals` で参照することができる．

```{r}
bb_data <- MASS::Animals |> rownames_to_column() |> as_tibble()
```

このデータは28種類の動物の平均的な体重と脳の重さを調べたもの[^4]で，データの内容を表にすると以下のようになる．

[^4]: データの詳細は `help(Animals, package = "MASS")` で確認することができる．

```{r}
#| label: tbl-bb
#| tbl-cap: 28種類の陸生動物の平均的な体重と脳の重さ．
bb_data |> gt()
```

横軸に体重，縦軸に脳の重さを対数変換して描画したものが以下の図になる．

```{r}
#| label: fig-bb-loglog
#| fig-cap: 体重と脳の重さの両対数グラフ．
#| fig-width: 7
#| fig-height: 7
bb_p <- bb_data |>
    ggplot(aes(body, brain)) +
    geom_point(colour = alpha("royalblue", 0.75)) + 
    ggrepel::geom_text_repel(aes(label = rowname), # 各点の名前を追加
                             size = 3) + 
    scale_x_log10() + scale_y_log10() + # log-log plot を指定
    labs(title = "Brain and Body Weights",
         x = "body [kg]", y="brain [g]")
print(bb_p) # グラフを表示
```

両対数変換したデータにおいてはほぼ線形の関係が成り立っていることが見て取れるので，その関係式(回帰式)を推定して図示すると以下のようになる．

```{r}
#| label: fig-bb-lm
#| fig-cap: 回帰式とその信頼区間．
#| fig-width: 7
#| fig-height: 7
bb_p <- bb_p + # 回帰式を追加
    geom_smooth(method = "lm", 
                colour = "dodgerblue",
                fill = "dodgerblue")
print(bb_p) # グラフを表示
```

実際の分析においては， 回帰式の妥当性や外れ値などを考察する必要があり， そのための統計的な方法を学ぶ．

### 重回帰分析

例として配布データ `wine.csv` を `readr::read_csv()` 用いて読み込み，以下 `bw_data` で参照する．

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

このデータはボルドー地区の何年かにわたるワインの価格と気候の関係を調べたもの[^5]で，データは以下のとおりである．各列は生産年(VINT)，対数変換した価格(LPRICE2)，冬の降雨量(WRAIN)，育成期平均気温(DEGREES)，収穫期降雨量(HRAIN)，評価までの経過年数(TIME_SV)を意味する．

[^5]: 以下の分析は俗に "Ashenfelter のワイン方程式" と呼ばれる問題で，詳細は[記事](http://www.liquidasset.com/orley.htm)および[データと分析結果](http://www.liquidasset.com/winedata.html)で確認することができる．

```{r}
#| label: tbl-wine
#| tbl-cap: ボルドーワインの価格と気候の関係．
bw_data |> gt()
```

変数の間の関係を散布図として描画すると以下のようになる．

```{r}
#| label: fig-wine-pairs
#| fig-cap: ボルドーワインの価格と気候の関係．
#| fig-width: 7
#| fig-height: 7
bw_data |> 
    ggpairs(columns = 2:6,
            lower = list(continuous = wrap("smooth_loess", colour = "blue"))) +
    theme(axis.title.x = element_text(size = 8), # 文字の大きさを調整
          axis.title.y = element_text(size = 8)) 
```

ワインの価格を冬の降雨量(WRAIN)，育成期平均気温(DEGREES)，収穫期降雨量(HRAIN)，年数(TIME_SV)で説明する回帰式は以下のように推定される．

```{r}
bw_fit <- lm(LPRICE2 ~ . - VINT, # VINTを除く全て
             data = bw_data)
#' 推定結果を表にまとめるためのパッケージ gtsummary を利用
library(gtsummary)
library(broom.helpers)  # gtsummary のいくつかの関数で利用(インストールされていれば不要)
bw_fit |>
    tbl_regression(estimate_fun = label_style_sigfig(digits = 4)) |>
    add_glance_source_note(include = c(r.squared, adj.r.squared))
```

重回帰による予測値と実際の価格を比較すると以下のようになる． 良好な予測が行われていることが視覚的に確認できる．

```{r}
#| label: fig-bw-prediction
#| fig-cap: 重回帰による予測値と実際の価格
#| fig-width: 7
#| fig-height: 7
bw_fit |>
    broom::augment() |>
    ggplot(aes(x = LPRICE2, y = .fitted, label = VINT)) +
    geom_text(na.rm = TRUE) + # text で表示
    geom_abline(slope = 1, colour = "red") +
    labs(title = "Bordeaux Wine Price",
         x = "Price (log)", y = "Prediction")
```

### 主成分分析

県別の生活環境に関するデータ `jpamenity.csv` [^6]を用いる． 補助的な情報なども含めて整理したデータフレームを `ja_data` とする．

[^6]: 配布したデータ `jpamenity.csv` には，[政府統計の総合窓口](https://www.e-stat.go.jp)から取得した以下の25項目が，各県ごとにまとめられている． 昼夜間人口比率（％）,年少人口割合［15歳未満人口］（％）,老年人口割合［65歳以上人口］（％）,人口増減率（％）,粗出生率（人口千人当たり）,粗死亡率（人口千人当たり）,婚姻率（人口千人当たり）,離婚率（人口千人当たり）,高等学校数（15〜17歳人口10万人当たり）（校）,高等学校数（可住地面積100k㎡当たり）（校）,高等学校卒業者の進学率（％）,労働力人口比率［男］（％）,労働力人口比率［女］（％）,完全失業率［男］（％）,完全失業率［女］（％）,交通事故発生件数（人口10万人当たり）（件）,刑法犯認知件数（人口千人当たり）（件）,食料費割合［二人以上の世帯］（％）,住居費割合［二人以上の世帯］（％）,教育費割合［二人以上の世帯］（％）,平均貯蓄率［勤労者世帯］（％）,趣味・娯楽の平均時間［有業者・男］（時間.分）,趣味・娯楽の平均時間［無業者・男］（時間.分）,趣味・娯楽の平均時間［有業者・女］（時間.分）,趣味・娯楽の平均時間［無業者・女］（時間.分）． また `jpamenityitem.csv` には簡略化した項目名が記載されている．

```{r}
## 元のデータに県名と地方名を付加してデータフレームを作成する
ja_data <- bind_cols(
    read_csv(file = "data/prefecture.csv",
             col_select = c(2,4)) |>
    set_names("県名", "地方名"), # 列の名称を"県名"と"地方名"に変更
    read_csv(file = "data/jpamenity.csv",
             col_select = !1:2) |> slice(-1) |>
    set_names(names(read_csv(file = "data/jpamenityitem.csv")))) # 簡略化した項目名に変更
```

データの一部を表形式で表示すると以下のようになる．

```{r}
#| label: tab-ja-data
ja_data |>
    select(1:10) |> slice(1:15) |> # 大きな表なので一部を選択する(列:select，行:slice)
    gt()
```

関連の強い項目ごとに3つのグループに分け， それぞれのグループで散布図を作成する．

```{r}
#| label: fig-ja-pairs1
#| fig-cap: 県別の生活環境(人口動態に関連する項目)
#| fig-width: 7
#| fig-height: 7
ja_data |>
    GGally::ggscatmat(columns = 3:10, color = "地方名", alpha = .5) +
    theme(text = element_text(size = 8)) 
```

```{r}
#| label: fig-ja-pairs2
#| fig-cap: 県別の生活環境(教育・労働などに関連する項目)
#| fig-width: 7
#| fig-height: 7
ja_data |>
    GGally::ggscatmat(columns = 11:19, color = "地方名", alpha = .5) +
    theme(text = element_text(size = 8))
```

```{r}
#| label: fig-ja-pairs3
#| fig-cap: 県別の生活環境(貯蓄・余暇などに関連する項目)
#| fig-width: 7
#| fig-height: 7
ja_data |>
    GGally::ggscatmat(columns = 20:27, color = "地方名", alpha = .5) +
    theme(text = element_text(size = 8))
```

25次元のデータを主成分分析を用いて2次元に縮約し， 各県の特性をバイプロットで図示すると以下のようになる．

```{r}
#| label: fig-ja-biplot
#| fig-cap: 県別の生活環境の主成分分析
#| fig-width: 7
#| fig-height: 7
ja_fit <- ja_data |> 
    column_to_rownames(var = "県名") |>
    select(where(is.double)) |>
    prcomp(scale. = TRUE) # 主成分分析の実行
autoplot(ja_fit, # バイプロット
         data = ja_data,
         shape = FALSE, label = TRUE, loadings = TRUE, loadings.label = TRUE,
         label.family = label_family, loadings.label.family = label_family,
         label.size = 3, loadings.label.size = 3.5) 
```

### 判別分析

例として `MASS::biopsy` [^7] を用いる． 不要な項目を削除し整理したものを `bio_data` とする．

[^7]: データの詳細は `help(biopsy, package = "MASS")` で確認することができる．

```{r}
bio_data <- MASS::biopsy |>
    as_tibble() |>
    na.omit() |> # NA を除く
    select(-ID)  # IDを除く
```

このデータは乳癌の良性・悪性と9種類の生検との関係を700弱の患者について調べたものである． 生検同士の関係と良性・悪性を散布図に図示すると以下のようになる．

```{r}
#| label: fig-bio-pairs
#| fig-cap: 9種類の生検の散布図．
#| fig-width: 7
#| fig-height: 7
bio_data |>
  GGally::ggpairs(diag = list(mapping = aes(colour = class)),
                  lower = list(mapping = aes(colour = class)))
```

以下の図は，9種類の生検の相互の関係を主成分分析を用いて2次元に縮約し，良性・悪性との関係を視覚的に捉えたものである．

```{r}
#| label: fig-bio-pca
#| fig-cap: 主成分分析を用いて2次元に縮約した生検と良性・悪性の関係．
#| fig-width: 7
#| fig-height: 7
autoplot(bio_data |> select(where(is.numeric)) |> prcomp(), # 主成分分析
         data = bio_data,
         colour = "class") +
    theme(legend.position=c(.9,.85)) 
```

線形判別分析では，回帰分析と同様に説明変数(この場合は9種類の生検の値)の重み付けた値(判別関数値)によりクラスの判別を行う． このデータにより作成された判別関数値のクラスごとの分布を表示すると以下のようになる． クラスごとに大きく異なる分布となっていることがわかる．

```{r}
#| label: fig-bio-discriminant
#| fig-cap: 線形判別関数による判別関数値の分布．
#| fig-width: 7
#| fig-height: 7
library(MASS) # パッケージ MASS に含まれる判別分析のための関数を利用
bio_fit <- lda(class ~ ., data = bio_data)
bio_data |>
    mutate(x = predict(bio_fit)$x) |>
    ggplot(aes(x = x)) +
    geom_histogram(aes(fill = class), show.legend = FALSE) +
    facet_grid(class ~ .)
```

### クラスタ分析

データ例として `omusubi.csv` [^8] を用いる． 補助的な情報を付加して `om_data` とする．

[^8]: このデータは [おむすびの日アンケート](https://gohan.gr.jp/result/18/) で公開されていたものの(2009年版)であるが，現在は公開されていないようである．

```{r}
om_data <- bind_cols( 
    read_csv(file = "data/omusubi.csv"), 
    read_csv(file = "data/prefecture.csv")) # 県名・地方名の情報を付加
```

県別の人気比率を図示すると以下のようになる．

```{r}
#| label: fig-om-barplot
#| fig-cap: おむすびの具の県別人気アンケート (2009年度実施)．
#| fig-width: 7
#| fig-height: 7
om_data |>
    select(ume:etc,jp) |>
    set_names(c("梅","鮭","昆布","鰹","明太子","鱈子","ツナ","その他","県名")) |>
    pivot_longer(-県名) |>
    mutate(県名 = fct_rev(as_factor(県名)),
           name = as_factor(name)) |> 
    ggplot(aes(y = 県名, x = value)) +
    geom_bar(aes(fill = name),
             stat = "identity",
             position = position_stack(reverse=TRUE)) +
    labs(x = "人気比率", fill = "具材") +
    theme(axis.text.y = element_text(size = 9))
```

比率は確率分布と見ることができるので，県同士の人気比率の距離を分布間の距離の一つである Hellinger 距離 [^9] を用いて測り，階層的クラスタリングを行う． 結果をデンドログラムとして表示すると以下のようになる．

[^9]: 比率を表すベクトルを1/2乗してEuclid距離を計算すればHellinger距離に比例した量が得られる．

```{r}
#| label: fig-om-dendrogram
#| fig-cap: おむすびの具人気アンケート結果のもとづく階層的クラスタリング．県別の人気比率のHellinger距離を群平均法により分析した．
#| fig-width: 7
#| fig-height: 7
library(cluster)  # クラスタ分析のためのパッケージ
library(ggdendro) # ggplot でデンドログラムを描くためのパッケージ
om_agnes <- 
	om_data |> 
    select(ume:etc,jp) |>
    set_names(c("梅","鮭","昆布","鰹","明太子","鱈子","ツナ","その他","県名")) |>
    column_to_rownames(var = "県名") |> 
    sqrt() |> 
    agnes()
om_agnes |> as.dendrogram() |>
    ggdendrogram(rotate = TRUE, theme_dendro = FALSE) +
    labs(x = "県名", y = "距離") +
    theme(axis.text.y = element_text(size = 9))
```

### 時系列解析

例として `datasets::AirPassengers` を用いる[^10]．

[^10]: `datasets` のデータは起動時に読み込まれでるので，特に準備はいらない．データの詳細は `help(AirPassengers)` で確認することができる．

```{r}
library(tsibble)    # 時系列を扱うためのパッケージ
library(feasts)     # ggplotで時系列を扱うための拡張パッケージ
library(fable)      # 
```

時系列のためのデータ形式になっており，Console 上ではそのまま適切に表示されるが，表形式にするには若干工夫が必要となる．

```{r}
#| label: tbl-ap
#| tbl-cap: 1949年から1960年における米国における月ごとの国際線乗客数の変遷．
ap_tsbl <- AirPassengers |> as_tsibble() # 時系列向きの形式に変換
ap_tsbl |> 
    mutate(Year = year(index),
           Month = month(index, label = TRUE)) |>
    as_tibble() |> select(!index) |> # 
    pivot_wider(names_from = Month) |>
    gt() |>
    tab_header(title = "Monthly Airline Passenger Numbers",
               subtitle = "1949-1960")
```

横軸に時間をとり，その変遷を図示すると以下のようになる．

```{r}
#| label: fig-ap-timeseries
#| fig-cap: 米国の国際線乗客数の変遷．
#| fig-width: 7
#| fig-height: 7
ap_tsbl |>
    autoplot(value) +
    labs(title = "AirPassengers",
         x = "Year", y = "Passengers/1000")
```

ほぼ線形のトレンドと増大する分散変動があることがわかる． 分散変動を安定化するために以下では対数変換したデータを扱う． また，時系列の予測の精度を評価するために，訓練データを試験データに分割する．

```{r}
ap_train <- ap_tsbl |> filter_index(~ "1958-12") # 訓練データ
ap_test  <- ap_tsbl |> filter_index("1959-01" ~ .) # 試験データ(2年分)
```

階差を取ることによって定常化できるかどうか検討する．

```{r}
#| label: fig-ap-difference
#| fig-cap: 対数変換と階差による時系列の定常化の検証．
#| fig-width: 7
#| fig-height: 7
ap_train |> 
    gg_tsdisplay(difference(log(value)),
                 plot_type = "partial") 
```

対数変換した訓練データでSARIMAモデルを推定し，試験データの予測を行う．

```{r}
#| label: fig-ap-prediction
#| fig-cap: SARIMAモデルの推定と予測．
#| fig-width: 7
#| fig-height: 7
ap_fit <- ap_train |> model(sarima = ARIMA(log(value))) # SARIMAモデルの推定(次数は自動推定)
ap_fit |> 
    forecast(h = nrow(ap_test)) |>
    autoplot(ap_train) +
    autolayer(ap_test, value, colour = "purple") +
    labs(title = "Forecast from SARIMA model",
         x = "Year", y = "Passengers/1000")
```

## RStudio の基本的な使い方

### コンソール上での計算

与えられた式の計算を実行してみる．

```{r}
#' 1 x 2 + 3^2 の計算
1 * 2 + 3^2
#' sin(2π) の計算
sin(2*pi)
#' √2 + |-0.6| の計算
sqrt(2) + abs(-0.6)
```

::: callout-note
演算子の前後に空白があっても良い(空白は可読性のために適宜利用する)．

2文字以上の演算子の途中には空白は入れられない(例: 羃乗の別記法の"\*\*")．

πは `pi` として定義されている．

`e-16` は10の-16乗の意味である．

計算結果は数値誤差のため厳密に0にならないことがある． 利用環境下で扱うことのできる数値については ".Machine" に保存されている． 特に上記の丸め誤差については ".Machine\$double.eps" を参照せよ．
:::

### 関数の実行

引数が1つの関数の計算例は以下のようになる．

```{r}
#' 正弦関数(引数が1つ)の計算例
sin(x = pi/2) 
sin(pi/2) # 引数名は省略でき，前の行とこの行は同じ結果になる
```

以下は擬似コードなので， `a, b` を適当な数値に置き換えて実行しなさい．

```{r}
#| eval: false
#' 対数関数(引数が2つ)の計算例
log(a, b) # 底を b とする a の対数 
log(x = a, base = b) #上と同値
log(base = b, x = a) #上と同値(引数名があれば順序は自由に変えられる)
log(b, a) # = log(x=b,base=a) (引数名がなければ規定の順序で解釈される)
log(a) # 自然対数 =log(a,base=exp(1))
```

### ヘルプの使い方

ヘルプを参照する例は以下のようになる．

```{r}
#| eval: false
#' 関数 log() に関するヘルプの例
help(log) # Helpタブに結果は表示される
?log # 上と同値
example(log) # ヘルプ内の例を実行
help.search("log") # "log"に関連する項目は？
??"log" # 上と同値
```

### オブジェクトの操作

オブジェクトの代入の例は以下のようになる．

```{r}
#' 数値を変数 foo に代入する
(foo <- 3) # foo <- 3; print(foo) と等価
#' 変数 foo を用いて計算し，結果を bar に代入する
bar <- sin(2/3*pi) + cos(foo * pi/4) # 計算結果は表示されない
#' 変数 bar の内容を表示する
print(bar)
```

### 自作関数の定義

関数を定義するには関数 `function()` を用いる．

```{r}
#' 縦の長さ a, 横の長さ b (既定値は1) の長方形の面積
foo <- function(a, b = 1){ 
  out <- a * b 
  return(out) # 計算結果を外に返却
}
```

上記で定義した関数の実行例は以下のようになる．

```{r}
foo(2, 3) # foo(a = 2, b = 3) と同義
foo(2) # foo(a = 2, b = 1) と同義
```

関数名を付けない関数(無名関数)を作ることもできる．

```{r}
#' 変数や関数を定義して計算する方法
(x <- 1:10/10)
foo <- function(x) sin(x)/x # 式の計算結果を返却(returnを省略可)
foo(x)
#' 変数や関数を定義せずに計算する方法
(function(x) sin(x)/x)(1:10/10)
(\(x) sin(x)/x)(1:10/10) # R 4.1 以降の短縮表現
```

無名関数の典型的な使い方を例示する．

```{r}
#' 関数 1/(1+x^2) を区間 [0,1] で積分 (関数オブジェクトを渡す書き方)
f <- function(x) 1/(1+x^2) # 関数を定義
integrate(f, 0, 1) # pi/4 
#' 関数 e^(-x^2/2) を実軸全体で積分(Gauss積分) (関数を定義しない書き方)
integrate(\(x)exp(-x^2/2), -Inf, Inf) # sqrt(2*pi)
```