Rで標準化効果量を計算する(Rewriting the code of SAS proc %stddiff to R)

どうしてもstandardized difference(以下stdiff)を計算
する必要が出た。

いろいろ調べたが、どうもSASがスタンダードになっている
みたい。

しかし、どうしても順序変数のstdiffの計算方法が分からない。
いきなりメールを送って良いものか迷いながら、知り合いとか、
一度も会ったことにない知り合いの知り合いとか、捨て身の覚悟で
失礼を承知で突撃したが、結局直接的な情報は得られなかった。

しかし、直接作者にメールしても良いんじゃね？という貴重な
アドバイスが頂けた。多分自分だけじゃ、作者にメールなどどいう
一歩は絶対に（確信）踏み出せなかったと思う。

紆余曲折は会ったが、アドバイスに従いSASコードを作成した
Yang D先生に直接メールを送ることにした。

勇気を出してメールを送るも２－３日音沙汰なし。

やはり突然のメールは失礼すぎたかと激しく凹んだ。

しかし、変な日本語名のアドレスなので迷惑メールに分類
されたかもしれないではないか！そこに一縷の望みをかける
ことにした。

今後も必要になるかと、英語の実名アドレスを新たに習得して、
再度送信した所、２０分以内（早っ）にSASコードを頂けた。
しかも、フリーで、改変もOKだよと言っていただいた。まさに神。

しかも調度よいタイミングでSAS University Editionが出ている
ではないか！
インストールに失敗し結局５ギガ以上のダウンロードを繰り返し
数日かけてやっと動くようになった。

メールが届きソースコードが手に入ったので早速、初SAS。
ググりながら、はじめて触るSASと格闘する。CSVの読み込方が
わからず、結局直接データを全部書き込むという原始的荒業を
必要としたものの、SASコードを読み解いていった。
この年で更に新たな言語を学ぶことになるとは思わなかった。
まあ、楽しいんですけどね。つか、SAS使えるなら、Rに翻訳したい
ことが山ほどありまくりんぐ。夢が広がるわ。SASも楽しいぞ。

計算する必要が出てから苦節１０日間(メールが来て１日ｗｗ)。
やっとSASと同じ結果をRで吐き出すことに成功した。

こんな苦労は他の人には味わってほしくないので、今回の成果を
ブログにて報告する。自己責任で関数は自由に使って頂いて
構わない。

なお、論文に使う際は、もとのSASコードのpaper
Yang D, Dalton JE: A unified approach to measuring the
effect size between two groups using sas. 2012.
を引用するのは必須と思う。うちのblogを書くのはacknowledgementで
十分じゃね。いずれgithubのことも勉強してインストールしやすい
ように登録したいな。

＃＃＃以下のコードは元論文とそのSASコードを参考に
http://support.sas.com/resources/papers/proceedings12/335-2012.pdf
作成したものである。元論文の記述と実際のSASコードに違いがあるので、
今回は元論文ではなくSASコード(proc %stddiff)の方を正しいとして実装している。
なぜなら、そっちの方が既報のJAMAとかの論文の数値と一致するからである。
したがって、元論文の記載とは一部異なっていることをご了承いただきたい。
重み付けは各自好きに改変したら良い。
なお、私は他の関数で重み付けmean/sdや割合を計算しておいてから、こいつに
放り込む使い方で対処している。

##標準化効果量 standardized difference
###Yang D, Dalton JE: A unified approach to measuring the 
###effect size between two groups using sas. 2012.
###Rewriting the code of SAS proc %stddiff to R
###by @rs_fun_jp October 14, 2014, Japan Standard Time

##standardized difference for categorical variables

p1 <- c(12/200, 88/200,100/200)  #group treat:   rate of event by categoly. sum of p1 = 1.00
p0 <- c(5/181,133/181,43/181)    #group control: rate of event by categoly. sum of p0 = 1.00


stddiff.cat <- function(p1, p0){   #standardized difference for categorical variables
S2 <- matrix(0, nrow = (length(p1) - 1), ncol = (length(p1) - 1))

for (i in seq_len(length(p0) - 1)){
      for (j in seq_len(length(p0) - 1)){
            if (i == j){
                  S2[i,j] <- (p1[i]*(1-p1[i]) + p0[i] * (1-p0[i])) / 2
            }else{
                  S2[i,j] <- - (p1[i]*p1[j] + p0[i]*p0[j]) / 2  #this line is different from original paper, but same as SAS proc %stddiff
            }}
}

sqrt(t(p1[-length(p1)] - p0[-length(p0)]) %*% solve(S2) %*% (p1[-length(p1)] - p0[-length(p0)]))
}

stddiff.cat(p1, p0)


##standardized difference for continuous variables  #for ordinal data, simply use rank() and calc. mean/SD
mu1 <- 76.26 #mean treat group
sd1 <- 7.26  #SD treat group
mu0 <- 77.12 #mean control group
sd0 <- 6.90  #SD control group

stddiff.cont <- function(mu1, sd1, mu0, sd0){
      (mu1 - mu0)/sqrt((sd1^2 + sd0^2)/2)
}

stddiff.cont(mu1, sd1, mu0, sd0)

####standardized difference for ordinal variables
group.t <- c(rep(1, 12), rep(2,88),rep(3,100))  #ordinal variable, group treat
group.c <- c(rep(1,5),rep(2,133),rep(3,43))     #ordinal variable, group control

stddiff.ordinal <- function(group.t, group.c){
temp <- c(group.t, group.c)
g = c(rep(1, length(group.t)), rep(2, length(group.c)))
r <- rank(temp, ties.method = "average")
data <- data.frame(temp, r, g)
data2 <- split(data, data$g)
(mean(data2[[1]]$r) - mean(data2[[2]]$r))/sqrt((var(data2[[1]]$r) + var(data2[[2]]$r))/2)
}

stddiff.ordinal(group.t, group.c)