N <- 1000 # 母集団サイズ
R <- 0.3 # 欲しい母相関係数の値
y <- rnorm(N, mean = 0, sd = 1)
x <- rnorm(N, mean = 0, sd = 1)
slope_xy <- cor(x, y, method = "pearson") * sd(y) / sd(x)
intercept_xy <- mean(y) - slope_xy * mean(x)
e <- y - (intercept_xy + slope_xy * x)
sd_x <- sqrt(sd(x)^2 * (N - 1)/N)
sd_e <- sqrt(sd(e)^2 * (N - 1)/N)
std_x <- (x - mean(x))/sd_x
std_e <- (e - mean(e))/sd_e
u <- R * std_x + sqrt(1 - R^2) * std_e
sd_u <- sqrt(sd(u)^2 * (N - 1)/N)
## サンプリング
n <- 50 # サンプルサイズの指定
times <- 100000 # 抽出回数
corrs <- vector(length = times) # 標本相関係数を格納するヴェクトル
for (i in 1:times) {
j <- sample(N, n)
corrs[i] <- cor(std_x[j], u[j], method = "pearson")
}『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』
第5章 2変数の関連の推定と検定
0 全体の章構成
1 積率相関係数についての推測統計
\[r_{xy}=\frac{s_{xy}}{s_{x}s_{y}}=\frac{1}{n}\sum_{i=1}^{n}(\frac{x_{i}-\bar{x}}{s_{x}})(\frac{y_{i}-\bar{y}}{s_{y}})\]
標本相関係数の標本抽出分布
hist(corrs, xlim = c(-1,1), col = "#99CC0080", border = "#99CC00",
main = sprintf("標本相関係数%d個のヒストグラム\n（母相関係数ρ=%.2f, 標本相関係数の平均%.3f）",
times, R, mean(corrs)),
sub = sprintf("母集団サイズ%d, 標本サイズ%.d", N, n), cex.main = 1.0,
xlab = "標本相関係数rの値", ylab = "標本の個数", family = "serif")
axis(side = 1, R, font = 2, family = "serif")
text(c(min(corrs), max(corrs)), 1000, round(c(min(corrs), max(corrs)),3), family = "serif")
1-1 積率相関係数のt検定(t-test)
積率相関係数の検定統計量 \[\frac{r}{\sqrt{1-r^2}}\sqrt{n-2}\;\sim\;t_{(n-2)}\]
標本サイズと標本相関係数から，ゼロ仮説「母相関係数=0」のt検定を行う（有意確率を表示する）スクリプトを示す。
片側検定と両側検定の両方を表示しているが，両側検定での有意確率は片側検定での有意確率の2倍である。
よって，両側検定で有意になる場合，同じ有意水準の片側検定では必ず有意になるが，片側検定で有意になったとしても，同じ有意水準の両側検定で有意になるとは限らない。
nやrに自由な値を設定して試してみよう。
rの絶対値はもちろん1を超えてはならない。
1-2 標本サイズと有意確率
本文の表のように，n，rの幾つかの組合せに対して，95%限界値，有意確率を計算してみよう。本文よりも，限界値に対応する標本相関係数の値など少し追加してある。
1-3 母相関係数の区間推定
サンプルサイズと標本相関係数が与えられた時に，信頼区間を計算する方法は，本文脚注にも書いてある通りである。
フィッシャーのZ変換を行ったものはそのままZと表記する事が多いが[南風原 2002: 116]，ここでは他の多くの標準正規変量と区別し易いようにあえてLと表記しておく。
\[L=\frac{1}{2}ln\frac{1+r}{1-r}\sim\;N(\mu_{L}, \sigma_{L}^2)=N\left(\frac{1}{2}ln\frac{1+\rho}{1-\rho}, \left(\frac{1}{\sqrt{n-3}}\right)^2\right)\] \[\frac{L-\mu_{L}}{\sigma_{L}}\sim\;N(0, 1^2)\] 95%信頼区間は， \[L-1.96\sigma_{L}<\mu_{L}<L+1.96\sigma_{L}\] 母相関係数\(\rho\)の信頼区間は，下記の関係式（逆関数）から計算する。
\[\rho=\frac{e^{2\mu_{L}}-1}{e^{2\mu_{L}}+1}\]
ここで，フィッシャーのZ変換を行ったものが実際に正規分布で近似できるかをシミュレイションで確かめてみよう【2024.10.19修正】。
上の1で使用した模擬母集団とそこからのランダムサンプリングの結果を利用する。
hist(corrs, xlim = c(-1, 1),
col = "#99CC0080", border = "#99CC00",
main = sprintf("標本相関係数%d個のヒストグラム\n（母相関係数ρ=%.2f, 標本相関係数の平均%.3f）",
times, R, mean(corrs)),
cex.main = 1.0,
sub = sprintf("母集団サイズ%d, 標本サイズ%.d", N, n),
xlab = "標本相関係数rの値", ylab = "標本の個数",
family = "serif")
axis(side = 1, R, font = 2, family = "serif")
text(c(min(corrs), max(corrs)), 1000,
round(c(min(corrs), max(corrs)),3), family = "serif")
L <- 1/2 * log((1 + corrs)/(1 - corrs)) # Lのシミュレイション値
muL <- 1/2 * log((1 + R)/(1 - R)) # Lの期待値
hist(L, xlim = c(muL - 6/sqrt(n - 3), muL + 6/sqrt(n - 3)),
breaks = seq(-6, 6, by = .05),
col = "#ffcc3399", border = "#ffcc33",
main = sprintf("フィッシャーZ変換後のL：%d個のヒストグラム\n（Lの期待値= %.2f, Lの平均 %.3f）",
times, muL, mean(L)),
cex.main = 1.0,
sub = sprintf("母集団サイズ%d, 標本サイズ%.d", N, n),
xlab = "Lの値", ylab = "標本の個数", family = "serif")
text(c(min(L), max(L)), 1000, round(c(min(L), max(L)),3),
family = "serif")
確かに，かなり正規分布に近付いている様子が分かる。
2 分割表の独立性についてのカイ二乗検定
カイ二乗統計量 \[\chi^2=\sum_{i=1}^{I}\sum_{j=1}^{J}\left(\frac{n_{ij}-f_{ij}}{\sqrt{f_{ij}}}\right)^2\;\sim\;\chi^2_{df=(I-1)(J-1)}\]
2-1 カイ二乗分布(chi-squared distribution)
カイ二乗分布の重ね描き
本文の，「自由度1,4,6,8,9,12のカイ二乗分布」のグラフを描くスクリプトは以下である。
x <- seq(0, end <- 20, by = 0.1)
df0 <- 1
# 最初に自由度1で基本のグラフを描く
plot(x, dchisq(x, df = df0), type = "l", ylim = c(0, .20),
bty = "n", las=1,
xlab = "カイ二乗変量", ylab = "確率密度",
main="自由度1,4,6,8,9,12のカイ二乗分布",
family="serif")
for (df1 in c(4, 6, 8, 9, 12)) {
par(new = T) # グラフを重ね描きする為のコマンド
plot(x, dchisq(x, df = df1), type = "l", ylim = c(0, .20),
bty = "n", axes = F, ylab = "", xlab = "",
lwd = 1 + df1/5,
lty = floor(df1/2))
}
棄却域付きカイ二乗分布グラフ
棄却域を塗り潰したカイ二乗分布のグラフは下記。自由度と有意水準は変更する事が出来る。
df0 <- 6 # 自由度を指定する
p0 <- .95 # 有意水準を指定する
x <- seq(0, end <- 20, by = 0.1) # endはx軸の右端
plot(x, dchisq(x, df = df0), type = "l", ylim = c(0, .20),
bty = "n", las = 1,
xlab = "カイ二乗変量", ylab = "確率密度",
main = paste("自由度", df0, "のカイ二乗分布の",
(1 - p0)*100, "%棄却域"),
family = "serif", cex.main = 1)
abline(h = 0, col = "gray")
rej <- seq(qchisq(p0, df0), end, by = .1) # 棄却域を設定
segments(rej, 0, rej, dchisq(rej, df0), col = "#CC330080")
segments(min(rej), 0, max(rej), 0, lwd = 3, col = "#CC3300")
axis(side = 1, at = round(min(rej),1), las = 2,
font = 2, family = "serif")
カイ二乗分布の定義： \[\chi_{(df=\nu)}^2=z_{1}^2+z_{2}^2+...+z_{\nu}^2=\left(\frac{x_{1}-\bar{x}_{1}}{s_{x_{1}}}\right)^2+\left(\frac{x_{2}-\bar{x}_{2}}{s_{x_{2}}}\right)^2+...+\left(\frac{x_{\nu}-\bar{x}_{\nu}}{s_{x_{\nu}}}\right)^2\]
2-2 カイ二乗検定の前提条件
ポワソン分布（計数分布）\(x_{ij}\;\sim\;Poisson(\lambda=f_{ij}),\;mean=\lambda,\;s.d.=\sqrt{\lambda}\)
\(\lambda\geq5\)ならば，\(x_{ij}\;\sim\;Poisson(\lambda=f_{ij})\;\sim\;N(\lambda,\lambda)\)
標準化すると， \[\frac{x_{ij}-f_{ij}}{\sqrt{f_{ij}}}\;\sim\;N(0.1^2)\]
分割表の独立性についてのカイ二乗検定において，検定統計量がカイ二乗分布で近似できる条件は，各セルの度数がポワソン分布に従い，そのポワソン分布が正規分布で近似できるということである。
本文中の二つのポワソン分布を重ね描きするスクリプト
x1 <- c(0:15)
lambda1 <- 2
lambda2 <- 5
plot(x1, dpois(x1, lambda1), type = "h", bty = "n",
lwd = 5, col = "#00000060", xlab = "", ylab = "",
las = 1, family = "serif", cex.main = 1.0, ylim = c(0, .3),
main = sprintf("平均%dと%dのポワソン分布", lambda1, lambda2))
par(new = T)
plot(x1, dpois(x1, lambda2), type = "b", bty = "n", axes = F,
lwd = 1, col = "#000000", xlab = "", ylab = "",
las = 1, family = "serif", ylim = c(0, .3))
少し改変したスクリプト
lambda <- c(3, 6) # 二つのポワソン分布の平均をヴェクトルで指定
x1 <- c(0: (max(lambda)*2.5))
plot(x1 - .1, dpois(x1, min(lambda)),
type = "h", bty = "n", xlim = range(x1),
lwd = 5, col = "#0000ff60",
xlab = "", ylab = "",
las = 1, family="serif",
cex.main = 1.0,
ylim = c(0, Y <- dpois(min(lambda), min(lambda))),
main = sprintf("平均%dと%dのポワソン分布",
min(lambda), max(lambda)))
par(new = T)
plot(x1 + .1, dpois(x1, max(lambda)),
type = "h", bty = "n", axes = F,
lwd = 5, col = "#ff000060",
xlab = "", ylab = "", xlim = range(x1),
las = 1, family = "serif", ylim = c(0, Y))
2-3 データに対するモデルの適合度検定(Goodness of Fit test)
社会学では，「カイ二乗検定」と云うと「独立性の検定」と同一視する傾向があるが，データに対して「2変数独立モデル」が適合するかどうかを検定していると考えれば「適合度検定」の一種である。
カイ二乗統計量は社会統計学のより発展的な場面で色々と登場してくる事になるが，いずれも「データとモデルのズレ」を表す統計量として統一的に考えるのが良い。
ここでは，chisq.test( )関数を使用した最も単純な適合度の検定を紹介しよう。
少人数の無作為標本に対して，「夫は外で働き，妻は家庭を守るべきである」という考え方について，そう思うかそう思わないか（性別役割分業意識）を，「そう思う／どちらかと言えばそう思う／どちらかと言えばそう思わない／そう思わない」の4件法で尋ねた結果，それぞれ5人，8人，15人，12人であったとする。
いずれのモデルも（そして他にも多くの可能なモデルが）データに適合する（適合しないとは言い切れない）と云う結果になる。
そのいずれかが正しいと証明される訳でも無く，カイ二乗統計量や有意確率の値からいずれかが「より正しい」と言えるかどうかにも議論がある。
適合度検定と考えられるものには，第6章の基礎1-1や第7章基礎1-5の等分散性の検定（各群の母分散が等しいモデルが適合するかどうか），第6章発展1の正規性の検定（変数の分布が正規分布に適合するかどうか）がある。
それ以外にも，第11章発展1-1の構造方程式モデリング(SEM; Structural Equation Modelling)では，設定したモデルの適合度がカイ二乗検定される（それ以外に，GFI[Goodness-of-Fit Index]を始めとして幾つかの適合度指標が算出されるが，GFIらについては統計的検定がされる訳では無い）。
第12章基礎1-4では，2項ロジスティック回帰分析に対する適合度検定としてHosmer-Lemeshow検定を紹介した。
第12章発展1のログリニア・モデル(log-linear model)も「モデルがデータに適合するか」を重視した分析手法であると言える。
適合度検定の多くは，統計的検定一般がそうである様に，サンプルサイズが大きいと有意になり易い。
大標本では「僅かなズレでも敏感に検出してしまう」のである。
適合度検定で統計的に有意になると云う事は，「モデルがデータに適合している」と仮定すると（ゼロ仮説），ズレの指標が偶然とは思えない程に大きくなるので，ゼロ仮説は棄却されると云う事であり，つまり「このモデルはデータには適合しない」と云う結論になるので，モデルが当てはまる事を期待していた分析者にとっては嬉しい結果ではない事も多いだろう。
発展1 分割表のさらなる検定
発展1-1 分割表の残差分析(analysis of residual)
上の2-1でデータフレイムdata01，その中の変数data01$sex，data01$q1601b，分割表table1などを既に準備済みであるので，それを利用して調整済み標準化残差（ハーバーマン残差）を表示したりモザイクプロットを表示したりしてみよう。
table1 # 分割表 機会平等なら格差やむをえない
性別 そう思う どち…思う どちらとも どち…思わない そう思わない
男性 53 57 32 12 8
女性 35 81 76 13 12ch # chisq.test(table1)の結果
Pearson's Chi-squared test
data: table1
X-squared = 19.041, df = 4, p-value = 0.0007715round(ch$stdres, 3) # 小数第三位までに四捨五入して表示 機会平等なら格差やむをえない
性別 そう思う どち…思う どちらとも どち…思わない そう思わない
男性 3.783 -0.429 -3.258 0.550 -0.255
女性 -3.783 0.429 3.258 -0.550 0.255mosaicplot(table1, type = "pearson", shade = T,
main = "モザイクプロット", las = 1)
発展1-2 イェーツの連続性修正(Yates’ correction for continuity)とフィッシャーの正確検定(Fisher’s exact test)
第5章の【練習問題】の解答
1) の答え
r <- .25 # 標本相関係数
n <- 30 # サンプルサイズ
p0 <- .95 # 信頼水準（信頼係数）
t0 <- r/sqrt(1-r^2)*sqrt(n-2)
t0 # t検定統計量[1] 1.36626qt((1-p0)/2, n-2); qt((1+p0)/2, n-2) # 100p%限界値[1] -2.048407[1] 2.048407pt(-1*t0, n-2)+(1-pt(t0, n-2)) # 両側検定での有意確率[1] 0.1827305有意確率は18.3%程度であり，有意水準を1%にしようが5%にしようが10%にしようが有意にはならない。つまりゼロ仮説は棄却出来ない。片側検定として考えるなら有意水準を半分にすればよいが，それでも5%では有意にならない。
2) の答え
r <- .25
n <- 70
p0 <- .95
t0 <- r/sqrt(1-r^2)*sqrt(n-2)
t0[1] 2.129163qt((1-p0)/2, n-2); qt((1+p0)/2, n-2)[1] -1.995469[1] 1.995469pt(-1*t0, n-2) + (1-pt(t0, n-2))[1] 0.03686488両側5%検定で有意になり，ゼロ仮説は棄却される。1)とは標本相関係数の値は同じであるが，標本サイズが2倍以上になっており，弱い関連に対してより敏感に検出出来る様になっている。
3) の答え
t01 <- matrix(c(49,47,84,37,37,17,17,14,6,25,25,17), ncol=4)
t01 [,1] [,2] [,3] [,4]
[1,] 49 37 17 25
[2,] 47 37 14 25
[3,] 84 17 6 17ch1 <- chisq.test(t01); ch1
Pearson's Chi-squared test
data: t01
X-squared = 30.377, df = 6, p-value = 3.333e-05高度に有意となり，母集団では独立であるとのゼロ仮説は棄却される。
独立性が否定されたので，行比率や列比率を計算して，どの様な偏りがあるのかを確認しておく。
addmargins(prop.table(t01, margin=1))[1:3,] Sum
0.3828125 0.2890625 0.1328125 0.1953125 1
0.3821138 0.3008130 0.1138211 0.2032520 1
0.6774194 0.1370968 0.0483871 0.1370968 1addmargins(prop.table(t01, margin=2))[,1:4]
0.2722222 0.4065934 0.4594595 0.3731343
0.2611111 0.4065934 0.3783784 0.3731343
0.4666667 0.1868132 0.1621622 0.2537313
Sum 1.0000000 1.0000000 1.0000000 1.00000004) の答え
カイ二乗検定のみならず正確検定も行っておこう。
t01s <- round(t01/3, 0); t01s [,1] [,2] [,3] [,4]
[1,] 16 12 6 8
[2,] 16 12 5 8
[3,] 28 6 2 6chisq.test(t01s)Warning in chisq.test(t01s): Chi-squared approximation may be incorrect
Pearson's Chi-squared test
data: t01s
X-squared = 9.5046, df = 6, p-value = 0.1471fisher.test(t01s)
Fisher's Exact Test for Count Data
data: t01s
p-value = 0.1437
alternative hypothesis: two.sidedaddmargins(prop.table(t01s, margin=1))[1:3,] Sum
0.3809524 0.2857143 0.14285714 0.1904762 1
0.3902439 0.2926829 0.12195122 0.1951220 1
0.6666667 0.1428571 0.04761905 0.1428571 1いずれにしても有意にはならず，独立モデルは棄却出来ない。3)と比較すると，分割表の行数・列数は等しく，したがって自由度も同じである。各セルの度数を一律に三分の一にしているので，行比率の表もほぼ同じである。
しかしカイ二乗統計量の値は，30.377と9.5046と3倍余り異なり，その結果検定結果は1%有意と10%でも有意でないと云う様に対照的である。サンプルサイズが小さいと，同じ様な関連が「偶然に過ぎないかも知れない」と見做されてしまうのである。
5) の答え
調整済み標準化残差（ハーバーマン残差）が標準正規分布に近似的に従う事を利用する。
調整済み標準化残差は，カイ二乗検定の結果を格納したオブジェクトのstdresで取り出せる。
ch1$stdres [,1] [,2] [,3] [,4]
[1,] -2.711816 1.508699 1.5961802 0.6057539
[2,] -2.650739 1.835024 0.6875128 0.8682866
[3,] 5.378449 -3.351823 -2.2948443 -1.4770222有意水準を両側5%とすると，絶対値が1.96を超えているセルが，有意に独立モデルから乖離しているセルであり，[3行, 1列]が期待値より有意に，かなり多い。[3, 2], [3, 3], [1, 1], [1, 2]は期待値より有意に少ない。つまり，3行1列のセルに著しく集中している事が分かる。
実はこの分割表は，演習用データの学歴3区分と従業上の地位4区分のものであり，高学歴者ほど常時雇用者が多いと云う結果を示しているのである。
addmargins(t01) Sum
49 37 17 25 128
47 37 14 25 123
84 17 6 17 124
Sum 180 91 37 67 375tbl <- table("学歴3区分"=data01$edu2, "従業上の地位"=data01$job)
dimnames(tbl)[[1]] <- c("1 初等", "2 中等", "3 高等")
dimnames(tbl)[[2]] <- c("1 常時雇用", "2 臨時雇用", "3 自営", "4 無職")
addmargins(tbl) 従業上の地位
学歴3区分 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 375round(prop.table(tbl, margin=1)*100, 1) 従業上の地位
学歴3区分 1 常時雇用 2 臨時雇用 3 自営 4 無職
1 初等 38.3 28.9 13.3 19.5
2 中等 38.2 30.1 11.4 20.3
3 高等 67.7 13.7 4.8 13.7ウェブ増補1 xtabs( )関数による分割表の作成
第2章と第5章を中心として，本書では分割表を作成する為に table( )関数を使用しているが，他によく使用される関数として，xtabs( )があるのでここで補っておく。
分割表の基本
まずは7ケース分のデータ・フレイムdata01を作成してから，table( )とxtabs( )で三元分割表を作成してみる。その後，ftable( )を利用して三元分割表を見易くする。
コマンドの文法と出力スタイルの両方を見比べてみよう。
data01 <- data.frame("sex"=c( 1, 1, 1, 1, 2, 2, 2),
"edu"=c(12,14,16,16,12,14,16),
"outcome"=c( 3, 4, 4, 5, 5, 4, 3))
data01 # 7ケース分のデータ・フレイム, NA無しtable(edu, outcome, sex, data=data01) # table( )ではこの書き方はNGError in eval(expr, envir, enclos): object 'edu' not foundtable(data01$edu, data01$outcome, data01$sex), , = 1
3 4 5
12 1 0 0
14 0 1 0
16 0 1 1
, , = 2
3 4 5
12 0 0 1
14 0 1 0
16 1 0 0# xtabs( )での書き方. table( )よりも若干親切な出力表示と，lm( )に似た構文.
xtabs(~ edu + outcome + sex, data=data01), , sex = 1
outcome
edu 3 4 5
12 1 0 0
14 0 1 0
16 0 1 1
, , sex = 2
outcome
edu 3 4 5
12 0 0 1
14 0 1 0
16 1 0 0# 三元分割表をftableで表示するには，変数の順序に注意
ftable(table(data01$edu, data01$outcome, data01$sex)) 1 2
12 3 1 0
4 0 0
5 0 1
14 3 0 0
4 1 1
5 0 0
16 3 0 1
4 1 0
5 1 0ftable(table(data01$sex, data01$edu, data01$outcome)) 3 4 5
1 12 1 0 0
14 0 1 0
16 0 1 1
2 12 0 0 1
14 0 1 0
16 1 0 0ftable(xtabs(~ edu + outcome + sex, data=data01)) sex 1 2
edu outcome
12 3 1 0
4 0 0
5 0 1
14 3 0 0
4 1 1
5 0 0
16 3 0 1
4 1 0
5 1 0 この限りでは，特に重回帰分析（一般線型モデル）などに慣れている人ほど，xtabs( )の方を好むだろう。出力の表記だけをとってもxtabs( )の方がやや親切である。
但しtable( )の場合は，以下の様に with( ) と組合せる事も出来る。
with(data01, table(edu, outcome, sex)), , sex = 1
outcome
edu 3 4 5
12 1 0 0
14 0 1 0
16 0 1 1
, , sex = 2
outcome
edu 3 4 5
12 0 0 1
14 0 1 0
16 1 0 0NAの処理
実際の調査データ分析では，欠損値（欠測missing）の処理が必須である。この点も違いがあるので比べておこう。
data02 <- data.frame("sex"=c( 1, 1, 1,NA, 2, 2, 2),
"edu"=c(12,14,NA,16,12,14,16),
"outcome"=c( 3, 4, 4, 5, 5,NA, 3))
data02 # 7ケース分のデータ・フレイム, NA有りtable(data02$edu, data02$outcome, data02$sex), , = 1
3 4 5
12 1 0 0
14 0 1 0
16 0 0 0
, , = 2
3 4 5
12 0 0 1
14 0 0 0
16 1 0 0xtabs(~ edu + outcome + sex, data=data02), , sex = 1
outcome
edu 3 4 5
12 1 0 0
14 0 1 0
16 0 0 0
, , sex = 2
outcome
edu 3 4 5
12 0 0 1
14 0 0 0
16 1 0 0table(data02$edu, data02$outcome, data02$sex, useNA="ifany"), , = 1
3 4 5 <NA>
12 1 0 0 0
14 0 1 0 0
16 0 0 0 0
<NA> 0 1 0 0
, , = 2
3 4 5 <NA>
12 0 0 1 0
14 0 0 0 1
16 1 0 0 0
<NA> 0 0 0 0
, , = NA
3 4 5 <NA>
12 0 0 0 0
14 0 0 0 0
16 0 0 1 0
<NA> 0 0 0 0xtabs(~ edu + outcome + sex, data=data02, na.action=na.pass, exclude=NULL), , sex = 1
outcome
edu 3 4 5 <NA>
12 1 0 0 0
14 0 1 0 0
16 0 0 0 0
<NA> 0 1 0 0
, , sex = 2
outcome
edu 3 4 5 <NA>
12 0 0 1 0
14 0 0 0 1
16 1 0 0 0
<NA> 0 0 0 0
, , sex = NA
outcome
edu 3 4 5 <NA>
12 0 0 0 0
14 0 0 0 0
16 0 0 1 0
<NA> 0 0 0 0ftable(table(data02$edu, data02$outcome, data02$sex, useNA="ifany")) 1 2 NA
12 3 1 0 0
4 0 0 0
5 0 1 0
NA 0 0 0
14 3 0 0 0
4 1 0 0
5 0 0 0
NA 0 1 0
16 3 0 1 0
4 0 0 0
5 0 0 1
NA 0 0 0
NA 3 0 0 0
4 1 0 0
5 0 0 0
NA 0 0 0ftable(xtabs(~ edu + outcome + sex, data=data02, na.action=na.pass, exclude=NULL)) sex 1 2 NA
edu outcome
12 3 1 0 0
4 0 0 0
5 0 1 0
NA 0 0 0
14 3 0 0 0
4 1 0 0
5 0 0 0
NA 0 1 0
16 3 0 1 0
4 0 0 0
5 0 0 1
NA 0 0 0
NA 3 0 0 0
4 1 0 0
5 0 0 0
NA 0 0 0NAの処理は，xtabs( )の方が何だか面倒である。
データ・フレイムへの（逆）変換
分割表からデータ・フレイムへの変換（逆変換？）を行う関数，as.data.frame.table( )があるが，これを使用する場合には，table( )の結果よりもxtabs( )の結果の方が便利である。上のdata01を使って実際に示してみるので，どう云う事か違いをよく見比べて読み取って欲しい。
t1 <- table(data01$edu, data01$outcome, data01$sex); t1, , = 1
3 4 5
12 1 0 0
14 0 1 0
16 0 1 1
, , = 2
3 4 5
12 0 0 1
14 0 1 0
16 1 0 0xt1 <- xtabs(~ edu + outcome + sex, data=data01); xt1, , sex = 1
outcome
edu 3 4 5
12 1 0 0
14 0 1 0
16 0 1 1
, , sex = 2
outcome
edu 3 4 5
12 0 0 1
14 0 1 0
16 1 0 0as.data.frame.table(t1) # データフレイムに逆変換as.data.frame.table(xt1) # データフレイムに逆変換ウェイト付け
上のas.data.frame.table( )関数の結果をよく見ると，出力されたデータ・フレイムにはFreqと云う変数が新たに加わっている。この例では，それぞれのパタンに該当するのはせいぜい1ケースしかなかったのでFreqは0か1の値しかないが，全く同じ変数値プロファイルを持つケースが複数あれば，Freqもその2以上のケース数を表す。
逆に，その様なFreq変数（変数名は別でも良い）を持つデータ・フレイムから，その度数によって重み付けを行った分割表を作成する事も xtabs( )によって可能になる。
# weightを付ける事が出来る。
data01$Freq <- c(1,2,3,4,5,6,7)
data01xtabs(~ edu + outcome + sex, data=data01), , sex = 1
outcome
edu 3 4 5
12 1 0 0
14 0 1 0
16 0 1 1
, , sex = 2
outcome
edu 3 4 5
12 0 0 1
14 0 1 0
16 1 0 0xtabs(Freq ~ edu + outcome + sex, data=data01), , sex = 1
outcome
edu 3 4 5
12 1 0 0
14 0 2 0
16 0 3 4
, , sex = 2
outcome
edu 3 4 5
12 0 0 5
14 0 6 0
16 7 0 0これを利用すると，単に度数で重みを付ける（頻度ウェイティング）だけでなく，サンプルの歪みの補正にしばしば使われる重み付け（確率ウェイティング）を行う事も出来る。原理的には，重回帰分析などのウェイト付けと同じである。ウェイト付けを行わない場合と行った場合では，そもそも分布が違ったものになるので，分析結果も異なり得る。
# 単純に度数ではなく，ケースへのweight付けを同様に行える。
data01$weight <- c(1.6, .4, 1.8 , 0.2, .4, .6, 2.0)
data01xtabs(weight ~ edu + outcome + sex, data=data01), , sex = 1
outcome
edu 3 4 5
12 1.6 0.0 0.0
14 0.0 0.4 0.0
16 0.0 1.8 0.2
, , sex = 2
outcome
edu 3 4 5
12 0.0 0.0 0.4
14 0.0 0.6 0.0
16 2.0 0.0 0.0# 線型モデル lm( )におけるweight付け
data01$fsex <- factor(data01$sex) # 性別は要因型（名義尺度）に
data01$cedu <- data01$edu - mean(data01$edu) # 教育年数は中心化
lm(outcome ~ fsex * cedu, data=data01) # weight付け無し
Call:
lm(formula = outcome ~ fsex * cedu, data = data01)
Coefficients:
(Intercept) fsex2 cedu fsex2:cedu
3.92208 -0.06494 0.36364 -0.86364 lm(outcome ~ fsex * cedu, data=data01, weights=weight) # weight付け有り
Call:
lm(formula = outcome ~ fsex * cedu, data = data01, weights = weight)
Coefficients:
(Intercept) fsex2 cedu fsex2:cedu
3.6734 0.1838 0.2725 -0.7725