群馬大学 | 医学部 | サイトトップ | 医学情報処理演習

医学情報処理演習：2011年度第7回の関数と補足テキスト

補足テキスト

テキストには【対応のある多群間の比較】を掲載していないので，ここで説明する。

【例1】（出典：丹後『医学への統計学』朝倉書店，例題8.6，p.163-4）
催眠に関する研究において，恐怖，幸福，落胆，沈着の感情が催眠状態にある8人の被検者に与えられた。その結果，得られた皮膚電位(mV)の測定値が表74に示されている。感情の違いにより皮膚電位に違いが認められるといえるか。
被検者恐怖幸福落胆沈着
123.122.722.522.6
257.653.253.753.1
310.5 9.710.8 8.3
423.619.621.121.6
511.913.813.713.3
654.647.139.237.0
721.013.613.714.8
820.323.616.314.8
個人差が非常に大きいので，感情の違いを要因とする一元配置分散分析では差が見えにくい。個人差を考慮する必要がある。分布も正規分布とは言い難いので，こういう場合はフリードマンの検定が適している。しかし二元配置分散分析もできる。
コードは以下の通り。
x <- matrix(c(
 1,23.1,22.7,22.5,22.6,
 2,57.6,53.2,53.7,53.1,
 3,10.5, 9.7,10.8, 8.3,
 4,23.6,19.6,21.1,21.6,
 5,11.9,13.8,13.7,13.3,
 6,54.6,47.1,39.2,37.0,
 7,21.0,13.6,13.7,14.8,
 8,20.3,23.6,16.3,14.8),8,5,byrow=TRUE) # データを行列として入力
colnames(x) <- c("subj","fear","happiness","despair","calmness") # 刺激の種類を定義
matplot(t(x[,-1]), type="l", ylab="ev", xlab="stimuli", xaxt="n") # ４種類の刺激に対する個人ごとの電位をプロット
axis(1, 1:4, colnames(x[,-1])) # グラフの横軸に刺激の種類を表示
mx <- colMeans(x[,-1]) # 各刺激についての電位の平均値。apply(x[,-1],2,mean)と同値。
sx <- apply(x[,-1],2,sd) # 各刺激に対する電位の標準偏差
plot(1:4,mx,type="l",ylim=c(min(mx-sx),max(mx+sx)),ylab="ev",xlab="stimuli",xaxt="n") # ４種類の刺激に対する平均電位をプロット
arrows(1:4,mx-sx,1:4,mx+sx,code=3,angle=90) # 標準偏差のひげをプロット
axis(1, 1:4, colnames(x[,-1])) # グラフの横軸に刺激の種類を表示
dat <- data.frame(
 subj=rep(x[,1],4),
 emst=gl(4,8,labels=c("fear","happiness","despair","calmness")),
 ev=c(x[,2],x[,3],x[,4],x[,5])) # 縦長データフレームの定義
friedman.test(ev ~ emst | subj, data=dat) # フリードマンの検定。このデータでは有意水準5％で有意ではない。
library(car) # Anova()関数を使うためcarライブラリをロード
stimuli <- factor(colnames(x[,-1]), ordered=FALSE) # 刺激種類を要因型変数として定義
rep.anova <- Anova(lm(x[,-1]~1), idata=data.frame(stimuli), idesign= ~stimuli, type="III") # 反復測定分散分析実行
summary(rep.anova, multivariate=FALSE) # 結果の表示。G-G補正やH-F補正をしたp値は0.05より大きい。
# avova(lm(ev ~ emst+subj, data=dat)) として二元配置分散分析するより切れ味がよい。
# oneway.test(ev ~ emst, data=dat) として個人差を無視してしまうと，刺激による差はマスクされる。
４種の刺激に対する個人別反応グラフ
４種の刺激に対する平均反応グラフ
【例2】（出典：高橋・大橋・芳賀『SASによる実験データの解析』東京大学出版会，第12章）
10mgの降圧剤経口投与後の血圧変化データ（出典には詳しい情報がないが，mmHg単位の収縮期血圧であろう）。被験者6名（おそらくhealthy volunteers）を投与1時間前から8時間後まで1時間ごとに測定した値が下表のようであった。
被験者1時間前投与直前1h2h3h4h5h6h7h8h
1112119113105114110115114110111
2116110115110112107116115120118
3122123126114111113119123119124
4124130127110100127130134120124
5126121115122124117124132128120
6129135125122115110114124133131
＊経時的な観察では，一定の反応パタンがあるかどうかを知りたい（上昇傾向，下降傾向，U字，逆U字など）ことが多いので，既知の関数（１次関数，２次関数など）を当てはめることもよく行われる。
＊パタンがよくわからないとき，経時的観察データをそのまま分析する。個人差を考慮して，何らかの経時的変化があるかどうかを調べることになる（帰無仮説は一定の経時的変化の傾向はない，ということ）。
データ定義をして経時的変化パタンをプロットするコードは以下の通り。
z <- matrix(c(
 1,112,119,113,105,114,110,115,114,110,111,
 2,116,110,115,110,112,107,116,115,120,118,
 3,122,123,126,114,111,113,119,123,119,124,
 4,124,130,127,110,100,127,130,134,120,124,
 5,126,121,115,122,124,117,124,132,128,120,
 6,129,135,125,122,115,110,114,124,133,131),6,11,byrow=TRUE) # データを行列として入力
colnames(z) <- c("subj","T.1","T0","T1","T2","T3","T4","T5","T6","T7","T8") # 変数名定義
matplot(t(z[,-1]), type="l",ylab="BP(mmHg)",xlab="time",xaxt="n") # 個人ごとの経時変化パタンをプロット
axis(1,1:10,labels=colnames(z[,-1])) # 横軸ラベル表示
mx <- colMeans(z[,-1]) # 各時点での血圧平均値。apply(z[,-1],2,mean)と同値。
sx <- apply(z[,-1],2,sd) # 各時点での血圧の標準偏差
plot(1:10,mx,type="l",ylim=c(min(mx-sx),max(mx+sx)),ylab="ev",xlab="time",xaxt="n") # 各時点での平均血圧をプロット
arrows(1:10,mx-sx,1:10,mx+sx,code=3,angle=90) # 標準偏差のひげをプロット
axis(1,1:10,labels=colnames(z[,-1])) # 横軸ラベル表示
dat2 <- data.frame(
 subj=factor(rep(z[,1],10)),
 tvars=gl(10,6,labels=colnames(z[,-1])),
 BP=as.vector(z[,-1])) # フリードマンの検定や分散分析のための縦長データフレームの定義
friedman.test(BP ~ tvars | subj, data=dat2) # フリードマンの検定。個人差を考慮すると時点の有意な効果がある。
# 本来は下記だが，時点数より人数が少ないと行列が計算できずエラーになる
# library(car) # Anova()関数を使うためcarライブラリをロード
# times <- factor(colnames(z[,-1])) # 時点変数を要因型変数として定義
# rep.anova.2 <- Anova(lm(z[,-1]~1), idata=data.frame(times), idesign=~times, type="II") # 分析実行（→エラー）
# summary(rep.anova.2, multivariate=FALSE) # 分散分析の結果を表示
# そこで，下記を実行
res1 <- aov(BP ~ subj + tvars, data=dat2) # 通常の二元配置分散分析の実行結果を得る
summary(res1) # 結果を表示
mlmfit1 <- lm(z[,-1]~1) # anova.mlmによる反復測定分散分析
mlmfit0 <- lm(z[,-1]~0) # まず，このように２つのオブジェクトを計算する
anova.mlm(mlmfit1, mlmfit0, X=~1, test="Spherical") # これでSASのPROC GLMのREPEATEDと同じ結果が得られる
個人別血圧の経時変化グラフ
下がって上がる傾向がありそうにも思えるが定かではない。個人差はかなりある。
平均血圧の経時変化グラフ
※G-G補正ではp>0.05だが，H-F補正ではp<0.05。微妙なところだが，時間経過とともに５％水準で統計的に有意な血圧変化があったといえる。

今回の関数まとめ

今回使った関数や文の主なものをまとめます (A selected summary of functions and statements used in the 7th practice is shown here.)

関数名(name)機能(effect)使い方(usage)
str()オブジェクトの構造をみる(check the structure of a given object.)
gl()グループ変数を定義するgl(g, p, labels=1:g)で，p個体ずつg個のグループがつながった形の要因型変数が定義される。例えば，gl(3, 4, labels=c("X","Y","Z"))は，as.factor(c("X","X","X","X","Y","Y","Y","Y","Z","Z","Z","Z"))と同値である。
summary()オブジェクトの詳細出力を得る
aov()分散分析の実行(1)単純にaov(values ~ groups)という形では分散分析表が出力されないので，必ずsummary(aov(values ~ groups))という形にする。データフレーム内の変数を指定する場合は，オプションとしてdata=を用いる。【対応のない多群間の平均値の差の検定に相当】
anova()分散分析の実行(2)引数は線形モデルの結果オブジェクト。anova(lm(values ~ groups))という形で一元配置分散分析が実行できる。【対応のない多群間の平均値の差の検定に相当】
参考までに，二元配置の場合はanova(lm(values ~ gvar1+gvar2))のようにするが，交互作用効果がある場合はanova(lm(values ~ gvar1*gvar2))とする。gvar1*gvar2と書くと，gvar1+gvar2+gvar1:gvar2という意味になる。なお，反復測定データの一元配置分散分析（次々項参照）は，要因の一つを個人の効果とし，もう一つを時間の効果とした，繰り返しのない二元配置分散分析と同じ形になる。【対応のある多群間の平均値の差の検定に相当】
データフレーム内の変数を指定する場合は，lm()の中のオプションとしてdata=を用いる。
lm()線形モデルを当てはめるlm(values ~ groups)でvaluesをgroupsで説明するモデルの当てはめが行われる。データフレーム内の変数を指定する場合は，オプションとしてdata=を用いる。
Anova()carライブラリの分散分析関数library(car)の後で実行する。通常のanova()と同じく，lm(values ~ groups)を引数にとって分散分析を実行することもできるが，idata=とidesign=オプションを使って，反復測定分散分析を実行することができる。例えば，すべての対象者に対して３時点での１種類の測定値があり，時点１での値が変数T1，時点２での値が変数T2，時点３での値が変数T3に入っていて，同じ人は同じ順番にでてくる場合，
tvars<-as.factor(c("T1","T2","T3"))
として時点変数を定義するオブジェクトtvarsを生成してから，
res <- Anova(lm(cbind(T1,T2,T3)~1), idata=data.frame(tvars), idesign=~tvars, type="III")
として反復測定分散分析を実行し，結果を，
summary(res, multivariate=FALSE)
によって表示できる。【対応のある多群間の平均値の差の検定に相当】
oneway.test()Welchの拡張による一元配置分散分析の実行summary(aov())やanova(lm())は各群の等分散性を仮定しているので，2群の場合と同じく，等分散性を仮定しないWelchの拡張の方が良い結果が得られる。oneway.test(values ~ groups)とする。データフレーム内の変数を指定する場合は，オプションとしてdata=を用いる。【対応のない多群間の平均値の差の検定に相当】
bartlett.test()バートレットの検定を実行する多群の母分散がすべて等しいという帰無仮説を検定する。bartlett.test(values ~ groups)とする。データフレーム内の変数を指定する場合は，オプションとしてdata=を用いる。
pairwise.t.test()多重比較の実行多群から2群ずつの組み合わせ全てについて標準誤差をプールしたt検定を繰り返し検定の多重性を調整した有意確率を計算する。文法は，pairwise.t.test(values,　groups,　p.adjust.method="holm")である。data=オプションは使えないので，データフレーム内の変数はx$valuesのように参照するかattach()を用いた後で実行する必要がある。p.adjust.method=オプションを省略するとHolmの方法になる（"holm"と指定するのと同等）。"bon"でBonferroniの方法になる。
attach()データフレームのアタッチattach(データフレーム名）を実行した後では，データフレーム名$変数名としなくてもデータフレーム内の変数が参照できるようになる。必要が無くなったら必ずdetach(データフレーム名）すること。
TukeyHSD()テューキーの多重比較の実行引数は分散分析の結果。つまり，TukeyHSD(aov(values~groups))とする。テューキーの方法で検定の多重性を調整した，多群からの2群ずつの組み合わせ全てについて，平均値に差が無いという帰無仮説の検定が実行できる。データフレーム内の変数を指定する場合は，aov()の中のオプションとしてdata=を用いる。
cut()連続量を区切ってカテゴリ化文法は，cut(values, intervals)で，intervalsを「〜以上〜未満」にしたいときはright=FALSEオプションをつける。例えば，10人の人のBMIが
BMI <- c(16.5, 20, 25, 40, 28.5, 22.5, 23.5, 17, 35, 28)
であるとき，18.5未満をやせ，18.5以上25未満を標準，25以上を肥満とする3つの水準からなるカテゴリ変数OBESITYに変換するには，
OBESITY <- cut(BMI, c(0,18.5,25,50), right=FALSE)
とすればOK。intervalsはデータ全体を含まないとエラーになるので，左端は0，右端は50としたが，もちろんギリギリでも構わない。さらに，18.5未満を"THIN"，18.5以上25未満を"NORMAL"，25以上を"OBESE"というカテゴリ名に変えたければ，
levels(OBESITY) <- c("THIN","NORMAL","OBESE")
とする。
kruskal.test()クラスカル＝ウォリスの検定の実行３つ以上の水準のある要因型変数が，１つの量的変数に有意な影響を与えているかどうかを調べるためのノンパラメトリックな検定であるクラスカル＝ウォリスの検定を実行する。【対応のない多群間の分布の位置の差の検定に相当】第13回参照。
friedman.test()フリードマンの検定の実行複数の対象者について，３つ以上の要因による同種の測定値（時間経過とともに繰り返されるか，別の刺激による同じ量であるなど）があるとき，対象者の個人差を考慮しても，要因の違いがその量の違いに影響しているかどうかを調べるためのノンパラメトリックな検定であるフリードマンの検定を実行する。【対応のある多群間の分布の位置の差の検定に相当】第13回参照。

リンクと引用について