統計データ解析

第1講 サンプルコード

Quartoの使い方

Quarto とは

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

Quarto の構造

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

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

---
title: "統計データ解析"
subtitle: "第1講 サンプルコード"
date: "2024-10-02"
format:
1    html:
      toc: true
2      html-math-method: katex
      self-contained: true
      grid: 
3        margin-width: 350px
execute: 
4  echo: fenced
5reference-location: margin
citation-location: margin
tbl-cap-location: margin
fig-cap-location: margin
6warning: false
editor: visual
---

1

HTML形式の文書を作成することを宣言する．

2

TeXのレンダリングにKaTeXを利用する．

3

右側に大き目のマージンを指定する．

4

チャンク内の情報を全て表示するように指定する．

5

脚注，文献，表や図の説明をマージンに表示する．

6

作成した文書中に警告を表示しない．

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

¹ 以下ではこの他に 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 ² を利用して処理される． また，文書全体を通して1つのプロセスが処理し，チャンク間で情報は共有されるため，前のチャンクでの計算結果を後のチャンクでも使うことができる．

² knitr のオプションと Quarto の関係は Code Cells: Knitr に詳しく記載されている．

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

数式の記述の例

正規分布は \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 を用いる³． 拡張したデータフレーム形式である tibble に変換して bb_data で参照する．

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

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

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

⁴ データの詳細は help(Animals, package = "MASS") で確認することができる．

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

	body	brain
Mountain beaver	1.350	8.1
Cow	465.000	423.0
Grey wolf	36.330	119.5
Goat	27.660	115.0
Guinea pig	1.040	5.5
Dipliodocus	11700.000	50.0
Asian elephant	2547.000	4603.0
Donkey	187.100	419.0
Horse	521.000	655.0
Potar monkey	10.000	115.0
Cat	3.300	25.6
Giraffe	529.000	680.0
Gorilla	207.000	406.0
Human	62.000	1320.0
African elephant	6654.000	5712.0
Triceratops	9400.000	70.0
Rhesus monkey	6.800	179.0
Kangaroo	35.000	56.0
Golden hamster	0.120	1.0
Mouse	0.023	0.4
Rabbit	2.500	12.1
Sheep	55.500	175.0
Jaguar	100.000	157.0
Chimpanzee	52.160	440.0
Rat	0.280	1.9
Brachiosaurus	87000.000	154.5
Mole	0.122	3.0
Pig	192.000	180.0

Table 1: 28種類の陸生動物の平均的な体重と脳の重さ．

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

```{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) # グラフを表示
```

Figure 1: 体重と脳の重さの両対数グラフ．

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

```{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) # グラフを表示
```

Figure 2: 回帰式とその信頼区間．

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

重回帰分析

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

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

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

⁵ 以下の分析は俗に “Ashenfelter のワイン方程式” と呼ばれる問題で，詳細は記事およびデータと分析結果で確認することができる．

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

VINT	LPRICE2	WRAIN	DEGREES	HRAIN	TIME_SV
1952	-0.99868	600	17.1167	160	31
1953	-0.45440	690	16.7333	80	30
1954	NA	430	15.3833	180	29
1955	-0.80796	502	17.1500	130	28
1956	NA	440	15.6500	140	27
1957	-1.50926	420	16.1333	110	26
1958	-1.71655	582	16.4167	187	25
1959	-0.41800	485	17.4833	187	24
1960	-1.97491	763	16.4167	290	23
1961	0.00000	830	17.3333	38	22
1962	-1.10572	697	16.3000	52	21
1963	-1.78098	608	15.7167	155	20
1964	-1.18435	402	17.2667	96	19
1965	-2.24194	602	15.3667	267	18
1966	-0.74943	819	16.5333	86	17
1967	-1.65388	714	16.2333	118	16
1968	-2.25018	610	16.2000	292	15
1969	-2.14784	575	16.5500	244	14
1970	-0.90544	622	16.6667	89	13
1971	-1.30031	551	16.7667	112	12
1972	-2.28879	536	14.9833	158	11
1973	-1.85700	376	17.0667	123	10
1974	-2.19958	574	16.3000	184	9
1975	-1.20168	572	16.9500	171	8
1976	-1.37264	418	17.6500	247	7
1977	-2.23503	821	15.5833	87	6
1978	-1.30769	763	15.8167	51	5
1979	-1.53960	717	16.1667	122	4
1980	-1.99582	578	16.0000	74	3
1981	NA	535	16.9667	111	2
1982	NA	712	17.4000	162	1
1983	NA	845	17.3833	119	0
1984	NA	591	16.5000	119	-1
1985	NA	744	16.8000	38	-2
1986	NA	563	16.2833	171	-3
1987	NA	452	16.9833	115	-4
1988	NA	808	17.1000	59	-5
1989	NA	443	NA	82	-6

Table 2: ボルドーワインの価格と気候の関係．

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

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

Figure 3: ボルドーワインの価格と気候の関係．

ワインの価格を冬の降雨量(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))
```

Characteristic	Beta	95% CI 1	p-value
WRAIN	0.0012	0.0002, 0.0022	0.024
DEGREES	0.6164	0.4190, 0.8138	<0.001
HRAIN	-0.0039	-0.0055, -0.0022	<0.001
TIME_SV	0.0238	0.0090, 0.0387	0.003
R² = 0.828; Adjusted R² = 0.796
1 CI = Confidence Interval

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

```{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")
```

Figure 4: 重回帰による予測値と実際の価格

主成分分析

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

⁶ 配布したデータ jpamenity.csv には，政府統計の総合窓口から取得した以下の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()
```

県名	地方名	昼夜人口比	年少人口比	老年人口比	人口増減率	粗出生率	粗死亡率	婚姻率	離婚率
北海道	北海道	100.0	11.7	26.0	-0.47	7.09	10.63	4.86	2.12
青森県	東北	100.0	12.1	27.0	-0.95	6.79	12.81	4.33	1.78
岩手県	東北	99.7	12.4	27.9	-0.84	7.12	12.33	4.32	1.52
宮城県	東北	100.2	13.0	22.9	-0.09	8.05	9.51	5.30	1.70
秋田県	東北	99.9	11.1	30.7	-1.12	6.16	13.98	3.78	1.41
山形県	東北	99.8	12.6	28.3	-0.78	7.13	12.81	4.24	1.46
福島県	東北	99.6	12.9	26.1	-1.41	7.02	11.94	4.73	1.64
茨城県	関東	97.2	13.2	23.8	-0.51	7.78	10.20	4.92	1.79
栃木県	関東	99.1	13.2	23.2	-0.40	8.02	10.43	5.13	1.85
群馬県	関東	99.9	13.4	24.9	-0.45	7.49	10.63	4.64	1.77
埼玉県	関東	88.6	13.0	22.0	0.07	7.90	8.20	5.10	1.86
千葉県	関東	89.5	12.8	23.2	-0.31	7.89	8.59	5.19	1.86
東京都	関東	118.4	11.3	21.3	0.26	8.12	8.25	6.75	1.91
神奈川県	関東	91.2	13.0	21.5	0.10	8.32	7.94	5.68	1.85
新潟県	中部	100.0	12.5	27.2	-0.64	7.45	11.97	4.35	1.37

関連の強い項目ごとに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)) 
```

Figure 5: 県別の生活環境(人口動態に関連する項目)

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

Figure 6: 県別の生活環境(教育・労働などに関連する項目)

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

Figure 7: 県別の生活環境(貯蓄・余暇などに関連する項目)

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) 
```

Figure 8: 県別の生活環境の主成分分析

判別分析

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

⁷ データの詳細は 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)))
```

Figure 9: 9種類の生検の散布図．

以下の図は，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)) 
```

Figure 10: 主成分分析を用いて2次元に縮約した生検と良性・悪性の関係．

線形判別分析では，回帰分析と同様に説明変数(この場合は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 ~ .)
```

Figure 11: 線形判別関数による判別関数値の分布．

クラスタ分析

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

⁸ このデータは おむすびの日アンケート で公開されていたものの(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))
```

Figure 12: おむすびの具の県別人気アンケート (2009年度実施)．

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

⁹ 比率を表すベクトルを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))
```

Figure 13: おむすびの具人気アンケート結果のもとづく階層的クラスタリング．県別の人気比率のHellinger距離を群平均法により分析した．

時系列解析

例として datasets::AirPassengers を用いる¹⁰．

¹⁰ 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")
```

Monthly Airline Passenger Numbers
1949-1960
Year	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
1949	112	118	132	129	121	135	148	148	136	119	104	118
1950	115	126	141	135	125	149	170	170	158	133	114	140
1951	145	150	178	163	172	178	199	199	184	162	146	166
1952	171	180	193	181	183	218	230	242	209	191	172	194
1953	196	196	236	235	229	243	264	272	237	211	180	201
1954	204	188	235	227	234	264	302	293	259	229	203	229
1955	242	233	267	269	270	315	364	347	312	274	237	278
1956	284	277	317	313	318	374	413	405	355	306	271	306
1957	315	301	356	348	355	422	465	467	404	347	305	336
1958	340	318	362	348	363	435	491	505	404	359	310	337
1959	360	342	406	396	420	472	548	559	463	407	362	405
1960	417	391	419	461	472	535	622	606	508	461	390	432

Table 3: 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")
```

Figure 14: 米国の国際線乗客数の変遷．

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

```{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") 
```

Figure 15: 対数変換と階差による時系列の定常化の検証．

対数変換した訓練データで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")
```

Figure 16: SARIMAモデルの推定と予測．