臨床データに良く見るplot図を自作してみた
手違いで、過去の記事を消してしまいました
キャッシュから持ってきてとりあえず復元してみる
はじめに
kaz_yosさまへ
元記事を削除してしまいました。
それにともなって、貴重なご指摘も消えてしまいました。
この場を借りて、再びお礼と、お詫びをさせてください。
ここから
久しぶりに、R関連のことを書く。
ちょっと臨床系の質問をされて、調べようと思ったのものの、グラフの名前を調べるより自作した方が早そうだったので、自作した。
どんなグラフかというと、次のようなグラフで群間の測定値の散布図なんですが、ある閾値で丸めて横にドットを並べる様なグラフです。
PMDAのHPから似た図を引用します。こんなグラフです。
確かに、臨床系の論文ではよく目にするグラフですが、改めて考えると不思議なグラフです。
箱ヒゲ図で十分な気もしますし、症例数が少なくて箱ヒゲ図に値しないなら、幹葉図でも十分伝わります。
領域特有のグラフなのかな？よく分かりません。
自作する前に幹葉図を確認
で、幹葉図はstem()なので、
> dat <- matrix(sort(rnorm(200, mean=100, sd=10)), ncol=4) > head(dat) [,1] [,2] [,3] [,4] [1,] 74.02473 94.13281 100.6842 107.0629 [2,] 76.10390 94.16766 100.7503 107.1371 [3,] 77.73353 94.18114 100.8742 107.1417 [4,] 79.37207 94.39920 101.6981 107.1970 [5,] 81.01395 94.41547 101.7421 107.4019 [6,] 81.54789 94.48055 102.1720 107.4612 > > stem(dat) The decimal point is 1 digit(s) to the right of the | 7 | 4 7 | 689 8 | 12333444 8 | 5677788999 9 | 0000001111111222222223444444444444 9 | 555555666666667777778888888999999999 10 | 00000011111222223333333333344444444 10 | 55555555556666667777777777777788888889999 11 | 0000000111122223344 11 | 556667789 12 | 002 12 | 13 | 4
で、stem()の返り値ですが、これは残念ながらplot()系関数で描画できないようです。
（方法があるなら教えてください。）
> plot(stem(dat)) The decimal point is 1 digit(s) to the right of the | 7 | 4 7 | 689 8 | 12333444 8 | 5677788999 9 | 0000001111111222222223444444444444 9 | 555555666666667777778888888999999999 10 | 00000011111222223333333333344444444 10 | 55555555556666667777777777777788888889999 11 | 0000000111122223344 11 | 556667789 12 | 002 12 | 13 | 4 以下にエラー plot.window(...) : 有限な 'xlim' の値が必要です 追加情報: 警告メッセージ： 1: In min(x) : min の引数に有限な値がありません：Inf を返します 2: In max(x) : max の引数に有限な値がありません：-Inf を返します 3: In min(x) : min の引数に有限な値がありません：Inf を返します 4: In max(x) : max の引数に有限な値がありません：-Inf を返します > str(stem(dat)) The decimal point is 1 digit(s) to the right of the | 7 | 4 7 | 689 8 | 12333444 8 | 5677788999 9 | 0000001111111222222223444444444444 9 | 555555666666667777778888888999999999 10 | 00000011111222223333333333344444444 10 | 55555555556666667777777777777788888889999 11 | 0000000111122223344 11 | 556667789 12 | 002 12 | 13 | 4 NULL
多分、オブジェクト構造を理解すれば回避できると思うのですが、.Cで呼び出す関数を見るのも面倒ですし、自作した方が早いと思いました。
> stem function (x, scale = 1, width = 80, atom = 1e-08) { if (!is.numeric(x)) stop("'x' must be numeric") x <- x[!is.na(x)] n <- as.integer(length(x)) if (is.na(n)) stop("invalid length(x)") if (n == 0) stop("no non-missing values") if (scale <= 0) stop("'scale' must be positive") .C("stemleaf", as.double(x), n, as.double(scale), as.integer(width), as.double(atom), PACKAGE = "base") invisible(NULL) } <bytecode: 0x9b28610> <environment: namespace:graphics>
箱ヒゲ図を少し改変
要するに、ある桁数で丸めてしまって、その桁数に入った値を横に並べたらいいだけ。
参考にしたは、箱ヒゲ図です。
> boxplot(dat)
これに生データを重ねてみると、こうなります。
> boxplot(dat, pch=NA) > lapply(1:ncol(dat), function(i) points(x=rep(i,nrow(dat)), y=dat[,i]))
x軸が同値なので、プロットが重なって見にくい。
プロットするときは、小数点以下の違いとかは視認できないので意味がないこともある。
ちょっと工夫すると、まだ見やすい。
> boxplot(dat, pch=NA) > lapply(1:ncol(dat), function(i) points(x=rep(i,nrow(dat))+rnorm(nrow(dat), sd=.1), y=trunc(dat[,i])))
このグラフをイメージして自作すれば、目的のグラフはできそうです。
> stem2 <- function(dat, digits=1, pch=20, ...){ + digits <- -log10(digits) + boxplot(dat, border=0, ...) + lapply(1:ncol(dat), function(i){ + y <- as.data.frame(table(round(dat[,i], digits)),stringsAsFactors=FALSE) + y$Var1 <- as.numeric(y$Var1) + apply(y, 1, function(j){ + x <- seq(j[2])*1/(max(y$Freq)+2) + x <- x + i - median(x) + points(x=x, y=rep(j[1], length(x)), pch=pch) + }) + }) + }
ループが二重になって見苦しいですが、これを使うと次の様になります。
stem2(dat, digits=1, pch=20)
と自作したものの、多分、既成のグラフ関数があると思うので、知っている人がいたら教えて下さい。
久しぶりに、Rの記事を書いたので、はてな記法を忘れていて難儀しました…
