1 t-検定を連続実行, 箱ひげ図描画, 有意水準

1.1 データ

x1 <- iris[c("Sepal.Length", "Species")]
x2 <- airquality[c("Ozone", "Month")]
x1 & x2
Sepal.Length Species
5.1 setosa
4.9 setosa
4.7 setosa
4.6 setosa
5.0 setosa
5.4 setosa
Ozone Month
41 5
36 5
12 5
18 5
NA 5
28 5

1.2 全組み合わせを描画

# 関数定義: t.testを連続実行,箱ひげ図, 有意水準を図示
plev <- function(x, plv=0.05, wseg=-1){
  v <- x[[1]] ; fct <- x[[2]]
  
  # factor水準の全組み合わせ
  cmbs <- combn(levels(factor(fct)), 2, simplify = F) 
  
  # t.test 及び、箱ひげ図のformula
  frm <- as.formula(paste(names(x), collapse = "~"))
  
  # 全組み合わせのt.test と p-valueを取得 
  tts <- lapply(seq_along(cmbs), function(i){
    subd <- subset(x, fct %in% cmbs[[i]])
    t.test(frm, data = subd)
  })
  names(tts) <- sapply(cmbs, function(i)paste(i, collapse = "_vs_"))
  pvs <- sapply(tts, function(x)x$p.value)
  
  # 箱ひげ図描画
  if(wseg <= 0 ) wseg <- (max(v, na.rm = T)-min(v, na.rm = T))/(length(cmbs)*2)
  b <-boxplot(formula = frm, data = x, 
              las = 2, xlab = "", 
              ylim = c(min(v, na.rm = T), max(v, na.rm = T) + (wseg*length(cmbs))))
  
  # セグメントの座標, テキストラベル
  xpos <- combn(1:nlevels(factor(fct)), 2, simplify = F)
  xpos0 <- sapply(xpos, "[", 1)
  xpos1 <- sapply(xpos, "[", 2)
  ypos0 <- seq(from = max(v, na.rm = T) + wseg, 
               to = max(v, na.rm = T) + wseg*length(cmbs), 
               by = wseg)
  ypos1 <- ypos0 - wseg/2
  tposy <- ypos0 + wseg/4
  tposx <- xpos0 + (xpos1-xpos0)/2
  
  labs <- ifelse(pvs<0.01, "**", ifelse(pvs<0.05, "*", "")) 
  
  # セグメントを追記
  for(i in seq_along(cmbs)){
    xp0 <- xpos0[i]; xp1 <- xpos1[i]
    yp0 <- ypos0[i]; yp1 <- ypos1[i]
    
    if (pvs[i] < plv){
      segments(x0 = xp0, y0 = yp0, x1 = xp1, y1 = yp0, lwd = 2)
      segments(xp0, yp0, xp0, yp1, lwd = 2)
      segments(xp1, yp0, xp1, yp1, lwd = 2)
      text(x = tposx[i], y=tposy[i], labels=labs[i], cex = 2, col = 2)    
    }
  }
  return(tts)
}

# データ
x1 <- iris[c("Sepal.Length", "Species")]
x2 <- airquality[c("Ozone", "Month")]

# 実行
par(mfcol=c(1,2))
tt1 <- plev(x1)
tt2 <- plev(x2, wseg = 20)

1.3 組み合わせごとに描画

# 関数定義
pair_plev <- function(x, plv=0.05, wseg=-1){
  v <- x[[1]] ; fct <- x[[2]]
  
  # factor水準の全組み合わせ
  cmbs <- combn(levels(factor(fct)), 2, simplify = F) 
  
  # t.test 及び、箱ひげ図のformula
  frm <- as.formula(paste(names(x), collapse = "~"))
  
  # 全組み合わせのt.test と p-valueを取得 
  tts <- lapply(seq_along(cmbs), function(i){
    subd <- subset(x, fct %in% cmbs[[i]])
    t.test(frm, data = subd)
  })
  names(tts) <- sapply(cmbs, function(i)paste(i, collapse = "_vs_"))
  pvs <- sapply(tts, function(x)x$p.value)
  
  # テキストラベルとセグメントを追記する領域の指定
  labs <- ifelse(pvs<0.01, "**", ifelse(pvs<0.05, "*", ""))
  if(wseg <= 0 ) wseg <- (max(v, na.rm = T)-min(v, na.rm = T))/(length(cmbs)*2)
  
  # 箱ひげ図描画、セグメントを追記
  for(i in seq_along(cmbs)){
    dat <- subset(x, fct %in% cmbs[[i]])
    v <- dat[[1]]; dat[[2]] <- factor(dat[[2]])
    boxplot(formula = frm, data = dat, ylim = c(min(v, na.rm = T), max(v, na.rm = T) + wseg*2))
    
    yp0 <- max(v, na.rm = T) + wseg
    yp1 <- yp0 - wseg/2
    tpy <- yp0 + wseg/2
    segments(1, yp0, 2, yp0, lwd = 2)
    segments(1, yp0, 1, yp1, lwd = 2)
    segments(2, yp0, 2, yp1, lwd = 2)
    text(x = 1+0.5, y=tpy, labels=labs[i], cex = 2, col = 2)
  }
  return(tts)

}

# データ
x1 <- iris[c("Sepal.Length", "Species")]
x2 <- airquality[c("Ozone", "Month")]

par(mfrow=c(1,3))
res <- pair_plev(x1)

par(mfrow=c(2,5))
res <- pair_plev(x2)

2 環境.

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.8.3      highr_0.9         pillar_1.8.0      bslib_0.3.0      
##  [5] compiler_4.0.3    jquerylib_0.1.4   plyr_1.8.7        tools_4.0.3      
##  [9] digest_0.6.29     viridisLite_0.4.0 jsonlite_1.7.2    evaluate_0.14    
## [13] lifecycle_1.0.3   tibble_3.1.6      pkgconfig_2.0.3   rlang_1.0.6      
## [17] cli_3.5.0         DBI_1.1.1         rstudioapi_0.13   yaml_2.2.1       
## [21] xfun_0.25         rsko_0.1.0        fastmap_1.1.0     kableExtra_1.3.4 
## [25] xml2_1.3.2        httr_1.4.2        dplyr_1.0.8       stringr_1.4.0    
## [29] knitr_1.33        systemfonts_1.0.2 generics_0.1.3    vctrs_0.5.1      
## [33] sass_0.4.0        webshot_0.5.2     tidyselect_1.1.2  svglite_2.0.0    
## [37] glue_1.6.2        R6_2.5.1          fansi_1.0.3       rmarkdown_2.10   
## [41] purrr_0.3.4       magrittr_2.0.3    scales_1.2.0      htmltools_0.5.2  
## [45] ellipsis_0.3.2    assertthat_0.2.1  rvest_1.0.1       colorspace_2.0-3 
## [49] utf8_1.2.2        stringi_1.7.6     munsell_0.5.0