MultiRMDにおけるテンプレートの反復利用

T. Yoshiizumi

2017/03/18

　当ページは MultiRMDの簡単なサンプルにつぐ シリーズ第３回です。

　統計ソフトRとmarkdownの組合せによってレポートを生成することが主題です。

　身長および体重を性別によって分類し、その分布をみます。

　サンプルのファイルは mtest02.zip に同梱されています。

　　　　

＜目　次＞

1.　素材データと今回のねらい

　500人分の身長と体重を乱数発生により生成し、データフレームにします。

　性別は女性・男性・無回答の３分類で、これも乱数により決めます。各々おおよそ３分の１ずつになります。

　最初の部分だけ掲げると下のとおり。

ＩＤ 性別 身長 体重
1 男性 174.3 62.6
2 無回答 166.8 60.9
3 女性 169.6 51.1

　　　　

　今回は身長と体重の分布状況をみる訳ですが、
身長のデータも体重のデータも、その扱い方は共通です。

　summaryを適用したり、標準偏差を算出したり、度数分布表を作ったりします。

　それらをレポートする markdownのひな形も身長と体重で共通です。

　一つの markdownテンプレートを反復して利用できる、というのが今回のサンプルで示したいことです。

　前回の「MultiRMDの簡単なサンプル」(サイコロの目の表示)でもテンプレートの反復利用は出てきましたが、今回は少し実践的です。テンプレートの中にchunkが出てきたりもします。

目次に戻る

2.　レポートに盛り込む情報

　素材データから次ぎのような表を作成してレポートに盛り込みます。

　「身長」の例を掲げます。

  最小値 1／4位 中央値 平均値 3／4位 最大値
女性 145.5 159.5 163.8 163.7 167.5 176.4
男性 150.6 162.7 166.4 166.6 170.3 178.6
無回答 151.3 160.1 165.1 164.5 168.2 177.5
全体 145.5 161.5 165.1 165.0 168.8 178.6

　　　　

  平均 標本標準偏差 不偏標準偏差 正規性
女性 163.68 5.70 5.72 Yes(p=0.463762)
男性 166.60 5.64 5.66 Yes(p=0.451549)
無回答 164.49 5.66 5.68 Yes(p=0.463955)
全体 164.97 5.81 5.81 Yes(p=0.284448)

　　　　

  女性 男性 無回答 計
145〜150 3(　1.7) 0(　0.0) 0(　0.0) 3(　0.6)
150〜155 5(　2.9) 4(　2.2) 8(　5.5) 17(　3.4)
…… …… …… …… ……
175〜180 3(　1.7) 14(　7.7) 5(　3.4) 22(　4.4)
合計 174(100.0) 181(100.0) 145(100.0) 500(100.0)

目次に戻る

3.　mtest02.rmd の内容

　chunk定義ファイル mtest02.r を読み込んで、最終的にhtmlを生成するためのrmdファイル(mtest02.rmd)は下のとおり。

−−−−　ここから
---
title: "Test#2 for MultiRMD 〜 身長と体重の分布"
author: "Black Lady"
date: "2017/03/18"
output:
  html_document:
    md_extensions: -ascii_identifiers
    toc: true
---

```{r include=FALSE}
# source("multirmd.r")
chunk.def <- "mtest02.r"
set.mrmd(chunk.def)
read_chunk(chunk.def)
```
```{r child="$$final.rmd"}
```
```{r include=FALSE}
file.remove("$$final.rmd")
```
−−−−　ここまで

　上は、「簡単なサンプル」のときの mtest01.rmd と基本的に同じです。

　前回のサイコロの目を表示するサンプルでは、chunk定義ファイルに含まれているchunkが一つだけでした。diceというchunkです。

　しかし、今回は複数あります(後述)。

　set.mrmd() は $$final.rmd というファイルを書き出しますが、その中では複数のchunkが、定義されている順番で呼び出される形になっています。

　そして、最後に mrmd.buf をテンポラリファイルに出力し、それを実行処理するという流れが $$final.rmd に書き込まれます。

目次に戻る

4.　chunk定義ファイル mtest02.r

　mtest02.r には、chunkの定義(Rプログラム)が書かれています。

　また、それぞれのchunkにおいて、関連のmarkdownテンプレートも書かれています。テンプレートなしのchunkもありますが。

　chunk定義には次のものがあります。

　　　　

　def_func で定義している関数(sd.distr)は、与えられたベクトルから平均, 標本標準偏差, 不偏標準偏差, 正規性を取得して返します。

　正規性の有無は shapiro.test によります。

　set_data というchunkでは、素材データを生成し、念のためそれをcsvファイルとして書き出します。後で検証しなおしたいときなどに利用して下さい。

　女性・男性・無回答で少し違いが出るように、乱数発生を工夫しています。

　以下、markdownテンプレートのあるchunkについて述べます。

　　　　

(1)　first 〜 kableに関する留意点

　firstというchunkでは、身長と体重の分布をみる前に、性別の人数がどんな内訳になっているかを確認します。

　それを「はじめに」として表示します。

　たとえば次のとおい。

−−−−　ここから

はじめに

身長と体重について、性別を手がかりに分類・考察する。

性別の内訳(人数)は下のとおり。

  女性 男性 無回答 総数
人数(人) 185 158 157 500
比率(%) 37.0 31.6 31.4 100.0

−−−−　ここまで

　上記を表示するためのmarkdownテンプレートの中に、表を出力する kable() が出てきます。kableは、Rのライブラリknitrが提供する関数です。

　このkableという関数の使い方についてつまずいたので、それについて記します。

　　　　

整数か小数点数かを揃える

　内訳の人数を把握するのは容易です。

　素材データがデータフレームとして変数dtfに代入されている場合
frq <- table(dtf[,"性別"]) とするだけで、女性、男性、無回答の人数がfrqに入ります。

　その構成比は pct <- prop.table(frq)*100.0 で取得できます。

　「合計」を付け加えるなら frq <- addmargins(frq) で大丈夫です。

　そして、mtx <- cbind(人数=frq, 比率=pct) とすればマトリックスが生成されます。

　この変数mtxをkableに引き渡すと、たとえば次のような表になります。

  人数 比率
女性 145 29.0
男性 173 34.6
無回答 182 36.4

　「人数」の列は整数、「比率」は小数点数になります。何も問題ありません。

　ところが、この縦長の表を横長にしたいとなると少々面倒です。

　kableは、縦方向、つまり列ごとに数値の書式を整えようとするようです。

　マトリックスの行と列を入れ替えること自体は容易です。

　mtx <- t(mtx) とすれば入れ替えることができます。

　しかし、これをkableに引き渡すと、整数であってほしいものが小数になったりします。

　それを回避するには数値を文字に変換します。

　frq <- sprintf("%d", frq); pct <- sprintf("%.1f", pct) とすれば、
見た目でいうと frqが整数、pctは小数になります。実態は文字です。

　　　　

左揃え・右揃えの調整(align)

　さて、今度は表のセル内における左揃えか右揃えかの問題が生じます。

　kableは、数値は右揃えにしてくれますが、文字だと左揃えにします。

　今回、やむなく数値を文字にした訳ですが、左揃えでなく右揃えにしたい場合は、kableのalignオプションを使います。

　kable(mtx, align="rrrr") とすれば、４列すべて右揃えになります。

　４列とは「女性、男性、無回答、総数」の四つです。

　表として出力すると、これら４列よりも左に「人数」　「比率」という行名(rownames)があります。しかし、この行名は、alignオプションの対象外です。

　変数mtxの実態は、あくまで４列です。

　なので、alignの指定は "lrrrr" ではなく "rrrr" です。

　言わずもがなですが、右揃えは r, 左揃えは l, 中央揃えなら c です。

　alignは、align=c(rep("r",4)) のように、ほんとはベクトルで指定するもののようです。ですが、"rrrr" でも大丈夫のようです。

　　　　

　参考まで、kableに引き渡すマトリックス(count.mtx)を生成するRプログラムを下に掲げておきます。

−−−−　ここから
frq <- table(dtf[,"性別"])
pct <- addmargins(prop.table(frq)*100.0)
frq <- addmargins(frq)
frq <- sprintf("%d", frq)
pct <- sprintf("%.1f", pct)
mtx <- cbind(frq, pct)
rownames(mtx) <- c(levels(dtf[,"性別"]), "総数")
colnames(mtx) <- c("人数(人)", "比率(%)")
count.mtx <- t(mtx)
−−−−　ここまで

目次に戻る

(2)　distribution 〜 chunk名の省略

　身長と体重の各々について、summary() により最小値、平均値、最大値などを検出します。

　また、summaryでは得ることができない標準偏差や正規性もみることにします。

　ここのchunkにはmarkdownテンプレートが二つあります。

　template[1] は、本論に入る前の導入の文章です。

　二つ目の template[2] は、summaryの結果を表形式で表示し、また、それとは別に標準偏差等も表形式で表示するフォーマットになっています。

　二つ目のテンプレートには %s という２文字が出てきます。何カ所かにあります。

　この２文字は、「身長」に置き換えたうえで mrmd.buf にストックします。

　そして、同じテンプレートの %s を今度は「体重」に置き換えてから mrmd.buf にストックします。

　つまり、テンプレートを二度使います。

　mrmd.buf の中身は、最後にrmdファイルとして書き出し、htmlに変換します。

　　　　

テンプレート内に書かれたchunkの名前に関する留意点

　テンプレートを反復利用する場合に注意してほしいのは、テンプレート内に出てくるchunkの扱いです。

　template[2] には ```{r ……} のchunk記述が出てきますが、この場合、chunk名を指定するとエラーになります。

　一連の処理の中で、同一のchunk名を持つchunk定義が重複して出てくると、エラーが発生するからです。

　chunk名を指定すると、テンプレートの一回目の利用のときはいいとしても、二回目の利用の際に「既にそのチャンク名は使われています」という趣旨のエラーが発生します。

　テンプレート内にchunk定義を記述する場合は、chunk名を省略するのが無難です。

　　　　

summaryの結果を表示するRプログラムの要点

　distributionのところにあるRプログラムでは、apply関数群を多用しています。

　処理手順を簡潔に記述できて便利ですが、文章で説明するのは筆者の力量を超えるので、ここでは、summaryの結果を表示するRプログラムの要点を列記するにとどめます。

　下に出てくる変数numsは c("身長", "体重")、
keyは "性別"、key.lvlは c("女性", "男性", "無回答") です。

　上の処理で得られた summary.tbl を使うと、
summary.tbl$身長, summary.tbl$体重 によって、表示したいマトリックスを取り出すことができます。それをkable関数に引き渡します。

　上記は中核的な部分を取り出しただけで、実際にはリストやマトリックスのnamesを付け替える処理などが必要になります。

　　　　

標準偏差や正規性を表示するRプログラムの要点

　distributionのところには標準偏差や正規性を取得するRプログラムもあります。

　手順は、前述のsummaryに関するRプログラムと同じです。

　summary関数の代わりにsd.distrという関数を用います。

　この関数は自前で設けました。def_func というchunkに書かれています。

　平均値、標本標準偏差、不偏標準偏差、正規性を取得します。

　正規性は shapiro.test による検定で、そのp値を返します。検定できないときは -1 を返します。

　関数sd.distrの戻り値は変数 sd.tbl に記録します。

　ただ、そのままだと小数点以下の桁数が多いのと、正規性の結果が分かりにくいので、表示用に調整し、変数sd.tbl2に記録します。数値でなく文字にします。

　正規性は、yes, no とp値を表示するようにします。

　kableには sd.tbl2$身長, sd.tbl2$体重 を引き渡します。

目次に戻る

(3)　histgram 〜 度数分布における区分け

　身長と体重の各々について、度数分布の表とヒストグラムを表示します。

　ここのchunkにはmarkdownテンプレートが二つあります。

　template[1] は、本論に入る前の導入の文章です。表に関する注意書きを記しています。

　二つ目の template[2] は、度数分布を表形式で表示し、また、そのヒストグラム(グラフ)も表示するフォーマットになっています。

　度数分布のマトリックスは、変数 frq.tbl$身長, frq.tbl$体重 で取り出せるようにし、それぞれをkableに引き渡します。

　グラフは、身長を image01.png、体重を image02.png に書き出します。

　それを画像表示するためのmarkdown記述は、各々 下のようになります。

![身長のヒストグラム](image01.png)
![体重のヒストグラム](image02.png)

　上の２行を変数 img.list に格納します。

　そして、img.list$身長, img.list$体重 それぞれを template[2] の中で取り込むようにしています。

　以下、Rプログラムについて少し説明します。

　　　　

度数分布の区分け

　グラフを作成するときだけでなく、度数分布表を作る際も hist() を用いていることを少し説明します。

　度数分布表を作る場合、区分けをどうするかが問題になります。

　身長でいえば「160〜165」　「165〜170」といった区分けです。全部で何等分するのがいいか迷います。

　hist() は、この区分けを自動的に行ってくれます。理論上で理想的といえる区分けの方法があるようですが、hist関数は、その方法で区分けします。バランスのとれた区分けということでしょうか。

h <- hist(dtf[,"身長"], include.lowest=T, right=F, plot=F)

　上のようにするとグラフは作成されず、変数hに区分けごとの人数などの情報がセットされます。

　h$breaks にそれぞれの区分けの区切りの値がセットされます。それをたどれば、区分けの下限値と上限値を確認できる訳です。

　「160〜165」とか「165〜170」といった一連の区分けの名前をベクトル化するには次のようにします。変数lblがそのベクトルになります。

brk <- h$breaks
lbl <- sapply(1:(length(brk)-1), function(i)
    sprintf("%s〜%s", brk[i], brk[i+1]))

　　　　

「性別ごと×区分けごと」の人数の算出

　h$breaks によって500人全員の理想的な区分けを得ました。

　次は、その区分けに従って「女性、男性、無回答」それぞれの人数を数え上げる処理が必要になります。

　いろいろなやり方があるとおもいますが、今回は下の手順を採りました。

○　ステップ１

　500人全員に区分けの名前を割り当てます。

　身長163の人なら「160〜165」という名前を割り当てます。

y <- cut(dtf[,"身長"], breaks=brk, labels=lbl, right=FALSE)

　上のようにすれば、区分けの名前 500人分のデータが変数yにセットされます。

○　ステップ２

　上記のyがあれば、table関数で「性別ごと×区分けごと」の人数を算出できます。

x <- dtf[,"性別"]
frq <- table(x, y)

　変数frqは、数値から構成される行列です。各カテゴリーの人数が格納されています。

　あとは frq を基にしてパーセントの値を算出したり、合計の値を追加します。

　詳細はプログラム本体を参照して下さい。

　プログラムの記述は少し煩雑ですが、基本は上に書いたことだけです。

　　　　

グラフの作成

　グラフの作成は、plot, linesといった関数で行いますが、女性、男性、無回答の３グループを一つのグラフの上にプロットする関係で色分けします。

　first というchunkの中で
key.color <- c("red", "blue", "green3", ……)
というのを定義しています。

　これは、女性をred、男性をbulue、無回答をgreen3にするためのものです。

　グラフの作成については下のサイトを参考にさせていただきました。

Rで複数のヒストグラムを比較する方法あれこれ - Rプログラミングの小ネタ

　詳細はプログラム本体を参照して下さい。ずいぶん端折るようで申し訳ありません。

　　　　

　以上でこのページは終了です。

　　次回は MultiRMDにおけるエラー処理および条件分岐 です。