第３回で行った演習の解答例です。必ずしもこのコードが最適という訳ではありませんので参考程度に利用してください。なお、勉強会当日説明とコードを変えてあるものもあります。
　
Rの関数はヘルプで記載されている引数の順番通りに引数を指定した場合は引数名を記述する必要はありませんが、あえて明示的に引数名を記述している場合もあります。また、関数の前に記述しているパッケージ名はパッケージが読み込まれている場合には記述必要はありませんが、追加パッケージに限り、どのパッケージの関数かが分かるようにあえて記述してあります。
　

演習１ メトリクス統計分析から

メトリクス統計分析入門の演習1.4からの問題です。生産性に加えて、工数予実割合も加味してプロジェクトの評価をしようと考えています。以下の条件、方法に従い、各プロジェクトをA～Eの5段階で評価してください（詳細はテンプレート参照方）。

  • データ分布の確認は省略してください（外れ値の判定も不要です）
  • 分析手順はテンプレートに記載してあります
    • 欠損値があるレコード（プロジェクト）は外してください
  • 分からない点は周りのメンバーと協力してチャレンジしてみてください

　

解答例

この問題のポイントはdplyrパッケージの基本的な処理方法を知る点にあります。

# データの読み込み(for Ubuntu)
x <- "../data/data_utf8.csv" %>%
  read.csv(encoding = "UTF-8") %>%
  # Rは日本語の変数名（列名）を扱うのが不得手なので変数名（列名）が日本語の
  # 場合には英数字に変換しておくことを推奨（または日本語は''で括る）
  dplyr::rename(pj = 'プロジェクト名', prod = '生産性', rate = '工数予実割合') %>% 
  # 欠損値（NA）があれば無条件に削除
  tidyr::drop_na()

# 平均値と標準偏差を計算
y <- x %>% 
  # 関数名を""で括ると変数名に入るので識別しやすくなる
  dplyr::summarise_if(is.numeric, c("mean", "sd"), na.rm = TRUE)

# 評価換算のための基礎数値を設定
breaks <- c(-Inf, -1.5, -0.5, 0.5, 1.5, Inf)    # 区切る幅
labels <- c("E", "D", "C", "B", "A")            # 因子の順番に注意

# Zスコアと偏差値、ランクの計算
x %>% 
  # scale関数の返り値はマトリクス型なのでベクトル型に変換しておく
  dplyr::mutate(prod_z = as.vector(scale(prod, y$prod_mean, y$prod_sd)),
                rate_z = as.vector(-scale(rate, y$rate_mean, y$rate_sd)),
  # mutate関数は関数内での計算結果をすぐに参照できる
                z = (prod_z + rate_z)/2,
                ss = z*10 + 50,
  # cut関数は階級計算のための関数
                rank = cut(z, breaks = breaks, labels = labels)) %>% 
  # 数字の丸めはどこで実行しても可
  dplyr::mutate(prod_z = round(prod_z, 2), rate_z = round(rate_z, 2),
                z = round(z, 2), ss = round(ss)) %>% 
  # 表示の際に変数名（列名）を日本語に戻すとわかりやすくなる
  dplyr::rename('プロジェクト' = pj, '生産性' = prod, '工数予実割合' = rate,
                '生産性Z値' = prod_z, '工数予実割合Z値' = rate_z,
                'Zスコア' = z, '偏差値' = ss, '評価' = rank)

　

Tips

summarise関数

summarize_if関数は条件に合致する変数（列）に対して指定の計算を行います。以下の例では、数値の列だけを対象にするため.predicateにis.numeric関数を指定しています。関数名が渡せれば良いので指定の際に()は必要ありません。.funsには計算したい関数または数式を指定します。ベクトル型にすることで複数の指定を一度にできます。なお、c関数はベクトル変数を作成するための関数です。

dplyr::summarise_if(x,
                    .predicate = is.numeric,
                    .funs = c("mean", "sd"), na.rm = TRUE)

.funsで複数の関数を指定する場合""で括らないくても指定は可能ですが、変数名（列名）が長くなる点に留意してください。

dplyr::summarise_if(x,
                    .predicate = is.numeric,
                    .funs = c(mean, sd), na.rm = TRUE)

関数指定に違和感がある場合は''で指定するとコードの記述に統一感が出ます。

dplyr::summarise_if(x,
                    .predicate = 'is.numeric',
                    .funs = c('mean', 'sd'), na.rm = TRUE)

　

scale関数

scale関数のcenterとscaleは自動的に計算されますので指定しなくても構いません。計算結果はマトリクス型の返り値の属性として格納されています。

scale(x$prod) %>% head()
##             [,1]
## [1,]  0.05133555
## [2,] -0.38567905
## [3,]  0.72946167
## [4,]  0.41730838
## [5,]  0.10946065
## [6,]  1.03730940
scale(x$prod, center = y$prod_mean, scale = y$prod_sd) %>% head()
##             [,1]
## [1,]  0.05133555
## [2,] -0.38567905
## [3,]  0.72946167
## [4,]  0.41730838
## [5,]  0.10946065
## [6,]  1.03730940

　

cut関数

cut関数はbreaksで指定する階級のどこに入るかを計算する関数です。labelsを指定しない場合にはどの階級なのかを示す(a, b]（\(a \lt x \le b\)の意味）で表示されます。

cut(x$rate, breaks = c(-Inf, -1.5, -0.5, 0.5, 1.5, Inf))
##  [1] (0.5,1.5]  (0.5,1.5]  (0.5,1.5]  (0.5,1.5]  (0.5,1.5]  (0.5,1.5] 
##  [7] (0.5,1.5]  (0.5,1.5]  (0.5,1.5]  (0.5,1.5]  (0.5,1.5]  (1.5, Inf]
## [13] (1.5, Inf] (0.5,1.5]  (0.5,1.5]  (1.5, Inf] (1.5, Inf] (1.5, Inf]
## [19] (0.5,1.5]  (0.5,1.5]  (1.5, Inf] (1.5, Inf] (0.5,1.5]  (1.5, Inf]
## [25] (0.5,1.5]  (1.5, Inf] (1.5, Inf] (0.5,1.5]  (1.5, Inf] (0.5,1.5] 
## [31] (0.5,1.5]  (0.5,1.5]  (1.5, Inf] (0.5,1.5]  (1.5, Inf] (0.5,1.5] 
## [37] (0.5,1.5]  (0.5,1.5]  (1.5, Inf] (1.5, Inf] (1.5, Inf] (1.5, Inf]
## [43] (0.5,1.5]  (1.5, Inf] (0.5,1.5]  (-0.5,0.5] (1.5, Inf] (1.5, Inf]
## [49] (0.5,1.5]  (1.5, Inf] (0.5,1.5]  (1.5, Inf]
## Levels: (-Inf,-1.5] (-1.5,-0.5] (-0.5,0.5] (0.5,1.5] (1.5, Inf]

この結果を集計すると度数分布表になります。
labelsを指定すると指定したラベル（文字型）で階級が置き換えられます。

cut(x$rate, breaks = c(-Inf, -1.5, -0.5, 0.5, 1.5, Inf),
    labels = c("不可", "可", "良", "優", "秀"))
##  [1] 優 優 優 優 優 優 優 優 優 優 優 秀 秀 優 優 秀 秀 秀 優 優 秀 秀 優
## [24] 秀 優 秀 秀 優 秀 優 優 優 秀 優 秀 優 優 優 秀 秀 秀 秀 優 秀 優 良
## [47] 秀 秀 優 秀 優 秀
## Levels: 不可 可 良 優 秀

　

演習２ アンスコムのデータ例

anscombeのデータ例を整然データに変換してください。

　

解答例

この演習のポイントは典型的な縦横変換（ロング-ワイド変換）の際に行を識別するための情報が必要であるという点です。

anscombe %>%
  # 行を識別できるようにIDを付与しておきます
  tibble::rownames_to_column("id") %>%
  # id以外の変数（列）をkeyとvalueに渡します
  tidyr::gather(key = "key", value = "value", -id) %>%
  # key例のxn, ynを文字と数値に分割します
  tidyr::separate(key, c("axis", "group"), 1) %>%
  # axisとvalueを変数（列）と展開します
  tidyr::spread(axis, value) %>% 
  # 付与したIDは不要なので削除しておきます
  dplyr::select(-id) %>% 
  # groupで並べ替えしておきます
  dplyr::arrange(group)

　

Tips

rownames_to_column関数

rownames_to_column関数は行名を変数（列）として扱えるようにします。変数（列）は文字型になりますので、数値型（整数型）として扱いたい場合にはrowid_to_column関数を用いてください。なお、varで指定するのは変数名（列名）です。

anscombe %>% tibble::rownames_to_column(var = "ID")
anscombe %>% tibble::rowid_to_column(var = "ID")

　

gather関数

gather関数は指定した変数（列）の名前をkeyで指定する変数（列）へ、その値をvalueで指定する列にまとめる関数です。変数（列）の指定には-を用いて「除外する」指定方法もあります。なお、keyとvalueを指定する際に新しくできる変数（列）は""で括ることが推奨されています。

anscombe %>% 
  tibble::rowid_to_column(var = "ID") %>% 
  tidyr::gather(key = "key", value = "value", x1, x2, x3, x4, y1, y2, y3, y4)

　

spread関数

spread関数はkeyで指定した変数（列）の値を変数名（列名）にvalueで指定した変数（列）の値を該当するkeyの変数名（列名）の中身に振り分ける関数です。なお、gather関数と異なりkeyとvalueを指定する際に既存の変数（列）は""で括らないことが推奨されています。

anscombe %>% 
  tibble::rowid_to_column(var = "ID") %>% 
  tidyr::gather(key = "key", value = "value", x1, x2, x3, x4, y1, y2, y3, y4) %>% 
  tidyr::spread(key = key, value = value)

　

演習３ データフレームの制限

irisデータセットを以下のような形に変換してください。変換後の行数に注意してください。

　

解答例

この演習のポイントはデータフレーム型の制限を知る点にあります。

(petal <- iris %>% 
  dplyr::select(Petal.Length, Petal.Width, Species) %>% 
  tidyr::gather(key = "Petal", value = "Petal.value", -Species) %>% 
  tibble::rowid_to_column("ID"))
(sepal <- iris %>% 
  dplyr::select(Sepal.Length, Sepal.Width, Species) %>% 
  tidyr::gather(key = "Sepal", value = "Sepal.value", -Species) %>% 
  tibble::rowid_to_column("ID"))
(petal %>% 
  dplyr::left_join(sepal, by = c("ID")) %>% 
  dplyr::select(-ID))

　

演習４ クロス集計

mtcarsデータセットを用いてgear（ギア数）とcyl（シリンダ数）でmpg（燃費）の平均値をクロス集計してください。

　

解答例

この演習のポイントはクロス集計（ピボットテーブル）の作り方を知る点にあります。

mtcars %>% 
  # gearとcylでグループ化します
  dplyr::group_by(cyl, gear) %>% 
  # グループ単位でmpgの平均値を求めます
  dplyr::summarise(mpg = mean(mpg, na.rm = TRUE)) %>% 
  # 横に広げてクロス集計表の形にします
  tidyr::spread(key = cyl, value = mpg)

Tips

group_by

group_by関数は指定の変数（列）でグループ化する関数です。任意の数の変数（列）を指定することが可能ですが、spread関数では一つの変数（列）しか展開できませんのでクロス集計表にする場合は以下のような形になります。

mtcars %>% 
  dplyr::group_by(cyl, gear, am) %>% 
  dplyr::summarise(mpg = mean(mpg, na.rm = TRUE)) %>% 
  dplyr::mutate(mpg = round(mpg, digits = 2)) %>% 
  tidyr::spread(key = cyl, value = mpg)

　

count関数

この演習では統計量のクロス集計を行いましたが度数集計を行う場合はcount関数を用います。count関数はgroup_by関数とsummarize関数の機能を持ちあわせた関数ですのでgroup_byする必要はありません。

mtcars %>% 
  dplyr::count(cyl, gear) %>% 
  tidyr::spread(key = cyl, value = n)
# 上記と等価のコード
mtcars %>% 
  dplyr::group_by(cyl, gear) %>% 
  dplyr::summarize(n = n()) %>% 
  tidyr::spread(key = cyl, value = n)

　

演習５ 時系列の集計

オープン・クローズチャートを作成するために2017年に起票されたRedmineのチケットを週単位でオープンチケットとクローズチケットに分けて集計してください。チケットデータは“data”フォルダ内にあります。
　
集計条件

  • 2017年に起票されたチケット（開始日が2017/1/1から2017/12/31）が対象
  • 集計は週次（2017年の第1週から第53週）
  • オープンチケットは開始日作成日で、クローズチケットは終了日で集計
  • 累計も計算する
  • データがない週も考慮する

変数名は日本語ですが英語に変更しておいた方が処理の記述が楽です。

　

解答例

この演習のポイントは今までの演習で学んだことを応用する点にあります。

(redmine <- "../data/redmine_data_utf8.csv" %>% 
  readr::read_csv() %>% 
  # 日本語変数名は英数文字に変換しておきます
  dplyr::select(ID = '#', status = 'ステータス', open = '作成日',
                close = '終了日') %>% 
  # 日時のデータを日付のデータに変換しておきます
  dplyr::mutate(open = lubridate::as_date(open),
                close = lubridate::as_date(close)))
# Openチケットの数を数えます
(open <- redmine %>% 
  # 集計対象の2017年のチケットのみに絞ります
  dplyr::filter(open >= "2017-1-1" & open <= "2017-12-31") %>% 
  # チケットがオープンになった日の週番号を求めます
  dplyr::mutate(week = lubridate::week(open)) %>% 
  dplyr::group_by(week) %>%
  # フラグの数を用いて集計します
  dplyr::summarise(open = n()))
# Closedチケットの数を数えます
(close <- redmine %>% 
  # ステータスが"Closed"のものだけを対象とする
  dplyr::filter(status == "Closed") %>% 
  # 集計対象の2017年にオープンしたチケットのみに絞ります
  dplyr::filter(open >= "2017-1-1" & open <= "2017-12-31") %>% 
  # 更にその中で2017年にクローズしたチケットのみに絞ります
  dplyr::filter(close >= "2017-1-1" & close <= "2017-12-31") %>% 
  # チケットがクローズになった日の週番号を求めます
  dplyr::mutate(week = lubridate::week(close)) %>% 
  # 週番号でグループ化して
  dplyr::group_by(week) %>%
  # フラグの数を用いて集計します
  dplyr::summarise(close = n()))
# 週次の集計
# 週番号のデータフレームを作ります
data.frame(week = seq(1:53)) %>% 
  # 週番号を元にチケットを集計したデータを結合します
  dplyr::left_join(open, by = "week") %>% 
  dplyr::left_join(close, by = "week") %>% 
  # NAは0(zero)に変換しておきます（この場合はチケットがない=0(zero)が成り立つので）
  dplyr::mutate(open = tidyr::replace_na(open, 0),
                close = tidyr::replace_na(close, 0)) %>% 
  # オープン数とクローズ数の累計を計算します
  dplyr::mutate(cumopen = cumsum(open), cumclose = cumsum(close)) %>% 
  # 表示のために変数名を日本語に変更します（selectは順番の変更ができる）
  dplyr::select('週番号' = week, 'オープン数' = open, 'オープン数累計' = cumopen,
                'クローズ数' = close, 'クローズ数累計' = cumclose)

　

Tips

cumulative関数

cumulative関数は累計（累和）、累積などを計算するウィンドウ関数と言われるものです。詳細はチートシート のVector Functions項を参照してください。

# 累計（累和）
cumsum(1:10)
##  [1]  1  3  6 10 15 21 28 36 45 55
# 累積
cumprod(1:10)
##  [1]       1       2       6      24     120     720    5040   40320
##  [9]  362880 3628800

　

演習６ 発想の転換

ggplot2::mpgデータセットをclassとtransでクロス集計し、classごとにtransの比率を求めてください。

　

解答例

この問題のポイントは複数列に同一操作を行いたい場合にはtidy data形式にすることで一括で処理できるようになる点にあります。

ggplot2::mpg %>% 
  # クロス集計して度数を求めます
  dplyr::count(class, trans) %>% 
  # 一度、展開して全項目にデータが入るようにします（NAで埋める）
  tidyr::spread(key = class, value = n) %>%
  # tidy data形式に変換します
  tidyr::gather(key = "key", value = "value", -trans) %>%
  # NAを一括で0に変換します
  dplyr::mutate(value = tidyr::replace_na(value, replace = 0L)) %>% 
  # 比率を求めるために再度、展開します
  tidyr::spread(key, value) %>% 
  # 変数（列）ごとに比率を求めます
  dplyr::mutate_if(is.numeric, prop.table) %>% 
  # 再びtidy data形式に変換し
  tidyr::gather(key = "key", value = "value", -trans) %>% 
  # 0をNAに一、その他は小数点以下2桁に丸めるます
  dplyr::mutate(value = dplyr::if_else(value == 0, NA_real_, round(value, 2))) %>% 
  # 再び展開
  tidyr::spread(key = key, value = value)

　

演習７ 演習１改

演習１において欠損値があるレコード（行）を削除せずにプロジェクトを評価してください。Zスコアは生産性のZスコアと工数予実割合のZスコアの平均値としますが、欠損値がある場合は片方のZスコアの値とします。
　

解答例

この演習のポイントはNAの扱い方を知る点にあります。

# データの読み込み(for Ubuntu)
x <- "../data/data_utf8.csv" %>%
  read.csv(encoding = "UTF-8") %>%
  dplyr::rename(pj = 'プロジェクト名', prod = '生産性', rate = '工数予実割合')

# 平均値と標準偏差の計算
y <- x %>% 
  dplyr::summarise_if(is.numeric, c("mean", "sd"), na.rm = TRUE)

# 評価換算のための基礎数値
breaks <- c(-Inf, -1.5, -0.5, 0.5, 1.5, Inf)
labels <- c("E", "D", "C", "B", "A")            # 因子の順番に注意

# Zスコアと偏差値、ランクの計算
x %>% 
  dplyr::mutate(prod_z = as.vector(scale(prod, y$prod_mean, y$prod_sd)),
                rate_z = as.vector(-scale(rate, y$rate_mean, y$rate_sd)),
                z = ifelse(is.na(prod_z), rate_z, (prod_z + rate_z)/2),
                z = ifelse(is.na(rate_z), prod_z, (prod_z + rate_z)/2),
                ss = z*10 + 50,
                rank = cut(z, breaks = breaks, labels = labels)) %>% 
  dplyr::mutate(prod_z = round(prod_z, 2), rate_z = round(rate_z, 2),
                z = round(z, 2), ss = round(ss)) %>% 
  dplyr::rename('プロジェクト' = pj, '生産性' = prod, '工数予実割合' = rate,
                '生産性Z値' = prod_z, '工数予実割合Z値' = rate_z,
                'Zスコア' = z, '偏差値' = ss, '評価' = rank)

　

注意 ここから先は

以降の演習はdplyrやtidyrの基本動作を理解した上で演習してください。
　

　

演習８ 層別の計算

irisデータセットを用いて以下のように品種、弁毎に四分位値を計算してください（ググればコードは分かると思いますが、かなりの難問。だけど覚えると便利）。

　

解答例

この演習のポイントはmutate関数では処理できない複数の返り値がある場合の処理方法を知る点にあります。

iris %>% 
  tidyr::gather(part, value, -Species) %>% 
  dplyr::group_by(Species, part) %>% 
  dplyr::do(qt = quantile(.$value)) %>%
  cbind(do.call(rbind, .$qt)) %>%
  dplyr::select(-qt)

　

演習９ 層別の計算２

irisデータセットを用いて回帰モデル（Sepal.Length ~ Sepal.Width）を層別に求めてみましょう。回帰モデルの返り値は四分位値よりも複雑なので要約にはbroomパッケージを用います。

　

解答例

この演習のポイントは演習８より複雑なリスト型の返り値を処理する方法を知る点にあります。

iris %>% 
  dplyr::group_by(Species) %>% 
  dplyr::do(lm_res = lm(Sepal.Length ~ Sepal.Width, data = .)) %>% 
# coefficientsを表示します
  broom::tidy(lm_res)

　

Tips

予測値や残差などを表示する

予測値（fitted value）や残差などを表示する場合にはaugument関数をもちいます。

iris %>% 
  dplyr::group_by(Species) %>% 
  dplyr::do(lm_res = lm(Sepal.Length ~ Sepal.Width, data = .)) %>% 
# 予測値、残差などを表示する場合
  broom::augment(lm_res)

　

モデル評価などを表示する

回帰モデル自体の評価結果を表示する場合にはglance関数をもちいます。

iris %>% 
  dplyr::group_by(Species) %>% 
  dplyr::do(lm_res = lm(Sepal.Length ~ Sepal.Width, data = .)) %>% 
# モデル評価値などを表示する場合
　 broom::glance(lm_res)

　

演習１０ purrr

演習９をpurrrパッケージを用いて解いてみましょう。purrrパッケージの使い方はググってみましょう （dplyr::do関数は将来的に廃止される可能性があるらしいので参考問題です）。
　
余談：「超難問」の超の使い方はあっている？
　　　可能性は「高い低い」、「大きい小さい」、「あるない」？

　

解答例

この演習のポイントは将来的にはシュリンクする可能性があると噂あれているdplyr::do関数の後継となるpurrrパッケージの有用性を知ることにあります。

# purrrを使う場合
iris %>%
  # 種別ごとにリスト化します
  split(.$Species) %>%
  # 種別ごとに回帰モデルを計算します
  purrr::map(~ lm(Sepal.Length ~ Sepal.Width, data = .x)) %>%
  # 結果をデータフレーム型にします
  purrr::map_dfr(broom::tidy, .id = "Species")

　

License

CC BY-NC-SA 4.0, Sampo Suzuki