ひゃまだのblog

ひゃまだ(id:hymd3a)の趣味のブログ

Rでtukeyの多重検定を行って有意差のアルファベットをつける

(2023-07-26 初稿 - 2023-12-04 修正)

このサイトでも、以前に以下のページをアップした。

最初に紹介したサイトは、筆者としてはそこそこの需要があるようで、皆さんに見ていただいているようである。

このページは、Rを使ってtukeyの多重比較を行い、有意差がある場合にアルファベットまで付ける方法を紹介する。

なお、このページは、以下のサイトを参照させていただいた。多謝。

Rとパッケージのインストール

まずは、Rを使えるようにする。

Debianの場合

$ sudo apt install r-base

さらに、多重検定のアルファベットを付けるために、multcompパッケージをインストールする。

$ sudo apt install r-cran-multcomp

Debianで無い場合、Rを起動してから以下のコマンドを入力。

install.packages("multcomp", dependencies = TRUE)

「yes」を押下すれば、依存関係のパッケージも含めてインストールしてくれる。

データファイルの準備

表計算ソフトで表を作って、csvでダウンロードまたはエクスポート。

例では r-test.csvとして保存。

id Value Method
1 2.0 M1
2 3.0 M1
3 2.0 M1
4 6.0 M2
5 7.0 M2
6 8.0 M2
7 10.0 M3
8 11.0 M3
9 12.0 M3

Tukeyの実行と有意差のアルファベットを付ける

$ R
> library(multcomp)
> data <- read.table("r-test.csv", sep=",", header=T)
> data$Method <- factor(data$Method)
> res <- aov(Value ~ Method, data)
> summary(res)
            Df Sum Sq Mean Sq F value   Pr(>F)    
Method       2 112.89   56.44   72.57 6.26e-05 ***
Residuals    6   4.67    0.78                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> tuk<-glht(res, linfct=mcp(Method="Tukey"))
> cld(tuk, decreasing=T)
 M1  M2  M3 
"c" "b" "a" 

無事に有意差のアルファベットも出力できた。

エラーメッセージが出る場合

筆者も最初に実行したとき、以下のエラーメッセージが出た。

mcp2matrix(model, linfct = linfct) でエラー: 
  Variable(s) ‘Method’ of class ‘character’ is/are not contained as a factor in ‘model’.

ネットで検索すると、Methodを要素(factor)にしないといけないとのこと。再度、以下の行があることを確認して。

> data$Method <- factor(data$Method)

スクリプトの作成と実行

毎回、上記のとおりコマンドを打つのは大変である。

以下のサイトに、スクリプトにして実行する方法の記述があったので早速参考にしてスクリプトを作ってみた。多謝。

作成したスクリプトは以下のとおりで、tukey-label.Rと命名した。

### tukey-label.R  
# ver0.01 2023-07-26 start
# ver0.02 2023-07-26 add commandArgs

#[Usage]
# $ Rscript tukey-label.R data-file.csv

library(multcomp)
args <- commandArgs(trailingOnly=T)
data <- read.table(args[1], sep=",", header=T)
data$Method <- factor(data$Method)
res <- aov(Value ~ Method, data)
summary(res)
tuk <- glht(res, linfct=mcp(Method="Tukey"))
cld(tuk, decreasing=T)

使い方は、スクリプトの中にも書いたけど、以下のように実行する。

$ Rscript tukey-label.R r-test.csv
            Df Sum Sq Mean Sq F value   Pr(>F)    
Method       2 112.89   56.44   72.57 6.26e-05 ***
Residuals    6   4.67    0.78                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 M1  M2  M3 
"c" "b" "a" 

ちなみに、筆者は平均値が大きい順にアルファベット文字を付けるけど、小さい順に付ける場合は、以下のように変更する。

cld(tuk)

とっても簡単にできた。

こんなことなら、もっと早くからやれば良かった。

関連ページ