Google
WWWを検索
サイト内を検索

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.namescol.namesTRUEにして読み込むと、そのままクロス表になります。

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\)列の頻度の合計を乗じたものを、全体の頻度の合計で二回割ったものが期待確率です。期待頻度は、それにさらに全体の頻度をかけたものになります。

観測値と期待頻度の差は二項分布に、近似的に正規分布に従います。よって標準化された差の二乗は、χ二乗分布に従うと看做すことができます。よってχ二乗検定を用いて、独立性の検定を行うことができます1。また、各セルの差も正規分布を用いて外れ値ではないか検定ができます。

(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)と大きなサンプリングサイズのシミュレーションで値を計算しても、話が変わるようなことはほとんどないと思います。

2.1 検定統計量のプロット

ほとんど利用されていない気がしますが、χ二乗検定の検定統計量の各セルごとの値の符号付平方根がプロットができます(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}}\)です。

2.2 クロス表のプロット

プロットでは無く表を画くことがほとんどですが、ヒートマップで視覚化することもできます。ただし、標準の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で同様のことができるようです。

3 信頼区間の計算

医療統計の教科書どおりに計算すれば出てきますが、コードを確認しておきましょう。

3.1 オッズ比

医療統計でよくみるオッズ比と信頼区間の計算をしましょう。

# 2×2表をつくる
M <- matrix(c(283, 145, 169, 130), 2, dimnames = list(c("men", "women"), c("negative", "positive")))
# オッズを計算
odds <- M[,2]/M[,1]
# オッズ比を計算
odds_ratio <- odds[2]/odds[1]
# 対数オッズ比の標準誤差を計算
se_log_odds <- sqrt(sum(1/M))
# 対数オッズ比の100*1-a%信頼区間を計算
a <- 0.05
interval <- c(qnorm(a/2), 0, qnorm(1-a/2))
names(interval) <- c("lower", "estimate", "upper")
log_odds_confint <- se_log_odds*interval + log(odds_ratio)
# オッズ比の信頼区間を計算
(odds_confint <- exp(log_odds_confint))
   lower estimate    upper 
1.107895 1.501326 2.034470 

対数オッズ比の分散はクロス表の4つのセルの観測数のそれぞれの逆数の和で近似することができるので、そのように計算しています。

なお、実用ではepitoolsパッケージなどを使う方がよいでしょう。2行以上の行列も処理してくれますし、独立性検定なども複数の手法で同時にやってくれます。

library(epitools)
(r_odds <- oddsratio.wald(M))$measure
                        NA
odds ratio with 95% C.I. estimate    lower   upper
                   men   1.000000       NA      NA
                   women 1.501326 1.107895 2.03447

3.2 リスク比

オッズとリスクの計算の分母が異なることに注意すれば、リスク比も同様に計算できます。

risk <- M[,2]/(M[,1] + M[,2])
risk_ratio <- risk[2]/risk[1]
se_log_risk <- sqrt(sum(1/M[1,2] + 1/M[2,2] - 1/(M[1,1] + M[1,2]) - 1/(M[2,1] + M[2,2])))
log_risk_confint <- se_log_risk*interval + log(risk_ratio)
names(log_risk_confint) <- c("lower", "estimate", "upper")
(risk_confint <- exp(log_risk_confint))
   lower estimate    upper 
1.063839 1.264336 1.502619 

これも実用では、epitoolsパッケージなどを使う方がよいでしょう。

library(epitools)
(r_risk <- riskratio.wald(M))$measure
                        NA
risk ratio with 95% C.I. estimate    lower    upper
                   men   1.000000       NA       NA
                   women 1.264336 1.063839 1.502619

3.3 絶対リスク差

あまり使われていないみたいですが、割り算が入らない分だけ容易に計算できます。

rd <- c(risk[2] - risk[1], use.names = FALSE)
se_rd <- sqrt(M[1,1]*M[1,2]/(M[1,1]+M[1,2])^3 + M[2,1]*M[2,2]/(M[2,1]+M[2,2])^3)
(rd_confint <- se_rd*interval + rd)
     lower   estimate      upper 
0.02486461 0.09883347 0.17280232 

これは実用では、ビルトインのprop.testを用いるのがよいでしょう。

prop.test(M, correct = FALSE)

    2-sample test for equality of proportions without continuity correction

data:  M
X-squared = 6.8976, df = 1, p-value = 0.008631
alternative hypothesis: two.sided
95 percent confidence interval:
 0.02486461 0.17280232
sample estimates:
   prop 1    prop 2 
0.6261062 0.5272727 

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