群馬大学 | 医学部 | サイトトップ | 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である】 | |
カテゴリ変数間の関連性:ファイ係数とクラメールのV | library(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) | |
正規分布する量がカテゴリ間で差がないか:一元配置分散分析+多重比較 | ||
正規分布しない量×カテゴリ変数:クラスカル・ウォリスの検定+ホルムの方法で調整した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.