『宇宙怪人しまりす　統計よりも重要なことを学ぶ』ヨクナール使用とかぜ症状の回復

shimarisu_pfun.R

#佐藤俊哉『宇宙怪人しまりす　統計よりも重要なことを学ぶ』（朝倉書店）より
X = matrix(c(58, 22,
            62, 38), byrow = TRUE, nrow = 2, ncol = 2)
print(X)

res_chisq = chisq.test(X, correct = FALSE)
print(res_chisq$p.value)
#[1] 0.1375639

x = X[,1]
n = rowSums(X)
res_prop = prop.test(x, n = n, correct = FALSE)
print(res_prop$p.value)
#[1] 0.1375639
####
#「差がない」とした推定値で求めた標準誤差
phat = x/n
p_pooled = sum(x)/sum(n)
delta = phat[1]-phat[2]
se0 = sqrt(p_pooled*(1-p_pooled)*sum(1/n))
print(pchisq(abs(delta/se0)^2, 1, lower.tail = FALSE))
#[1] 0.1375639
print(se0)
#[1] 0.07071068
###
#差を許す推定値で求めた標準誤差
se = sqrt(sum(phat*(1-phat)/n))
print(se)
#[1] 0.06962893
pvalfun_wald <- function(mu0, delta, se){
  pchisq(((delta-mu0)/se)^2, df = 1, lower.tail = FALSE)  
}

mu_v <- seq(-0.2, 0.4, length.out=501)
res_p = sapply(mu_v, pvalfun_wald, delta=delta, se=se)
#png("pfun.png")
plot(mu_v, res_p, type="l", xlab = "risk difference", ylab = "p-value")
#dev.off()

cat("H0:diff=0",round(pvalfun_wald(0, delta=delta, se=se)*100, 1),"%")
#H0:diff=0 13.2 %

cat("H0:diff=0.2",round(pvalfun_wald(0.2, delta=delta, se=se)*100, 1),"%")
#H0:diff=0.2 17.2 %