【R言語】Steel検定のコードを写経して自分なりのコードを作ってみる

こんにちは。R言語の勉強をしている生命科学系研究者のえいこです。

前回の記事で、Steel検定の良いパッケージが無いので作ってしまおう！という宣言をしました。

※Steel検定は簡単に言うと、Dunnett検定のノンパラメトリック版。（と説明されているサイトが多い）処理群と対照群のデータに順位をつけて検定する方法でMann-WhitneyのU検定に似ています。

前回は、function関数の使い方について、簡単にまとめました。

今回は、こちらのコードを写経しながら、自分の扱いやすいコードを書いていきたいと思います。

私は手を動かして覚えていくタイプなので、まずは参考にしたコードをぜんぶノートに書き出しました。

そのあと、それぞれのコードで何をやっているのかを確認しながら自分でRStudioにコードを書いて一つ一つ確認していきます。

元々ネットに載っているコードをコピペでも良いのでは？

読んでいる方は「なんで、こんな面倒なことをしているのか？」って疑問を持つのでは無いかと思います。

勉強熱心だなーって思うかもしれませんが、実際は...参考にしたコードを実行したところ、あんまりうまく動か無かったからです。

関数の使い勝手があまりよくわからなかったので、どんなことが書いてあるのか確認してみたい！というのが出発点です。

コピペでうまく動いたら、そのままやっていたのにーとも思いつつ。

初心者はコードをシンプルに！

参考にしたコードは、if文とか入り乱れていてプログラミングを始めたばっかりの私からしたらすごく複雑なコードでした。

実際は、エラー表示や関数がうまく動かないのを防ぐために色々書かれているだけ。

もっとシンプルに！わかりやすく！ということで、if文をごっそり削除して必要最低限だけのコードにしていきます。

群馬大学の青木先生が作成されたコードをベースに、p値も計算するコードを書き足していくイメージで書いていきました。

実際に書いてみよう！計算に必要な変数を定義していく

まずは、p値の計算に必要なパッケージをインストールしておきます。

install.packages(“mvtnorm”)
library(mvtnorm)

次に、関数の名前の定義していきましょう。

SteelTest <-function(x, g, control=NULL, alternative=c("two-side", "less", "greater"))

※青木先生のコードでは必然的に最初に書かれているグループがコントロールとして認識されてしまいますが、control群が指定できるコードを採用しました。

関数内の計算を定義していきます。

{
alternative <- match.arg(alternative)
 OK<- complete.cases(x, g)
 x<- x[OK]
 g <- g[OK]
 g<- factor(g)
 k <- nlevels(g)​
if (is.null(control)){
   control <- levels(g)[1]
 }​

一行目の”alternative<- match.arg(alternative)”は何気なく書かれていますが、かなり便利なので残しておきました。

alternativeはp値を出すときの計算方法を指定する引数なのですが、"two-side", "less", "greater"と3つあります。

two-sideってわざわざ書くの、めんどくさい！っていう時に、twoでもOK！tでもOK！ってしてくれる関数です。

“match.arg”は一部でも一致するものがあれば、その引数を渡すという関数。

使い方はこんな感じで使えます。

princess <- function(name=c("シンデレラ", "ラプンツェル", "ベル", "アリエル")){
 name <- match.arg(name)
 words<- paste("やっぱりディズニープリンセスといえば、",name, "だよね！")
 return(words)
}

今、娘が絶賛ディズニープリンセスにどハマり中なのでこんな関数を作ってみました。

で、match.argを使うと一文字だけでも引数を渡してくれるようになるってわけです。

二行目の”complete.cases”は欠損値（NA）が入っている行を除去する関数です。これを入れておくと「NAがあって計算できない！」ってことがないので、一応データのクオリティチェックのために入れておきます。

群をfactro型に変換して、nlevelsで群の個数を変数kに渡しています。

最後に、controlがNULLだった場合には1番目のグループをコントロールにしますよーっていうコードを書き足しています。

検定用の計算コードを書いていく（というか写経）

これで、準備はOK！ρ値の計算をしていきます。青木先生のコードをほぼコピペです。

get.rho <- function(ni) # ρを計算する
 {
   l <- length(ni)
   rho <- outer(ni, ni, function(x, y) { sqrt(x/(x+ni[1])*y/(y+ni[1])) })
   diag(rho) <- 0
   return(sum(rho[-1, -1])/(l-2)/(l-1))
 }
 
 ni <- table(g) #number of data
 a<- length(ni) #number of groups
 xc <- x[g==control] #control data
 nc<- length(xc) #number of control data
 rho<- ifelse(sum(nc==ni)==a, 0.5, get.rho(ni)) #rho value

次に、コントロール群との順位検定をしていきます。for文を回すので、最初にからのベクターを用意しておきます。

#prepare the empty vector
 vc <- c()
 vt <- c()
 vp <- c()
 vs <- c()
 lg <- levels(g)

順位検定の計算は、青木先生の計算方法に従いました。

 for (i in lg){
   if (control == i){
     next
   }
   r<-rank(c(xc, x[g==i]))
   R <- sum(r[1 : nc])
   N <- nc + ni[i]
   E <- nc*(N+1)/2
   V<- nc*ni[i]/N/(N-1)*(sum(r^2)-N*(N+1)^2/4)
   t<- abs(R-E)/sqrt(V)

Rのrank関数は値が小さいものから1位、2位と順位をつけていきます。

変数rのところで、コントロール群と処理群の2群に関して順位づけをしているのがわかりますね。

いよいよ、p値の計算をしていきます。トライフーズさんのコードをそのままお借りしました。

#calculate p value
   corr <- diag(a-1)
   corr[lower.tri(corr)]<- rho
   pmvt.lower<- Inf
   pmvt.uppwer <- Inf
   if (alternative == "less"){
     pmvt.lower <- -t
     pmvt.upper <- Inf
   } else if (alternative == "greater"){
     pmvt.lower <- t
     pmvt.uppwer <- Inf
   } else {
     t<- abs(t)
     pmvt.lower <- -t
     pmvt.uppwer <- t
   }
   p <- 1- pmvt(lower=pmvt.lower,
                upper = pmvt.uppwer,
                delta= numeric(a-1),
                df = 0, corr = corr,
                abseps = 0.0001)

最後に、有意差がぱっと見でわかるようにstarという変数を作ってみました。

if (p < 0.001){
     star <- "***"
   } else if (p< 0.01){
     star <- "**"
   } else if (p<0.05){
     star <- "*"
   } else {
     star <- ""
   }

計算結果をデータフレームで出力する

最後に、それぞれベクトルでまとめていた結果を一つのデータフレームにして、出力します。

   vc <- c(vc, c(paste(i, control, sep=":")))
   vt <- c(vt,c(t))
   vp <- c(vp,c(p))
   vs <- c(vs, c(star))
 }
 df<- data.frame(comparison=vc,
                 t.value = vt,
                 rho=rho, p.value=vp, star=vs)
 rownames(df) <- NULL
 return(df)
}

これでおしまいです。

実際に作ったデータで、関数を回してみましょう。

できました！Aとcontrolで有意に差があることが一目瞭然でわかりますね。

次回は、群馬大学の青木先生のSteel-Dwass検定のコードを少し変えて結果をデータフレームに出力するコードを書いてみようと思います。

それでは、また！