『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』

第2章　2変数の関連の記述統計

0 全体の章構成

• 第 1章　1変量の記述統計の基礎
• 第 2章 2変数の関連の記述統計
• 第 3章　推測統計の基礎
• 第 4章　統計的検定の一般型(general form of statistical testing)
• 第 5章　2変数の関連の推定と検定
• 第 6章　2群の母平均の差のt検定
• 第 7章　平均値の差の分散分析(ANOVA)
• 第 8章　単回帰分析(SRA)
• 第 9章　重回帰分析（Ⅰ）(MRA)
• 第10章　重回帰分析（Ⅱ）(MRA)
• 第11章　主成分分析と因子分析
• 第12章　ロジスティック回帰分析
• 第13章　階層構造のあるデータの分析

1 ピアソンの積率相関係数

グラフ

script

#### 不偏分散から素の標準偏差を求める関数を定義
sd2 <- function(x) sqrt(sd(x)^2*(length(x)-1)/length(x))

#### 標本分散を用いて標準化する関数を定義
std <- function(x) (x-mean(x))/sd2(x)

N <- 500 # 母集団サイズ
y <- rnorm(N, mean=0, sd=1)
x <- rnorm(N, mean=0, sd=1)

e <- lm(y ~ x)$residuals

R <- 0.4 # 欲しい相関係数の値

u <- R*std(x) + sqrt(1-R^2)*std(e) # 相関がR, 標本分散が1になる。

w <- std(x)*10+50
v <- u*10+50

# cor(x, u)

col_a <- "#ff9966"
par(family= "serif")
plot(w, v, col="darkgreen", pch=1, bty="n", 
     xlim=c(0, 100), ylim=c(0, 100), las=1,
     main=paste0("相関係数", sprintf("%.2f", R), "の散布図(n=", N, ")"),
     xlab="2変数とも算術平均50，標準偏差10", ylab="")
segments(mean(w), 5, mean(w), 95, lty=3, lwd=1, col="purple")
segments(5, mean(v), 95, mean(v), lty=3, lwd=1, col="purple")
text(95, 95, "第Ⅰ象限: 偏差積＋", adj=c(1, 1), col="#ff0000")
text( 5, 95, "第Ⅱ象限: 偏差積−", adj=c(0, 1), col="#0000ff")
text( 5,  5, "第Ⅲ象限: 偏差積＋", adj=c(0, 0), col="#ff0000")
text(95,  5, "第Ⅳ象限: 偏差積−", adj=c(1, 0), col="#0000ff")

1-1 散布図(scatter plot)と相関係数(correlation coefficient)

　本文のような散布図を描くスクリプトを示す。一度に全て実行するよりも，少しずつ実行して結果を確認しながら進むのが良い。散布図も二つ描画されるので，ひとつずつ実行して比べるのが良い。

x0 <- rnorm(n <- 100, mean=0, sd=1) # 100個の標準正規変量
y0 <- rnorm(n, mean=0, sd=1)

lm0 <- lm(y0 ~ x0); names(lm0) # 単回帰分析。格納された情報を確認。

 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        

r0 <- lm0$residuals # 残差のヴェクトル
x1 <- (x0-mean(x0))/sqrt(mean((x0-mean(x0))^2)) # x0を標準化
r1 <- (r0-mean(r0))/sqrt(mean((r0-mean(r0))^2)) # 残差を標準化

rxy <- 0.6 # ここで相関係数を任意に指定する。
z1 <- rxy*x1 + sqrt(1-rxy^2)*r1 # ここが肝心
x2 <- x1*10+50 # x1を「偏差値」風に変換。本質的ではない。
z2 <- z1*10+50
cor(x2, z2) # 積率相関が指定したものと一致する事を確認。

[1] 0.6

基本の散布図

plot(x2, z2, xlim=c(10,90), ylim=c(10,90), bty="n",
     main=sprintf("相関係数%.2fの散布図(n=%d)", cor(x2, z2), length(x2)),
     xlab="", ylab="", family="serif", las=1)

xの偏差

\[x_{i}-\bar{x}\;\;\;,\;\;i=1,2,3,...,n\]

# x軸とy軸の描画範囲を基本の散布図より広げている。
plot(x2, z2, xlim=c(0,100), ylim=c(0,100), bty="n",
     main="xの偏差の符号",
     xlab="", ylab="", family="serif", las=1)
abline(v=mean(x2), lty=2)
text(80,90,"xの偏差＋")
text(20,90,"xの偏差－")

yの偏差

\[y_{i}-\bar{y}\;\;\;,\;\;i=1,2,3,...,n\]

plot(x2, z2, xlim=c(0,100), ylim=c(0,100), bty="n",
     main="yの偏差の符号",
     xlab="", ylab="", family="serif", las=1)
abline(h=mean(z2), lty=2)
text(20,90,"yの偏差＋")
text(20,10,"yの偏差－")

4つの象限に分割した散布図

偏差積 \[(x_{i}-\bar{x})(y_{i}-\bar{y})\;\;\;,\;\;i=1,2,3,...,n\]

plot(x2, z2, xlim=c(0,100), ylim=c(0,100), bty="n",
     main=sprintf("相関係数%.2fの散布図(n=%d)", cor(x2, z2), length(x2)),
     xlab="", ylab="", family="serif", las=1)
abline(v=mean(x2), lty=2); abline(h=mean(z2), lty=2)
text(80,90,"第Ⅰ象限：\n偏差積＋")
text(20,90,"第Ⅱ象限：\n偏差積－")
text(20,10,"第Ⅲ象限：\n偏差積＋")
text(80,10,"第Ⅳ象限：\n偏差積－")

1-2 偏差積和と共分散(covariance)

偏差積 \[(x_{i}-\bar{x})(y_{i}-\bar{y}),\;\;\;\;i=1,2,3,...,n\] 偏差積和 \[\sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\bar{y})\] 共分散 \[s_{xy}=\color{red}{\frac{1}{n}\sum_{i=1}^{n}}(x_{i}-\bar{x})(y_{i}-\bar{y})\]

　上の式の色分けで示した様に，共分散は「偏差積の平均」である。
　共分散を幾つかの方法で計算して結果を比較してみよう。多くの統計ソフトの出力する共分散や分散が手計算と一致しない理由も合わせて説明している。

x <- rnorm(n <- 20, mean=0, sd=1)
y <- rnorm(n, mean=0, sd=1)

x # xの値自体を表示

 [1] -0.08383040 -0.07673258  1.48005025  0.58480852 -0.22706396 -0.20303145
 [7]  0.91240535  0.52202371  1.58490391  0.54029290 -0.67725767 -0.43695608
[13]  1.17936907 -1.10323971 -0.98884270  0.95186320  0.52679188  0.75386165
[19]  1.04791747 -0.69285921

x-mean(x) # xの（算術平均からの）偏差を表示

 [1] -0.3635541 -0.3564563  1.2003265  0.3050848 -0.5067877 -0.4827552
 [7]  0.6326816  0.2423000  1.3051802  0.2605692 -0.9569814 -0.7166798
[13]  0.8996454 -1.3829634 -1.2685664  0.6721395  0.2470682  0.4741379
[19]  0.7681938 -0.9725829

(x-mean(x))*(y-mean(y)) # xの偏差とyの偏差の偏差積を表示

 [1] -0.06447075  0.17831587  1.11994889 -0.60699099  0.16987186 -0.65509818
 [7] -0.56531956  0.03412091  0.26316239 -0.42237535 -0.14918153 -0.13416101
[13]  0.51757186 -0.73326587 -0.81181221 -0.02663410 -0.20752053  0.55333127
[19]  0.17899558  0.07715016

sum((x-mean(x))*(y-mean(y))) # 偏差積和を表示

[1] -1.284361

sum((x-mean(x))*(y-mean(y)))/length(x) # 共分散を表示・方法1

[1] -0.06421806

mean((x-mean(x))*(y-mean(y))) # 共分散を表示・方法2

[1] -0.06421806

# Rの標準関数で共分散を計算
cov(x, y)

[1] -0.06759796

# 値が一致しない！何故か？
# 統計ソフトは，分散・共分散を，nで割った値ではなく，n-1で割った値（不偏推定値）で示す。
# 以下の様に，n=length(x)ではなくn-1で割ると値が一致する。
sum((x-mean(x))*(y-mean(y)))/(length(x)-1) # 共分散を表示・方法3

[1] -0.06759796

1-3 ピアソンの積率相関係数(Pearson’s product-moment correlation coefficient)

積率相関係数 \[r_{xy}=\frac{s_{xy}}{s_{x}s_{y}}=\frac{\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})^2}\sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\bar{y})^2}}=\frac{1}{n}\sum_{i=1}^{n}(\frac{x_{i}-\bar{x}}{s_{x}})(\frac{y_{i}-\bar{y}}{s_{y}})\]

　3つの標準正規変量を発生させて演習を行おう。また，実際の社会調査データに近付ける為に，敢えて欠損値(NA，項目無回答item nonresponseと言う)を混ぜる操作をする。

x <- rnorm(n <- 20, mean=0, sd=1)
y <- rnorm(n, mean=0, sd=1)
z <- rnorm(n, mean=0, sd=1)

x[c( 1, 8,15)] <- NA # ヴェクトルxの 1, 8, 15番目の値をNAに置換
y[c( 3,10,17)] <- NA
z[c( 5,12,19)] <- NA

cor(x, y) # オプションを指定しないと，NAを含む計算は結果もNAになる。

[1] NA

cor(x, y, use="pairwise") # ペアのいずれにもNAを含まないものに限定

[1] 0.07128685

cor.test(x, y) # cor.test( )関数では自動でpairwise処理

    Pearson's product-moment correlation

data:  x and y
t = 0.24757, df = 12, p-value = 0.8086
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4773475  0.5799314
sample estimates:
       cor 
0.07128685 

# cor( )は3つ以上のヴェクトルの相関係数行列を計算できる。
cor(cbind(x, y, z), use="pairwise") # ペアでNAが無ければ有効

           x          y         z
x 1.00000000 0.07128685 0.1382841
y 0.07128685 1.00000000 0.5540965
z 0.13828412 0.55409650 1.0000000

cor(cbind(x, y, z), use="complete") 　# 全てのヴェクトルでNAを含まないケースだけが有効

           x          y          z
x 1.00000000 0.04422053 0.07294946
y 0.04422053 1.00000000 0.51002013
z 0.07294946 0.51002013 1.00000000

cor.test(y, z) # cor( )のpairwiseの方と結果が一致する事を確認。

    Pearson's product-moment correlation

data:  y and z
t = 2.3058, df = 12, p-value = 0.03978
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.03330983 0.83823996
sample estimates:
      cor 
0.5540965 

cor.test(z, x)

    Pearson's product-moment correlation

data:  z and x
t = 0.48368, df = 12, p-value = 0.6373
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4233574  0.6231432
sample estimates:
      cor 
0.1382841 

（補足）曲線的関係と積率相関係数

グラフ

script

set.seed(1968) # 結果を再現する為の乱数種の指定
cx <- rnorm(N <- 200, mean=50, sd=10)
cy <- 80*(cx-15)*(85-cx)/(15*85) + sample(c(-1:1), length(cx), replace=T) 

plot(cx, cy, xlim=c(20, 80), ylim=c(30, 90),
     main=paste("相関係数r=",round(cor(cx, cy, method="pearson"),2),", n=", N),
     xlab=paste("横軸変数の平均",round(mean(cx),1),", 標準偏差", round(sd(cx),1)), 
     ylab=paste("縦軸変数の平均",round(mean(cy),1),", 標準偏差", round(sd(cy),1)),
     col="#ff996680", pch=16, bty="n", family="serif")
segments(mean(cx), 30, mean(cx), 90, lty=2, lwd=2, col="#0000ff60")
segments(20, mean(cy), 90, mean(cy), lty=4, lwd=2, col="#0000ff60")
text( 80, 90, "第1象限;1st quadrant", adj=c(1, 1), col="#0000ff", family="serif", cex=0.7)
text( 20, 90, "第2象限;2nd quadrant", adj=c(0, 1), col="#0000ff", family="serif", cex=0.7)
text( 20, 30, "第3象限;3rd quadrant", adj=c(0, 0), col="#0000ff", family="serif", cex=0.7)
text( 80, 30, "第4象限;4th quadrant", adj=c(1, 0), col="#0000ff", family="serif", cex=0.7)

1-4 相関係数と線形変換(linear transformation)

\[v_{i}=ax_{i}+b\] \[\bar{v}=\frac{1}{n}\sum_{i=1}^{n}v_{i}=\frac{1}{n}\sum_{i=1}^{n}(ax_{i}+b)=\frac{1}{n}\sum_{i=1}^{n}(ax_{i})+\frac{1}{n}\sum_{i=1}^{n}b=a\frac{1}{n}\sum_{i=1}^{n}x_{i}+b=a\bar{x}+b\] \[\begin{alignat}{1} s_{v}^2&=\frac{1}{n}\sum_{i=1}^{n}(v_{i}-\bar{v})^2=\frac{1}{n}\sum_{i=1}^{n}(ax_{i}+b-a\bar{x}-b)^2\\ &=\frac{1}{n}\sum_{i=1}^{n}(ax_{i}-a\bar{x})^2=a^2\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})^2=a^2s_{x}^2\\ \end{alignat}\]

\[s_{v}=|a|s_{x}\]

　乱数を利用して，任意の相関関係にある2変数を発生させ，それを線形変換すると各種の統計量がどうなるかを簡単に示してみる。 　相関係数rは-1から1の範囲で好きな値を指定出来る。そうすると，厳密にその相関関係になる変数x1とy1が乱数発生される。 　次に，x1とy1から，一部予想不能な乱数を含んで線形変換したx2とy2が生成される。元のx1，y1と線形変換後のx2，y2で，平均，標準偏差，共分散などは全て変化してしまうが，相関係数は不変である事が確認出来る。

\[\begin{alignat}{2} v_{i}&=ax_{i}&&+b\\ w_{i}&=cy_{i}&&+d\\ \end{alignat}\]

\[\begin{alignat}{1} s_{vw}&=\frac{1}{n}\sum_{i=1}^{n}(v_{i}-\bar{v})(w_{i}-\bar{w})=\frac{1}{n}\sum_{i=1}^{n}(ax_{i}-a\bar{x})(cy_{i}-c\bar{y})\\ &=ac\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\bar{y})=ac\cdot s_{xy}\\ \end{alignat}\]

\[r_{vw}=\frac{s_{vw}}{s_{v}s_{w}}=\frac{ac\cdot s_{xy}}{|a|s_{x}\cdot|c|s_{y}}=\frac{ac}{|ac|}r_{xy}\]

r <- 0.7 # 最初に相関係数の値を任意に設定
n <- 400 # ケース数を指定
u <- rnorm(n, mean=0, sd=1) # 材料としてのヴェクトルu
x0 <- rnorm(n, mean=0, sd=1); x1 <- scale(x0) # このxは後で使用
e0 <- lm(u ~ x1)$residuals; e <- scale(e0) # 残差を取り出す
y1 <- sqrt(1-r^2)*e + r*x1 # x1とy1の相関が厳密にrになる
cor(x1, y1) # 確かに設定した相関になっている事を確認

     [,1]
[1,]  0.7

plot(x1, y1)　# 簡単に散布図を描いてみる

# xとyを，乱数を交えて線型変換してみる（結果が予め決まっていない）
x2 <- runif(1, min=1.1, max=5)*x1 + runif(1, min=-2, max=2)
y2 <- runif(1, min=0.1, max=0.9)*y1 + runif(1, min=-2, max=2)

mean(x1); mean(x2) # 変換前後で平均が変わる

[1] -1.940153e-17

[1] 0.003289395

sd(x1); sd(x2) # 変換前後で標準偏差も変わる

[1] 1

[1] 4.526269

mean(y1); mean(y2)

[1] -3.172159e-17

[1] 1.405561

sd(y1); sd(y2)

[1] 1

[1] 0.487787

cov(x1, y1); cov(x2, y2) # 変換前後で共分散も変わる

     [,1]
[1,]  0.7

         [,1]
[1,] 1.545499

cor(x1, y1); cor(x2, y2) # 変換前後で相関係数は不変

     [,1]
[1,]  0.7

     [,1]
[1,]  0.7

plot(x2, y2)　# 簡単に散布図を描いてみる

線形変換の一例としての標準化

既に学んだ標準化(standardization)も線形変換（一次変換）である。\(\bar{x}\)と\(s_{x}\)が定数である事に注意すれば分かる。

\[\color{red}{z_{i}}=\frac{\color{red}{x_{i}}-\color{blue}{\bar{x}}}{\color{blue}{s_{x}}}=\frac{1}{\color{blue}{s_{x}}}\color{red}{x_{i}}-\frac{\color{blue}{\bar{x}}}{\color{blue}{s_{x}}},\;\;\;\;\;i=1,2,3,...,n\]

1-5 相関係数行列と散布図行列

　上の1-3で発生させたx, y, zの三つのヴェクトルで例示してみよう。既に相関係数行列は上で示した。有効ケースの判断として“pairwise”と“complete”の二通りがあり，結果が異なる事に注意しなければならない。

# cbind( )は「列columnを結合する」と云う関数
# ヴェクトルを列結合した結果を u と云うオブジェクトに代入している。
# 以下のコマンドを実行した後で，R Consoleで u と打つとuが何かが分かる。
cor(u <- cbind(x, y, z), use="pairwise") 

           x          y         z
x 1.00000000 0.07128685 0.1382841
y 0.07128685 1.00000000 0.5540965
z 0.13828412 0.55409650 1.0000000

cor(u, use="complete")

           x          y          z
x 1.00000000 0.04422053 0.07294946
y 0.04422053 1.00000000 0.51002013
z 0.07294946 0.51002013 1.00000000

# ここで3つのヴェクトルを含むuに，散布図を描く plot( )関数を適用しても，ひとつの散布図しか描かれない。
# plot(x, y)と見比べると同じである事が分かる。
plot(u)
plot(x, y)

# 散布図行列を描くには，単に列結合するだけでなく，それを「データ・フレイム」に変換しなければならない。
w <- data.frame(u) # 3つのヴェクトルを列結合したuをデータ・フレイム形式にする。
plot(w) # 散布図行列を描く

　3つのヴェクトルのうちの2つの組合せで散布図が描かれているが，横軸と縦軸が入れ替わっているだけのものがあるので3組6枚の散布図になっている。次の様に2つのヴェクトルでの散布図を3組描いて，よく見比べるとどれがどれなのか分かるだろう。

par(mfrow=c(1,3)) # グラフ画面を1行3列に分割する
plot(x, y)
plot(x, z)
plot(y, z)

1-5-1 “psych”パッケイジのpairs.panels( )関数

　追加パッケイジのpsychを使用すると，もっと情報量の多い散布図行列が描ける。下の例では，対角線に各変数のヒストグラム，右上三角に2変数ごとの積率相関係数，左下三角に2変数ごとの散布図と回帰直線が描かれている。

# install.packages("psych", repos="https://cran.ism.ac.jp")
library(psych)

Attaching package: 'psych'

The following objects are masked from 'package:ggplot2':

    %+%, alpha

pairs.panels(w, scale=F, lm=T, ellipses=F)

2 分割表(contingency table)（クロス表cross-tabulation）

2-1 分割表(contingency table)の構成と読み取り

　データ・フレイム“d01”内の変数“sex”と“q1900”（大文字と小文字も区別されるので注意！）の分割表（クロス表）を作成するには，with(d01, table(sex, q1900))とする。これをオブジェクト（“以下ではt01”とする）に代入して表示すると良い。

# まずはデータファイルから読み込んでデータフレイムを作成する。
d01 <- read.csv("practice.csv")

with(d01, table(sex, q1900)) # そのまま表示して終わり

   q1900
sex  1  2  3  4
  1  6 25 26 33
  2  6 27 37 48

t01 <- with(d01, table(sex, q1900)) # オブジェクトに保存してから表示
t01

   q1900
sex  1  2  3  4
  1  6 25 26 33
  2  6 27 37 48

　代入文だけでは結果を表示しないことには注意が必要である。代入した結果も表示したい場合には，上のように単にオブジェクト名を指定するか，もしくは代入文全体を括弧 ( ) に入れて，(t01 <- table(d01\(sex, d01\)q1900)) などとする。 　上の例で分かるように，table( )は表の本体を表示するだけの非常に素っ気ない関数である。これに周辺度数を付加する関数が addmargins( )である。
　比率を計算するのは prop.table( )である。全体比率，行比率，列比率の3つがある事に注意しよう。

addmargins(t01)

     q1900
sex     1   2   3   4 Sum
  1     6  25  26  33  90
  2     6  27  37  48 118
  Sum  12  52  63  81 208

prop.table(t01) # 全体比率

   q1900
sex          1          2          3          4
  1 0.02884615 0.12019231 0.12500000 0.15865385
  2 0.02884615 0.12980769 0.17788462 0.23076923

prop.table(t01, margin=1) # 行比率

   q1900
sex          1          2          3          4
  1 0.06666667 0.27777778 0.28888889 0.36666667
  2 0.05084746 0.22881356 0.31355932 0.40677966

prop.table(t01, margin=2) # 列比率

   q1900
sex         1         2         3         4
  1 0.5000000 0.4807692 0.4126984 0.4074074
  2 0.5000000 0.5192308 0.5873016 0.5925926

　全体比率の表に周辺度数を付加すると次のようになる。

addmargins(prop.table(t01))

     q1900
sex            1          2          3          4        Sum
  1   0.02884615 0.12019231 0.12500000 0.15865385 0.43269231
  2   0.02884615 0.12980769 0.17788462 0.23076923 0.56730769
  Sum 0.05769231 0.25000000 0.30288462 0.38942308 1.00000000

　「小数ではなく％にして，その小数第一位までだけ表示したい」という場合には，次の様にすると良い。round( )は四捨五入の関数であり，最初の引数に数値オブジェクト（この場合は分割表）を指定し，二番目の引数に小数点以下の桁数を指定する。小数第一位までとすると，小数第二位を四捨五入する。上の表とよく見比べてみよう。

t02 <- addmargins(prop.table(t01))
round(t02*100, 1)

     q1900
sex       1     2     3     4   Sum
  1     2.9  12.0  12.5  15.9  43.3
  2     2.9  13.0  17.8  23.1  56.7
  Sum   5.8  25.0  30.3  38.9 100.0

　関数 table( ) は，欠損値NAは自動的に除外して集計する設定となっている。しばしばこれでは不十分であるので，NAを含んで全てのケースについての分割表を作成するオプション useNA=“no|ifany|always” を紹介する。このオプションを省略した場合（デフォルト）では，useNA=“no” となる。

table(d01$sex, d01$q1900) # NA除外

   
     1  2  3  4
  1  6 25 26 33
  2  6 27 37 48

table(d01$sex, d01$q1900, useNA="ifany") # NAがあれば表示

   
      1   2   3   4 <NA>
  1   6  25  26  33   75
  2   6  27  37  48  101

table(d01$sex, d01$q1900, useNA="always") # NAを常に表示

      
         1   2   3   4 <NA>
  1      6  25  26  33   75
  2      6  27  37  48  101
  <NA>   0   0   0   0    0

　ついでに，分割表にラベルを貼るやり方も紹介しておこう。 分割表を t01 とする時，names(dimnames(t01)) は行変数の名前と列変数のラベル，dimnames(t01)[[1]] は行変数のカテゴリのラベル，dimnames(t01)[[2]] は列変数のカテゴリのラベルを意味する。 　以下のスクリプトと結果をよく見て，ラベルを付ける前と付けた後を比較して欲しい。列のラベルと数値がずれている様に見えるが，コピーしてみると大丈夫である。

(t01 <- table(d01$sex, d01$q1900)) # 代入して表示

   
     1  2  3  4
  1  6 25 26 33
  2  6 27 37 48

names(dimnames(t01)) <- list("性別", "性別役割分業")
dimnames(t01)[[1]] <- c("1 男性", "2 女性")
dimnames(t01)[[2]] <- c("1 そう思う", "2 どち…思う", "3 どち…思わない", "4 思わない")
t01

        性別役割分業
性別     1 そう思う 2 どち…思う 3 どち…思わない 4 思わない
  1 男性          6           25               26         33
  2 女性          6           27               37         48

NAを含む場合は以下の様にする。

(t03 <- table(d01$sex, d01$q1900, useNA="always"))

      
         1   2   3   4 <NA>
  1      6  25  26  33   75
  2      6  27  37  48  101
  <NA>   0   0   0   0    0

names(dimnames(t03)) <- list("性別", "性別役割分業")
dimnames(t03)[[1]] <- c("1 男性", "2 女性", "NA")
dimnames(t03)[[2]] <- c("1 そう思う", "2 どち…思う", "3 どち…思わない", "4 思わない", "NA")
t03

        性別役割分業
性別     1 そう思う 2 どち…思う 3 どち…思わない 4 思わない  NA
  1 男性          6           25               26         33  75
  2 女性          6           27               37         48 101
  NA              0            0                0          0   0

2-1-1 3元分割表（3重クロス表）

　ここでは2変数の連関を説明しているが，3変数以上の分割表（連関表，クロス表）についても紹介しておこう。 　以下は，演習用データ・ファイルをdata01と云う名前のデータ・フレイムに読み込んで，性別（層変数）×学歴（行変数）×職業（列変数）の3元分割表を作成する例である。 　まずは，table( )関数で単純に分割表を作成，周辺度数を付加した分割表を示す。

d01 <- read.csv("practice.csv") # csvファイルをデータフレイムd01として読み込む
t04 <- table(d01$edu2, d01$job, d01$sex)
names(dimnames(t04)) <- list("学歴", "職業", "性別")
dimnames(t04)[[3]] <- c("1 男性", "2 女性")
dimnames(t04)[[2]] <- c("1 正規", "2 非正規", "3 自営", "4 無職")
dimnames(t04)[[1]] <- c("1 高校", "2 短大", "3 四大")
t04

, , 性別 = 1 男性

        職業
学歴     1 正規 2 非正規 3 自営 4 無職
  1 高校     37        6     10      3
  2 短大     25        1      4      1
  3 四大     64        2      4      3

, , 性別 = 2 女性

        職業
学歴     1 正規 2 非正規 3 自営 4 無職
  1 高校     12       31      7     22
  2 短大     22       36     10     24
  3 四大     20       15      2     14

m04 <- addmargins(t04)
m04

, , 性別 = 1 男性

        職業
学歴     1 正規 2 非正規 3 自営 4 無職 Sum
  1 高校     37        6     10      3  56
  2 短大     25        1      4      1  31
  3 四大     64        2      4      3  73
  Sum       126        9     18      7 160

, , 性別 = 2 女性

        職業
学歴     1 正規 2 非正規 3 自営 4 無職 Sum
  1 高校     12       31      7     22  72
  2 短大     22       36     10     24  92
  3 四大     20       15      2     14  51
  Sum        54       82     19     60 215

, , 性別 = Sum

        職業
学歴     1 正規 2 非正規 3 自営 4 無職 Sum
  1 高校     49       37     17     25 128
  2 短大     47       37     14     25 123
  3 四大     84       17      6     17 124
  Sum       180       91     37     67 375

　この3元分割表を2次元に縮約して表示する関数として ftable( ) がある。上の3元分割表はそれぞれ以下の様になる。

ft04 <- ftable(t04); ft04

                性別 1 男性 2 女性
学歴   職業                       
1 高校 1 正規            37     12
       2 非正規           6     31
       3 自営            10      7
       4 無職             3     22
2 短大 1 正規            25     22
       2 非正規           1     36
       3 自営             4     10
       4 無職             1     24
3 四大 1 正規            64     20
       2 非正規           2     15
       3 自営             4      2
       4 無職             3     14

fm04 <- ftable(m04); fm04

                性別 1 男性 2 女性 Sum
学歴   職業                           
1 高校 1 正規            37     12  49
       2 非正規           6     31  37
       3 自営            10      7  17
       4 無職             3     22  25
       Sum               56     72 128
2 短大 1 正規            25     22  47
       2 非正規           1     36  37
       3 自営             4     10  14
       4 無職             1     24  25
       Sum               31     92 123
3 四大 1 正規            64     20  84
       2 非正規           2     15  17
       3 自営             4      2   6
       4 無職             3     14  17
       Sum               73     51 124
Sum    1 正規           126     54 180
       2 非正規           9     82  91
       3 自営            18     19  37
       4 無職             7     60  67
       Sum              160    215 375

　この様に簡潔な表形式で3元分割表の情報を表示する事が出来るが，度数だけではやや読み取りにくい。しかし比率を表示させるにはやや工夫が必要になる。
　これだと周辺度数に相当する部分は表示されないものの，男性と女性の違いはかなり分かり易く読み取れる。

p04 <- t04 # 比率の表の入れ物を用意する
p04[,,1] <- prop.table(t04[,,1], margin=1)
p04[,,2] <- prop.table(t04[,,2], margin=1)
p04

, , 性別 = 1 男性

        職業
学歴         1 正規   2 非正規     3 自営     4 無職
  1 高校 0.66071429 0.10714286 0.17857143 0.05357143
  2 短大 0.80645161 0.03225806 0.12903226 0.03225806
  3 四大 0.87671233 0.02739726 0.05479452 0.04109589

, , 性別 = 2 女性

        職業
学歴         1 正規   2 非正規     3 自営     4 無職
  1 高校 0.16666667 0.43055556 0.09722222 0.30555556
  2 短大 0.23913043 0.39130435 0.10869565 0.26086957
  3 四大 0.39215686 0.29411765 0.03921569 0.27450980

fp04 <- round(ftable(p04)*100, 1)
fp04

                性別 1 男性 2 女性
学歴   職業                       
1 高校 1 正規          66.1   16.7
       2 非正規        10.7   43.1
       3 自営          17.9    9.7
       4 無職           5.4   30.6
2 短大 1 正規          80.6   23.9
       2 非正規         3.2   39.1
       3 自営          12.9   10.9
       4 無職           3.2   26.1
3 四大 1 正規          87.7   39.2
       2 非正規         2.7   29.4
       3 自営           5.5    3.9
       4 無職           4.1   27.5

2-2 分割表を表現するグラフ

　本文で紹介した，棒グラフ，帯グラフ，モザイクプロットのスクリプトと実例を紹介する。 まずは分割表を準備する。

t05 <- table(d01$sex, d01$q1900) # オブジェクト名tbl
names(dimnames(t05)) <- list("性別", "性別役割分業")
dimnames(t05)[[1]] <- c("1 男性", "2 女性")
dimnames(t05)[[2]] <- c("1 そう思う", "2 どち…思う", "3 どち…思わない", "4 思わない")
t05 # まず分割表（特に行変数と列変数）を確認する。

        性別役割分業
性別     1 そう思う 2 どち…思う 3 どち…思わない 4 思わない
  1 男性          6           25               26         33
  2 女性          6           27               37         48

　まずは棒グラフ barplot( ) を使用する。初めての場合はとにかく最も単純なものを実行して結果を解読するのが良いだろう。その後，必要に応じて色々なオプションを指定する事を学んでいく。

barplot(t05) # 積み上げ棒グラフ

　分割表に対してそのまま棒グラフを描くと，この様に積み上げ棒グラフになる。列変数ごとに，行変数の値を積み上げている事が分かる。横に並んだ単純な棒グラフを描きたい場合は，beside=T というオプションが必要になる。

barplot(t05, beside=T) #単純な棒グラフ

　積み上げは積み上げでも，男女別に積み上げたい場合は，転置行列の関数 t( ) を使用する。

barplot(t(t05)) # 行と列を入れ替えたいので，転置行列の関数 t( ) を使用。

　凡例が無いと見づらいので，legend.text= オプションで凡例を指定する。

barplot(t(t05), legend.text=dimnames(t05)[[2]]) # 凡例を入れる

　パーセンテージで表現する「帯グラフ」を描くには，分割表を比率にして積み上げ棒グラフを描けば良い。まずは prop.table(t05, margin=1) によって行比率の分割表に変換し，それを上と同様に転置する t(prop.table(t05, margin=1)) 。それを凡例付きで積み上げ棒グラフにする。序に色も付けておこう。

barplot(t(prop.table(t05, margin=1)), legend.text=dimnames(t05)[[2]], col=terrain.colors(4)) # 帯グラフ（比率にして積み上げ）

　最後に，モザイクプロットを描く。

mosaicplot(t05) # モザイクプロット

　これでは余りに素っ気ないので，幾つかオプションを指定して装飾しよう。最後の二行の text( ) は度数を書き込んでいるものだが，中身はやや難しいので分からなくて良い。

mosaicplot(t05,
     main="男女ごとの性別役割意識モザイクプロット",
     col=terrain.colors(4))
text(sum(t05[1,]*.5)/sum(t05), 1-cumsum(t05[1,])/sum(t05[1,]), t05[1,])
text(sum(t05[1,]+t05[2,]*.5)/sum(t05), 1-cumsum(t05[2,])/sum(t05[2,]), t05[2,])

2-3 名義尺度(nominal scale)の連関係数(association coefficient)

　まず演習用に2行2列の分割表を用意しよう。後の式の見易さの為にオブジェクト名は一文字でxとしておく。行列を作成する matrix( ) 関数の書式をよく確認しておくと良い。

x <- matrix(c(10, 8, 5, 6), ncol=2, nrow=2)
x

     [,1] [,2]
[1,]   10    5
[2,]    8    6

　テクストの定義式に従ってユールのQを計算すると次のようになる。 \[Q=\frac{ad-bc}{ad+bc}\]

Q <- (x[1,1]*x[2,2]-x[1,2]*x[2,1])/(x[1,1]*x[2,2]+x[1,2]*x[2,1])
Q

[1] 0.2

　ユールのQやφ係数の分子 ad-bc は行列式と呼ばれるもので，det( )で求められる。

det(x)/(x[1,1]*x[2,2]+x[1,2]*x[2,1])

[1] 0.2

　別の2×2分割表でQを計算したい時，行列 x の中身だけ変えれば，これらの計算式がそのまま利用出来る。

※ table( )関数で作成した分割表には det( ) は使えないので，その場合には table( ) から matrix( ) に変換する。

# data.frameから分割表を作成する
t06 <- table("sex"=d01$sex, "univ"=d01$edu1)
t06

   univ
sex   0   1
  1  87  74
  2 164  51

# tableオブジェクトにはdet()関数は使えない
det(t06)

Error in UseMethod("determinant"):  'determinant' をクラス "table" のオブジェクトに適用できるようなメソッドがありません 

## クラスtableのオブジェクトをクラスmatrixに変換
x <- matrix(t06, ncol=2)
det(x)

[1] -7699

次はφ係数を計算してみよう。色々な計算の仕方が出来る。

\[\phi=\frac{ad-bc}{\sqrt{(a+b)(c+d)(a+c)(b+d)}}\]

phi <- det(x)/sqrt(sum(x[1,])*sum(x[2,])*sum(x[,1])*sum(x[,2]))
phi

[1] -0.2336201

x2 <- addmargins(x); x2

            Sum
     87  74 161
    164  51 215
Sum 251 125 376

det(x)/sqrt(x2[3,1]*x2[3,2]*x2[1,3]*x2[2,3])

[1] -0.2336201

本文026頁に間違いを発見しました。

最後はオッズ比である。度数が0のセルが無ければ単純である。

\[odds\;ratio=\frac{\frac{a}{b}}{\frac{c}{d}}=\frac{ad}{bc}\]

OR <- x[1,1]*x[2,2]/(x[1,2]*x[2,1])
OR

[1] 0.3656065

　0セルがあるかどうかは，すべてのセルの積を prod(x) で計算し，これが0になるかどうかで判断出来るので，prod(x)が0であれば0.5，正の値であれば0となる変数dを作成して，これをすべてのセルに加えてオッズ比を計算すれば良い。

x <- matrix(c(10, 8, 1, 6), ncol=2, nrow=2)
d <- 0; if (prod(x)==0) d <- 0.5
OR <- (x[1,1]+d)*(x[2,2]+d)/((x[1,2]+d)*(x[2,1]+d))
x; d; OR

     [,1] [,2]
[1,]   10    1
[2,]    8    6

[1] 0

[1] 7.5

x <- matrix(c(10, 8, 0, 6), ncol=2, nrow=2)
d <- 0; if (prod(x)==0) d <- 0.5
OR <- (x[1,1]+d)*(x[2,2]+d)/((x[1,2]+d)*(x[2,1]+d))
x; d; OR

     [,1] [,2]
[1,]   10    0
[2,]    8    6

[1] 0.5

[1] 16.05882

　以下のようにすれば，判別条件を含んで一つのコマンドで書く事も出来るが，単に好みの問題であり，分からなくても支障ない［多少洗練させた.2020.08.30］。

x <- matrix(c(10, 8, 1, 6), ncol=2, nrow=2)
OR <- (x[1,1]+(d <- 0.5*(prod(x)==0)))*(x[2,2]+d)/((x[1,2]+d)*(x[2,1]+d))
OR

[1] 7.5

x <- matrix(c(10, 8, 0, 6), ncol=2, nrow=2)
OR <- (x[1,1]+(d <- 0.5*(prod(x)==0)))*(x[2,2]+d)/((x[1,2]+d)*(x[2,1]+d))
OR

[1] 16.05882

2-3-1 対数オッズ比

　のちのロジスティック回帰分析などでは，オッズ比そのものではなく，オッズ比の自然対数である「対数オッズ比」が登場する事になる。オッズ比は正の実数で，1を境に非対称であるが，対数オッズ比は0を境に対称的に変化する。

\[log(odds\;ratio)=log\left(\frac{ad}{bc}\right)\]

x <- matrix(c(6, 3, 8, 2), ncol=2, nrow=2)
OR <- (x[1,1]+(d <- 0.5*(1 - ceiling(prod(x)/(prod(x)+.1)))))*(x[2,2]+d)/((x[1,2]+d)*(x[2,1]+d))
LogOR <- log(OR)
x; OR; LogOR

     [,1] [,2]
[1,]    6    8
[2,]    3    2

[1] 0.5

[1] -0.6931472

x <- matrix(c(8, 6, 2, 3), ncol=2, nrow=2)
OR <- (x[1,1]+(d <- 0.5*(1 - ceiling(prod(x)/(prod(x)+.1)))))*(x[2,2]+d)/((x[1,2]+d)*(x[2,1]+d))
LogOR <- log(OR)
x; OR; LogOR

     [,1] [,2]
[1,]    8    2
[2,]    6    3

[1] 2

[1] 0.6931472

　0～+∞の値を取りうるオッズ比と，-∞～+∞の値を取りうる対数オッズ比の関係をグラフで確認してみよう。描画の都合上，x軸のオッズ比は0.1から10の範囲を指定する。オッズ比が0.1なら対数オッズ比は-2.302585，オッズ比が10なら対数オッズ比は2.302585である。

curve(log(x), from=0.1, to=10, bty="n", col="#d8f255",
     main="オッズ比と対数オッズ比の関係", xlab="OR; Odd Ratio", ylab="log(OR)")
abline(h=0, lty=2)
abline(v=1, lty=2)

2-4 独立性(independence)とカイ二乗統計量(chi-square statistic)

分割表(contingency table)；クロス表

観測度数（実現度数）の表

	そう思う	やや思う	あまり思わない	思わない	行周辺度数
高校まで	n11	n12	n13	n14	n1・
短大・高専	n21	n22	n23	n24	n2・
四大・院	n31	n32	n33	n34	n3・
列周辺度数	n・1	n・2	n・3	n・4	n

期待度数（理論度数）の表: counterfactual

	そう思う	やや思う	あまり思わない	思わない	行周辺度数
高校まで	f11	f12	f13	f14	n1・
短大・高専	f21	f22	f23	f24	n2・
四大・院	f31	f32	f33	f34	n3・
列周辺度数	n・1	n・2	n・3	n・4	n

期待度数 \(f_{ij}=\frac{n_{i\cdot}\times n_{\cdot j}}{n}\)

　分割表について，カイ二乗統計量などを計算する方法を説明しよう。Rのカイ二乗検定の関数 chisq.test( ) を使えば簡単に求められるが，より手計算に近い方法も紹介しよう。

\[\chi^2=\sum_{i=1}^{I}\sum_{j=1}^{J}\left(\frac{x_{ij}-f_{ij}}{\sqrt{f_{ij}}}\right)^2=\sum_{i=1}^{I}\sum_{j=1}^{J}\frac{(x_{ij}-f_{ij})^2}{f_{ij}}\]

# 分割表を準備
x <- matrix(c(6, 3, 3, 17, 15, 19, 26, 21, 15, 25, 22, 31), nrow=3)
# 周辺度数を付加
addx <- addmargins(x); addx

                Sum
     6 17 26 25  74
     3 15 21 22  61
     3 19 15 31  68
Sum 12 51 62 78 203

# 期待度数の表を，ヴェクトル同士の積 %*% を使用して作成する。
f <- addx[1:3, 5] %*% t(addx[4, 1:4])/sum(x)
f

                                        
[1,] 4.374384 18.59113 22.60099 28.43350
[2,] 3.605911 15.32512 18.63054 23.43842
[3,] 4.019704 17.08374 20.76847 26.12808

# 残差の表
x-f

                                              
[1,]  1.6256158 -1.5911330  3.399015 -3.433498
[2,] -0.6059113 -0.3251232  2.369458 -1.438424
[3,] -1.0197044  1.9162562 -5.768473  4.871921

# 標準化残差（ピアソン残差）の表
stdres <- (x-f)/sqrt(f); stdres

                                                
[1,]  0.7772477 -0.3690231  0.7149725 -0.6439047
[2,] -0.3190814 -0.0830512  0.5489544 -0.2971137
[3,] -0.5086010  0.4636199 -1.2657810  0.9531177

# 標準化残差の二乗和，即ちカイ二乗統計量
sum(stdres^2)

[1] 5.148682

　chisq.test( )関数を使用し，そこから情報を取り出す方法を以下で紹介する。カイ二乗検定の結果の要点(c0)が表示され，カイ二乗値が上で求めた値と（四捨五入の範囲内で）一致している事が確認出来る。検定の自由度は6，有意確率は.5249となっている。

# カイ二乗検定の関数
c0 <- chisq.test(x); c0

Warning in chisq.test(x): Chi-squared approximation may be incorrect

    Pearson's Chi-squared test

data:  x
X-squared = 5.1487, df = 6, p-value = 0.5249

names(c0)

[1] "statistic" "parameter" "p.value"   "method"    "data.name" "observed" 
[7] "expected"  "residuals" "stdres"   

　names( )でc0に格納されている情報名の一覧が表示される。情報の取り出し方は以下を参考にして欲しい。また，上の手計算の結果と照合して欲しい。

c0$statistic # カイ二乗統計量

X-squared 
 5.148682 

c0$observed # 観測度数の表

     [,1] [,2] [,3] [,4]
[1,]    6   17   26   25
[2,]    3   15   21   22
[3,]    3   19   15   31

c0$expected # 期待度数の表

         [,1]     [,2]     [,3]     [,4]
[1,] 4.374384 18.59113 22.60099 28.43350
[2,] 3.605911 15.32512 18.63054 23.43842
[3,] 4.019704 17.08374 20.76847 26.12808

c0$residuals # 標準化残差の表（ただの残差ではない！）

           [,1]       [,2]       [,3]       [,4]
[1,]  0.7772477 -0.3690231  0.7149725 -0.6439047
[2,] -0.3190814 -0.0830512  0.5489544 -0.2971137
[3,] -0.5086010  0.4636199 -1.2657810  0.9531177

c0$stdres # 調整済み標準化残差（Haberman残差）の表

           [,1]       [,2]       [,3]       [,4]
[1,]  1.0051802 -0.5349743  1.0761703 -1.0293612
[2,] -0.3933116 -0.1147562  0.7875508 -0.4527094
[3,] -0.6429685  0.6570061 -1.8624223  1.4894319

sum(c0$residuals^2) # 標準化残差の二乗和＝カイ二乗値

[1] 5.148682

　クラメール（クラマー）のVの計算には，カイ二乗統計量を利用する。 \[V=\sqrt{\frac{\chi^2}{n\times min(I-1, J-1)}}\]

V <- sqrt(c0$statistic/(min(dim(x)-1)*sum(x)))
names(V) <- "Cramer's V" # 出力のラベル表示を修正するだけなので省略可
V

Cramer's V 
 0.1126121 

2-5 順序尺度(ordinal scale)の連関係数(rank correlation coefficient)

　任意の分割表から，本文にあるペアのタイプのA，B，C，D，Eの数を全て数え上げ，それを元にグッドマンとクラスカルのγや，ケンドールのτ3種類を計算するスクリプトを，前半と後半に分けて紹介する。

# 分割表を準備
x0 <- matrix(c(6, 3, 3, 17, 15, 19, 26, 21, 15, 25, 22, 31), nrow=3)

## rank Correlation
table01 <- x0     # 元になる分割表を x0 とし，table01にコピーして以下で使用。
n.col <- ncol(x0) # x0の列数
n.row <- nrow(x0) # x0の行数

sum01 <- sum(table01) # 全度数
sum.pair <- sum01*(sum01-1)/2　 # ペアの総数
# それぞれのタイプのペアの数の初期値を全て0として定義しておく。
con.pair <- dis.pair <- r.tie <- c.tie <- all.tie <- 0

# 各タイプのペアの数を数え上げる
for (i in 1:n.row) {
 for (i2 in 1:n.row) {
  for (j in 1:n.col)  {
   for (j2 in 1:n.col) {
    if (i>i2 & j>j2) con.pair <- con.pair+table01[i,j]*table01[i2,j2]
    if (i>i2 & j<j2) dis.pair <- dis.pair+table01[i,j]*table01[i2,j2]
    if (i==i2 & j>j2) r.tie <- r.tie+table01[i,j]*table01[i2,j2]
    if (i>i2 & j==j2) c.tie <- c.tie+table01[i,j]*table01[i2,j2]
    if (i==i2 & j==j2) all.tie <- all.tie+table01[i,j]*(table01[i2,j2]-1)/2
 }}}}

# それぞれのペアの数を表示
con.pair; dis.pair; r.tie; c.tie; all.tie

[1] 5165

[1] 4363

[1] 4740

[1] 4166

[1] 2069

# ペアの総数を確認
all.tie+r.tie+c.tie+con.pair+dis.pair

[1] 20503

sum.pair

[1] 20503

　これらからつぎに，4種類の順序連関係数を計算する。

\[Goodman\;\&\;Kruskal's\; \gamma=\frac{n(A)-n(B)}{n(A)+n(B)}\] \[Kendall's\;\tau_{b}=\frac{n(A)-n(B)}{\sqrt{n(A)+n(B)+n(C)}\sqrt{n(A)+n(B)+n(D)}}\]

## Goodman and Kruskal's Gamma
gamma <- (con.pair-dis.pair)/(con.pair+dis.pair)
## Kendall's tau-a
taua <- (con.pair-dis.pair)/sum.pair
## Kendall's tau-b
taub <- (con.pair-dis.pair)/sqrt((con.pair+dis.pair+r.tie)*(con.pair+dis.pair+c.tie))
## Kendall's tau-c
tauc <- (con.pair-dis.pair)*2*min(dim(table01))/(sum01^2*(min(dim(table01))-1))
gamma; taua; taub; tauc

[1] 0.08417296

[1] 0.03911623

[1] 0.05737565

[1] 0.0583853

相関係数の関数でKendallのtau-bを計算

Kendallの\(\tau_{b}\)だけはRの標準の相関係数関数でローデータから計算できる。
上記の分割表は演習用データpractice.csvの学歴3区分edu2と性別役割分業意識q1900の分割表である。

# 上の分割表と同じ分割表
with(d01, table(edu2, q1900)) # 分割表の確認

    q1900
edu2  1  2  3  4
   1  6 17 26 25
   2  3 15 21 22
   3  3 19 15 31

with(d01,      cor(edu2, q1900, method="kendall", use="pairwise"))

[1] 0.05737565

with(d01, cor.test(edu2, q1900, method="kendall", use="pairwise"))

    Kendall's rank correlation tau

data:  edu2 and q1900
z = 0.9294, p-value = 0.3527
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
0.05737565 

発展1 関連の有無について注意すべきこと

発展1-1　生態学的誤謬(ecological fallacy)

　少し煩雑だが，生態学的誤謬を視覚的に確認するスクリプトを作ってみたので，#### ここまでで一旦実行してグラフを確認 ####の指示に従って，スクリプトを部分的に実行しながらグラフを確認して行って欲しい。

rxy <- -0.7 # グループ内での個人レヴェルでの相関
x0 <- rnorm(30, mean=0, sd=1); y0 <- rnorm(30, mean=0, sd=1) # 10人×3グループ分のタネ

## 第1グループ
z1 <- rxy*(x1 <- scale(x0[1:10])) + sqrt(1-rxy^2)*scale(lm(y0[1:10] ~ x0[1:10])$residuals)
x11 <- x1*10+58; z11 <- z1*10+68

## 第2グループ
z1 <- rxy*(x1 <- scale(x0[11:20])) + sqrt(1-rxy^2)*scale(lm(y0[11:20] ~ x0[11:20])$residuals)
x12 <- x1*10+50; z12 <- z1*10+55

## 第3グループ
z1 <- rxy*(x1 <- scale(x0[21:30])) + sqrt(1-rxy^2)*scale(lm(y0[21:30] ~ x0[21:30])$residuals)
x13 <- x1*10+45; z13 <- z1*10+40


# 3つのグループごとのプロット
plot(x11,z11, xlim=c(0,100), ylim=c(0,100), bty="n", pch=1,
     main="生態学的相関", xlab="", ylab="", family="serif", las=1)
text(50, 95, "それぞれのグループ内でのピアソン相関は全て-0.7")
#### ここまでで一旦実行してグラフを確認 ####

par(new=T); plot(x12,z12, xlim=c(0,100), ylim=c(0,100), bty="n",
     main="", xlab="", ylab="", family="serif", las=1, pch=17)
#### ここまでで一旦実行してグラフを確認 ####

par(new=T); plot(x13,z13, xlim=c(0,100), ylim=c(0,100), bty="n",
     main="", xlab="", ylab="", family="serif", las=1, pch=2)
#### ここまでで一旦実行してグラフを確認 ####

# 全体での個人レヴェルでの相関係数
x_ind <- c(x11, x12, x13)
z_ind <- c(z11, z12, z13)
text(50, 15, paste("全体での個人レヴェルでのピアソン相関", round(cor(x_ind, z_ind), 3))) # 個人レヴェルでの相関
#### ここまでで一旦実行してグラフを確認 ####

# グループ平均のプロット
x_agg <- c(mean(x11), mean(x12), mean(x13))
z_agg <- c(mean(z11), mean(z12), mean(z13))
par(new=T);
plot(x_agg, z_agg, xlim=c(0,100), ylim=c(0,100), bty="n",
     main="", xlab="", ylab="", family="serif", las=1, pch=19, col="red", cex=2)
text(50, 10, # アグリゲイトレヴェルでの相関
     paste("グループレヴェルでのピアソン相関", round(cor(x_agg, z_agg), 3)), col="red")

発展1-2　シンプソンのパラドクス(Simpson’s paradox)

　本文の例について，それぞれクラメールのVを計算して示す。

A <- matrix(c(550, 450, 80, 20), ncol=2) # A大学
B <- matrix(c( 50, 150,120,180), ncol=2) # B大学

VA <- sqrt(chisq.test(A)$statistic/(min(dim(A)-1)*sum(A)))
VB <- sqrt(chisq.test(B)$statistic/(min(dim(B)-1)*sum(B)))
Z <- A + B # A大学とB大学のデータの合併
VZ <- sqrt(chisq.test(Z)$statistic/(min(dim(Z)-1)*sum(Z)))
names(VA) <- names(VB) <- names(VZ) <- "Cramer's V"

A; VA

     [,1] [,2]
[1,]  550   80
[2,]  450   20

Cramer's V 
 0.1420887 

B; VB

     [,1] [,2]
[1,]   50  120
[2,]  150  180

Cramer's V 
 0.1508172 

Z; VZ

     [,1] [,2]
[1,]  600  200
[2,]  600  200

Cramer's V 
         0 

発展1-3　スピアマンの順位相関係数(Spearman’s rank-order coefficient)

　2変数のうち一方でも，尺度水準が間隔尺度とは言えず順序尺度でしかない場合には，ピアソンの積率相関係数を計算する事が出来ず，ケンドールのτの様な順序連関係数やスピアマンの順位相関係数を計算するしかない。取り得る値の種類が少ないか同順位(tie)が多い場合には順序連関係数を用いれば良いと考えられるので，順位相関係数は，取り得る値の種類は多くて同順位になる事は少ないが，しかし間隔尺度とは見做せない場合に適用される事になる。社会調査データでは余りこうした変数は多くは無いだろう。よって以下では，間隔尺度以上の変数に対して，順位相関係数を用いる事の意義に注目する。 　本文のように，外れ値が無い場合と有る場合とで，積率相関と順位相関がどのように変わるかをシミュレイションするスクリプトを挙げておく。実行するたびに形が変わる。

x0 <- rnorm(n=100, mean=0, sd=1)
y0 <- rnorm(n=100, mean=0, sd=1)

r0 <- lm(y0 ~ x0)$residuals
z1 <- (rxy <- 0.2)*scale(x0) + sqrt(1-rxy^2)*scale(r0)

x3 <- x2 <- x0*10+40
z3 <- z2 <- z1*10+40
xc <- 50; zc <- 50
x3[x2 > xc & z2 > zc] <- x2[x2 > xc & z2 > zc] + 20
z3[x2 > xc & z2 > zc] <- z2[x2 > xc & z2 > zc] + 20

plot(x2[x2 <= xc | z2 <= zc],z2[x2 <= xc | z2 <= zc], 
     xlim=c(0,100), ylim=c(0,100), bty="n",
     main=sprintf("積率r=%.2f，順位ρ=%.2f．(n=%d)",
     cor(x2, z2), cor(x2, z2, method="spearman"), length(x2)),
     xlab="外れ値無し", ylab="", family="serif", las=1)
par(new=T)
plot(x2[x2 > xc & z2 > zc],z2[x2 > xc & z2 > zc], 
     xlim=c(0,100), ylim=c(0,100), bty="n", col="#f6aa00",
     main="", xlab="", ylab="", family="serif", las=1, pch=16)

plot(x3[x2 <= xc | z2 <= zc],z3[x2 <= xc | z2 <= zc], 
     xlim=c(0,100), ylim=c(0,100), bty="n",
     main=sprintf("積率r=%.2f，順位ρ=%.2f．(n=%d)",
     cor(x3, z3), cor(x3, z3, method="spearman"), length(x2)),
     xlab="外れ値有り", ylab="", family="serif", las=1)
par(new=T)
plot(x3[x2 > xc & z2 > zc],z3[x2 > xc & z2 > zc],
     xlim=c(0,100), ylim=c(0,100), bty="n", col="#f6aa00",
     main="", xlab="", ylab="", family="serif", las=1, pch=16)

第2章の【練習問題】の解答

1) の答え

# データフレイムの準備の復習
getwd() # Working Directoryを確認して，そこにデータのcsvファイルを置く

[1] "E:/科研費基盤A申請/社会統計学/Statistics"

data <- read.csv("practice.csv") # csvファイルをデータフレイムとして読み込む

# 散布図――散布図のスクリプトは一つずつ実行して比較しよう
plot(data$income, data$q1700)

# 幾つかオプションを付ける
plot(data$income, data$q1700,
     bty="n", main="本人年収と幸福度")

# 更にオプションを付ける
plot(data$income, data$q1700,
     bty="n", main="本人年収と幸福度",
     xlab="本人年収income（万円）",
     ylab="幸福度（低0←　→10高）")

# 相関係数の値と回帰直線を書き込む
r <- cor(data$income, data$q1700, use="complete")
plot(data$income, data$q1700,
     bty="n", main="本人年収と幸福度",
     xlab="本人年収income（万円）",
     ylab="幸福度（低0←　→10高）",
     sub=paste("ピアソン相関", round(r, 4)))
abline(lm(data$q1700 ~ data$income), lty=2)

　因みに，本人年収incomeを世帯年収fincomeに置き換えると，ピアソン相関は約0.30に上がる。自分で試してみよう。

2) の答え

# 最も単純な分割表
t01 <- table(data$sex, data$edu1)
t01

   
      0   1
  1  87  74
  2 164  51

# 分割表にラベルを付ける
names(dimnames(t01)) <- c("性別", "大卒／非大卒")
dimnames(t01)[[1]] <- c("男性", "女性")
dimnames(t01)[[2]] <- c("非大卒", "大卒")
t01

      大卒／非大卒
性別   非大卒 大卒
  男性     87   74
  女性    164   51

# 周辺度数付きの分割表
t02 <- addmargins(t01)
t02

      大卒／非大卒
性別   非大卒 大卒 Sum
  男性     87   74 161
  女性    164   51 215
  Sum     251  125 376

# 行％の表
prop.table(t01, margin=1)

      大卒／非大卒
性別      非大卒      大卒
  男性 0.5403727 0.4596273
  女性 0.7627907 0.2372093

# 列％の表
prop.table(t01, margin=2)

      大卒／非大卒
性別      非大卒      大卒
  男性 0.3466135 0.5920000
  女性 0.6533865 0.4080000

### ユールのQ
# 上で作成したスクリプトを再利用する為に，xにt01を行列にして代入
x <- matrix(t01, ncol=2)
Q <- det(x)/(x[1,1]*x[2,2]+x[1,2]*x[2,1]); Q

[1] -0.4645508

# 分割表を行列化しない場合は，下記のスクリプトを使用
x <- t01
Q <- (x[1,1]*x[2,2]-x[1,2]*x[2,1])/(x[1,1]*x[2,2]+x[1,2]*x[2,1])
Q

[1] -0.4645508

### φ係数
x <- matrix(t01, ncol=2)
phi <- det(x)/sqrt(sum(x[1,])*sum(x[2,])*sum(x[,1])*sum(x[,2]))
phi

[1] -0.2336201

# 行列化しない方法
x <- t01
phi <- (x[1,1]*x[2,2]-x[1,2]*x[2,1])/sqrt(sum(x[1,])*sum(x[2,])*sum(x[,1])*sum(x[,2]))
phi

[1] -0.2336201

### オッズ比
x <- t01
OR <- (x[1,1]+(d <- 0.5*(1 - ceiling(prod(x)/(prod(x)+.1)))))*(x[2,2]+d)/((x[1,2]+d)*(x[2,1]+d))
OR

[1] 0.3656065

3) の答え

table(data$sex, data$edu2)

   
     1  2  3
  1 56 31 74
  2 72 92 51

table(data$sex, data$edu2, useNA="ifany")

   
     1  2  3 <NA>
  1 56 31 74    4
  2 72 92 51    4

table(data$sex, data$edu2, useNA="always")

      
        1  2  3 <NA>
  1    56 31 74    4
  2    72 92 51    4
  <NA>  0  0  0    0

# 分割表にラベルを付けるには下記のように。
t10 <- table(data$sex, data$edu2, useNA="always")
names(dimnames(t10)) <- c("性別", "学歴") 
dimnames(t10)[[1]] <- c("男性", "女性", "NA")
dimnames(t10)[[2]] <- c("初等学歴", "中等学歴", "高等学歴", "NA")
t10

      学歴
性別   初等学歴 中等学歴 高等学歴 NA
  男性       56       31       74  4
  女性       72       92       51  4
  NA          0        0        0  0

4) の答え

x <- table(data$edu2, data$q0101)
c0 <- chisq.test(x)
x; c0

   
     1  2  3  4
  1 50 51 17 10
  2 28 50 30 14
  3 16 65 36  8

    Pearson's Chi-squared test

data:  x
X-squared = 29.632, df = 6, p-value = 4.617e-05

V <- sqrt(c0$statistic/(min(dim(x)-1)*sum(x)))
names(V) <- "Cramer's V" # 出力のラベル表示を修正
V

Cramer's V 
 0.1987694 

cor(data$edu2, data$q0101, method="kendall", use="complete")

[1] 0.177631

5) の答え

t21 <- table(data$edu2, data$job)
t22 <- table(data$edu2, data$job, data$sex)

# これでは見にくいので，ラベルを付ける。
names(dimnames(t21)) <- c("学歴3区分", "従業上の地位")
names(dimnames(t22)) <- c("学歴3区分", "従業上の地位", "性別")
dimnames(t21)[[1]] <- dimnames(t22)[[1]] <- c("1 初等", "2 中等", "3 高等")
dimnames(t21)[[2]] <- dimnames(t22)[[2]] <- c("1 常雇", "2 臨時", "3 自営", "4 無職")
dimnames(t22)[[3]] <- c("1 男性", "2 女性")
t21

         従業上の地位
学歴3区分 1 常雇 2 臨時 3 自営 4 無職
   1 初等     49     37     17     25
   2 中等     47     37     14     25
   3 高等     84     17      6     17

t22

, , 性別 = 1 男性

         従業上の地位
学歴3区分 1 常雇 2 臨時 3 自営 4 無職
   1 初等     37      6     10      3
   2 中等     25      1      4      1
   3 高等     64      2      4      3

, , 性別 = 2 女性

         従業上の地位
学歴3区分 1 常雇 2 臨時 3 自営 4 無職
   1 初等     12     31      7     22
   2 中等     22     36     10     24
   3 高等     20     15      2     14

# 3元分割表の「第1層」と「第2層」を足せば2元分割表に一致する。
t22[,,1] + t22[,,2]

         従業上の地位
学歴3区分 1 常雇 2 臨時 3 自営 4 無職
   1 初等     49     37     17     25
   2 中等     47     37     14     25
   3 高等     84     17      6     17

# 3元（3次元）分割表，3重クロス表の変数を入れ替える。
# 変数は順に，行変数，列変数，層変数となっている。
# 行変数と列変数の入れ替えは余り重要ではない。
# 以下の結果をよく見比べて解読しよう。

table(data$edu2, data$job, data$sex)

, ,  = 1

   
     1  2  3  4
  1 37  6 10  3
  2 25  1  4  1
  3 64  2  4  3

, ,  = 2

   
     1  2  3  4
  1 12 31  7 22
  2 22 36 10 24
  3 20 15  2 14

table(data$sex, data$job, data$edu2)

, ,  = 1

   
     1  2  3  4
  1 37  6 10  3
  2 12 31  7 22

, ,  = 2

   
     1  2  3  4
  1 25  1  4  1
  2 22 36 10 24

, ,  = 3

   
     1  2  3  4
  1 64  2  4  3
  2 20 15  2 14

table(data$sex, data$edu2, data$job)

, ,  = 1

   
     1  2  3
  1 37 25 64
  2 12 22 20

, ,  = 2

   
     1  2  3
  1  6  1  2
  2 31 36 15

, ,  = 3

   
     1  2  3
  1 10  4  4
  2  7 10  2

, ,  = 4

   
     1  2  3
  1  3  1  3
  2 22 24 14