Go to Qiita Advent Calendar Top
LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

Day 1

サイバーフィジカルシステムに使える技術 Advent Calendar 2025

@3084LAB(澤橋 龍之介)

【R】有意差バーを自動的に配置

Posted at

1. はじめに

Rではggsignifパッケージを使うことにより，有意差バーをグラフに数値的な位置で配置することができます．
しかし，配置する位置の定義は自分自身で行う必要があります．
本記事ではデータの最大値や組み合わせに合わせて自動的に配置する手法を共有します．
※群数が多いデータではまだ試しておりません．
※綺麗にする前の汚いコードです．

2. 対応する統計手法

今回は統計手法として，
Wilcox signed rank test + hochberg法でp値補正した多重比較をしました．
ひとまずdunn.testとpairwise.wilcox.testの結果が扱えます．

3. コード

### +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
### 箱ひげ図を作成するプログラム(R)
###                                         Ryunosuke Sawahashi
###
### version:
###   R: 4.2.1    ggplot2: 3.5.1
### Q1からQ9まで対応 2桁になるとsubstring()でバグ
### +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

### ggplotの座標定義 (箱ひげ図の場合)
###  y
###     |
### 100 |           *
###     | (x,y)> -------- <-(xend,yend)
###  80 |          
###     |        |      
###  60 |       +-+     |
###     |       | |    +-+
###  40 |       | |    | |
###     |       +-+    +-+
###  20 |        |      |
###     |
###   0 |-------------------------
###            "None" "Base"   
###
###      x       1.0    2.0 ...


### -----------------------------------------------------------
### 初期インストール
### -----------------------------------------------------------
### 依存関係
### ggplot2: 拡張したグラフ描画
#install.packages("ggplot2")
### usethisをロード中にrlang周りでエラーが出たら
### remove.packages("rlang")実行とrlang.dllを一度削除
#install.packages("devtools")
### ggpattern:  ggplotのパターン塗りつぶしに使える
#install_github("coolbutuseless/ggpattern")
### dunn.text:  dunn検定
#install.packages("dunn.test")
### ggsignif:   有意差のアスタリスクを簡潔に表示
#install.packages("ggsignif")
### tidyverse:  文字列を扱いやすくする関数
#install.packages("tidyverse")

### -----------------------------------------------------------
### ライブラリ読み込み
### -----------------------------------------------------------
library("ggplot2")
library("dunn.test")
library("exactRankTests")
library("ggsignif")
library("tidyverse")

# 有意水準
alpha <- 0.05

# ラベルの順序(グラフ描画や有意差バー作成用)
#TODO: CSVのcondition列から自動抽出
labelOrder <- c("None","Device","Normal","Triple")

### -----------------------------------------------------------
### 箱ひげ図描画関数
### {augument}dat(list): 数値データ
### {augument}pvalue(list): p値
### -----------------------------------------------------------
drawGraph <- function(dat, plotPoints, questN) {
  
  # フォントをwindows内のフォントに変更(Windows限定)
  # ggsave()で反映されないなど問題あり
  #windowsFonts(TNR = windowsFont("Times New Roman"))

  # RGui内にwindowを作成(横, 縦)
  win.graph(4,5)

  # 箱ひげ図
  g <- ggplot(dat) + 
  geom_boxplot(
    aes(x=condition,y=score),
    fill=c("#e78d1a", "#1a592f", "#3599D5", "#5f4165"),
    width=0.5,
  ) +
  #theme_minimal(base_family="serif") +
  theme(
    # フォント
    axis.text.x = element_text(size=16, family="serif"),
    axis.text.y = element_text(size=16, family="serif"),
    axis.title.x = element_text(size=16, family="serif"),
    axis.title.y = element_text(size=16, family="serif"),
    aspect.ratio = 1.5,  						# プロットエリアの縦/横比
    #プロットエリア
    panel.grid.major.y = element_line(color="black", linetype=4),	# 水平Majorグリッド
    panel.background = element_rect(fill=NA),			# 背景をclear
    panel.grid.major.x = element_blank(),				  # 垂直グリッドOFF
    panel.grid.minor.y = element_blank(),         # 水平Minorグリッド OFF
    panel.border = element_rect(fill=NA, linewidth=1.0)		# プロットエリアの枠線
  ) +
  #scale_y_continuous(limits=c(0, 100)) +
  # y軸メモリは
  scale_y_continuous(breaks=seq(0,4,1), limits=c(0,4)) +	# y軸メモリ
  xlab("") +              	# x軸ラベル
  ylab("Task time [s]") + 	# y軸ラベル
  # 有意差バー
  geom_text(
        aes(x=(plotPoints$x+plotPoints$xend)/2,
        y=plotPoints$y+1, label=plotPoints$annotation), plotPoints, size=8, family="serif") +
  geom_segment(aes(x=plotPoints$x, xend=plotPoints$xend, 
      y=plotPoints$y, yend=plotPoints$yend), plotPoints)

  ### 画像出力
  # pngで出力
  ggsave(paste("q", questN,".png", sep=""), dpi=600)
}


### -----------------------------------------------------------
### 有意差のあるラベルのリストを作成する関数
### {augument}pvalues(list): p値
### -----------------------------------------------------------
makePValueBar <- function (pvalues, comparisons, maxvalues) {
  numP <- 0       # 有意差のペアの数
  pList <- c()    # p値
  cmpList <- c()  # 比較群　ex) "Base - Mass"
  # 有意差のあるペアだけ抽出
  index <- 0
  for (p in pvalues) {
    index <- index + 1
    if (p < alpha) {
      numP = numP + 1
      pList[numP] <- p
      cmpList[numP] <- comparisons[index]
    }
  }
  print(numP)
  print(pList)
  print(cmpList)
  x <- c()
  xend <- c()
  xIndex <- c()     # オフセット前の数値を使うため
  xendIndex <- c()  # オフセット前の数値を使うため
  y <- c()
  yend <- c()
  annotation <- c()

  if (numP > 0) {
    # (1) x軸
    # ラベルを数値に変換
    # 空白とハイフン( - )を除去しxとxendに代入
    counter = 0
    for (c in cmpList) {
      counter = counter + 1
      # 分割
      # ex) "Base - None" -> [1]"Base" [2]"None"
      xCoodinates <- unlist(strsplit(c, " - "))
      # ラベルの順序リスト(labelOrder)を参照して要素番号の数値として変換
      # ex) x1:"Base"->2 x2:"None"->1
      # NOTE: str_which()だと部分一致 match()が完全一致
      x1 <- str_which(labelOrder, xCoodinates[1])
      x2 <- str_which(labelOrder, xCoodinates[2])
      # x, xendを小,大になるようにする
      if (x1 < x2) {
        xIndex[counter] <- x1
        xendIndex[counter] <- x2
      } else {
        xIndex[counter] <- x2
        xendIndex[counter] <- x1
      }
      # 隣り合うエラーバーと線が重複しないようにoffset
      x[counter] <- xIndex[counter] + 0.1
      xend[counter] <- xendIndex[counter] - 0.1
    }

    # (2) y軸
    # 隣り合うラベル以外なら1段y軸を高くする
    # 同じ高さにするとエラー？
    counter = 0
    for (c in cmpList) {
      counter = counter + 1
      # ラベルが隣合う場合
      if (abs(xendIndex[counter] - xIndex[counter]) == 1) {
        # 2つの群の最大値を比較
        baseline <- max(maxvalues[xIndex[counter]], maxvalues[xendIndex[counter]])
        y[counter] <- baseline + 5 * counter
        yend[counter] <- baseline + 5 * counter
      }
      # 1つ飛ばしの場合
      else if (abs(xendIndex[counter] - xIndex[counter]) == 2) {
        baseline <- max(maxvalues[xIndex[counter]], maxvalues[xIndex[counter]+1], maxvalues[xIndex[counter]+2])
        y[counter] <- baseline + 5 + 5 * counter
        yend[counter] <- baseline + 5 + 5 * counter
      }
      # 両端の場合 (2つ飛ばし以上の場合)
      else {
        y[counter] <- 115
        yend[counter] <- 115
      }
    }

    # (3) 記号
    counter = 0
    for (p in pList) {
      counter = counter + 1
      if (p < 0.001) {
        annotation[counter] <- "***"
      }
      else if (p < 0.01) {
        annotation[counter] <- "**"
      }
      else if (p < 0.05) {
        annotation[counter] <- "*"
      }
      else {
        annotation[counter] <- ""
      }
    }
    print(annotation)
    
    # データフレーム作成
    #print(length(x))
    #print(length(y))
    #print(length(xend))
    #print(length(yend))
    #print(length(annotation))
    frame <- data.frame(x, y, xend, yend, annotation)
    # y軸が低い方から並べ替え (ggsignifが正しく描画されない?)
    frame <- frame %>% arrange(y, x)

  } else {
    frame <- data.frame(NULL)
  }

  return(frame)
}

### -----------------------------------------------------------
### 検定
### -----------------------------------------------------------
analyzeData <- function(dat, questN) {
  print("-----------------------------------------------------")
  print(paste("------------------------- Q", questN, " ------------------------", sep=""))
  print("-----------------------------------------------------")

  # 各群の正規性の検定
  print("shapiro-wilk test")
  i = 1
  isNormality = TRUE
  for (label in labelOrder) {
    shapiro.res <- shapiro.test(x=dat$score[dat$condition==labelOrder[i]])
    # 正規性検定の帰無仮説が棄却された場合
    if (shapiro.res$p < alpha) {
      isNormality = FALSE
    }
    print(paste0(labelOrder[i], ": ", shapiro.res$p))
    i = i + 1
  }
  print(paste0("正規性: ",isNormality))

  #print("--------------------------------------")
  #print("[method1] dunn+bonferroni(1)")
  #dunn.res <- dunn.test(dat$score, dat$condition, method="bonferroni")
  #print(dunn.res$P.adjust)
  #print(dunn.res$comparison)
  

  #print("--------------------------------------")
  #print("dunn+bonferroni(2)")
  #dunn.res1 <- dunn.test(data.q1$score, data.q1$condition, method="holm")
  #wilcox.res <- wilcox.test(data.q1$score, data.q1$condition, paired=T, p.adjust.method="bonferroni")
  #wilcox.res <- wilcox.test(x=data.q1$score, g=data.q1$condition, paired=T, p.adjust.method="bonferroni")

  print("--------------------------------------")
  print("[method2] wilcox+BH")
  #wilcox.res <- wilcox.exact(data.q1$score[data.q1$condition=="None"],data.q1$score[data.q1$condition=="Mass"], paired=T, )
  wilcox_bonfe.res <- pairwise.wilcox.test(
    x=dat$score, g=dat$condition, paired=T, p.adjust.method="hochberg")
  print(wilcox_bonfe.res)
  summary(wilcox_bonfe.res)
  # 調整後の p 値を抽出
  adjusted_p_values <- wilcox_bonfe.res$p.value
  # p値の結果の表示
  print(adjusted_p_values)
  #count <- 1
  #for (p in adjusted_p_values) {
  #  print(adjusted_p_values[count])
  #  count <- count + 1
  #}

  print("--------------------------------------")
  print("[method3] wilcox+holm")
  wilcox.res <- pairwise.wilcox.test(
    x=dat$score, 
    g=dat$condition, 
    paired=T, 
    p.adjust.method="holm"
  )
  print(wilcox.res)
  #summary(wilcox.res)

  return(wilcox_bonfe.res)
}

### -----------------------------------------------------------
### グラフ描画準備関数
### -----------------------------------------------------------
preGraph <- function(dat, res) {
  # 各群の最大値
  maxvalues <- c()
  index = 0
  for (lbl in labelOrder) {
    index = index + 1
    maxvalues[index] <- max(dat$score[dat$condition==lbl])
  }
  print(maxvalues)
  # p値
  # (1)dunn.testの場合
  #pvalues <- res$P.adjusted
  # (2) wilcoxの場合
  pvalues <- array(dim=c(1,6))
  pvalues[1] <- res$p.value[5]
  pvalues[2] <- res$p.value[6]
  pvalues[3] <- res$p.value[9]
  pvalues[4] <- res$p.value[1]
  pvalues[5] <- res$p.value[2]
  pvalues[6] <- res$p.value[3]
  
  # 比較群
  # (1)dunn.testの場合
  #comparisons <- res$comparisons
  # (2) wilcoxの場合
  # 調整後の p 値を抽出
  adjusted_p_values <- res$p.value
  # 比較の組み合わせを取得する
  comparisons1 <- as.data.frame(as.table(adjusted_p_values))
  #print(comparisons1)  # 比較ペアと p 値をデータフレームとして表示
  # Var1 と Var2 を " - " で結合して 'comparison' 列を作成
  comparisons1$comparison <- paste(comparisons1$Var1, comparisons1$Var2, sep = " - ")
  # 不要な列を削除 (任意)
  comparisons1 <- comparisons1[, c("comparison", "Freq")]
  # 列名を分かりやすく変更
  colnames(comparisons1) <- c("comparison", "Adjusted_p_value")
  comparisons <- array(dim=c(1,6))
  comparisons[1] <- comparisons1$comparison[5]
  comparisons[2] <- comparisons1$comparison[6]
  comparisons[3] <- comparisons1$comparison[9]
  comparisons[4] <- comparisons1$comparison[1]
  comparisons[5] <- comparisons1$comparison[2]
  comparisons[6] <- comparisons1$comparison[3]

  # エラーバーの定義作成
  plotPoints <- makePValueBar(pvalues, comparisons, maxvalues)
  print(plotPoints)
  return(plotPoints)
}

# TODO: フォルダ内を対称のcsvと本スクリプトのみにして全部のcsvに同一処理
filename <- c(
  #"Processed_Q1_Data.csv",
  #"Processed_Q2_Data.csv",
  #"Processed_Q3_Data.csv",
  #"Processed_Q4_Data.csv",
  #"Processed_Q5_Data.csv",
  "Processed_Task_Time_Data.csv"
)

for (i in 1:length(filename)) {
  ### 前処理
  # CSVデータ読み込み
  dat <- read.csv(filename[i], header=T)
  # "id"列削除
  dat <- dat[, colnames(dat)!= "id"]
  # ラベル順序入れ替え(for graph)
  dat <- transform(dat, condition=factor(condition, levels=labelOrder))
  # 検定
  res <- analyzeData(dat, i)
  # グラフ描画準備
  plotPoints <- preGraph(dat, res)
  # グラフ描画
  drawGraph(dat, plotPoints, i)
}

4. 解説

有意差バー以外の箇所も含め時間があるときに解説を付けます．

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?