なぜグラフを描くのかと問われれば要約統計量を見やすい形で表現したいからであり、ヒトの目という不完全なデバイスを通してデータを認識するためにはグラフとよばれる表現手段を介したほうが効率が良いからである。

と考えるならば、分析のためにグラフを描くという行為は統計解析と切り離せないものであり、画像を描くためだけのライブラリを選択することはあまりよろしいことではないと思う。

そして望ましくは、適切に前処理されたマトリックスから望むグラフが作られるべきであり、グラフを描くためにデータを（ハッシュなんかに）加工するというのはあまりいいやり方とは言えないんじゃないかなぁと思っている。

これが、Rでグラフを描く大きな理由。データの前処理はRでやることはほとんどなくて、大体はperl,pythonを使うので、rpy2は重宝する。

import rpy2.robjects.lib.ggplot2 as ggplot2
from rpy2.robjects import r
from rpy2.robjects.packages import importr
base = importr('base')

datasets = importr('datasets')
mtcars = datasets.mtcars

pp = ggplot2.ggplot(mtcars) + \
     ggplot2.aes_string(x='wt', y='mpg', col='factor(cyl)') + \
     ggplot2.geom_point() + \
     ggplot2.geom_smooth(ggplot2.aes_string(group='cyl'), method='lm')

#pp.plot()
r['ggsave'](file="mypng.png", plot=pp)
#r['dev.off']()

これで、重量(x)と燃費(y)のプロットに対して、気筒毎に回帰直線と信頼区間を描くことができる。

分析手法をテンプレート化して自動的にレポートを出力したい時にはプログラミングできる手段があるといいのでrpy2は良いライブラリだと思う。

グラフィックスのためのRプログラミング
H. ウィッカム
丸善出版 / 4200円 ( 2012-02-29 )

富士山をプロットする関数をみかけたので、Rで描いてみた。

fx <- function(x){x^4-x^2+6}
sfx <- function(x){12/(abs(x)+1)}
gx <- function(x){1/2*cos(6*x)+7/2}

plot(gx,-2,2,xlim=c(-7,7),ylim=c(0,7),xlab="",ylab="")
par(new=T)
plot(sfx,-7,-1,xlim=c(-7,7),ylim=c(0,7),axes=F,xlab="",ylab="")
par(new=T)
plot(sfx,1,7,xlim=c(-7,7),ylim=c(0,7),axes=F,xlab="",ylab="")
par(new=T)
plot(fx,-1,1,xlim=c(-7,7),ylim=c(0,7),axes=F,xlab="",ylab="")

ggplotで関数だけ描く方法わからんなぁとつぶやいたら、教えてもらったのでggplot2でも描いた。

gx <- function(x){
x[-2>x | x>2] <-NA
1/2*cos(6*x)+7/2
}

fx <- function(x){
x[x<(-1)] <- NA
x[x>1] <- NA
x^4-x^2+6
}

fx <- function(x){
x[-1<x & x<1] <-NA
12/(abs(x)+1)
}

ggplot(data.frame(x=c(-7,7),y=0), aes(x,y))+stat_function(fun=gx) \
+stat_function(fun=fx)+stat_function(fun=sx)

ggplot2だと線がうまくつながらなかった。

Ggplot2: Elegant Graphics for Data Analysis (Use R!)
Hadley Wickham
Springer-Verlag New York Inc (C) / ￥ 5,259 ( 2009-08-30 )

昨日の化合物の類似性ネットワークのやつは何がメリットなのかなぁと考えてみたけど、一つにはedgeの数を自由に決められるあたりなんじゃないかと。あんまり低い閾値だとなんでもかんでもつながって何見てるかわからんし、逆に厳しい閾値にしちゃうとシングルトンばかりになっちゃうし。

というわけで、プロット欠かせてそれ見て閾値決める必要はあるんだろうなと。最終的にはeddgesの数とか、ノードあたりに生えているedgeの数の平均とかを統計的に処理して自動で閾値を決めるのが楽なんだろうけど、、、（一般的な方法あるのかな？）。X-meansに倣って、BICとかAIC使ってみようかなぁ。

で、プロット描いた。最初に書いたコードがナイーブすぎて遅くていらっときたので書き直した。まだかなり遅いけど、耐えられるのでこれでよしとした。

import pybel

mols = list(pybel.readfile("sdf", "pc_sample.sdf"))
fps = [x.calcfp() for x in mols] 

similarity_matrix = []
for i in range(len(fps)):
    for j in range(i+1,len(fps)):
        similarity = fps[i] | fps[j]
        similarity_matrix.append(similarity)

threshold = 0.0
print "threshold\tnodes"
while threshold < 100:
    t = threshold / 100.0
    new_matrix = [e for e in similarity_matrix if e > t]
    nodes = len(new_matrix)
    print "%2.1f\t%d" % (threshold, nodes)

    similarity_matrix = new_matrix
    threshold += 0.1

あとはRで

library(ggplot2)
nodes <- read.table("xxxx",header=T)
qplot(threshold,edgess, data=nodes, log="y")
dev.copy(device=png,filename="xxxxxx")
dev.off()

node数が1158なのでedgeの数が10000から100あたりの間でちょうどいい閾値を探したいんだけど、これだけじゃわからん。やっぱシングルトンの数とかノードあたりに生えているエッジ数の平均とか数えないとだめかな。

自分の血圧のデータが3年くらい溜まっているので、収縮期血圧と測定日時のプロットをして、スムージングもしてみた。冬場は血圧が上昇する傾向にある。測定始める前は180超えて190近くあったから、かなり落ち着いた感じ。

library(ggplot2)
data <- read.csv("/Users/kzfm/mybp.csv",header=FALSE)
time <- strptime(data$V1, "%Y-%m-%d")
high <- data$V2
low  <- data$V3
mybp <- data.frame(time=time,high=high,low=low)
p <- qplot(data=mybp,time,high)
p +  stat_smooth() + stat_smooth(colour="green", method="lm", formula=y~poly(x,9), size=1,se=FALSE)

2008年から2009年の頭のあたりで急激に落ちているのは環境要因ですな。そして、ああいう働き方は二度としないしできないだろうなぁ。