r-statistics-fanの日記

統計好きの現場の臨床医の覚書のようなもの

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

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

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

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

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

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

勇気を出してメールを送るも2-3日音沙汰なし。

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

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

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

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

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

計算する必要が出てから苦節10日間(メールが来て1日ww)。
やっと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)