🐯

Rによる計数データ（カテゴリカルデータ）解析

2022/03/04に公開

はじめに

本記事では、医療分野における計数データの基本的な解析のRでの実装を紹介する。計数データは人数やイベント回数のような整数値をとるデータで、カウントデータやカテゴリカルデータと呼ばれることもある。本記事では、2項分布の解析、2標本におけるオッズ比、リスク比、リスク差の推定、また多群での解析など幅広く扱う。

生物統計の初学者よりは一定の知識がある者を対象とし数式やRの使用方法に関する記述を省いたためわかりづらい箇所があるかもしれないが、ご容赦いただければ幸いである。

【参考文献】

医学統計学ハンドブック

Package
仮想的なランダム化比較試験のデータ生成
データの要約
- カテゴリ化した変数の要約
二項分布の解析
- Wald法
- Score法
- Exact法
2*2表の解析
- 独立性の検定
- 粗解析によるリスク比、オッズ比、リスク差
- 層別解析による共通リスク比、オッズ比、リスク差（Mantel-Haenzel推定量）
- 回帰分析によるリスク比、オッズ比、リスク差
- 対応のある検定（McNemar test）
- （補足）McNemar testのI*J表への拡張
I*J表の解析
- 2*J 表
- I*2 表
- I*J 表（名義尺度*名義尺度）
- I*J 表（順序尺度*順序尺度）
サンプルサイズ計算
- 単群のサンプルサイズ計算
- 対応のない2標本のサンプルサイズ計算
- McNemar検定のサンプルサイズ計算
まとめ
- （メモ）今後、追加する可能性のある解析

Package

利用するPackageは以下の通りである。

#package
library(tidyverse)
library(gtsummary)
library(epiR)
library(geepack)
library(broom)
library(coin)
library(TrialSize)

仮想的なランダム化比較試験のデータ生成

本記事では、仮想的なランダム化比較試験（RCT）のデータを生成し解析した。仮想的なRCTとして、200名の2型糖尿病患者もしくは予備群を対象にある薬剤の効果を検討する試験を考える。割付群としては、Placebo群、2.0 mg, 5.0mg, 10.0mg群の4群に均等に割り付け、アウトカムはHbA1cの変化量とした。データ生成時には、アウトカムは連続量として与えているが、本記事内では計数データの解析を扱うためカテゴリ化した変数のみを解析対象とする。

#データ生成
set.seed(1234)
n <- 200
th1 <- 
  data.frame(ID = c(1:n), 
             pre_drug = sample(c(0,1), n, replace=TRUE, prob = c(.6, .4)) ) %>%
  mutate(TRT = case_when(
    ID <= n/4 ~ "placebo",
    ID <= 2*n/4 ~ "2.0 mg",
    ID <= 3*n/4 ~ "5.0 mg",
    ID <= n ~ "10.0 mg")) %>% 
  mutate(TRT_c = case_when(
    ID <= n/4 ~ "0",
    ID <= 2*n/4 ~ "1",
    ID <= 3*n/4 ~ "2",
    ID <= n ~ "3")) %>% 
  mutate(TRT_c2 = ifelse(ID <= n/4,"0","1")) %>% 
  mutate(dose = case_when(
    ID <= n/4 ~ 0,
    ID <= 2*n/4 ~ 2,
    ID <= 3*n/4 ~ 5,
    ID <= n ~ 10)) %>% 
  mutate(pre = round(runif(n=n, min=6.5, max=7.4),2)) %>% 
  mutate(pre = ifelse(pre_drug == 1, pre-0.2, pre)) %>% 
  mutate(post = case_when(
    ID <= n/4 ~ round(pre + rnorm(n=n/4, mean=0.2, sd=0.4),2),
    ID <= 2*n/4 ~ round(pre + rnorm(n=n/4, mean=0, sd=0.4),2),
    ID <= 3*n/4 ~ round(pre + rnorm(n=n/4, mean=-0.2, sd=0.4),2),
    ID <= n ~ round(pre + rnorm(n=n/4, mean=-0.4, sd=0.4),2))) %>% 
  mutate(post = ifelse(pre_drug == 1, post-0.1, post)) %>% 
  mutate(change = post - pre)

print(head(th1))

変数名	概要
ID	ID
pre_drug	試験前の薬の服用有無
TRT	割付群
TRT_c	Placebo群、2.0 mg, 5.0mg, 10.0mgをそれぞれ(0, 1, 2, 3)で表記
TRT_c2	Placebo群（0）と治療群（1）の2値
pre	ベースラインでのHbA1c [%]
post	試験終了時のHbA1c [%]
change	post - pre

##   ID pre_drug     TRT TRT_c TRT_c2 dose  pre post change
## 1  1        0 placebo     0      0    0 7.09 7.48   0.39
## 2  2        1 placebo     0      0    0 6.78 7.16   0.38
## 3  3        1 placebo     0      0    0 6.59 6.76   0.17
## 4  4        1 placebo     0      0    0 6.99 7.37   0.38
## 5  5        1 placebo     0      0    0 6.77 6.99   0.22
## 6  6        1 placebo     0      0    0 6.96 7.36   0.40

th1のpre, post, changeをもとにカテゴリ化した変数をいくつか作成した。

th2 <- th1 %>% 
  mutate(pre2 = ifelse(pre < 7.0,1,0)) %>% 
  mutate(post2 = ifelse(post < 7.0,1,0)) %>% 
  mutate(change2 = ifelse(change < 0,1,0)) %>% 
  mutate(change3 = case_when(
    change < -0.5 ~ 2,
    change < 0 ~ 1,
    change < 0.5 ~ -1,
    TRUE ~ -2 ))

print(head(th2))

変数名	概要
pre2	preのHbA1cが7.0%未満なら1、7.0%以上なら0の指示変数
post2	postのHbA1cが7.0%未満なら1、7.0%以上なら0の指示変数
change2	HbA1cが改善していれば1，していなければ0の指示変数
change3	HbA1cが0.5%以上減少していれば2，改善で1, 増加で-1, 0.5%以上の増加で-2 をとる変数

##   ID pre_drug     TRT TRT_c TRT_c2 dose  pre post change pre2 post2 change2
## 1  1        0 placebo     0      0    0 7.09 7.48   0.39    0     0       0
## 2  2        1 placebo     0      0    0 6.78 7.16   0.38    1     0       0
## 3  3        1 placebo     0      0    0 6.59 6.76   0.17    1     1       0
## 4  4        1 placebo     0      0    0 6.99 7.37   0.38    1     0       0
## 5  5        1 placebo     0      0    0 6.77 6.99   0.22    1     1       0
## 6  6        1 placebo     0      0    0 6.96 7.36   0.40    1     0       0
##   change3
## 1      -1
## 2      -1
## 3      -1
## 4      -1
## 5      -1
## 6      -1

データの要約

各群ごとに、カテゴリカル変数については人数を集計し割合を算出し、連続量に対しては平均（標準偏差）、中央値（IQR）、範囲を算出した。

#データの要約
table1 <- th1 %>% 
  dplyr::select(pre_drug,pre,post,change,TRT_c) %>%
  tbl_summary(
    by = TRT_c,
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c("{mean} ({sd})","{median} ({p25}, {p75})", "{min}, {max}"),
    missing = "no") %>% 
  add_n() %>% 
  modify_header(label = "**Variable**")

カテゴリ化した変数の要約

th2で新たに定義したカテゴリ化した変数について、いくつかの分類のもと表を作成した。

TRT_c(Placebo群、2.0 mg, 5.0mg, 10.0mg群)とchange3(4 levels)の4*4のクロス表を作成した。

table2 <- th2 %>% 
  tbl_cross(row = TRT_c, col = change3, percent = "row")

TRT_c(Placebo群、2.0 mg, 5.0mg, 10.0mg群)とchange2(2 levels)の4*2のクロス表を作成した。

table3 <- th2 %>% 
  tbl_cross(row = TRT_c, col = change2, percent = "row")

TRT_c2(Placebo群、投薬群)とchange2(2 levels)の2*2のクロス表を作成した。

table4 <- th2 %>% 
  tbl_cross(row = TRT_c2, col = change2, percent = "row")

適宜、紹介する手法に合わせたクロス表を利用する。

二項分布の解析

二項分布の解析を行った。単群の計数データに対する解析になる。そのため、投薬群のみを用いて解析した。

#2項分布の解析
th3 <- th2 %>% 
  filter(TRT_c2 == "1") %>% 
  summarise(
    n = n(),
    x = sum(change2))

一般に、2項分布における $p$ の推定量 $\hat{p}$ は $\hat{p} = x/n$ で与えられる。しかし分散 $V=p(1-p)$ の $p$ が明らかでないために信頼区間の推定やp値の算出でいくつかの手法に分かれる。Score法かExact法が多くの場面で適当であるため、Wald法、Score法、Exact法について紹介する。

Wald法

関数による算出がわからなかったため、式から算出した。

##Wald
x <- th3$x
n <- th3$n

pi <- x/n
v <- x*(n-x)/(n^3)
p0 <- 0.5

z <- (pi-p0)/sqrt(v)
p_val <- 2*(1-pnorm(abs(z)))
wald_cl <- round(c(max(0,pi-1.96*sqrt(v)),min(pi+1.96*sqrt(v),1)),2)

print(p_val)

## [1] 1.490234e-05

点推定値と95％信頼区間についてWald法では0.67[0.59-0.74]となった。

print(c(round(pi,2),wald_cl))

## [1] 0.67 0.59 0.74

Score法

点推定値と95％信頼区間についてScore法では連続補正なしの場合には、0.67[0.59-0.74]となった。連続補正ありの場合でも0.67[0.59-0.74]となった。

##Score
prop.test(x,n, p=p0, correct=F)　#連続補正なし

## 
##  1-sample proportions test without continuity correction
## 
## data:  x out of n, null probability p0
## X-squared = 16.667, df = 1, p-value = 4.456e-05
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.5878976 0.7371124
## sample estimates:
##         p 
## 0.6666667

prop.test(x,n, p=p0, correct=T)　#連続補正あり

## 
##  1-sample proportions test with continuity correction
## 
## data:  x out of n, null probability p0
## X-squared = 16.007, df = 1, p-value = 6.312e-05
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.5844682 0.7401791
## sample estimates:
##         p 
## 0.6666667

Exact法

点推定値と95％信頼区間についてExact法では、0.67[0.59-0.74]となった。

##Exact
binom.test(x,n, p=p0, alternative="two.sided")

## 
##  Exact binomial test
## 
## data:  x and n
## number of successes = 100, number of trials = 150, p-value = 5.448e-05
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.5851570 0.7414436
## sample estimates:
## probability of success 
##              0.6666667

2*2表の解析

2*2表の解析を紹介するために、TRT_c2とchange2によるクロス表を解析した。


x <- xtabs(formula = ~TRT_c2+change2, data = th2)
x

##       change2
## TRT_c2   0   1
##      0  41   9
##      1  50 100

独立性の検定

Yatesの連続補正を加えた $\chi^2$ 統計量は以下のように算出できる。

## chisq
chisq.test(x)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  x
## X-squared = 33.881, df = 1, p-value = 5.858e-09

Exact検定のp値は以下のように算出できる。

## exact
fisher.test(x)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x
## p-value = 2.451e-09
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   3.925819 22.794715
## sample estimates:
## odds ratio 
##    9.00197

粗解析によるリスク比、オッズ比、リスク差

リスク比、オッズ比、リスク差の推定について紹介する。epi.2by2を用いれば簡単に推定できる。リスク比は"Inc risk ratio"、オッズ比は"Odds ratio"、リスク差[%]は"Attrib risk in the exposed *"を見れば良い。

##RR, OR, RD(Attrib risk in the exposed|ARisk)
x2 <- matrix(rev(as.vector(x)), ncol=2, byrow = F)
print(x2)

##      [,1] [,2]
## [1,]  100   50
## [2,]    9   41

epi.2by2(x2, method = "cohort.count", conf.level = .95)

##              Outcome +    Outcome -      Total        Inc risk *        Odds
## Exposed +          100           50        150              66.7        2.00
## Exposed -            9           41         50              18.0        0.22
## Total              109           91        200              54.5        1.20
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 3.70 (2.03, 6.76)
## Odds ratio                                     9.11 (4.10, 20.22)
## Attrib risk in the exposed *                   48.67 (35.62, 61.72)
## Attrib fraction in the exposed (%)            73.00 (50.69, 85.22)
## Attrib risk in the population *                36.50 (23.81, 49.19)
## Attrib fraction in the population (%)         66.97 (42.46, 81.04)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 35.817 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

層別解析による共通リスク比、オッズ比、リスク差（Mantel-Haenzel推定量）

計数データの交絡因子の対処としては、大きく層別解析と回帰分析がある。回帰分析も後述するが、ここでは層別解析を取り上げる。層別解析のアイデアは、交絡因子により複数の層に分類し、各層ごとに推定を行いその結果を統合することで共通効果を推定するというものである。交絡因子がカテゴリカル変数であればよいが、連続量の場合にはカテゴリカル変数に変換するか回帰分析により調整することになる。

各層ごとの結果の統合方法としては、重み付き最小二乗推定量と、Mantel-Haenzel推定量がしばしば用いられるが、ここではMantel-Haenzel推定量(MH推定量)を紹介する。

試験前の投薬有無のpre_drugを調整する層別解析を行った。先ほど同様に、リスク比は"Inc risk ratio"、オッズ比は"Odds ratio"、リスク差[%]は"Attrib risk in the exposed *"を見れば良いが、共通効果を見る場合は（M-H）を確認する。粗解析（crude）は層別しない場合の解析であるため、先ほどの結果と一致する。今回はRCTを想定しており、データ上で交絡が生じていないため粗解析と層別解析の結果はほぼ一致する。

##staratified RR, OR, RD
th3 <- th2 %>% 
  dplyr::select(TRT_c2, change2, pre_drug) %>% 
  mutate(TRT_c3 = ifelse(TRT_c2==1,-1,0)) %>% 
  mutate(change3 = ifelse(change2==1,-1,0)) 
  
stratified_x <- table(th3$TRT_c3, th3$change3, th3$pre_drug,
                      dnn = c("TRT", "Outcome", "Pre_Drug"))

print(stratified_x)

## , , Pre_Drug = 0
## 
##     Outcome
## TRT  -1  0
##   -1 60 35
##   0   5 26
## 
## , , Pre_Drug = 1
## 
##     Outcome
## TRT  -1  0
##   -1 40 15
##   0   4 15

epi.2by2(stratified_x, method = "cohort.count", conf.level = 0.95)

##              Outcome +    Outcome -      Total        Inc risk *        Odds
## Exposed +          100           50        150              66.7        2.00
## Exposed -            9           41         50              18.0        0.22
## Total              109           91        200              54.5        1.20
## 
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio (crude)                         3.70 (2.03, 6.76)
## Inc risk ratio (M-H)                           3.71 (2.03, 6.78)
## Inc risk ratio (crude:M-H)                     1.00
## Odds ratio (crude)                             9.11 (4.10, 20.22)
## Odds ratio (M-H)                               9.31 (4.17, 20.81)
## Odds ratio (crude:M-H)                         0.98
## Attrib risk in the exposed (crude) *           48.67 (35.62, 61.72)
## Attrib risk in the exposed (M-H) *             48.78 (35.05, 62.50)
## Attrib risk (crude:M-H)                        1.00
## -------------------------------------------------------------------
##  M-H test of homogeneity of PRs: chi2(1) = 0.042 Pr>chi2 = 0.838
##  M-H test of homogeneity of ORs: chi2(1) = 0.013 Pr>chi2 = 0.909
##  Test that M-H adjusted OR = 1:  chi2(1) = 35.834 Pr>chi2 = <0.001
##  Wald confidence limits
##  M-H: Mantel-Haenszel; CI: confidence interval
##  * Outcomes per 100 population units

回帰分析によるリスク比、オッズ比、リスク差

リスク比、オッズ比、リスク差は回帰分析からも推定することができ、共変量の調整を行える。層別解析とは異なり連続量の調整も行えるが、その分モデルの仮定は追加されることになる。基本的には一般化線形モデル（GLM）^[1]の枠組みで解析できる。オッズ比はLink関数にLogit、リスク比はLink関数にLog、リスク差はLink関数に恒等関数を用いればよい。

ただし、対数二項分布だとうまく計算できなかったので、今回はリスク比の推定には修正ポアソン回帰を用いた。修正ポアソン回帰はポアソン回帰の過分散に対処するために分散にロバスト分散を利用する手法であり、ロバスト分散を使うためにGEEのpackageを利用している。

層別解析同様に、pre_drugで調整して推定した。

##REG
OR_model <- glm(change2 ~ TRT_c2 + pre_drug, data=th2, binomial(link="logit"))
# RR <- glm(change2 ~ TRT_c2 + pre_drug, data=th2, binomial(link = "log"))
RR2_model <- geeglm(change2 ~ TRT_c2 + pre_drug, data=th2, poisson(link = "log"), id=ID)
RD_model <- glm(change2 ~ TRT_c2 + pre_drug, data=th2, binomial(link = "identity"))

log_OR <- OR_model[["coefficients"]][["TRT_c21"]]
log_OR_cl <- confint(OR_model)[2, ]

OR <- c(exp(log_OR),exp(log_OR_cl))

RR <- tidy(RR2_model, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE)[2,c(2,6,7)] %>% as.matrix
RD <- c(RD_model[["coefficients"]][["TRT_c21"]],confint(RD_model)[2, ])

names(OR) <- c("Estimate","95%cl_low","95%cl_high")
names(RD) <- c("Estimate","95%cl_low","95%cl_high")

print(OR)

オッズ比

調整オッズ比の点推定値と95％信頼区間は9.35[4.36-22.1]となった。

##   Estimate  95%cl_low 95%cl_high 
##   9.350588   4.360122  22.069783

print(RR)

リスク比（修正ポアソン回帰）

調整リスク比の点推定値と95％信頼区間は3.71[2.03-6.77]となった。

##      estimate conf.low conf.high
## [1,] 3.711331 2.034354  6.770687

print(RD)

リスク差

調整リスク差の点推定値と95％信頼区間は0.49[0.35-0.60]となった。

##   Estimate  95%cl_low 95%cl_high 
##  0.4858404  0.3453432  0.6036405

対応のある検定（McNemar test）

対応のあるデータを比べる場合には、McNemar testが用いられる。投薬群に限定し、HbA1cが7.0%の閾値で2値化した変数を試験前後で比較する。McNemar testは、 $pre=0 \rightarrow post=1$ と $pre=1 \rightarrow post=0$ の人数、つまりクロス表の非対角成分が情報を持つ検定である。そのためサンプルサイズがどれほど大きくても、前後で変化した人が少なければ検出力はあまり上がらない。

##McNemar test
th4 <- filter(th2, TRT_c2 ==1 )
mcnemar <- xtabs(formula = ~pre2 + post2, data = th4)
print(mcnemar)

##     post2
## pre2  0  1
##    0 27 28
##    1 14 81

mcnemar.test(mcnemar, correct = FALSE)

## 
##  McNemar's Chi-squared test
## 
## data:  mcnemar
## McNemar's chi-squared = 4.6667, df = 1, p-value = 0.03075

（補足）McNemar testのI*J表への拡張

McNemar testは2*2表で利用されるが、さらにI*J表拡張した手法としてMcNemar-Bowker Testが存在する。プログラムはMcNemar testと同様だが、対称性の検定を行っていることに注意する。

##(補足)McNemar-Bowker Test
tmp <- matrix(c(19,12,12,9,34,12,0,12,38), ncol=3, byrow=TRUE)
print(tmp)

##      [,1] [,2] [,3]
## [1,]   19   12   12
## [2,]    9   34   12
## [3,]    0   12   38

mcnemar.test(tmp, correct = FALSE)

## 
##  McNemar's Chi-squared test
## 
## data:  tmp
## McNemar's chi-squared = 12.429, df = 3, p-value = 0.00605

I*J表の解析

対応のないI*J表のデータ解析を紹介する。本記事内では、行が治療等の割付群、列がアウトカムを表しており、I*2表、2*J表、I*J表に分けてコードを示す。しかし基本的な考え方は同じであるため、はじめに基本的な考え方を提示する。

2*2表からI*J表の解析を考えるにあたり、最も考慮すべきは尺度が名義尺度か順序尺度かの区別をする必要がある点である。

名義尺度*名義尺度であればピアソンの $\chi^2$ 統計量と似た統計量で解析でき、直感的に理解しやすい。具体的には、I*J表の多重超幾何分布を当てはめることで検定統計量として条件付き $\chi^2$ 統計量を利用できる。これは、ピアソンの $\chi^2$ 統計量を $(N-1)/N$ 倍した値と等しくなる。条件付き $\chi^2$ 統計量によるI*J表の解析は一般化Cochran-Mantel-Haenszel検定と呼ばれ、cmh_testで実行できる。自由度は、 $df=(I-1)*(J-1)$ である。

順序尺度の場合は順序を考慮するためのスコアを割り当てることを考える。例えば、順序のある $m$ 水準の各群を $c_1,c_2,...,c_m$ とすると、それぞれの群に $x_1,x_2,...,x_m$ というスコアを割り当てる。このとき $x_1 \leq x_2 \leq x_3,..., \leq x_{m-1} \leq x_m$ で少なくとも1つは $<$ が成立するような $x_i$ であれば何でも良い。1, 2, 3, 4 ...,と順に割り当てていくだけではなく、今回の仮想的な試験の薬剤の量としてPlacebo群、2.0 mg, 5.0mg, 10.0mg群を反映させた(0, 2, 5, 10)のようなスコアを割り当てても良い。

クロス表の $\chi^2$ 統計量の自由度は、 $(行数-1)*(列数-1)$ と行と列に応じて定まるがスコアを割り当てた行もしくは列は自由度を求めるときに計算に含めない。つまり、順序尺度*順序尺度であれば自由度は1になり、名義尺度(I 水準)*順序尺度の場合は自由度は $df=I-1$ となる。

"xxx"と呼ばれるなどという記載がある場合には、数理的に親しい概念が存在することが多い。例えば、I*J表で行と列ともに順序尺度の場合の解析はMantelの相関検定と呼ばれる。確かにMantelの相関検定の検定統計量 $\chi^2$ を書き下していくと、行のスコアと列のスコアをそれぞれXとYの値としたときの相関係数 $r$ と $\chi^2=(N-1)*r^2$ という関係にあることがわかる。しかしこのような関係の導出は本記事の範囲外とする。

以下、Rのコードを紹介する。chisq_testでsocresを指定すればスコアを割り当てて解析することができる。

2*J 表

2*J 表で列が順序尺度のとき、特にCochran-Armitage検定と呼ばれる。

## 2*J 表
x_2_j <- xtabs(formula = ~TRT_c2 + change3, data = th2)
x_2_j

##       change3
## TRT_c2 -2 -1  1  2
##      0  8 33  7  2
##      1  6 44 67 33

cmh_test(x_2_j)

## 
##  Asymptotic Generalized Cochran-Mantel-Haenszel Test
## 
## data:  change3 by TRT_c2 (0, 1)
## chi-squared = 37.097, df = 3, p-value = 4.388e-08

chisq_test(x_2_j, scores = list("change3"= c(1,2,3,4)))

## 
##  Asymptotic Linear-by-Linear Association Test
## 
## data:  change3 (ordered) by TRT_c2 (0, 1)
## Z = -5.6872, p-value = 1.292e-08
## alternative hypothesis: two.sided

I*2 表

I*2 表で行が順序尺度でスコアに単に1, 2, 3,...と割り当てたとき、Wilcoxon検定と呼んだりする。

## I*2 表
x_i_2 <- xtabs(formula = ~TRT_c + change2, data = th2)
x_i_2

##      change2
## TRT_c  0  1
##     0 41  9
##     1 28 22
##     2 14 36
##     3  8 42

cmh_test(x_i_2)

## 
##  Asymptotic Generalized Cochran-Mantel-Haenszel Test
## 
## data:  change2 by TRT_c (0, 1, 2, 3)
## chi-squared = 52.544, df = 3, p-value = 2.294e-11

chisq_test(x_i_2, scores = list("TRT_c"= c(1,2,3,4)))

## 
##  Asymptotic Linear-by-Linear Association Test
## 
## data:  change2 by TRT_c (0 < 1 < 2 < 3)
## Z = -7.1759, p-value = 7.185e-13
## alternative hypothesis: two.sided

chisq_test(x_i_2, scores = list("TRT_c"= c(0,2,5,10)))

## 
##  Asymptotic Linear-by-Linear Association Test
## 
## data:  change2 by TRT_c (0 < 1 < 2 < 3)
## Z = -6.8141, p-value = 9.488e-12
## alternative hypothesis: two.sided

IJ 表（名義尺度名義尺度）

条件付き $\chi^2$ 統計量による一般化Cochran-Mantel-Haenszel検定である。

## I*J 表
x_i_j <- xtabs(formula = ~TRT_c + change3, data = th2)
x_i_j

##      change3
## TRT_c -2 -1  1  2
##     0  8 33  7  2
##     1  5 23 18  4
##     2  0 14 28  8
##     3  1  7 21 21

cmh_test(x_i_j)

## 
##  Asymptotic Generalized Cochran-Mantel-Haenszel Test
## 
## data:  change3 by TRT_c (0, 1, 2, 3)
## chi-squared = 68.528, df = 9, p-value = 2.958e-11

IJ 表（順序尺度順序尺度）

I*J 表でともに順序尺度のときはMantelの相関検定と呼ばれる。

chisq_test(x_i_j, scores = list("TRT_c"= c(1,2,3,4),
                                "change3"= c(1,2,3,4)))

## 
##  Asymptotic Linear-by-Linear Association Test
## 
## data:  change3 (ordered) by TRT_c (0 < 1 < 2 < 3)
## Z = 7.4665, p-value = 8.238e-14
## alternative hypothesis: two.sided

IJ 表（順序尺度名義尺度）

片側の尺度のみにスコアを割り当てることも可能である。

chisq_test(x_i_j, scores = list("TRT_c"= c(1,2,3,4)))

## 
##  Asymptotic Generalized Pearson Chi-Squared Test
## 
## data:  change3 by TRT_c (0 < 1 < 2 < 3)
## chi-squared = 57.888, df = 3, p-value = 1.661e-12

サンプルサイズ計算

単群、対応のない2*2 表と対応のある2*2 表での割合の検定のサンプルサイズ計算のコード例を示す。

単群のサンプルサイズ計算

##単群のサンプルサイズ計算
OneSampleProportion.Equality(alpha=0.05, beta=0.2, p=0.8, differ=0.05)

## [1] 502.3283

対応のない2標本のサンプルサイズ計算

##対応のない2標本のサンプルサイズ計算
2 * TwoSampleProportion.Equality(alpha=.05, beta=0.2, p1=0.25, p2=0.15, k=1)
# 2 * power.prop.test(p1 = 0.15, p2 = 0.25, sig.level = 0.05, power = 0.8)[["n"]]

## [1] 494.4794

McNemar検定のサンプルサイズ計算

##McNemar検定（正規近似）のサンプルサイズ計算
p01 <- 0.15
p10 <- 0.45
n <- McNemar.Test(alpha=0.05, beta=0.2, p01/p10, p01+p10)
n

## [1] 49.90101

まとめ

計数データ（カテゴリカルデータ）について、二項分布、2*2表、I*J表の解析とサンプルサイズ計算を紹介した。特に2標本でのリスク比、オッズ比、リスク差は粗解析だけではなく、層別解析（Mantel-Haenzel推定量）、回帰分析による解析も行った。またI*J表の解析では、一般化Cochran-Mantel-Haenszel検定にとどまらず順序尺度へのスコアの割当も言及した。対応のあるデータの解析としては、McNemar検定とMcNemar-Bowker検定のコードを示した。もし本記事内で誤りを見つけられた方は、ご連絡頂ければ幸いである。

（メモ）今後、追加する可能性のある解析

拡張Mantel検定
カッパ係数

サポートして頂けるとモチベーションに繋がりますのでぜひ宜しくお願いします。データ解析や臨床研究でのご相談があれば、お気軽にTwiiterもしくはメールにてご連絡下さい。

作成者：Masahiro Kondo
作成日：2022/3/4
連絡先：m.kondo1042(at)gmail.com

脚注

GLMがよくわからない場合は、一般化線形モデル入門を読むとよい。 ↩︎

Rによる計数データ（カテゴリカルデータ）解析

はじめに

目次

Package

仮想的なランダム化比較試験のデータ生成

データの要約

カテゴリ化した変数の要約

二項分布の解析

Wald法

Score法

Exact法

2*2表の解析

独立性の検定

粗解析によるリスク比、オッズ比、リスク差

層別解析による共通リスク比、オッズ比、リスク差（Mantel-Haenzel推定量）

回帰分析によるリスク比、オッズ比、リスク差

対応のある検定（McNemar test）

（補足）McNemar testのI*J表への拡張

I*J表の解析

2*J 表

I*2 表

IJ 表（名義尺度名義尺度）

IJ 表（順序尺度順序尺度）

IJ 表（順序尺度名義尺度）

サンプルサイズ計算

単群のサンプルサイズ計算

対応のない2標本のサンプルサイズ計算

McNemar検定のサンプルサイズ計算

まとめ

（メモ）今後、追加する可能性のある解析

Discussion

はじめに

目次

Package

仮想的なランダム化比較試験のデータ生成

データの要約

カテゴリ化した変数の要約

二項分布の解析

Wald法

Score法

Exact法

2*2表の解析

独立性の検定

粗解析によるリスク比、オッズ比、リスク差

層別解析による共通リスク比、オッズ比、リスク差（Mantel-Haenzel推定量）

回帰分析によるリスク比、オッズ比、リスク差

対応のある検定（McNemar test）

（補足）McNemar testのI*J表への拡張

I*J表の解析

2*J 表

I*2 表

I*J 表（名義尺度*名義尺度）

I*J 表（順序尺度*順序尺度）

I*J 表（順序尺度*名義尺度）

サンプルサイズ計算

単群のサンプルサイズ計算

対応のない2標本のサンプルサイズ計算

McNemar検定のサンプルサイズ計算

まとめ

（メモ）今後、追加する可能性のある解析

Discussion

IJ 表（名義尺度名義尺度）

IJ 表（順序尺度順序尺度）

IJ 表（順序尺度名義尺度）