---
title: "統計データ解析I"
subtitle: "第4講 練習問題 解答例"
date: "`r Sys.time()`"
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
editor_options: 
  chunk_output_type: console
---

## 準備

以下で利用する共通パッケージを読み込む．

```{r}
library(conflicted)  # 関数名の衝突を警告
conflicts_prefer(    # 優先的に使う関数を指定
    dplyr::filter(),
    dplyr::select(),
    dplyr::lag(),
    )
library(tidyverse)
```

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

### 問題

`datasets::airquality` に対して以下の条件を満たすデータを取り出しなさい．

-   7月のオゾン濃度 (`Ozone`)

-   風速 (`Wind`) が時速10マイル以上で，かつ気温 (`Temp`) が華氏80度以上の日のデータ

-   オゾン (`Ozone`) も日射量 (`Solar.R`) も欠測 (`NA`) でないデータの月 (`Month`) と日 (`Day`)

### 解答例

`datasets::airquality` を `tibble` 形式に変換する．

```{r}
aq_tbl <- as_tibble(airquality) # tibble 形式へ変換
```

7月のオゾン濃度のデータを取り出す．

```{r}
aq_tbl |>
    filter(Month == 7) |> # Month %in% 7 などとしても良い
    select(Ozone)
```

風速時速10マイル以上かつ気温が華氏80度以上のデータを取り出す．

```{r}
aq_tbl |>
    filter(Wind >= 10 & Temp >= 80)
```

日射量が欠測でないデータの月と日を取り出す．

```{r}
aq_tbl |>
    filter(!is.na(Ozone) & !is.na(Solar.R)) |>
    select(Month,Day) # 書いた順に並ぶ
aq_tbl |>
    filter(!is.na(Ozone) & !is.na(Solar.R)) |>
    select(Month:Day) # もともと並んでいるので c(Month,Day) と同じ
```

## データの読み込みと操作

### 問題

サンプルデータ (`jpdata`) を読み込んで操作しなさい．

::: callout-note
<https://www.e-stat.go.jp> より取得したデータ (文字コード : utf8)

-   jpdata1.csv : 県別の対象データ
-   jpdata2.csv : 対象データの内容説明
-   jpdata3.csv : 県と地域の対応関係

(地域から探す / 全県を選択 / 項目を選択してダウンロード)
:::

-   作業ディレクトリ内にサブディレクトリ `data` を作成し，その中に置いて読み込む．

    ```{r}
    #| eval: false
    jp_data <- read_csv(file = "data/jpdata1.csv")
    jp_item <- read_csv(file = "data/jpdata2.csv")
    jp_area <- read_csv(file = "data/jpdata3.csv")
    ```

    日本語に問題がある場合は英語版を読み込む．

    ```{r}
    #| eval: false
    jp_data_en <- read_csv(file = "data/jpdata1-en.csv")
    jp_area_en <- read_csv(file = "data/jpdata3-en.csv")
    ```

### 解答例

ファイルからデータを読み込む．

```{r}
jp_data <- read_csv(file = "data/jpdata1.csv")
jp_item <- read_csv(file = "data/jpdata2.csv")
jp_area <- read_csv(file = "data/jpdata3.csv")
```

::: callout-note
文字コードは関数 `readr::guess_encoding()` により推測される．

```{r}
#| eval: false
guess_encoding("data/jpdata1.csv")
```
:::

以下にデータ操作の例を示す．

```{r}
#' 項目名の内容を確認
jp_item
#' 最初の10県を表示
head(jp_data, n = 10) 
#' 最後の10県を表示
tail(jp_data, n = 10)
#' 地方名を確認
View(jp_area) # 左上ペインに表として表示される
#' 人口の最大値を見る
with(jp_data, max(人口))
#' どの都道府県の人口多いか調べる ('?which.max' を参照)
with(jp_data, which.max(人口)) # 行番号が返る
jp_data[with(jp_data, which.max(人口)),]
#' 行を抽出するには関数 dplyr::filter() を使うのが簡単
jp_data |> filter(人口 == max(人口))
#' 人口の多い順に都道府県を並べる ('?order' を参照)
with(jp_data, order(人口, decreasing = TRUE)) # 行番号を取得
jp_data[with(jp_data, order(人口, decreasing = TRUE)),] # 並べ替え
with(jp_data, 県名[order(人口, decreasing = TRUE)]) # 県名のみを表示
#' 行を並べ替えるには関数 dplyr::arrange() を使うのが簡単 
jp_data |> arrange(desc(人口)) |> print(n = Inf) # ('?dplyr::desc' を参照)
```

上記と同じ内容を英語版のデータで実行する場合は以下のようになる．

```{r}
jp_data_en <- read_csv(file="data/jpdata1-en.csv")
jp_area_en <- read_csv(file="data/jpdata3-en.csv")
#' データを確認する
head(jp_data_en, n = 10) 
tail(jp_data_en, n = 10)
View(jp_area_en)
#' 人口の最大値を見る
with(jp_data_en, max(population))
#' 人口の最大値を与える県の番号を見る
with(jp_data_en, which.max(population))
#' どの都道府県の人口多いか調べる
jp_data_en |> filter(population == max(population))
#' 人口の多い順に都道府県を並べる
jp_data_en |> arrange(desc(population)) |> View()
```

## データの集計

### 問題

以下の指示に従ってサンプルデータ (`jpdata`) を整理しなさい．

-   県別の人口密度を求めよ．

-   地方別の人口密度を求めよ．

    ::: callout-warning
    県ごとに人口が異なるので単純に人口密度を平均してはいけない．
    :::

-   地方別の1000人当たりの婚姻・離婚数を概算せよ． なお，データの記述では「人口1000人当たり」とあるが，この「人口」とは若年層は婚姻不可として除いた「婚姻可能な人口1000人当たり」と考えて計算しなさい．

### 解答例

県別の人口密度 (人口/面積) を求める．

```{r}
with(jp_data, 人口/面積) # 値のみ返す
```

都道府県別の人口密度を求める．

```{r}
jp_data |>
    select(県名, 人口, 面積) |>  # 県名，人口，面積を選択
    mutate(人口密度 = 人口/面積) # 人口密度を計算
```

地方別の人口密度 (地方の総人口/地方の総面積) を求める．

```{r}
jp_data |>
    mutate(地方 = as_factor(jp_area[["地方"]])) |> # 地方の情報を付加
    group_by(地方) |>                              # 地方ごとにグループ化
    summarize((across(c(人口, 面積), sum))) |>     # 地方ごとに人口と面積を集計
    mutate(人口密度 = 人口/面積)                   # 地方別の人口密度を計算
```

地方別の婚姻・離婚数(1000人あたり)を求める．

```{r}
jp_data |>
    mutate(婚姻可能 = 人口-若年,                   # 婚姻可能な人口を推計 
           婚姻数 = 婚姻可能*婚姻/1000,            # 婚姻数/1000人から人数を推計
           離婚数 = 婚姻可能*離婚/1000,            # 離婚数/1000人から人数を推計
           地方 = as_factor(jp_area[["地方"]])) |> # 地方の情報を付加
    group_by(地方) |>                              # 地方ごとにグループ化
    summarize((across(婚姻可能:離婚数, sum))) |>   # 地方別の合計を計算
    mutate(婚姻 = 婚姻数/婚姻可能*1000,            # 1000人あたりの婚姻数を計算
           離婚 = 離婚数/婚姻可能*1000)            # 1000人あたりの離婚数を計算
```

::: callout-tip
日本語の扱いで不具合がある場合は

-   jp_data -\> jp_data_en
-   人口 -\> population
-   若年 -\> young_population
-   婚姻 -\> marriage
-   離婚 -\> divorce

など書き換えて試してみて下さい．
:::