日々のつれづれ

不惑をむかえ戸惑いを隠せない男性の独り言

par(new=TRUE)と軸の話

Rで2軸PLOT(y軸の左右に値を表記)する方法 - My Life as a Mock Quantで良く分からないコメントをしてしまった…反省。

いつまでたっても、相手の話していることの本質を理解できない。
そんな自分がいつも嫌になる。

で、ちょっと頭を整理してみた。

  • xlimが共通でylimが違うプロット。

スケールの違う二変量の動態をxlimの変遷とともに比較するときに便利。

> x <- seq(0,13,0.01)
> 
> jpeg("par_new1.jpg")
> plot(x, dexp(x), type = "l", ylab = "", col = 'red')
> par(new =T)
> plot(x, dnorm(x -10), type = "l", axes = FALSE, ylab = "", col = 'blue')
> axis(4)
> dev.off()

で、こうなる。

  • で、ここで使うpar(new=TRUE)は前のプロットの上に新しいグラフを上書きする。

この場合、

    • plot(x, dexp(x), type = "l", ylab = "", col = 'red')
    • plot(x, dnorm(x -10), type = "l", axes = FALSE, ylab = "", col = 'blue')

のxlimとylimは別物になる。でも、x軸の引数が同じ値を使っているので問題ない。

それは、こうすれば確認できる。

> x <- seq(0,13,by=.01)
> y <- seq(10,13,by=.01)
> 
> jpeg("par_new2.jpg")
> plot(x,dexp(x),type = "l",ann=FALSE,col=2)
> par(new =T)
> plot(y,dexp(y),type = "l",ann=FALSE,col=4,axes=FALSE)
> axis(1,line=3)
> axis(4)
> legend("topright",legend=c("seq(0,13,by=.01)","seq(10,13,by=.01)")
+ ,col=1:2,pch=NA,lty=1)
> dev.off()

で、こうなる。

  • で、xlimとylimの両方を共有できるときはlines関数などを利用する。
> x <- density(rnorm(1000,mean=5,sd=3))
> y <- density(rnorm(1000,mean=2,sd=2))
> 
> jpeg("par_new3.jpg")
> plot(x,type = "l",ann=FALSE,col=2,xlim=range(x$x,y$y),ylim=range(x$y,y$y))
> lines(y,col=4)
> legend("topright",legend=c("rnorm(1000,mean=5,sd=3)","rnorm(1000,mean=2,sd=2)")
+ ,col=1:2,pch=NA,lty=1)
> dev.off()

で、こうなる。

  • じゃあ、グラフを重ねてみたいときに、xlimもylimも違うシチュエーションってあったっけ?

というのが、今日の反省。
これまで何してきたっけ??

例えば、同じサンプルを別の測定法で測定することがあったとする。
できれば同じデータセットにして解析したい。
本来同じものを測定しているのだから、同じ値が出るはず。でもそうならない。

そもそも測定系が違うのだからダイナミックレンジが違う。だから当たり前。
でも、別の測定法のデータだから一緒に解析できないのはもったいない。

そういうとき、ブリッジングデータをとって両者の計量性を比較したりする。
これにはデータの関係は単純にplot関数を使ったりする。

仮にyeastデータセットのeightyカラムとninetyカラムが同じサンプルを別の測定系から得たデータとしてみる。

> library(som)
> data(yeast)
> dat <- log2(subset(yeast,subset=eighty>0&ninety>0,select=c(eighty,ninety)))
> 
> jpeg("par_new4.jpg")
> plot(dat, pch=20)
> dev.off()

散布図は概要を知るには良い方法と思います。
でも、yeastデータセットはデータ点が多すぎて分布が良く分かりません。
そういうときはhist関数やdensity関数を使ったりする。

> jpeg("par_new5.jpg")
> layout(1:2)
> hist(dat$eighty,col=rgb(r=0,g=0,b=1,alpha=.5),xlim=c(0,16))
> hist(dat$ninety,col=rgb(r=1,g=0,b=0,alpha=.5),add=TRUE)
> 
> x <- density(dat$eighty)
> y <- density(dat$ninety,bw=x$bw)
> plot(x,col=1,xlim=c(0,16))
> lines(y,col=2)
> dev.off()

ただ、xlimを共通にできるのはxの1とyの1が同じ、つまり単位量が同じときだけ。
xとyが別の測定系だった場合はx軸の単位は意味を持たない。
こういうときは、両方の最大と最小を揃えて、両方の測定系のデータの特徴を見たりする。

> jpeg("par_new6.jpg")
> plot(x,col=1,xlab="",main="x:black, y:red")
> par(new=TRUE)
> plot(y,col=2,axes=FALSE,ann=FALSE)
> axis(1,line=3,col=2)
> axis(4,col=2)
> dev.off()

こんな感じになって、測定系xはダイナミックレンジが測定系yよりも狭く、得に高シグナルを検出できていないと考えられると思います。

2つの測定系を比較する方法は他にもありますが、簡易的に計量値の違いを見るときにはこれぐらいで良いと思ってます。