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

統計処理の羅針盤

最終更新: 2008年 9月 5日 (金曜日) 12時42分 :平均値の差の検定について加筆修正

このページは,適切な統計手法を既に知っているけれども,それに対応するRの関数がわからない人のために,一覧表を提供する目的で作成した。しかし,統計手法を網羅することは事実上不可能なので,基礎的な解析に絞ることにした。どうせ基礎的な解析に絞るなら使う順番に出ている方が便利だろうと判断し,一般的な手順を同時に示すことにした(多少煩雑なのでまとめ関数も作ってみた)。ただし,グラフィカルな表示は必須であるにしても,デザインベースな解析では常に正確な確率を出すというのも一つの考え方だし,モデルベースな解析ならglm()やnlm()をきわめるとか,関数定義してoptim()をきわめるというのも一つの考え方であろうから,この手順が常にマニュアル的に適用可能とは限らない。

なお,ここでは慣例に従って(たぶん現状ではその方がニーズが大きいと思うので)仮説検定を中心にまとめたが,Kenneth J. RothmanがEpidemiology: An Introduction (Oxford University Press)の中で繰り返し強調しているように,仮説検定よりも区間推定(もちろん,それが帰無仮説に対応する値を含んでいるかどうかだけに着目するのではなくて,区間そのものを見るべきだということは,言うまでもない)の方がずっと情報量が多いし,値の意味は大きい。そもそも,統計解析結果から因果推論をするときには先行研究や生物学的なメカニズムとの整合性など,さまざまな知識を総合することが必要であって,検定結果を機械的に使うべきではない。少なくとも,ありがちな間違った解釈をしないようにしたいものである。例えば,有意確率が小さいことは,帰無仮説の下でそのような観測値が得られる可能性が小さいことを示すに過ぎず(サンプルサイズが大きくなるだけで達成されてしまう),対立仮説における差の大きさや関連の強さを意味するものでは決してないことは,肝に銘じておくべきであろう。

より包括的な手法の一覧(と使用例)は,群馬大学の青木先生のサイトにRによる統計処理としてまとまっているので,是非参考にされたい。

解析する前に

データファイルができてから

できた後でも外れ値(とくに二次元以上の)がある場合など,必要があれば,いつでも再コーディングやデータ修正,対象からの削除に戻ってやり直す。その意味では,統計処理としてやった手続きを保存しておける方がいい。Excelはマクロに覚えさせるにも限界があるので,Rのような統計ソフトを使った方がいいことになる。以下,Rでの処理を一覧表の形で示す。代表的な変数の型として,データフレーム,ベクトル,数値,文字があるが,変数名は予約語以外任意であるが,下の表内ではデータフレームをdat,量的な変数ベクトルをXまたはYで示し,カテゴリ変数ベクトルをC,二値変数(事象生起を値1で,事象が起こらないのを値0で表すものとする)ベクトルをB,順序変数ベクトルをIで示す。

Windowsでは,Rguiの起動は,デスクトップのRのアイコンをダブルクリックするだけでいい。>という記号が表示されて入力待ちになる。この記号をプロンプトと呼ぶ。Rへの対話的なコマンド入力は,基本的にプロンプトに対して行う。そうやって入力した手続きは,FileメニューのSave History ...で保存でき,後でFileのSourceで呼び出せば再現できる。プロンプトに対してsource("プログラムファイル名")としても同じことになる(但し,Windowsではファイルパス中,ディレクトリ(フォルダ)の区切りは/または\\で表すことに注意)。また,上向き矢印キーで既に入力したコマンドを呼び戻すことができる。

段階用途Rの命令
基本操作終了q()
付値(計算結果を任意の変数に保存する)<-
関数定義(例:平均と標準偏差をリストとして返す関数)meansd <- function(X) { list(mean(X),sd(X)) }
ライブラリインストール(例:CRANからvcdをダウンロード)install.packages("vcd")
読み込み1行目が変数名のタブ区切りテキストファイルのデータフレームへの読み込みdat <- read.delim("ファイル名")
1行目が変数名のコンマ区切りテキストファイルのデータフレームへの読み込みdat <- read.csv("ファイル名")
カテゴリ変数宣言(データが文字列の場合は自動的にfactorになるので不要)dat$C <- as.factor(dat$C)
順序変数宣言dat$I <- as.ordered(dat$I)
量的変数宣言dat$X <- as.single(dat$X)
確認データ構造表示str(dat)
データフレーム内変数名一覧names(dat)
一変数グラフ表示度数分布図barplot(table(C))
正規確率プロットqqnorm(X)
ヒストグラムhist(X)
箱ヒゲ図boxplot(X)
一変数基本集計度数分布表table(C)
五数要約値[最小,Q1,中央値,Q3,最大]fivenum(X)
サンプルサイズlength(X)
合計sum(X)
平均mean(X)【注:sum(X)/length(X)と同値】
不偏分散var(X)【注:sum((X-mean(X))^2)/(length(X)-1)と同値】
不偏標準偏差sd(X)【注:sqrt(var(X))と同値】
一変数の検定分布の正規性(シャピロ=ウィルクの検定)shapiro.test(X)
母平均の検定t.test(X,mu=母平均)
母比率の検定binom.test(table(B)[2],length(B),p=母比率)
二変数グラフ表示カテゴリ×カテゴリ=モザイクプロットmosaicplot(table(C1,C2))
量×カテゴリ=層別箱ヒゲ図boxplot(X~C)
量×カテゴリ=平均値とエラーバーを重ね書きしたストリップチャートstripchart(X~C); points(IX<-c(1.15,2.15),MX<-tapply(X,C,mean),pch=18); SX<-tapply(X,C,sd); arrows(IX,MX-SX,IX,MX+SX,code=3,angle=90,length=0.1)
量×量=散布図plot(Y~X)【plot(X,Y)と同値】
散布図に回帰直線を重ね書きabline(lm(Y~X))
二変数基本集計カテゴリ×カテゴリ=クロス集計表table(C1,C2)
量×量=ピアソンの相関係数(量が正規分布するとき)cor(X,Y)
量×量=スピアマンの順位相関係数(量が正規分布しないとき)cor(X,Y,method="spearman")
二変数の検定カテゴリ変数間の独立性のカイ二乗検定chisq.test(table(C1,C2))
カテゴリ変数間の独立性のFisherの直接確率fisher.test(table(C1,C2))【まとめ関数
オッズ比とその信頼区間library(vcd); summary(oddsratio(table(B1,B2),log=F))【注:対数オッズ比を計算するにはlog=Fをつけなければ良い。また,オッズ比やリスク比を計算する場合は,有意性を検定するよりも,信頼区間を計算してp値関数を想像しながら解釈する方がずっとinformativeである】
カテゴリ変数間の関連性:ファイ係数とクラメールのVlibrary(vcd); assoc.stats(table(C1,C2))
2回の繰り返しの一致度:カッパ係数library(vcd); Kappa(table(C1,C2))
順序変数×カテゴリ変数の出現頻度の傾向=Cochran-Armitageの検定prop.trend.test(table(B,I)[2,],table(I))
2つのカテゴリ間で正規分布する量の分散に差があるか:等分散性の検定(いわゆるF検定)var.test(X~B)
2つのカテゴリ間で正規分布する量の平均に差があるか:平均値の差のWelchの検定Welchの方法で自由度を調整するため,t.test(X~B)。【位置母数の比較に至るまとめ関数
2つのカテゴリ間で正規分布しないけれども分布の形は違わない量に差があるか:Wilcoxonの順位和検定wilcox.test(X~B)
対応のある2つの正規分布する量の差の検定:paired-t検定t.test(X,Y,paired=T)【内容的にt.test(X-Y,mu=0)と同等】
正規分布する量の分散がカテゴリ間で差がないか:バートレットの検定bartlett.test(X~C)
正規分布する量がカテゴリ間で差がないか:一元配置分散分析+多重比較等分散ならaov(X~C)で一元配置分散分析し,Cの主効果が有意ならoneway.test(X~C,var.equal=F)で,Welchの方法の拡張による一元配置分散分析を行い,多重比較はTukeyHSD(aov(X~C))またはpairwise.t.test(X,C)により行う(前者はTukeyの方法,後者はHolmの方法)。ただし一元配置分散分析の結果が有意でない場合の解釈には注意が必要。library(multcomp)すれば,simtest(X~C,type="Dunnett")でダネットの多重比較(対照群との比較),simtest(X~C,type="Williams")でウィリアムズの多重比較もできる。不等分散でもやってしまう場合もあるが,クラスカル・ウォリスの検定を使うこともある。
正規分布しない量×カテゴリ変数:クラスカル・ウォリスの検定+ホルムの方法で調整したWilcoxonの順位和検定を多重実行kruskal.test(X~C)とpairwise.wilcox.test(X,C)
量×量:無相関の検定(量が正規分布するとき)cor.test(X,Y)
量×量:無相関の検定(量が正規分布しないとき)cor.test(X,Y,method="spearman")
測定誤差がない量(または原因と考えられる量)Xによって,結果と考えられる量Yのばらつきを説明:回帰分析summary(lm(Y~X))

3つ以上の変数の関係をみたいときには無限に近いパタンがあるので,ここでは論じない。2つの変数の回帰関係に,カテゴリ変数による群間で差があるかどうかを見たい場合は共分散分析をする(例えばsummary(glm(Y~B+X+B:X))でB:Xの係数が有意にゼロと差があれば傾きが2群間で異なっていると判断でき,傾きに差がないとき,summary(glm(Y~B+X))でBをダミー変数化した変数の係数が有意にゼロと差があれば修正平均に差があると判断できる)。1つ以上の交絡要因の影響を調整して,2つのカテゴリ変数間の独立性をみたいときは,mantelhaen.test(C1,C2,C3)という関数(または,TMP <- table(C1,C2,C3)として3次元のクロス集計表TMPを作ってから,mantelhaen.test(TMP)としてもよい)が使える(複数階層で2つのカテゴリ変数間の関連性,例えばオッズ比の均質性を検定するためのWoolfの検定は,vcdライブラリのwoolf.test(table(B1,B2,C))で実行できる)。複数の独立変数で一つの変数のばらつきを説明したいときは,glm(Y~X1+X2+...)を使って重回帰分析やロジスティック回帰分析をしてもいい(独立変数群は線型和なら+で,交互作用だけなら:で,両方入れるときは*で結ぶ。a*bはa+b+a:bと同値)。そのときAICによって変数選択させたければstep()で囲めばいい(より細かく指定したいときはMASSライブラリのstepAIC()を用いる)。


リンクと引用について

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