intro_to_tidyverse_ppt_pdf.Rmd

---
title: "Rによる探索的データ分析入門"
author: "発電基盤開発課 高津一誠"
date: "2018年10月5日"
header-includes:               # LaTeXの設定
  - \usepackage{xltxtra}       # Unicode処理するために必要なもの
  - \usepackage{zxjatype}      # XeLaTeX で日本語の標準的な組版を行う
  - \usepackage[ipa]{zxjafont} # IPAフォントを使う
output:
  beamer_presentation: 
    latex_engine: xelatex      # TeX種類を指定
    keep_tex: yes              # TeXソースを残す
    toc: true                  # 目次作成
    fig_caption: true          # 図キャプション作成
    # number_sections: false     # 章番号作成
documentclass: beamer
classoption: xelatex,ja=standard
geometry: no
mainfont: "IPAPGothic"
monofont: "VL Gothic"
#mainfont: "VL PGothic"
#monofont: "VL Gothic"
---

```{r setup, include=FALSE}
library(knitr)
library(ggpubr)
library(cowplot)
library(tidyverse)
library(broom)

# グラフ中の日本語の文字化け対策
#   出力フォーマットがPDF(TeX)の場合のみ実行される
if (knitr::opts_knit$get("rmarkdown.pandoc.to") %in% c("beamer", "latex")) {
  options(device = function(file, width = 7, height = 7, ...) {
    cairo_pdf(tempfile(), width = width, height = height, ...)
  })
}
knitr::opts_chunk$set(dev="cairo_pdf", dev.args=list(family="VL Gothic"))

# Rコード部の設定
knitr::opts_chunk$set(echo = TRUE, fig.align="center", warning = F)

# ggplotグラフで日本語を使えるようにする
old <- theme_set(theme_bw(base_family = "VL Gothic"))

# baseグラフで日本語を使えるようにする
par(family = "VL Gothic")

# 表出力用に、数値の有効桁数を限定する
options(digits = 5)
```

## データ処理は探索的プロセス

実験データ分析は、探索的なプロセスです。
1回データ処理して終了することはあまりなく、多くの場合は仮説→検証を繰り返して適切な結果を得ることができます。
ここではRの可視化処理を使って、その簡単な例を示します。

![探索的プロセス](data_analysis_process.png)


## 実験データを読み込む

ここでは読み込みの手順は省略し、Rに組み込みのテスト用データを使うことにします。
使用するデータは生態学の計測データで、アヤメの花弁（petal）とガク（sepal）の長さと幅を、3つの種について50個体ずつ計測したデータです。

```{r, echo=FALSE}
iris <- iris %>% as_tibble()
print(iris, n = 5)
```

## 分析しやすいよう整形する

元々のデータでは部位ごとの計測項目が変数（列）になっていましたが、部位を変数にした方が扱いやすいと思います。
そこで、以下のように整形します。

```{r tidyr, echo=TRUE}
iris_long <- iris %>% 
  rownames_to_column("id") %>% 
  mutate(id = as.integer(id)) %>% 
  gather(key, value, matches("Length|Width")) %>% 
  separate(key, into = c("Part", "amount")) %>% 
  spread(amount, value)
```
```{r, echo=FALSE}
print(iris_long, n = 5)
```

## 可視化する

それでは、データを確認するために可視化してみましょう。
長さと幅の関係を種ごと部位ごとに確認してみることにします。

```{r 1st_plot, echo=TRUE, eval=FALSE}
iris_long %>% 
	ggplot(aes(x = Width, y = Length)) + 
	geom_point(aes(color = Part)) + 
	stat_smooth(method = "lm", color = "gray 40") + 
	facet_grid(Part ~ Species) +
	labs(x = "幅(cm)", y = "長さ(cm)", color = "部位")
```

## 可視化する

比較的よい相関があるようですが、よく見てみると、どの種でもガクの方が大きいようです。  

```{r 1st_plot, echo=FALSE, eval=TRUE, fig.height=3, fig.width=5, out.height='60%', fig.cap="最初の可視化"}
```

## 仮説その1

もしかすると、長さと幅の関係は花弁とガクで共通と考えたほうがいいのかもしれません。
確認してみましょう。

```{r theory_1, echo=TRUE, eval=FALSE}
iris_long %>% 
	ggplot(aes(x = Width, y = Length)) + 
	geom_point(aes(color = Part)) + 
	stat_smooth(method = "lm", color = "gray 40") + 
	facet_wrap(~ Species) +
	labs(x = "Width(cm)", y = "Length(cm)")
```

## 仮説その1

よい相関があるので、仮説は妥当だったようです。
更に見てみると、近似直線の傾きがどれも似ているようです。

```{r theory_1, echo=FALSE, eval=TRUE, fig.height=3, fig.width=5, out.height='60%', fig.cap="長さと幅の関係"}
```

## 仮説その2

長さと幅の傾きは種が異なっても共通かもしれません。
確認してみましょう。

```{r theory_2, echo=TRUE, eval=FALSE}
iris_long %>% 
  ggplot(aes(x = Width, y = Length, color = Species)) + 
	geom_point(aes(shape = Part), size = 3) + 
	stat_smooth(method = "lm") +
	labs(x = "Width(cm)", y = "Length(cm)")
```

## 仮説その2

95%信頼区間（グレーの領域）を考慮すると、同じ傾きである可能性は高そうです。

```{r theory_2, echo=FALSE, eval=TRUE, fig.height=3, fig.width=5, out.height='60%', fig.cap="傾きの比較"}
```

## 線形モデルの結果を数値で取得する

いままではグラフで確認していましたが、モデル化の結果を数値で取得することもできます。
係数は推定値と標準誤差がestimateとstd.errに、無相関のt検定の結果はstaticとp.valueに示されています。

```{r model, echo=TRUE}
lm_coef <- iris_long %>% 
  group_by(Species) %>% 
  summarise(list(lm(Length ~ Width) %>% tidy())) %>% 
  unnest()
```
```{r, echo=FALSE, fig.height=2, fig.width=5.5, out.height='30%', fig.cap="線形モデル係数"}
lm_coef %>% ggtexttable(rows = NULL, theme = ggpubr::ttheme("mBlue")) %>% 
  table_cell_bg(color = "red", linewidth = 2, column = 3, row = 3) %>%
  table_cell_bg(color = "red", linewidth = 2, column = 3, row = 5) %>% 
  table_cell_bg(color = "red", linewidth = 2, column = 3, row = 7)
```

## 線形モデルの結果を数値で取得する

ここから95%信頼区間を求めるには、以下のように直接計算することもできますし、

```{r, echo=TRUE}
ci_lwr <- function(coef, err, n, p = 0.95) {
  coef -1 * err * qt(df = n - 2, p = 1 - (1 - p)/2)
}

ci_upr <- function(coef, err, n, p = 0.95) {
  coef +1 * err * qt(df = n - 2, p = 1 - (1 - p)/2)
}

lm_coef_ci <- lm_coef %>% 
  filter(term == "Width") %>% 
  group_by(Species) %>% 
  summarise(ci_lwr = ci_lwr(estimate, std.error, 100),
            ci_upr = ci_upr(estimate, std.error, 100))
```

## 線形モデルの結果を数値で取得する

以下のように信頼区間を求める関数を使うこともできます。

```{r, echo=TRUE}
lm_coef_ci <- iris_long %>% 
  group_by(Species) %>% 
  summarise(list(lm(Length ~ Width) %>% 
                   confint(., "Width") %>% as_tibble())) %>% 
  unnest()
```

その結果は以下のようになり、傾きが同一である可能性があるといえます。

```{r, echo=FALSE, fig.height=1.2, fig.width=2.2, out.height='40%', fig.cap="傾きの信頼区間"}
lm_coef_ci %>% ggtexttable(rows = NULL, theme = ggpubr::ttheme("mBlue"))
```

## 線形モデルの妥当性確認

最後に、直線回帰が適切に行えているのか、残差の分布を確認してみましょう。

```{r, echo=TRUE, eval=TRUE}
lm_resid <- iris_long %>% 
  group_by(Species) %>% 
  summarise(list(lm(Length ~ Width) %>% augment())) %>% 
  unnest()

p1 <- lm_resid %>% 
  ggplot(aes(x = Width, y = .resid)) + 
  geom_point(alpha = 0.4) + 
  facet_wrap(~ Species, scales = "free_x")

p2 <- lm_resid %>% 
  ggplot(aes(x = .resid)) + 
  geom_histogram(bins = 8) + 
  facet_wrap(~ Species, scales = "free_x")
```

## 線形モデルの妥当性確認

説明変数（Width）に対して残差のバラツキはほぼ均等のようですし、残差の分布も偏っていないようなので、問題ないでしょう。

```{r check_lm, echo=FALSE, eval=TRUE, fig.height=3, fig.width=4, out.height='60%', fig.cap="残差プロット"}
plot_grid(p1, p2, ncol = 1)
```