度数とクロス集計
カテゴリカル変数のクロス集計表 (分割表) 、独立性の検定、連関、グラフ化について説明する。
度数表をつくる
度数表 (frequency table, 分割表/contingency table 、クロス集計表/cross-tab ともいう) を作成する方法はいくつかある。以下にその例を示す。また、青木先生のfrequency 関数も便利。
table 関数
table( ) 関数で度数表が作成できる。割合を算出したいときはe prop.table( ) 関数, 合計欄には margin.table( ) 関数が使える。
# 2元の度数表
mytable <- table(gear=mtcars$gear, carb=mtcars$carb) # "gear=..."としないと名前のない表になるので見づらい。
mytable
# 集計欄
margin.table(mytable, 1) # 行の合計
margin.table(mytable, 2) # 列の合計)
# 表に集計行を追加
addmargins(mytable)
# パーセンテージ
prop.table(mytable) # セルのパーセンテージ
prop.table(mytable, 1) # 行のパーセンテージ
prop.table(mytable, 2) # 列のパーセンテージ
# パーセンテージ (+ 合計欄) の関数
cprp <- function(x, p=1) {
xn <- addmargins(x)
cp <- t(t(xn)/xn[nrow(xn),]*p)[-nrow(xn), -ncol(xn)]
rp <- (xn/xn[, ncol(xn)]*p)[-nrow(xn), -ncol(xn)]
pt <- prop.table(x)
return(list(adm=xn, cp=cp, rp=rp, pt=pt))
}
cprp(mytable)
table( ) 関数は3元以上の集計表も作成できる。ただし、この場合は ftable( ) 関数を使って"flat"な表を作ったほうが見やすい。
# 3元の集計表
attach(mtcars)
mytable <- table(am, gear, carb)
mytable
ftable(mytable)
table 関数は欠損値を除いて集計する。 NAを含めたいときはtable関数の中でexclude=NULLと指定する。
x <- c(1,1,2,2,NA)
fx <- factor(c("a", "a", "b", "b", NA))
x
fx
table(x) # NAは表示されない。
table(fx)
table(x, exclude=NULL) # NAも集計される
table(fx, exclude=NULL)
nfx <- factor(fx, exclude=NULL)
levels(nfx) # 新しい因子変数を作るとNAが水準として認識される。
度数が0の項目を集計する（集計する項目を特定する）には、levelsオプションを指定してfactor変数にする。
x <- factor(1:2, levels=1:3)
y <- factor(c(1,3), levels=1:3)
xy <- data.frame(x,y)
table(xy)
# levelsを指定しない場合度数0がカウントされない。
x <- factor(1:2)
y <- factor(c(1,3))
xy <- data.frame(x,y)
table(xy)
xtabs関数
xtabs( ) 関数は数式形式で記述し、クロス集計表を作成する。
# 3元の集計表
mytable <- xtabs(~gear+carb+am, data=mtcars)
mytable # amで分割し、gearを行に、carbを列にする
ftable(mytable) # 見やすく表示
summary(mytable) # 独立性のカイ二乗検定
数式左辺に数値変数を記述すると、その変数は度数を表すものとして解釈され、集計が行われる。集計表をさらに計算するときなどに使う。
# 例示用データ
cls <- sample(c("red", "blue", "yellow"), 20, replace=T)
freq <- sample(1:5, 20, replace=T)
data.frame(cls, freq)
# これだと"red", "blue", "yellow" と "1", "2", .. "5" のクロス集計になる
xtabs(~freq+cls)
# こうするとfreqを度数として解釈し、合計を出す。
xtabs(freq~cls)
aggregate(freq, list(cls), sum) # つまりこれと同じ。
Crosstable 関数
gmodels パッケージのCrossTable( ) 関数はいろいろオプションが選べて便利。出力が昔のSASっぽい。
# 2元のクロス集計
library(gmodels)
CrossTable(mtcars$gear, mtcars$carb)
# 中身を見てみる
CTres <- CrossTable(mtcars$gear, mtcars$carb)
CTres # リストになってる
cprp(table(mtcars$gear, mtcars$carb)) # これと同じだった。やはり自分が欲しい関数はどっかでだれかが作っているのだ。
オプションで行列パーセントの表示を選べる。また、カイ二乗検定、フィッシャーの正確確率検定、マクネマー検定を実行できる。残差も出力できる。欠損値をケースに含めることもできる。詳しくはhelp(CrossTable) 。
独立性の検定
カイ二乗検定
2元の分割表はchisq.test( ) 関数で独立性の検定を行うことができる。デフォルトのp値は近似カイ二乗分布に基づくものである。オプションでモンテカルロシミュレーションでのp値を求めることができる。
# 心理統計学の基礎 p186
x <- c(12,10,4,14)
tx <- as.table(matrix(x, nr=2, dimnames=list(A=c("aru", "nai"), B=c("aru", "nai"))))
chires <- chisq.test(tx, correct=F) # 教科書と同じ結果、補正なし
chires$expected # 期待度数
str(chires) # その他の計算されるもの
chisq.test(tx) # デフォルト
chisq.test(tx, simulate.p.value=T) # モンテカルロシミュレーションでのp値
Fisher の直接確率検定
fisher.test(x) 関数で行う。x には2元のクロス集計表をオブジェクトとしてあたえる。
fisher.test(tx)
# spssの実行例。UCLAのサイトから
library(foreign)
dat <- read.spss("http://www.ats.ucla.edu/stat/spss/webbooks/reg/hsb2.sav", to.data.frame=T)
names(dat) <- tolower(names(dat))
tb <- table(dat$schtyp, dat$female)
tb
fisher.test(tb)
Mantel-Haenszel 検定
mantelhaen.test(x) 関数を用いる。xには3元のクロス集計表を与える。このうち、2つの要因の独立性を検定する。
# spssの実行例と一致。 uclaのサイトから。
tb <-array(c(
8,21, 2,35,
11,10, 2,13,
4, 7, 1,22,
19, 7, 2, 4),
dim = c(2,2,4),
dimnames = list(cat=1:0, iccpr=1:0, area=1:4))
tb
mantelhaen.test(tb)
対数線型モデル
MASSパッケージの loglm( )関数で実行できる。以下の例は立命館大学のspssのものから。結果は一致する。
# データの準備。3元のクロス集計表
tb <-array(c(
42,23,79,65,32,17,
4,11,12,41,8,24),
dim = c(2,3,2),
dimnames = list(
fault=c("low", "high"),
moral=c("high", "neutral", "low"),
verdict=c("guilty", "not guilty")))
tb
# 対数線形モデルの実行
library(MASS)
# 交互作用をすべて含むモデル
fit1 <- loglm(~fault+moral+verdict+fault*moral+fault*verdict+moral*verdict, data=tb)
fit1
# fault*verdictを除外したモデル
fit2 <- loglm(~fault+moral+verdict+fault*moral+moral*verdict, data=tb)
fit2
# 両モデルの比較
anova(fit1, fit2)
モザイクプロットの例および対数線形モデルの解説：Visualizing Loglinear Models
glmを使った例：analyzing loglinear models
連関
vcdパッケージのassocstats(tb) 関数を使うとファイ係数 (Phi-Coefficient) , 一致係数 (Contingency Coeff) 、クラメールの連関係数 (Cramer's V) がわかる。参考：
また、kappa(tb)関数でCohenのカッパ係数、重みづきカッパ係数がわかる。
参考：
名義尺度間の相関
Measures of Association in Crosstab Tables
# 心理統計学の基礎 p186
x <- c(12,10,4,14)
tb <- as.table(matrix(x, nr=2, dimnames=list(A=c("aru", "nai"), B=c("aru", "nai"))))
library(vcd)
assocstats(tb)
kappa(tb)
結果の視覚化
1次元のデータであれば 棒グラフや and 円グラフを使う。
複数のカテゴリカルデータの場合は vcd パッケージが利用できる (モザイクや連関プロット) 。
対応分析を行うca パッケージも使える。
集計表の変換
データの準備。3元のクロス集計表
tb <-array(c(
42,23,79,65,32,17,
4,11,12,41,8,24),
dim = c(2,3,2),
dimnames = list(
fault=c("low", "high"),
moral=c("high", "neutral", "low"),
verdict=c("guilty", "not guilty")))
tb
# フラットな表にする
ftable(tb)
data.frame(ftable(tb))
data.frame(tb) # 見づらい
# 元のデータを生成する。
x <- data.frame(ftable(tb))
x2 <- data.frame( lapply (x, function(i) rep(i, x[,"Freq"]))[-4])
x2
table(x2) # tbと同じ