層別解析：佐藤俊哉『宇宙怪人しまりす　統計よりも重要なことを学ぶ』（朝倉書店）第4話より

shimarisu_ch4_strata.R

library(DiagrammeR)
library(DiagrammeRsvg)
library(rsvg)
library(dplyr)
library(gt)

## 佐藤俊哉『宇宙怪人しまりす　統計よりも重要なことを学ぶ』（朝倉書店）より
g <- grViz("digraph{
                     graph[rankdir = TB]
                     node[shape = rectangle]
                     Age[label='年齢']
                     Stage[label = 'かぜの\n重症度']
                     Treat[label = 'ヨクナール\n使用']
                     Recover[label = 'かぜ症状の\n回復']
                     RFever[label = '解熱作用']
                     edge[color = black]
                     Age -> Stage
                     Age -> Recover
                     Stage -> Treat
                     Treat -> RFever
                     Stage -> Recover
                     RFever -> Recover
                     }")

print(g)
g %>% export_svg() %>% 
  charToRaw() %>% 
  rsvg_png("graph.png", width = 900, height = 900)
###
#リスク比
RR = function(X,level=0.95, confint=TRUE){
  x = X[,1]
  n = rowSums(X)
  phat = x/n
  beta = unname(log(phat[1]) - log(phat[2]))
  if(confint){
    se = sqrt(sum(1/x)-sum(1/n))
    alpha = (1-level)*0.5
    z = qnorm(1-alpha, 0, 1)
    lower = exp(beta - se*z)
    upper = exp(beta + se*z)
    return(list(estimates=exp(beta), 
                SE=se,
                CI=unname(c(lower,upper))))    
  }else{
    return(exp(beta))
  }
}
#オッズ比
OR = function(X,level=0.95, confint=TRUE){
  tau = unname((X[1,1]/X[1,2])/(X[2,1]/X[2,2]))
  if(confint){
    logor <- log(tau)
    se <- sqrt(sum(1/X))
    alpha = (1-level)*0.5
    z = qnorm(1-alpha, 0,1)
    lower = exp(logor - se*z)
    upper = exp(logor + se*z)
    return(list(estimates=tau, SE=se, CI=unname(c(lower,upper))))    
  }else{
    return(tau)
  }
}
#リスク差
RD = function(X, level=0.95, confint=TRUE){
  x = X[,1]
  n = rowSums(X)
  phat = x/n
  delta = unname(phat[1] - phat[2])
  if(confint){
    alpha = (1-level)*0.5
    z = qnorm(1-alpha,0,1)
    se = sqrt(sum(phat*(1-phat)/n))
    lower = delta - se*z
    upper = delta + se*z
    return(list(estimates=delta,
                SE=se,
                CI=unname(c(lower,upper))))    
  }else{
    return(delta)
  }
}
###
lname = list("ヨクナール" = c("使用","未使用"),
     "回復" = c("回復","未回復"),
     "重症度"= c("重症","軽症"))
dat = array(c(40,20,20,20,
              18,42,2,18), dim=c(2,2,2), dimnames = lname) 

# サブグループ解析（しまりす君が層別解析と間違えてやったもの）
res_or = apply(dat, 3, OR)
res_rr = apply(dat, 3, RR)
res_rd = apply(dat, 3, RD)

extract_ci <- function(res){
  cbind(
    sapply(res,function(x)x$estimates), 
    t(sapply(res,function(x)x$CI))
  )
}

tab_rd <- extract_ci(res_rd)
tab_rr <- extract_ci(res_rr)
tab_or <- extract_ci(res_or)


#差(％)
t1 = data.frame(round(tab_rd*100,1)) %>% 
  rename("点推定値"=X1,
         "下限（95％信頼区間）"=X2,
         "上限（95％信頼区間）"=X3) %>% 
  tibble::rownames_to_column("重症度") 
#t1

#比
t2 = data.frame(round(tab_rd,2)) %>% 
  rename("点推定値"=X1,
         "下限（95％信頼区間）"=X2,
         "上限（95％信頼区間）"=X3) %>% 
  tibble::rownames_to_column("重症度")
#t2

t3 = data.frame(round(tab_or,2)) %>% 
  rename("点推定値"=X1,
         "下限（95％信頼区間）"=X2,
         "上限（95％信頼区間）"=X3) %>% 
  tibble::rownames_to_column("重症度") 
#t3
tab1 = gt(bind_rows(t1,t2,t3)) %>% 
  tab_row_group(
    group = "回復割合のオッズ比",
    rows = 5:6
  ) %>% 
  tab_row_group(
    group = "回復割合の比",
    rows = 3:4
  )%>% 
  tab_row_group(
    group = "回復割合の差(％)",
    rows = 1:2
  ) 

print(tab1)
gtsave(tab1, filename = "table1.png")
###

#層化統合解析（a.k.a.層別解析）
M2 = apply(dat[1,,],2,sum) #ヨクナールグループの重症者数・軽症者数
r_s = dat[2,1,]/apply(dat[2,,], 2,sum) #経過観察のときの回復割合
expected = c(sum(r_s*M2), sum((1-r_s)*M2)) #もしヨクナールグループが経過観察だったら
M1 = apply(dat[1,,],1,sum) #ヨクナールグループの回復者数・未回復者数
tab_st = rbind("ヨクナール"=M1, "ヨクナール*"=expected)

tab2 = data.frame(tab_st) %>% 
  tibble::rownames_to_column("グループ") %>% 
  gt()

print(tab2)
gtsave(tab2, filename = "tab2.png")

cat("Risk ratio:", round(RR(tab_st, confint = FALSE),2), "\n")
#Risk ratio: 1.32 

cat("Risk difference:", round(RD(tab_st, confint = FALSE)*100,1), "%\n")
#Risk difference: 17.5 %