Rの関数:addz2ci {PropCIs}

Rの関数から addz2ci {PropCIs} を確認します。

関数 `addz2ci` とは

関数 addz2ci は、二項分布に基づく母割合（例：成功確率）の信頼区間を算出する手法の一つである、Agresti-Coull（アグレスティ・クール）法による近似信頼区間（Add-z² 信頼区間） を計算する関数です。

「Wald（ワルド）法」の信頼区間は、計算が容易である反面、成功回数または失敗回数が極端に少ない場合（特にゼロの場合）、標準誤差がゼロとして計算されてしまい、信頼区間が「幅を持たない一点」に潰れてしまいます。

全試行が失敗に終わったという観測結果のみをもって「真の成功確率は0%である」と判断することは、推論として不適切です。

Agresti-Coull 法は、観測された成功回数および失敗回数に一定の仮想的な回数(\(z^2 / 2\))を加算した上で Wald法を適用する、という補正を行います。

関数 `addz2ci` の引数

library(PropCIs)
args(addz2ci)

function (x, n, conf.level) 
NULL

x
- 観測された成功回数を表す数値です。
- 負ではない整数を指定します。
n
- 試行の総数を表す数値です。
- 正の整数を指定します。当然ながら x 以上の値である必要があります。
conf.level
- 推定する信頼区間の信頼水準を表す数値です。
- 0より大きく1より小さい数値を指定します。

シミュレーションコード

以下に、addz2ci 関数の機能を確認するためのサンプルデータを用いたシミュレーションコードを示します。

本シミュレーションでは、試行回数を20回に固定し、成功回数が0回から20回まで変動した場合の信頼区間を計算します。

従来型のWald法を手動で計算し、addz2ci 関数の出力結果と比較することで、極端なデータにおける本関数との差異を視覚的に確認します。

サンプルデータの設定

# パッケージの読み込み
library(ggplot2)
library(dplyr)

# 1. サンプルデータの設定
# 試行回数を20回に固定し、成功回数が極端なケースを含む配列を準備します。
n_trials <- 20
x_successes <- c(0, 1, 2, 5, 10, 15, 18, 19, 20)

cat(sprintf("設定: 試行回数 = %d 回\n", n_trials))
cat("検証する成功回数:", paste(x_successes, collapse = ", "), "\n\n")

# 結果を格納するための空のデータフレームを作成
df_results <- data.frame()

設定: 試行回数 = 20 回
検証する成功回数: 0, 1, 2, 5, 10, 15, 18, 19, 20

各成功回数に対する信頼区間の計算

# 2. 各成功回数に対する信頼区間の計算
for (x_val in x_successes) {
  
  # 手法A: addz2ci を用いた Agresti-Coull 法の計算
  res_addz2 <- addz2ci(x = x_val, n = n_trials, conf.level = 0.95)
  
  # 手法B: 標準的な Wald 法の計算 (比較用として手動計算)
  p_hat <- x_val / n_trials
  z_val <- qnorm(1 - (1 - 0.95) / 2)
  # Wald法の標準誤差
  se_wald <- sqrt(p_hat * (1 - p_hat) / n_trials)
  # 0から1の範囲に収まるように制限
  wald_lower <- max(0, p_hat - z_val * se_wald)
  wald_upper <- min(1, p_hat + z_val * se_wald)
  
  # 計算結果をデータフレームに追加
  df_results <- rbind(df_results, data.frame(
    Success = x_val,
    Method = "1. Add-z^2法 (本関数)",
    Estimate = res_addz2$estimate,
    Lower = res_addz2$conf.int[1],
    Upper = res_addz2$conf.int[2]
  ))
  
  df_results <- rbind(df_results, data.frame(
    Success = x_val,
    Method = "2. Wald法 (従来手法)",
    Estimate = p_hat,
    Lower = wald_lower,
    Upper = wald_upper
  ))
}

cat("【成功回数が 0回 の場合の結果】\n")
print(df_results %>% filter(Success == 0))
cat("\n")

cat("【成功回数が 10回 の場合の結果】\n")
print(df_results %>% filter(Success == 10))
cat("\n")

【成功回数が 0回 の場合の結果】
  Success                Method   Estimate Lower     Upper
1       0 1. Add-z^2法 (本関数) 0.08056258     0 0.1898096
2       0  2. Wald法 (従来手法) 0.00000000     0 0.0000000

【成功回数が 10回 の場合の結果】
  Success                Method Estimate     Lower     Upper
1      10 1. Add-z^2法 (本関数)      0.5 0.2992980 0.7007020
2      10  2. Wald法 (従来手法)      0.5 0.2808694 0.7191306

成功回数が 0回の場合

Wald法: 推定値 0.000、信頼区間 [0.000, 0.000]
Add-z²法: 推定値 0.081、信頼区間[0.000, 0.190]

Wald法では、観測された成功回数が0回であるという事実のみに基づき、推定値も信頼区間の幅も完全にゼロとなってしまっています。

つまり「20回の試行で1回も成功しなかった」という観測結果から、「真の成功確率は0%である」と不適切な推定を導いています。

一方、Add-z²法（本関数）では、推定値が約8.1%へと引き上げられ、上限値として約19.0%が算出されています。

仮想的な成功と失敗を加算する補正が働いた結果、「たまたま20回連続で失敗しただけであり、真の成功確率は最大で19%程度存在する可能性がある」という、不確実性を考慮した推論結果に至っています。

成功回数が 10回の場合

Wald法: 推定値 0.500、信頼区間 [0.281, 0.719]
Add-z²法: 推定値 0.500、信頼区間[0.299, 0.701]

成功回数が半分の10回（割合50%）となるケースでは、両手法の推定値は 0.5 で完全に一致しています。

また、信頼区間の下限および上限を見比べても、双方の手法でほぼ同等の区間幅が算出されています。

本結果は、割合が0%や100%といった極端な領域から離れるにつれて、Add-z²法による補正の影響が相対的に小さくなり、従来手法と同等の結果へと収束していく性質を表しています。

可視化プロットの作成

# 3. 可視化プロットの作成
p1 <- ggplot(df_results, aes(x = factor(Success), y = Estimate, color = Method)) +
  # 推定値を点で描画
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  # 信頼区間をエラーバーで描画
  geom_errorbar(aes(ymin = Lower, ymax = Upper), 
                position = position_dodge(width = 0.5), width = 0.4, linewidth = 1.2) +
  # 確率の理論的境界線 (0 と 1)
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50") +
  
  scale_color_manual(name = "計算手法", values = c("1. Add-z^2法 (本関数)" = "blue", "2. Wald法 (従来手法)" = "red")) +
  labs(
    title = "成功回数の変動に伴う信頼区間の推移 (試行数 n = 20)",
    x = "成功回数",
    y = "母割合の推定値および 95%信頼区間"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(p1)