佐藤俊哉『宇宙怪人しまりす統計よりも重要なことを学ぶ』第1話の最初の表の数値例 (ver.2)

shimarisu_ch1_reproduce.R

########
## 佐藤俊哉『宇宙怪人しまりす統計よりも重要なことを学ぶ』（朝倉書店）
## Chapter 1: 再現性
########
library(dplyr)
library(ggplot2)
TPratio = function(true_ratio,
                           beta = 0.8,
                           alpha = 0.05){
  study1 =cbind(true_ratio*c(beta,1-beta), 
                (1-true_ratio)*c(alpha*0.5, 1-alpha*0.5))
  study1[1,1]/sum(study1[1,])  
}

reproducibility = function(true_ratio,
                           beta = 0.8,
                           alpha = 0.05){
  study1 =cbind(true_ratio*c(beta,1-beta), 
                (1-true_ratio)*c(alpha*0.5, 1-alpha*0.5))
  (study1[1,1]*beta+study1[1,2]*alpha*0.5)/sum(study1[1,])  
}

z = seq(0, 0.5, length.out=101)
TPdf = expand.grid(alpha = c(0.01, 0.05, 0.1), 
                    beta = c(0.7,0.8)) %>% 
  group_by(alpha, beta) %>% 
  reframe(z=z,
          value = sapply(z, TPratio, alpha=alpha, beta=beta))

ggplot(TPdf, aes(x=z, y=value,
                  colour = alpha,
                  group = interaction(alpha, beta)))+
  facet_grid(col=vars(beta), labeller = label_both) +
  geom_line() + ylim(c(0,1)) + 
  labs(color="alpha",
       x = "positive ratio in the population",
       y = "true positive ratio among the test positive")+
  theme_bw(16)+scale_color_viridis_c()
ggsave("TP.png", width=8, height=7)

####
RPdf = expand.grid(alpha = c(0.01, 0.05, 0.1), 
                    beta = c(0.7,0.8)) %>% 
  group_by(alpha, beta) %>% 
  reframe(z=z,
          value = sapply(z, reproducibility, alpha=alpha, beta=beta))


ggplot(RPdf, aes(x=z, y=value,
                  colour = alpha,
                  group = interaction(alpha, beta)))+
  facet_grid(col=vars(beta), labeller = label_both) +
  geom_line() + ylim(c(0,1)) + 
  labs(color="alpha",
       x = "positive ratio in the population",
       y = "reproducibility")+
  theme_bw(16)+scale_color_viridis_c()
ggsave("RP.png", width=8, height=7)