差がゼロという帰無仮説以外の仮説の下でのp値の計算

『宇宙怪人しまりす　統計よりも重要なことを学ぶ』朝倉書店（以下しまりす本３）のpp.12「回復割合の差は20%」のp値を計算しようとしたが合わない。何かぼくのミスだと思うが原因がわからない。

x <- c(58, 62); n <- c(80, 100); X <- cbind(x, n-x)
colnames(X) <- c("回復", "未回復")
rownames(X) <- c("ヨクナール", "経過観察")
print(X)
s <- c(colSums(X)[1]/sum(X), colSums(X)[2]/sum(X))
# 帰無仮説：回復割合が処置と独立の場合の期待度数E
print(E <- rbind(ヨクナール=rowSums(X)[1]*s, 経過観察=rowSums(X)[2]*s))
E/rowSums(E) # 確認：回復割合にヨクナールと経過観察で差はない
1-pchisq(sum((X-E)^2/E), 1) # 帰無仮説の下でのp値
Z <- c(0.1, -0.1)
# ヨクナール群の回復割合が経過観察群より20%良い場合の期待度数EE
print(EE <- rbind(ヨクナール=rowSums(X)[1]*(s+Z), 経過観察=rowSums(X)[2]*(s-Z)))
EE/rowSums(EE) # 確認：回復割合がヨクナールの方が20％良い
1-pchisq(sum((X-EE)^2/EE), 1) # ヨクナールの回復割合が20%高い仮説の下でのp値

というRコードを実行すると、 差がないという帰無仮説（しまりす本３では「ゼロ仮説」）の下でのp値は本と同じく13.8%になるのだけれども、最後に表示される「回復割合の差が20%」という仮説の下でのp値は16.4%になってしまって、しまりす本３に書かれている17.2%とは合わない。しまりす本３ではpp.13-14に17.2%の算出方法が書かれていて、上のコードで使ったカイ二乗適合度検定とは違う方法（下記）なのだけれども、どちらでも一致するはずなんだがなあ。自分の頭の回転が鈍くて嫌になる。

rdp <- function(a, b, N1, N0, np=0) {
 RD <- a/N1-b/N0
 seRD <- sqrt(a*(N1-a)/N1^3+b*(N0-b)/N0^3)
 Z <- (RD-np)/seRD
 return(data.frame(est=RD, se=seRD, p=2*pnorm(-abs(Z))))
} 

とRで関数定義すると（注：当初、返り値をlistにしていたが、見やすさのためdata.frameに変更した）、確かにrdp(58, 62, 80, 100, np=0.2)は17.2%になるのだが、今度はrdp(58, 62, 80, 100)が13.2%になってしまって本と合わない。

だが実はこの計算はfmsbパッケージでriskdifference()に使っている式と同じなので、p値関数を描くpp.15は、以下の関数定義で実現する。 

# 割合の差についてのp値関数を描く関数定義
pvp2 <- function(a, b, N1, N0, ...) {
 ranges <- -1000:1000/1000
 RVs <- rdp(a, b, N1, N0, ranges)
 plot(RVs$p ~ ranges, type="l", ...)
}
# シーマ、『プロットして』
par(mar=c(5,4,4,4))
pvp2(58, 62, 80, 100, xlim=c(-0.2, 0.4), main="p値関数、信頼区間関数",
 xlab="設定した回復割合の差", ylab="", axes=FALSE)
grid()
axis(1, -2:4/10, sprintf("%d%%", -2:4*10))
axis(2, 0:5/5, sprintf("%d%%", 0:5*20), las=1)
mtext("両側p値", 2, line=3)
axis(4, 0:5/5, sprintf("%d%%", 5:0*20), las=1)
mtext("信頼係数", 4, line=3)
# シーマ、『グラフに80%信頼区間と95%信頼区間を追加して』
library(fmsb)
CI80 <- riskdifference(58, 62, 80, 100, conf.level=0.80)$conf.int
arrows(CI80[1], 0.2, CI80[2], 0.2, code=3)
text(58/80-62/100, 0.2, "80%信頼区間", pos=3)
CI95 <- riskdifference(58, 62, 80, 100, conf.level=0.95)$conf.int
arrows(CI95[1], 0.05, CI95[2], 0.05, code=3)
text(58/80-62/100, 0.05, "95%信頼区間", pos=3)