1 ファイルの読み込みと集計

回帰分析に使うデータセットと異なり変化があるので、それぞれの場合の読み込み方を説明します。

1.1 未集計データの場合

df01 <- read.table("https://wh.anlyznews.com/R/dataset/for_xtt.csv", 
   header = TRUE, sep = ",", 
   colClasses = c("character", "factor", "factor"))

2列目と3列名が因子型になるように指定しています。なお、colClassesを指定する代わりに、stringsAsFactors = TRUE, as.is = "name"を入れても同じ結果になります。

head(df01)

     name d1 d2
1    Noah  a  x
2 William  b  z
3  Alfred  b  z
4    Carl  c  z
5   Aksel  b  z
6    Emil  c  z

中身は見ての通り、3種類の文字列のデータセットです。

summary(df01)

     name           d1     d2    
 Length:100         a:32   x:14  
 Class :character   b:35   y:15  
 Mode  :character   c:33   z:71

因子型にしておくと、summaryで要素ごとの頻度が出ます。因子型の2列目と3列目が、今回の分析の対象となります。

ペンギンで覚えるRの集計関数で練習したように、集計してクロス表を作成します。

(xt <- xtabs(~ d1 + d2, df01))

   d2
d1   x  y  z
  a 10  0 22
  b  2  2 31
  c  2 13 18

1.2 集計済データの場合

集計済データで保存されている場合もあります。

df02 <- read.table("https://wh.anlyznews.com/R/dataset/for_xtt_summarized.csv", header = TRUE, sep = ",")
head(df02)

  d1 d2 Freq
1  a  x   10
2  b  x    2
3  c  x    2
4  a  y    0
5  b  y    4
6  c  y   10

xtabsのモデル式の左辺を指定すれば、集計済データを使ったクロス表に変換できます。

xt <- xtabs(Freq ~ d1 + d2, df02)

1.3 保存データが既にクロス表の場合

ファイルにクロス表が書かれている場合もあります。その場合はrow.namesとcol.namesをTRUEにして読み込むと、そのままクロス表になります。

xt <- as.matrix(read.table("https://wh.anlyznews.com/R/dataset/xtt.csv", header = TRUE, row.names = 1, sep = ","))

1行目の要素数が、2行目以降の要素数よりひとつ少ない形式であれば、row.names = 1は無くても読み込めます。

2 独立性の検定

十分な統制ではなく他に計量分析をかけることになるので計算しないことも多いと思いますが、表全体とセルごとの統計的仮説検定をかけるのが教科書的な方法です。

どの行かで表される因子とどの列かで表される因子が独立だとします。\(i\)行かつ\(j\)列である期待確率は、\(i\)行になる確率と\(i\)列になる確率を乗じたものになります。つまり、\(i\)行の頻度の合計と\(j\)列の頻度の合計を乗じたものを、全体の頻度の合計で二回割ったものが期待確率です。期待頻度は、それにさらに全体の頻度をかけたものになります。

観測値と期待頻度の差は二項分布に、近似的に正規分布に従います。よって差の二乗は、χ二乗分布に従うと看做すことができます。よってχ二乗検定を用いて、独立性の検定を行うことができます¹。

(xst <- chisq.test(xt, correct = FALSE)) # 2×2表でもイェーツの連続性補正は行わない

Warning in chisq.test(xt, correct = FALSE): カイ自乗近似は不正確かもしれません


    Pearson's Chi-squared test

data:  xt
X-squared = 24.373, df = 4, p-value = 6.722e-05

xt_p_value <- 2*pnorm(abs(xst$stdres), lower.tail = FALSE)
(xt_adj_p_value <- matrix(p.adjust(xt_p_value), nrow(xt_p_value), dimnames = dimnames(xt_p_value))) # 多重比較補正/Holm法

           x            y         z
a 0.01647524 0.0214869316 0.9624113
b 0.23439541 0.9624112772 0.2343954
c 0.65342761 0.0008567694 0.2343954

警告で近似が不正確かも知れないと脅されますが、chisq.test(xt, simulate.p.value = TRUE, B = 1000000)と大きなサンプリングサイズのシミュレーションで値を計算しても、話が変わるようなことはほとんどないと思います。

3 検定統計量のプロット

ほとんど利用されていない気がしますが、χ二乗検定の検定統計量の各セルごとの値の符号付平方根がプロットができます（Cohen-Friendly association plot）。

assocplot(xt)

行\(i\)列\(j\)のセルの観測値を\(f_{ij}\)、期待値を\(e_{ij}\)置いて、高さが\((f_{ij} - e_{ij}) / \sqrt{e_{ij}}\)、幅が\(\sqrt{e_{ij}}\)です。χ二乗検定の検定統計量は、\(\sum_{ij} (f_{ij} - e_{ij}) / \sqrt{e_{ij}}\)です。

4 クロス表のプロット

プロットでは無く表を画くことがほとんどですが、ヒートマップで視覚化することもできます。ただし、標準のheatmapはセルの値が表示できず、アスタリスクも入れられません。今回はimageを使って描画します。

plotXTT <- function(xt, xt.p.value = 1+1e06, pt.cex = 1.25){
   M <- xt
   for(i in 1:nrow(xt)){
      M[i, ] <- xt[nrow(xt) - i + 1, ]
   }
   rownames(M) <- rev(rownames(M))
   image(t(M), axes = FALSE, 
      col = colorRampPalette(c("white", "pink"))(max(M) + 1))
   at_x <- seq(0, 1, length.out = ncol(M))
   at_y <- seq(0, 1, length.out = nrow(M))
   axis(1, at = at_x, labels = colnames(M), lwd = 0, lwd.ticks = 1)
   axis(2, at = at_y, labels = rownames(M), lwd = 0, lwd.ticks = 1)
   A <- matrix(ifelse(xt.p.value <= 0.01, "^'***'", 
      ifelse(xt.p.value <= 0.05, "^'**'", 
        ifelse(xt.p.value <= 0.10, "^'*'", "")
      )
   ), nrow(M), ncol(M))
   for(i in 1:nrow(M)){
      for(j in 1:ncol(M)){
         text(at_x[j], at_y[i], 
            parse(text = sprintf("%.0f%s", M[i, j], A[i, j])), cex = pt.cex)
      }
   }
}

par(mar = c(7, 4, 1, 1))
plotXTT(xt, xt_adj_p_value)
mtext("Note) ***,  ** and * are statistically significant at the 1%, 5% and 10% level, respectively.", side = 1, line = par()$mgp[1] + 1, adj = 1, xpd = TRUE)
expr <- substitute(paste(chi^2, " statistic: ", s, ", DF: ", df, ", p-value: ", p), with(xst, list(s = statistic, df = parameter[["df"]], p = p.value)))
mtext(expr, side = 1, line = par()$mgp[1] + 2.5, adj = 1, xpd = TRUE)

ComplexHeatmap::Heatmapで同様のことができるようです。

正規分布を仮定しないFisherの正確検定をかける場合は、fisher.test(xt)です。↩︎