質的変数

myd（読みんこんだデータセット） には、女性を表す famale 変数と、内閣に対する支持を表す support 変数がある。これらの変数では、数字自体には特に意味がない。このような変数は、質的変数 (qualitative variables) と呼ばれる。また、female のように、ある特性を持っているかどうかを0 と 1 で表現する変数をダミー変数 (dummy variable) と呼ぶ。female の場合、変数名が female（女性）なので、女性という特性をもっているかどうかを表している。その特性を持っていれば（つまり、女性は）1、持っていなければ（つまり男性は）0で表される。

データセットのうち、何人が女性で何人が男性化を調べるために、表を作ってみよう。表は table() で作れる。

table(myd$female)

## 
##   0   1 
## 500 500

男性 (0) が500人、女性 (1) が500人いることがわかる。

このままでも分析はできるが、0と1だと何を表しているのかわかりにくい。そこで、0と1にラベル付けをしよう。質的変数は、factor という型に変形して、ラベルをつけると使いやすい。そのために factor() を使う。

myd$female <- factor(myd$female, levels = c(0, 1), labels = c("男性", "女性"))

もう一度、先ほどと同じように表を作ってみよう。

table(myd$female)

## 
## 男性 女性 
##  500  500

男性と女性というラベルが付いたことがわかる。

女性の人数を直接数えるには、次のようにすればよい。

sum(myd$female == "女性")

## [1] 500

一度ラベルをつけてしまうと、元の値は使えなくなるので注意が必要である。たとえば、女性ダミーで値が1の個数を数えると、

sum(myd$female == 1)

## [1] 0

となってしまう。

もう一つの質的変数である内閣への支持も表にしてみよう。

table(myd$support)

## 
##   0   1 
## 450 550

450人が内閣を不支持、550人が支持していることがわかる。やはり、0と1ではわかりにくいので、ラベルをつけよう。

myd$support <- factor(myd$support, levels = c(0, 1), labels = c("不支持", "支持"))

もう一度、先ほどと同じように表を作ってみよう。

table(myd$support)

## 
## 不支持   支持 
##    450    550

ラベルをつけることができた。

クロス表

上で見た2つの質的変数同士には、何か関係があるのだろうか。これを確かめるために、2つの質的変数の関係を表にしてみよう。複数の変数を使った表をクロス表 (corss table, contingency table)と呼ぶ。クロス表も table() で作れる。

table(myd$female, myd$support) # 行の変数を先に、列の変数を後に書く

##       
##        不支持 支持
##   男性    200  300
##   女性    250  250

myd$ をすべての変数につけるのが面倒なときは、with() を使って次のように書いてもいい。

with(myd, table(female, support))

##       support
## female 不支持 支持
##   男性    200  300
##   女性    250  250

変数名も表示されるので、こちらの方が便利かもしれない。

この表に名前をつけて後で使えるようにしよう。

tbl_fem_sup <- with(myd, table(female, support))

男性で内閣不支持が200人、男性で支持するのが300人、女性では不支持も支持も250人ずついることがわかる。

この表には、行 (row) や列 (column) ごとの合計値が示されていない。合計値は、あったほうが便利だが、表の中身がわかっていれば計算することができるので、必要不可欠な情報ではない。つまり、合計値は、データにとっては周辺的な情報であるし、物理的に見ても表の周辺に記載されることになるので、周辺度数 (margins) と呼ばれる。表を読む際には周辺度数があったほうがわかりやすいので、addmargins() で周辺度数を加えよう。

addmargins(tbl_fem_sup)

##       support
## female 不支持 支持  Sum
##   男性    200  300  500
##   女性    250  250  500
##   Sum     450  550 1000

表の周辺に合計 (Sum) が加えられた。

列の周辺度数のみを加えたいときは、margin = 1 を指定する。

addmargins(tbl_fem_sup, margin = 1)

##       support
## female 不支持 支持
##   男性    200  300
##   女性    250  250
##   Sum     450  550

行の周辺度数のみを加えたいときは、margin = 2 を指定する。

addmargins(tbl_fem_sup, margin = 2)

##       support
## female 不支持 支持 Sum
##   男性    200  300 500
##   女性    250  250 500

この表から、男性の方が内閣支持の割合が多いことがわかる。

割合 (proportion) そのものを表示したいときは、作ったテーブルに対して prop.table () を使う。 prop.table() を使うときに注意すべき点が2つある。 まず、何を100パーセントにするのか決める必要がある。 私たちのデータの場合、(1) 男性は男性だけで100パーセント、女性は女性だけで100パーセントにする、(2) 内閣支持者と内閣不支持者をそれぞれ100パーセントにする、(3) 全体（1,000人）を100パーセントにするという3通りが考えられる。私たちの表では、性別が行に、内閣の支持が列に表示されているので、(1) の場合を行パーセント、(2) を列パーセント、 (3) を全体パーセント と呼ぶ。ここでは、「性別によって内閣支持に違いがあるか」を調べたいとしよう。行にある変数である性別ごとの違いを知りたいので、行パーセントを指定する。行パーセント計算するときは、margin = 1 を指定する（列パーセントは margin = 2、全体パーセントは margin を指定しない）。addmargin() と margin 指定の仕方が異なるので、注意が必要である。

第2に、周辺度数を加える場合は、行と列のうち、原因として注目していないほうの周辺度数は、prop.table() を使う前に、注目している方の周辺度数は prop.table() より後に加えた方がよい。私たちは性別（行変数）による違いに注目しているので、

1. 列（内閣支持）の周辺度数を加える
2. 割合を計算する
3. 行（性別）の周辺度数を加える

という順番で表を作った方がよい。そのためには、次のようにすればよい。

tbl_fem_sup_m1 <- addmargins(tbl_fem_sup, margin = 1) # 列周辺度数を加える
tbl_fem_sup_p <- prop.table(tbl_fem_sup_m1, margin = 1) # 行パーセントで割合に変換する 
tbl_fem_sup_2 <- addmargins(tbl_fem_sup_p, margin = 2)  # 行周辺度数を加える

パーセント表示にしたいときは、100をかければよい。

(tbl_fem_sup_3 <- tbl_fem_sup_2 * 100)

##       support
## female 不支持 支持 Sum
##   男性     40   60 100
##   女性     50   50 100
##   Sum      45   55 100

男性の内閣支持率は60%、女性の内閣支持率は50%であることがわかる。

独立性の検定

私たちの標本（データ）では、女性の内閣支持理より、男性の内閣支持率の方が高い。これはたまたま得られた結果だろうか。それとも、母集団でも男性の内閣支持率の方が高いと考えられるだろうか。

これを確かめるために、統計的検定を行う。ここで検証する仮説は、以下のものである。

• 帰無仮説：性別と内閣支持率は独立である（つまり、男性の内閣支持率 = 女性の内閣支持率）
• 対立仮説：性別と内閣支持率には関連がある [独立ではない] （つまり、男性の内閣支持率 \(\neq\) 女性の内閣支持率）

この検定は、2つの変数か独立である（帰無仮説）か、独立でない（対立仮説）かを確かめるので、独立性の検定 と呼ばれる。また、検定に \(\chi^2\) （カイ二乗）分布を使うので、\(\chi^2\) 検定 と呼ばれることもある（英語では、\(\chi^2\) [chi-square] test of independence と呼ばれる）。

検定で使うカイ二乗分布の自由度は、分析する表の \((行数 - 1) \times (列数 - 1)\) である。周辺度数を除くと、私たちの表は2行（男性または女性）\(\times\) 2列（不支持または支持）なので、自由度 \((2-1)(2-1) = 1\) のカイ二乗分布を利用する。

有意水準を7% に設定すると、検定に使う臨界値は、qchisq() を使って、

qchisq(p = 0.07, df = 1, lower.tail = FALSE)

## [1] 3.28302

ということがわかる。検定統計量がこの臨界値より大きいとき、私たちは帰無仮説を棄却する。 反対に、検定統計量がこの臨界値以下のとき、私たちは帰無仮説を棄却しない。

検定統計量は、chisq.test() で計算できる（\(\ast\) イェーツの連続性補正は行わないので correct = FALSE とする）。

chisq.test(myd$female, myd$support, correct = FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  myd$female and myd$support
## X-squared = 10.101, df = 1, p-value = 0.001482

X-squared が検定統計量である。ここでは、10.101 という値が得られた。この値は、有意水準7%での臨界値である3.28302 より大きいので、帰無仮説は棄却される（\(p\)値の意味がわかるなら、\(p\)値を使って判断してもよい)。

したがって、内閣支持率は性別よって異なり、男性の方が女性よりも内閣を支持するという判断を下す。

量的変数

私たちのデータセットには、height（身長）と faheigh（父親の身長） という2つの量的変数（quantitative variables）がある。量的変数とは、簡単にいうと、数値自体に意味がある変数である。

量的変数を調べるときは、ます、基本的な統計量とヒストグラムを確認する。 まず、height について確認する。

mean(myd$height)

## [1] 165.8727

median(myd$height)

## [1] 166.15

var(myd$height)

## [1] 53.41806

sd(myd$height)

## [1] 7.308766

hist(myd$height, col = "skyblue", xlab = "身長 (cm)", ylab = "度数", main = "")

ヒストグラムに山が2つある。なぜだろうか？ このデータが身長のデータであり、男女共データに含まれていることを考えると、性別によって分布の山が変わりそうである。 男女別にヒストグラムを作ってみよう。

par(mfrow = c(2, 1))  # 図を2行1列に配置する
hist(myd$height[myd$female == "女性"], col = "dodgerblue", xlim = c(150, 180),
     main = "女性", xlab = "身長 (cm)", ylab = "度数")
hist(myd$height[myd$female == "男性"], col = "forestgreen", xlim = c(150, 180),
     main = "男性", xlab = "身長 (cm)", ylab = "度数")

やはり、性別によって分布が異なることがわかる。このようなことをわかるので、データ分析の前に必ずヒストグラムを確認する習慣をつけてほしい。

同様に、faheight に着いて確認する。

mean(myd$faheight)

## [1] 169.7533

median(myd$faheight)

## [1] 169.6

var(myd$faheight)

## [1] 31.34405

sd(myd$faheight)

## [1] 5.598576

par(mfrow = c(1, 1))  # 図を1つだけ表示する状態に戻す
hist(myd$faheight, col = "gray", xlab = "父親の身長 (cm)", ylab = "度数", main = "")

父親の身長は、（当たり前だが）子の性別によって変わることはないので、分布の山は1つである。

質的変数と量的変数の関係

身長のヒストグラムから明らかになったように、性別（質的変数）と身長（量的変数）の間に関係がありそうである。 質的変数と量的変数の関係は、図示すると理解しやすい。箱ひげ図 (box-and-whisker plot または単に box plot) と呼ばれる図を作って確認してみよう。箱ひげ図は、boxplot(量的変数 ~ 質的変数, data = データ名) で作る。

boxplot(height ~ female, data = myd, col = "gray",
        xlab = "性別", ylab = "身長 (cm)")

この図は、箱がデータの中心50%分を表している。箱の中にある太い線は、中央値（50パーセンタイル、第2四分位）である。箱の下端が第1四分位（25パーセンタイル）、箱の上端が第3四分位（75パーセンタイル）を示している。また、箱から上下に伸びている線がひげである。ひげの下端が最小値（外れ値を除く）であり、上端が最大値（外れ値を除く）である。ひげよりも外側にある点は外れ値（outlier）である。

\(\ast\)（第3四分位（箱の上端）より、1.5\(\times\)IQR以上大きい値、または第1四分位（箱の下端）より、1.5\(\times\)IQR 以上小さい値が外れ値と見なされる。ただし、IQR（interquartile range, 四分位範囲）とは、第3四分位と第1四分位の差（図で考えると、箱の高さ）である。）

この図を見ると、男性の身長の分布と女性の身長の分布は異なる分布であり、男性の身長のほうが高いようである。

しかし、これはあくまで標本の分布であり、母集団でも同じことが言えるかどうかは検定で確かめる必要がある。 これを確かめるためには、前回学習したウェルチ (Welch) の \(t\) 検定 を行えばよい（各自で実行してみること）。

量的変数同士の関係

次に量的変数同士（height とfaheight）の関係を確かめてみよう。

まず、2つの関係を図示する。2つの量的変数同士の関係は、散布図 (scatter plot) を使って図示する。plot() で散布図が描ける。

plot(height ~ faheight, data = myd, col = "slategray4",
     pch = 16,  # pch で点の形を指定する。
     xlab = "父親の身長 (cm)", ylab = "身長 (cm)", main = "")

## plot(myd$faheight, myd$height, col = "blue", pch = 16,
##     xlab = "父親の身長 (cm)", ylab = "身長 (cm)", main = "")  # これでも同じ

（pch については、ココ を参照。）

この図を見ると、父親の身長が高いほど、子の身長が高いという関係がありそうだ。ヒストグラムで確認したように、データには男女が含まれているが、性別によって身長の分布が異なり、点が2つのグループに分かれていることがわかる。

次に、2変数の直線的な関係の強さを測るために、相関係数 (correlation coefficient) を計算してみよう。cor() を使う。

with(myd, cor(height, faheight))

## [1] 0.2394457

相関係数は約0.24であり、弱い正の相関があるようである。

図から得られた情報と、相関係数から得られた情報は、整合的だろうか。 相関係数だけ見ると、父親の身長と子の身長にはあまり強い関係はなさそうだという結論になる。 しかし、散布図を見ると、2変数の間には非常に強い関係がありそうである。

図からわかる通り、散布図の点は2つのグループに分かれている。男女を色分けして散布図を作り直してみよう。

plot(height ~ faheight, data = myd, pch = 16,
     col = ifelse(female == "女性", "dodgerblue", "forestgreen"), # 性別で色分け
     xlab = "父親の身長 (cm)", ylab = "身長 (cm)", main = "")
## text(x = c(180, 155), y = c(155, 175), labels = c("女", "男"),
##      col = c("dodgerblue", "forestgreen"))                     # 以下の2行を一度に実行する
text(x = 180, y = 155, labels = "女", col = "dodgerblue")   # 指定した位置に「女」と書く
text(x = 155, y = 175, labels = "男", col = "forestgreen")  # 指定した位置に「男」と書く

\(\ast\)本当は、色の区別が難しい人に配慮して、色だけでなく、形も変えたほうがよい。plot() を使ってそれを実現するのはやや面倒なので、ggplot2 というパッケージを使う。

## install.packages("ggplot2") # インストール済みでない場合はまずインストールする
library("ggplot2")  # パッケージを読み込む
## theme_set(theme_gray(base_size = 12, base_family = "HiraKakuProN-W3")) # Macユーザはこの行も実行する

次のようにして散布図を作る（ggplot の使い方は、『計量経済学』で解説する）。

plt <- ggplot(data = myd, aes(x = faheight, y = height,
                              color = female, shape = female)) +
  geom_point() +
  scale_color_discrete(name = "性別") +
  scale_shape_discrete(name = "性別") +
  labs(x = "父親の身長 (cm)", y = "身長 (cm)") 
print(plt)

このように、性別ごとにグループができている。

男女別にして相関係数を求め直そう。 まず、女性の場合について求める。

with(myd[myd$female == "女性", ], cor(height, faheight))

## [1] 0.8147494

強い正の相関がある。

男性はどうだろうか。

with(myd[myd$female == "男性", ], cor(height, faheight))

## [1] 0.8121794

やはり強い正の相関がある。

この例からわかる通り、単に相関係数だけを求めると、2変数の関係を見誤る可能性がある。反対に、散布図だけに頼ると、ありもしないパタンを「見つけて」しまうことがある（偶然できた壁のシミが顔に見えることがあるように、人間は意味のないパタンを認識し易い）。 したがって、2つ（以上）の量的変数の関係について調べるときは、散布図と相関係数をセットで使うことが重要である。