このエントリーは統計学勉強会用の下書きである。
ここでは、Rでクロス集計表、相関係数、分散分析をする方法を述べる。 どんな難しい分析をするときにも、最初は記述統計、そして、この３つの分析をして、多変量解析に進む。 その意味で、基礎分析と言っても、初学者にとってのもの、という意味ではなく、すべての分析の基礎にあるという意味で捉えた方が適当である。
また、初学者にとっても、この３つの分析は重要である。 統計解析を普段しない人と話すときに、「統計学が苦手」「やろうとは思っているが難しい」といったことがでてくるのだが、びっくりするほど、難しいものではないように思う。難しいものは難しいし、簡単なものは簡単だ、というくらいであって、他のことと大差はない。基礎的な３つの分析であれば、たれでも十分にできるはずである。
おそらく、統計学の授業で最初に手計算をする先生にあたるのが、ほとんどの場合、統計学嫌いになる原因である。 うまく、授業をすれば、嫌いになるほど大変なものではないはずだ。
分析の決め方
変数は尺度水準によっていくつかに分かれる。 細分化された尺度水準を考慮する分析は存在するが、その分析をする頃には尺度水準についても詳しくなっているはずである。 とりあえず、分析をしたいという人であれば、２つの水準に分けるのがよいし、多くの場合２つの水準で十分である。
連続変数 2つの値の間の無限の数を持つ数値型の変数
例：身長、体重、車の走行距離、年齢、年収など。カテゴリカル変数 有限な数のカテゴリーによって構成される変数 例：生物学的性(男女)、人種、職種、大卒など。
２つの変数をみて、連続なのかカテゴリカルなのかを判断し、下の表で分析を選択する。
| 連続変数 | カテゴリカル変数 | |
|---|---|---|
| 連続変数 | 連続変数 | 分散分析 |
| カテゴリカル変数 | 分散分析 | クロス集計表 |
今回はAERパッケージのCPS1985というデータを例にして説明したい。
library(AER) #パッケージの呼び出し。CPS1985というデータがこのパッケージに含まれている head(CPS1985) #CPS1985データの行頭のみ表示 ?CPS1985 #データについて説明の表示 data(CPS1985) #データの呼び出し
クロス集計表
データの説明を読むとどのような変数があるかがわかる。クロス集計表はカテゴリカル変数である、gender(性別, 1=Female, 0=Male)とMARR(婚姻状態, 0=Unmarried, 1=Married)の２つの変数の関連について調べてみよう。
ct1<-table(CPS1985$gender,CPS1985$married) ct1 res1<- chisq.test(ct1) res1
変数名は「データ名$変数名」とドルマークを挟んで指定する。
結果
no yes male 101 188 female 83 162 X-squared = 0.028233, df = 1, p-value = 0.8666
有意差があるか否かは、カイ二乗検定の結果で判断する。 結果はp-valueと表示されている値である。
フィッシャーの正確確率検定
セルのケース数が一桁になる場合は、フィッシャーの正確確率検定を使った方がよい場合がある。Rではカイ二乗検定をしたときにフィッシャーテストをした方がよいというメッセージが出るので、そのメッセージを参考にするとよいだろう。
fisher.test(ct1)
chisq.testの頭の部分をfisherに変えればよい。 ちなみにchiはカイs、qは二乗のsquareの頭２文字である。testは検定のことである。
結果
p-value = 0.8551 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.7213704 1.5258105 sample estimates: odds ratio 1.048508
フィッシャーテストの方は95%信頼区間、オッズ比まで出してくれるので便利である。
調整済み標準化残差
調整済み標準化残差を調べたいときにはカイ二乗検定のあとに下記のスクリプトを追加する。
res1$stdres
結果
no yes male 0.259397 -0.259397 female -0.259397 0.259397
調整済み標準化残差は絶対値で、1％水準は2.58、0.1％水準は3.29、0.01%水準は3.89以上になる。 2×2のクロス集計表であれば、カイ二乗検定の結果と調整済み標準化残差の結果はほとんどの場合同一であるが、一つ変数に３カテゴリ以上ある場合には、調整済み標準化残差だけが有意を示す場合がある。このようなケースでは、リコードをどこで行えばよいかがわかるため、クロス集計表の分析をするときには、調整済み標準化残差の値は出しておいた方がよいだろう。
オッズ比
オッズ比の計算はRの標準のスクリプトにはない(ただし、フィッシャーテストをすれば出力できる)。 epitoolsパッケージを使うのが最も簡単であろう。
library(epitools) #epitoolsパッケージの読み込み oddsratio.wald(ct1)$measure
結果
odds ratio with 95% C.I.
estimate lower upper
male 1.000000 NA NA
female 1.048577 0.7327264 1.50058
クロス集計表が見にくいとき
データが数値の場合(CPS1985はカテゴリカル変数は文字列で入っているので問題ない)はexpressパッケージを使うとわかりやすくなるかもしない。
library("expss") #expssパッケージの呼び出し
cro(CPS1985$gender,CPS1985$married)
croの中に入れるのはtableの中に入れるものと同じ。 結果
| | | CPS1985$married | | | | | no | yes | | -------------- | ------------ | --------------- | --- | | CPS1985$gender | male | 101 | 188 | | | female | 83 | 162 | | | #Total cases | 184 | 350 |
この表は見やすくなっている以上に、Markdown形式であることがおそらく重要である。
Markdown形式では次のように表示される。
| CPS1985$married | |||
|---|---|---|---|
| no | yes | ||
| CPS1985$gender | male | 101 | 188 |
| female | 83 | 162 | |
| #Total cases | 184 | 350 |
例えば、Atomエディタに貼り付けて、Markdown preview / Markdown preview Enhancedで表示させて、プレビューされた内容をコピーしてExcelに張り付けるといった作業も可能である。 memiscパッケージを使うのが王道だが、手数が少ない分、このルートも使いどころはあるかもしれない。
相関係数
相関係数は連続変数同士の関連を調べる分析である。 age年齢とwage賃金について検討してみよう。
cor(CPS1985$age,CPS1985$wage)
corは相関のcorrelationの頭３文字である。 オプションを何も指定しないとピアソンの積率相関係数が求められる。 オプションは以下のようにつける。
cor(CPS1985$age,CPS1985$wage, , method="pearson") # ピアソンの積率相関係数 cor(CPS1985$age,CPS1985$wage, , method="kendall") # ケンドールの順位相関係数 cor(CPS1985$age,CPS1985$wage, , method="spearman") # スピアマンの順位相関係数
無相関か否か(関連が統計学的に有意か否か)を求めるには、testと打てばよい。
cor.test(CPS1985$age,CPS1985$wage)
結果
Pearson's product-moment correlation data: CPS1985$age and CPS1985$wage t = 4.1472, df = 532, p-value = 3.917e-05 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.09352049 0.25794433 sample estimates: cor 0.1769669
こちらもp-valueをみて有意か否かを判断する。
分散分析
カテゴリカル変数と連続変数の分析をする場合、分散分析を行うことになる。ここでは、エスニシティethnicityと賃金wageの関係について検討する。 エスニシティはカテゴリカル変数であり、賃金は連続変数である。
コマンド
res<- aov(wage ~ ethnicity, data =CPS1985) # 分散分析 summary(res) # 結果表示 with(CPS1985, tapply(wage, ethnicity, mean, na.rm=TRUE)) # 平均値の計算
aovは分散分析の英語であるAnalysis of Varianceの頭文字をとったコマンドである。ANOVAと略すことが多いがRではaovである。
連続変数 ~ カテゴリカル変数, データの順番で記述する。
aovコマンドでは平均値の比較ができない。 そこで、withから始まる計算式を使う。 meanは平均値、na.rm=TRUEは欠損値を含むケースを除外するコマンドである。
結果
Df Sum Sq Mean Sq F value Pr(>F) ethnicity 2 173 86.33 3.297 0.0378 * Residuals 531 13904 26.18 cauc hispanic other 9.277932 7.283333 8.058358
Pr(>F) がp-valueと同じものである。 下２段は平均値の比較であり、白人を意味するcauc(Caucasian)が最も賃金が高いという結果になっている。
追記：psychパッケージを利用する
library(psych) describe.by(CPS1985$wage, list(CPS1985$ethnicity))
