当ページは MultiRMDにおけるテンプレートの反復利用につぐ シリーズ第４回です。
統計ソフトRとmarkdownの組合せによってレポートを生成することが主題です。
身長および体重を男女別に分類し、その有意差を検証します。
その前に、有効回答がどれくらいあるかなどをチェックして、男性のデータしかない(女性のデータがない)といった処理困難な状況なら、Rプログラムを中断します。
サンプルのファイルは mtest03.zip に同梱されています。
400人分のＩＤ、性別、身長、体重を記録した data.csv を読み込んで処理します。
僅かですが、記入のない空欄があります。つまり、欠損値ありです。
今回のRプログラムでは欠損値があることを前提にして、有効回答を抽出するところから始めます。
有効回答がなければ、その旨を掲示するhtmlを出力しますが、その先の分析には進みません(条件分岐)。
なお、付属の dataEmpty.csv には有効回答がありません。
エラー処理については、
読み込むべきcsvファイルが存在しなかったり、あるいは、
性別として男性または女性のデータの一方しかなければ、
プログラムを中断します。男女の比較ができないためです。
ちなみに dataMaleOnly.csv には男性のデータしかありません。
プログラムを中断するときは、error.txt に「なぜ中断するか」を書き出します。
有効回答がある場合は統計的な処理に入ります。
数値データの分布をみる部分は、前回「MultiRMDにおけるテンプレートの反復利用」の mtest02.r と同じです。
それに加えて今回は、身長と体重の各々について男女間で有意差があるかどうかを検証します。
分布に正規性が認められるときは t検定を行い、
そうでなければウィルコクソンの順位和検定を行います。
条件分岐により、t検定用とウィルコクソン検定用のmarkdownテンプレートを使い分けます。
有効回答がどれくらいあるかをチェックした結果は、次のような箇条書きと表にして出力します。
−−−− ここから
| 身長 | 体重 | 総計 | |
|---|---|---|---|
| 男性 | 198( 98.0) | 201( 99.5) | 202(100.0) |
| 女性 | 192( 99.5) | 190( 98.4) | 193(100.0) |
| 無回答 | 2(100.0) | 2(100.0) | 2(100.0) |
| 合計 | 392( 98.7) | 393( 99.0) | 397(100.0) |
−−−− ここまで
有効回答がない場合は、その旨を告げる表示になります。上記とは異なります。
身長と体重の分布状況を表示する部分は、前回の「MultiRMDにおけるテンプレートの反復利用」と同じです。
summary関数の結果、標準偏差と正規性、度数分布表とヒストグラムです。
有意差の検定部分は、正規性が認められるか否か、有意差が認められるか否かを箇条書きで示した後で、t.test または wilcox.test の結果をそのまま記載します。
たとえば下のとおりです。
−−−− ここから
t検定の結果は下のとおり。
##
## Two Sample t-test
##
## data: $身長$男性 and $身長$女性
## t = 10.378, df = 388, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 5.362702 7.869653
## sample estimates:
## mean of x mean of y
## 166.7157 160.0995
−−−− ここまで
rmdファイルは、YAMLヘッダと二つのchunkから構成されます。
一つ目のchunkでは、Rプログラムで共通に用いる変数の値を設定します。
また、関数の定義もあります。
二つ目のchunkでは、まず mtest03_1.r の読み込みと実行を行い、その結果、有効回答があれば mtest03_2.r の読み込みと実行を行います。
今回は chunk定義ファイルが mtest03_1.r, mtest03_2.r の二つです。
rmdファイルの中身は下のとおり。
−−−− ここから
---
title: "Test＃3 for MultiRMD 〜 エラー処理および条件分岐"
author: "`r Sys.getenv('USERNAME')`"
date: "`r format(Sys.time(), '%Y/%m/%d')`"
output:
html_document:
md_extensions: -ascii_identifiers
toc: true
---
```{r include=FALSE}
# source("multirmd.r")
csv.file <- "data.csv"
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("身長", "体重") # 処理対象の数値欄の列名
key.item <- c("男性", "女性") # keyの各要素名。その並び順。
key.former <- c("m", "f") # 素材でのkeyの各要素名。上と同一なら不必要。
key.color <- c("blue", "red", "green3", "cyan", "magenta", "yellow", "gray")
call_chunk <- function(chunk.def) {
set.mrmd(chunk.def)
read_chunk(chunk.def)
child.file <- "$$final.rmd"
eval(call.child)
dummy <- file.remove(child.file)
}
```
```{r echo=FALSE, results='asis', include=TRUE}
call_chunk("mtest03_1.r")
if (valid_c > 0) # 有効回答がある場合
call_chunk("mtest03_2.r")
```
−−−− ここまで
素材データ data.csv では性別の欄が m, f となっています。
男性を m、女性を f とアルファベットで記入しています。
そのままレポートにしたのでは分かりにくいので、それぞれ「男性, 女性」に変更します。Rプログラムの factor() でそれを行います。
key.former が素材データ上での記入文字(m, f)、
key.item は「男性, 女性」です。
dtf[,"性別"] <- factor(dtf[,"性別"], levels=key.former, labels=key.item)
上のようにすれば「m, f」が「男性, 女性」になります。
もし素材データに最初から「男性, 女性」と記入されているなら、key.former という変数の設定は不要です。
性別の欄が空欄だと NA という欠損値になりますが、それには「無回答」を割り当てることになります。
関数 call_chunk() は、multirmdの典型的な処理手順を関数にしたものです。
chunk定義ファイル(chunk.def)を材料として、
set.mrmd(chunk.def) により $$final.rmd というファイルを書き出します。
次に read_chunk() によってchunk定義ファイルを読み込みます。
そして、$$final.rmd をchildとして取り込み(chunkの実行)、
最後に、$$final.rmd を削除します。
eval(call.child) の実態は、
cat(knit_child(child.file, quiet=TRUE), sep='') です。
この call_chunk() を他のchunk内で用いる場合、そのchunkオプションを
echo=FALSE, results='asis', include=TRUE とします。
そうしないと childファイルの取り込みが意図したようにいきません。
rmdファイルの２番目のchunkの中で、call_chunk() を２度呼び出しています。
まず mtest03_1.r を呼び出し(有効回答のチェック)、
有効回答があれば(つまり 変数 valid_c が０より大きければ)
mtest03_2.r をよびだします(分布や有意差のチェック)。
mtest03_1.r には二つのchunkがあります。setup, valid です。
validでは有効回答がどれくらいあるか、裏をかえせば、欠損値がどれくらいあるかをチェックします。
その前に、setupの方でいくつかエラーチェックも行っているのでそれについて触れます。
エラーチェックの要点をあげると下のとおりです。
file.exists("data.csv") が TRUE を返すなら data.csv が存在します。
FALLSE だとファイルがないのでRプログラムを中断させます。
q("no") によって、作業用ファイルを保存することなく中断できます。
単に中断しただけでは「なぜ中断してしまったのか」がユーザーに分からないので、error.txt というファイルにメッセージを書き込みます。
multirmd.r では変数 mrmd.error に "error.txt" を代入しているので、ファイル名として mrmd.error を使えます。
具体的なRプログラムとして掲げると下のとおり。
if (!file.exists(csv.file)) { # csvファイルが存在しない
cat(sprintf("'%s' がみつかりません．\n", csv.file), file=mrmd.error)
eval(call.clear) # 内部処理用のファイル・変数を消去
q("no") # Rプログラムを終了させる
}
eval(call.clear) の call.clear は、multirmd.r の中で定義されています。
$$MrmdTmp.rmd, $$final.rmd というファイルが存在していれば消去し、
mrmd.buf, mrmd.template という変数を空にします。
今回は「性別、身長、体重」の三つの列に注目する訳ですが、
この三つを取り出して別のデータフレームに記録するよう試みます。
その試みが失敗したら、Rプログラムを中断します。
dtf.org <- read.csv(csv.file, header=T, na.strings="")
dtf <- tryCatch({dtf.org[,c(key,nums)]},
error = function(e) {
cat(sprintf("%s: 列名「%s」は妥当ですか？\n", csv.file,
paste(c(key,nums), collapse="，")), file=mrmd.error)
return(NULL)})
if (is.null(dtf)) {
eval(call.clear)
q("no")
}
上のように「わざわざ別のデータフレームに記録する」という方法は、大規模データを取り扱うときは、やらない方がいいのかもしれません。
変数 key.item, key.former が定義されていない場合への対処を盛り込んでいます。
この二つが定義されていなくても、Rプログラムを中断することはしません。
levels() などで対応します。
key.item <- tryCatch({key.item},
error = function(e) {levels(dtf[,key])})
key.former <- tryCatch({key.former},
error = function(e) {key.item})
「性別」の要素として「男性, 女性」の２種類を想定している訳ですが、
有効回答中に「男性」しかない、あるいは「女性」しかないということだと、
この先の分析ができません。なので Rプログラムを中断します。
org_item <- levels(dtf[,key]) # 性別の要素(男性，女性など)を取得
if (length(key.item) != length(org_item)) { # 要素の個数が想定と違う
fmt <- paste0("%s\n「%s」の列の予定された要素 '%s' と\n",
"素材データの実際の要素 '%s' の個数が不一致です．\n")
cat(sprintf(fmt, csv.file, key, paste(key.item, collapse="，"),
paste(org_item, collapse="，")), file=mrmd.error)
eval(call.clear)
q("no")
}
上は、「性別」の要素の個数をチェックしているだけです。中身までは確認していません。
エラーチェックは、おおよそ以上のとおりです。
欠損値の取扱いは、multirmd.r の仕様と直接の関係はありませんが、メモとして記してみます。
素材データであるデータフレーム中で欠損値は NA になっています。
そこから有効回答を取り出す方法の要点を掲げてみます。
「身長，体重」の両方とも空欄だと、今回は分析の対象にならないのでデータフレームから消去します。
具体的には次のようにします。
for (i in nrow(dtf):1) {
if (length(dtf[i,nums][!is.na(dtf[i,nums])]) == 0)
dtf <- dtf[-i,]
}
１行ずつチェックして、二つとも空欄だと消去しています。
こうした１行ずつの処理よりも効率的な方法があるような気がしますが、よく分かりません。
性別の欄には m, f が記入されています。
これを「男性, 女性」に変更するのに factor() をつかいます。
空欄がないなら下のようにします。
dtf[,"性別"] <- factor(dtf[,"性別"], levels=c("m", "f"),
labels=c("男性", "女性"))
空欄があるときは次のとおり。
dtf[,"性別"] <- factor(dtf[,"性別"], levels=c("m", "f", NA),
labels=c("男性", "女性", "無回答"), exclude=NULL)
factor() は、普通は NA を対象外としますが、オプションの exclude=NULL を指定すると、「対象外とするものなし」になって NA も対象に加えるようです。
dtf[,"性別"] に NA があるかどうかは anyNA(dtf[,"性別"]) で分かります。
anyNA() が TRUE を返すなら NA が含まれています。
NA の扱いは関数によって違います。
ベクトルの要素数を返す length() は、NA も一つの要素として数えます。
もし NA を除いた数を得たいなら length(na.omit(x)) のようにします。
平均値を求める mean() は、NA が含まれているベクトルを与えると、NA を返します。
オプションとして na.rm=T を与えると数値を返しますが、同じベクトルから標準偏差や分散も求めるなら、予め NA を除去しておく方がやりやすいとおもいます。
一方、集計をするのに便利な table() は、デフォルトでは NA を対象外とします。factor() と同様です。
table() の場合、useNA(“no”｜”yes”) あるいは exclude というオプションで調整できます。
欠損値をどのように取り扱ったらいいか、実践ではいろいろ気をつかうところです。
validというchunkのところには二つのテンプレートがあります。
有効回答がある場合のものと、ない場合のものです。
有効回答がある場合、その件数や比率を表形式に取りまとめて変数 valid.tbl にセットします。
テンプレートでは、その valid.tbl を kable() に引き渡します。
有効回答がない場合は、読み込んだ件数(行数)が全部で何件か、そのうち有効回答でないものが何件かを箇条書きで表示します。表を示すことはしません。
mtest03_2.r には四つのchunkがあります。
def_func, distribution, histgram, difference の四つです。
difference以外は、前回「MultiRMDにおけるテンプレートの反復利用」の mtest02.r のものと同じです。
differenceでは、身長および体重について男女間で有意差があるかどうかを検証します。
有意差の検証は、t検定または wilcoxの検定で行います。
身長、体重の分布に正規性が認められるときは t検定、そうでなければ wilcoxの検定を用います。
テンプレートには、前置きに当たる文章、t検定の結果を表示するもの、および wilcox検定の結果を表示するものの三つがあります。
正規性の有無は、distributionというchunkで算出した、変数 sd.tbl をみれば分かります。それに応じて t検定用、wilcox検定用のテンプレートを使い分けます。
Rプログラム部分の大まかな流れを記しておきます。
まず memos というリストを設け、データを男女別に分類します。
memos$身長$男性, memos$身長$女性
memos$体重$男性, memos$体重$女性 を設定します。
男女双方とも正規性が認められる場合、
Ｆ検定で等分散か否かを確認した上で t検定を行い、
memos$身長$t.test にその検定結果を代入します。
正規性にこだわりすぎの感じがありますが、そこは目をつぶってください。
正規性に確証がないときは wilcox検定の結果を保存します。
memos$体重$wilcox.test に検定結果をまるごと代入します。
そして、テンプレートの方で memos を参照します。
以上でこのページは終了です。
multirmd.r について説明できていない部分があるので、あと１回くらいは解説を書きたいとおもいます。