【R】Steel-Dwass多重比較検定の結果にアルファベットの添え字をつける

R

Rでデータ解析を行う際、複数のグループ間の比較を行うことがよくあります。ここでは、Tukey-Kramer検定のノンパラメトリック版であるSteel-Dwass検定を用いて多重比較を行い、その結果をアルファベットの添え字形式で表示するカスタム関数 steel_dwass_test_cld を制作したので紹介します。

この関数は、NSM3pSDCFlig()関数をベースにSteel-Dwass検定を実行し、その結果からアルファベットの添え字を付けるように設計しています。ソースコードとして読み込むことで関数として使えるようになっています。

関数の詳細

関数と引数オプション

steel_dwass_test_cld(x, y, method = NA, n.mc = 10000, alphabets = letters, decreasing = TRUE, alpha = 0.05)
  • x: 数値ベクトル。比較対象のデータ。
  • y: 因子ベクトル。グループラベル。
  • method: Steel-Dwass検定の方法。"Exact", "Monte Carlo", "Asymptotic"のいずれか選択可能。
  • n.mc: モンテカルロ法を使用する場合のシミュレーション数。
  • alphabets: 結果表示用の文字。LETTERSで大文字、lettersで小文字を選択可能。
  • decreasing: 平均値を降順に並べるかどうかの論理値。
  • alpha: 有意水準。デフォルトは0.05

この関数では、以下のパッケージを使用しています

  • NSM3: 非パラメトリック多重比較検定のパッケージ。
  • tidyverse: データ操作や可視化のためのパッケージ群。
  • サンプルサイズが大きい場合にmethod = "Exact"を選択すると、警告されることがあります。
  • method = "Monte Carlo"を選択して実行すると、計算に時間がかかることがあります。

使用例と結果の出力

以下のコードは、steel_dwass_test_cld 関数を使ってSteel-Dwass検定を実行するサンプルです。

# サンプルデータ
set.seed(123)
group1 <- rnorm(20, mean = 1, sd = 2)
group2 <- rnorm(20, mean = 2, sd = 2)
group3 <- rnorm(20, mean = 3, sd = 2)
group4 <- rnorm(20, mean = 4, sd = 2)
group5 <- rnorm(20, mean = 5, sd = 2)

x = c(group1, group2, group3, group4, group5)
y = factor(rep(c("group1", "group2", "group3", "group4", "group5"), each = 20))

# Steel-Dwass検定の実行
steel_dwass_test_cld(x, y, method = "Asymptotic", n.mc = 10000, alphabets = letters, decreasing = TRUE, alpha = 0.05)
# 結果
group1 group2 group3 group4 group5
 "d"   "cd"   "bc"    "b"    "a" 

steel_dwass_test_cld()関数のソースコード

コード内容

steel_dwass_test_cld <- function(x,y,
                            method = NA,
                            n.mc = 10000,
                            alphabets = letters,
                            decreasing = TRUE,
                            alpha = 0.05) {
  
  # パッケージをインストールし、読み込む
  suppressWarnings({
    packages <- c("NSM3",
                  "tidyverse")
    package.check <- lapply(
      packages,
      FUN = function(x) {
        if (!require(x, character.only = TRUE)) {
          install.packages(x, dependencies = TRUE)
          library(x, character.only = TRUE)
        }
      }
    )
  })
  # データをdataframeにまとめる
  data <- data.frame(value = x, group = y)
  # Steel-Dwass検定の実行
  steel_dwass <- pSDCFlig(data$value, data$group, method = method, n.mc = n.mc)
  multcomp.df <- data.frame(labels = steel_dwass$labels,
                            p.val = steel_dwass$p.val)
  if (decreasing == TRUE) {
    #平均値の計算と並べ替え
    summarise <- data[1] %>% 
      group_by(data[2]) %>% 
      summarise_all(list(mean = mean)) %>% 
      arrange(-mean)
  } else {
    #平均値の計算と並べ替え
    summarise <- data[1] %>% 
      group_by(data[2]) %>% 
      summarise_all(list(mean = mean)) %>% 
      arrange(mean)
  }
  # p値表
  p.table <- matrix(nrow = length(attr(y, "levels")),
                    ncol = length(attr(y, "levels")))
  colnames(p.table) <- summarise[[1]]
  rownames(p.table) <- summarise[[1]]
  # 同じ列と行が一致する部分に空文字を入れる
  diag(p.table) <- ""
  # 右上半分を"-"に設定
  for (i in 1:nrow(p.table)) {
    for (j in 1:ncol(p.table)) {
      if (i < j) {
        p.table[i, j] <- "-"
      }
    }
  }
  # 特定のグループペアに対してp値または"*"を取得する関数
  get_p_val <- function(df, group1, group2) {
    # ラベルの組み合わせを作成
    label1 <- group1
    label2 <- group2
    
    # データフレーム内でラベルを検索
    p_val <- df$p.val[df$labels == label1 | df$labels == label2]
    
    # p値がalpha以下の場合は"*"を返し、それ以外は空文字を返す
    if (length(p_val) == 0) {
      return(NA)
    } else if (p_val <= alpha) {
      return("*")
    } else {
      return("")
    }
  }
  # NAのセルに記号を入れる
  for (i in 1:nrow(p.table)) {
    for (j in 1:ncol(p.table)) {
      if (p.table[i,j] %>% is.na() == TRUE) {
        p.table[i,j] <- get_p_val(
          multcomp.df,
          group1 = paste(rownames(p.table)[i],"-",colnames(p.table)[j]),
          group2 = paste(rownames(p.table)[j],"-",colnames(p.table)[i]))
      }
    }
  }
  # 各列に含まれる"*"の個数を数える
  star_counts <- apply(p.table, 2, function(x) sum(x == "*"))
  # アルファベットを追加する関数
  add_alphabets <- function(column, current_alphabet) {
    for (i in seq_along(column)) {
      if (column[i] == "") {
        column[i] <- current_alphabet
      }
    }
    return(column)
  }
  # 各列にアルファベットを追加
  new_matrix_data <- p.table
  previous_count <- NULL
  alphabet_index <- 1
  
  for (j in seq_along(colnames(p.table))) {
    current_count <- star_counts[j]
    if (is.null(previous_count) || current_count != previous_count) {
      current_alphabet <- alphabets[alphabet_index]
      alphabet_index <- alphabet_index + 1
    }
    new_matrix_data[, j] <- add_alphabets(p.table[, j], current_alphabet)
    previous_count <- current_count
  }
  # 行ごとに含まれる文字を集計
  cld <- apply(new_matrix_data, 1, function(row) {
    # 行内の"-"と"*"を除いた文字をユニークに集めてソート
    unique_chars <- unique(row[row != "-" & row != "*"])
    paste(sort(unique_chars), collapse = "")
  })
  return(cld[attr(y, "levels")])
}

ファイルのダウンロード

steel_dwass_test_cld()関数のソースコードはこちらからダウンロードしていただけます。

.txtファイルをsource()関数で読み込むとsteel_dwass_test_cld()関数が読み込まれます。


本記事では、RでSteel-Dwass検定の結果にアルファベットの添え字を付けるコードを紹介しました。

RでSteel-Dwass検定の結果に、アルファベットの添え字を付ける方法をネット上で探しても見つけられなかったので関数を自作しましたが、もっと簡単な方法などありましたらコメントやお問い合わせなどで教えていただけると嬉しいです。

この記事がお役に立てたら幸いです。

プロフィール
この記事を書いた人

農学の博士前期課程を修了した研究者。
植物生理と環境調節をテーマに研究しています。

Masaをフォローする
R統計
Masaをフォローする
タイトルとURLをコピーしました