群馬大学 | 医学部 | サイトトップ | Software Tips | RについてのTips

層別に分析するための丁寧な説明

このページは，メールでいただいた下記質問への回答として，Rで層別に分析するための方法と，ついでに関数定義の方法を，できるだけ丁寧に説明することを目的とする。なお，ここでいう層別とは，性別，年齢階級別のように，何らかの条件によって層ごとに分けて別々に分析するという意味である。

なお，データ入力から読み込みまでについては「データ入力から読み込みまでのできるだけ丁寧な説明」をご覧いただきたい。

最終更新： 2004年 4月 15日 (木曜日) 07時40分

質問（一部書き換えています）

例えば，男性だけの身長の平均と標準偏差，女性だけの身長の平均と標準偏差を計算させたいときは，Ｒでは，どのようにそれらを分割する命令を入力すればよいのでしょうか？　あるいは，年齢で４０才以上と未満にわけるとか？

それと，よく使用する命令を，何らかのコマンドで置き換えて記憶させておくことはできるのでしょうか？

データの例

10人の対象者についての身長と体重のデータが次の表のように得られているとする。

対象者ID	身長(cm)	体重(kg)	性別	年齢
1	170	70	M	54
2	162	50	F	34
3	166	72	M	62
4	170	75	M	41
5	164	55	F	37
6	159	62	F	55
7	168	80	F	67
8	183	78	M	47
9	157	47	F	49
10	185	100	M	45

これを１行目の変数名をPID, HT, WT, SEX, AGEとしてタブ区切りテキスト形式で保存したものがstsample.txtである（ちなみにOpenOffice.orgのcalc形式のデータはstsample.sxc [6.6 KB]）。これを，例えば"d:¥stsample.txt"としてダウンロードした場合，まずRのプロンプトに対して，

dat <- read.delim("d:/stsample.txt")

と打って，データをdatというデータフレームに付値することで，分析の準備が完了する。以下の手順を実行するためには，予めここまでやっておかねばならない。

データフレームの一部だけの解析をする

Rでは，変数名の後に"[ ]"で条件設定をすることで，変数の一部だけを分析することが可能である。質問にあった，男性だけの身長の平均と標準偏差（言うまでもないが念のために書いておくと，もちろん不偏標準偏差である）を出したければ，

cat("mean=",mean(dat$HT[dat$SEX=='M']),"sd=",sd(dat$HT[dat$SEX=='M']),"¥n")

とすればよい。しかし，同じ条件でたくさんの変数の一部だけの解析をしたいときに，いちいち[dat$SEX=='M']とつけるのは面倒だろう。そういう場合は，省力化のために関数定義をしてしまおう。

cmeansd <- function(X,C) {
cat("mean=",mean(X[C]),"sd=",sd(X[C]),"¥n") }

としておけば，次からは，

cmeansd(dat$HT,dat$SEX=='M')
cmeansd(dat$HT,dat$SEX=='F')

などととするだけで，男女別に身長の平均と標準偏差を表示することができる。表示するのでなく，次のように平均と標準偏差の値を返すような関数定義にすれば，

cmeansd.noprint <- function(X,C) {
list(mean=mean(X[C]),sd=sd(X[C])) }

得られた結果を別の関数で使うこともできるし，次のようにprintに渡せば，そのまま表示することもできる。

print(cmeansd.noprint(dat$HT,dat$SEX=='M'))

条件設定は一致することを意味する"=="だけでなく，不等号も使えるし，is.na()などの関数も使えるので，40歳以上だけについて身長の平均値と標準偏差を計算したければ，さっき定義したcmeansd()関数を使えばいいのだが，どうせだからNも表示するようにcmeansd()関数を再定義することにすると，

cmeansd <- function(X,C) {
cat("N=",length(X[C]),"¥t mean=",mean(X[C]),"¥t sd=",sd(X[C]),"¥n") }
cmeansd(dat$HT,dat$AGE>=40)

とできるし，"&"（かつ）や"|"（または）を使ってこれらの条件を組み合わせることもできるので，40歳以上の男性または30歳未満の女性についてとしたければ，やはり再定義後のcmeansd()関数を使えば，

cmeansd(dat$HT,(((dat$AGE>=40)&(dat$SEX=='M'))|((dat$AGE<30)&(dat$SEX=='F'))))

とすればよい。条件式も変数に付値できるので，例えば40歳以上と未満をそれぞれ出したければ（注："!"は論理式の否定を意味する），

overforty <- (dat$AGE>=40)
cmeansd(dat$HT,overforty)
cmeansd(dat$HT,!overforty)

とすればよい。

（注）なお，"[ ]"の中に数字を入れると，その順番のオブジェクトを参照することもできる。forループなどで用いる。

層別の分析をする

Rには，実は分類変数によって層別に任意の関数を適用する関数が用意されている。tapply()とかby()である。平均値と標準偏差を返す関数meansd()を定義してから，性別にそれを適用させるには，

meansd <- function(X) { list(mean=mean(X),sd=sd(X)) }
tapply(dat$HT,dat$SEX,meansd)
by(dat$HT,dat$SEX,meansd)

とすればよい。以上のプログラムをstratify.R，データファイル[stsample.txt]と同じディレクトリにおいてバッチモードで実行した結果のテキストファイルをstratify.outとしてリンクしておく。

（おまけ）t検定をカプセル化する

性別にサンプルサイズと平均値と標準偏差を一気に表示することができるなら，性差をみるためにt検定まで関数に含めてしまいたいのが人情だろう。２群間で分散に差が無ければ通常のt検定でいいが，差がある場合はWelchの検定をしなくてはいけないので，等分散性の検定の結果から場合分けして，有意水準を５％とするなら，次の関数で一気にt検定の結果として表に載せるのに必要な情報が表示される（欠損値がある場合も動作するようにis.na()関数を使って条件設定している）。ついでなので層別箱ヒゲ図も表示するように関数定義した。関数定義と，上のデータで身長の性差を検定するプログラムがautottest.R，データファイル[stsample.txt]と同じディレクトリにおいてバッチモードで実行した出力結果テキストファイルがautottest.outである（下図は同時に描かれる層別箱ヒゲ図）。

リンクと引用について

Correspondence to: minato@phi.med.gunma-u.ac.jp.

対象者ID	身長(cm)	体重(kg)	性別	年齢
1	170	70	M	54
2	162	50	F	34
3	166	72	M	62
4	170	75	M	41
5	164	55	F	37
6	159	62	F	55
7	168	80	F	67
8	183	78	M	47
9	157	47	F	49
10	185	100	M	45

対象者ID	身長(cm)	体重(kg)	性別	年齢
1	170	70	M	54
2	162	50	F	34
3	166	72	M	62
4	170	75	M	41
5	164	55	F	37
6	159	62	F	55
7	168	80	F	67
8	183	78	M	47
9	157	47	F	49
10	185	100	M	45

対象者ID	身長(cm)	体重(kg)	性別	年齢
1	170	70	M	54
2	162	50	F	34
3	166	72	M	62
4	170	75	M	41
5	164	55	F	37
6	159	62	F	55
7	168	80	F	67
8	183	78	M	47
9	157	47	F	49
10	185	100	M	45