低水準関数

MikuHatsune2014-05-10

Rでプロット画面を充実させたいとき、

par(new = TRUE)

で上書きしていくスタイルと、低水準関数と呼ばれるお絵かき関数群でプロットに付け足していくスタイルがある。
ダメ出し:par(new=TRUE) は使わないという話があって、これはどういうことかと質問されたのだが、上書きを重ねるとプロット軸が合っているのになぜか線が太くなっていく。

for(i in seq(1000)){
	plot(1)
	par(new=TRUE)
}


 
複雑な図を描こうと思ったら、白いキャンバスをとりあえずプロット関数で出しておいて、低水準関数でひたすら重ねていくスタイルが汎用性があっていいと思う。たぶん流行りのggplot系はこの手法でやっているんじゃないかな(適当
エクセルみたいな、GUIでゴリゴリやるタイプのプロットなら、自分の手で微調節できるが、低水準関数でゴリ押しすると座標を厳密に取得できるので、入力が変わった時に柔軟性が高い。
観察した値と、推定した値をプロットして、推定した折れ線の面積を図示して、各々数学的説明をつけつつ、軸のスケールやラベルも変更するということを地道にやるとこうなる。

# Isotonic regression
# データ
a1 <- list(
	observation = cbind(c(99,103,97,106,88,106,100,104,39,34,6,-5,29,-3,-12),
	c(93,60,119,91,86,84,33,41,29,10,17,14,1,-18,15)), 
	modified = cbind(
	c(100,100,100,100,99,99,99,99,39,34,12,12,12,0,0),
	c(93,90,90,90,86,84,37,37,29,12,12,12,1,0,0))
)
# 濃度点
iconc <- seq(nrow(a1$observation))

par(mar=c(4,4.5,2,2))
matplot(iconc, a1$observation, type="n", xlab="", ylab="viability (%)", xaxt="n", yaxt="n", cex.lab=2)
axis(1, at=iconc, labels=NA)
mtext("concentration (log scale)", s=1, line=2.5, cex=1.8)
axis(2, at=c(0,50,100),labels=c(0,50,100), cex.axis=1.6, las=1)
polygon(c(iconc, rev(iconc)), c(a1$modified[,1], rev(a1$modified[,2])), col=paste(grey(0.8), 95, sep=""), border="white")
for(i in seq(2)){lines(iconc, a1$modified[,i], type="l", lwd=3, cex=1.2, col=1, lty=c(1, 3)[i])}
for(i in seq(2)){points(iconc, a1$observation[,i], pch=c(16,22)[i], cex=1.2, col=1)} 
tx <- c(expression(italic(f["w"](x))),expression(italic(f["m"](x))),expression(italic(g["w"](x))),expression(italic(g["m"](x))))
legend("topright", legend=tx, bty="n", lwd=3, pch=c(16,22,NA,NA), lty=c(0,0,1,3), col=1, cex=2, merge=TRUE)
text(7.8, 50, substitute(Delta*italic(S[IR])), cex=1.6, adj=c(0.5, 0.5))
text(par()$usr[1], par()$usr[4], "A", cex=3, xpd=TRUE, adj=c(2.5,0))
xi <- c(10,11)
segments(iconc[xi], a1$modified[xi, 1], y1=-100, lty=3, col=grey(0), lwd=1)
for(i in seq(2)){text(iconc[xi[i]], par()$usr[3], substitute(italic(x[y]), list(y=xi[i])), xpd=TRUE, adj=c(0.5,1.6), cex=1.6)}
# 濃度制約
text(iconc[1]-0.3, 20, substitute(italic(group("(",1,")"))), pos=4, cex=1.5)
text(iconc[2]-0.2, 17, substitute(bgroup("{", italic(atop(g[w](x[i])>=g[w](x[j]),g[m](x[i])>=g[m](x[j]))), ""), list(i=xi[1], j=xi[2])), pos=4, cex=1.5)
# 細胞株制約
text(iconc[1]-0.3, -5, substitute(italic(group("(",2,")"))), pos=4, cex=1.5)
text(iconc[2]-0.2, -8, substitute(bgroup("{", italic(atop(g[w](x[i])>=g[m](x[i]),g[w](x[j])>=g[m](x[j]))), ""), list(i=xi[1], j=xi[2])), pos=4, cex=1.5)
# 逐一
for(k in seq(2)){
	arrows(iconc[xi[k]+1]+c(0.3,0.6)[k], a1$modified[xi[k],2]+c(35,23)[k], iconc[xi[k]], a1$modified[xi[k],2], lwd=2, length=0.12, angle=12)
	if(k==1){text(iconc[xi[k]+1]+c(0.3,0.6)[k], a1$modified[xi[k],2]+c(35,23)[k],
	 substitute(italic(g[z](x[y])), list(z=c("m")[k],y=c(10,11)[k])), cex=1.2, adj=c(0.5,-0.2))
	} else {text(iconc[xi[k]+1]+c(0.3,0.5)[k], a1$modified[xi[k],2]+c(35,23)[k],
	 substitute(italic(group("",list(g[w](x[11]), g[m](x[11])), ""))), cex=1.2, adj=c(0.33,-0.2))
	}
}
arrows(iconc[xi[k]-1]+0.2, a1$modified[xi[k]-1,1]+25, iconc[xi[k]-1], a1$modified[xi[k]-1,1], lwd=2, length=0.12, angle=12)
text(iconc[xi[k]-1]+0.2, a1$modified[xi[k]-1,1]+25, substitute(italic(g[w](x[10]))), adj=c(0.5, -0.2), cex=1.2)