Rで基礎分析(クロス集計表、相関係数、分散分析)

このエントリーは統計学勉強会用の下書きである。

ここでは、Rでクロス集計表、相関係数、分散分析をする方法を述べる。 どんな難しい分析をするときにも、最初は記述統計、そして、この3つの分析をして、多変量解析に進む。 その意味で、基礎分析と言っても、初学者にとってのもの、という意味ではなく、すべての分析の基礎にあるという意味で捉えた方が適当である。

また、初学者にとっても、この3つの分析は重要である。 統計解析を普段しない人と話すときに、「統計学が苦手」「やろうとは思っているが難しい」といったことがでてくるのだが、びっくりするほど、難しいものではないように思う。難しいものは難しいし、簡単なものは簡単だ、というくらいであって、他のことと大差はない。基礎的な3つの分析であれば、たれでも十分にできるはずである。

おそらく、統計学の授業で最初に手計算をする先生にあたるのが、ほとんどの場合、統計学嫌いになる原因である。 うまく、授業をすれば、嫌いになるほど大変なものではないはずだ。

分析の決め方

変数は尺度水準によっていくつかに分かれる。 細分化された尺度水準を考慮する分析は存在するが、その分析をする頃には尺度水準についても詳しくなっているはずである。 とりあえず、分析をしたいという人であれば、2つの水準に分けるのがよいし、多くの場合2つの水準で十分である。

  1. 連続変数 2つの値の間の無限の数を持つ数値型の変数
    例:身長、体重、車の走行距離、年齢、年収など。

  2. カテゴリカル変数 有限な数のカテゴリーによって構成される変数 例:生物学的性(男女)、人種、職種、大卒など。

2つの変数をみて、連続なのかカテゴリカルなのかを判断し、下の表で分析を選択する。

連続変数 カテゴリカル変数
連続変数 連続変数 分散分析
カテゴリカル変数 分散分析 クロス集計表

今回はAERパッケージのCPS1985というデータを例にして説明したい。

library(AER) #パッケージの呼び出し。CPS1985というデータがこのパッケージに含まれている
head(CPS1985) #CPS1985データの行頭のみ表示
?CPS1985 #データについて説明の表示
data(CPS1985) #データの呼び出し

クロス集計表

データの説明を読むとどのような変数があるかがわかる。クロス集計表はカテゴリカル変数である、gender(性別, 1=Female, 0=Male)とMARR(婚姻状態, 0=Unmarried, 1=Married)の2つの変数の関連について調べてみよう。

スクリプト

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の頭2文字である。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のクロス集計表であれば、カイ二乗検定の結果と調整済み標準化残差の結果はほとんどの場合同一であるが、一つ変数に3カテゴリ以上ある場合には、調整済み標準化残差だけが有意を示す場合がある。このようなケースでは、リコードをどこで行えばよいかがわかるため、クロス集計表の分析をするときには、調整済み標準化残差の値は出しておいた方がよいだろう。

オッズ比

オッズ比の計算はRの標準のスクリプトにはない(ただし、フィッシャーテストをすれば出力できる)。 epitoolsパッケージを使うのが最も簡単であろう。

スクリプト

library(epitools) #epitoolsパッケージの読み込み
oddsratio(ct1)$measure

結果

odds ratio with 95% C.I.
       estimate     lower    upper
male    1.00000        NA       NA
female  1.04829 0.7324027 1.502553

参照:Rでオッズ比と調整済み残差を出す

クロス集計表が見にくいとき

データが数値の場合(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の頭3文字である。 オプションを何も指定しないとピアソンの積率相関係数が求められる。 オプションは以下のようにつける。

 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と同じものである。 下2段は平均値の比較であり、白人を意味するcauc(Caucasian)が最も賃金が高いという結果になっている。