sankey plotと声優の誕生日

sankey plotというものを見た。
RでもrCharts, googleVis, riverplotパッケージあたりでできるらしいがこれを知る前に2時間くらいかかって作った。
つなぐのはロジスティクス関数で曲線を補完する。
4つの県と6つの高校でやってみたらこんな感じ。

　
声優の誕生月と誕生日の分布でやってみたけど長くなっただけできもかった。

　
barplotしてみた結果。各日で11人くらいの誕生日がいる。
最大は12月12日の21人。

a1 <- c(6, 12, 18, 24)
a2 <- rbind(c(6, 1, 0, 3), c(0, 2, 7, 1),
            c(0, 0, 0, 10), c(0, 5, 4, 1),
            c(0, 2, 6, 2), c(0, 2, 1, 7))


rownames(a2) <- paste(apply(replicate(nrow(a2), sample(LETTERS, 3, replace=TRUE)), 2, paste, collapse=""), "高校")
colnames(a2) <- paste(apply(replicate(ncol(a2), sample(LETTERS, 3, replace=TRUE)), 2, paste, collapse=""), "県")

b <- 1/100
col_r1 <- rainbow(ncol(a2))
col_r2 <- topo.colors(nrow(a2))
col_l <- grey(0.8)

bx <- 0.1
lcex <- 2
tcex <- 2
y1 <- c(0, cumsum(colSums(a2)/sum(a2)))
y2 <- c(0, cumsum(rowSums(a2)/sum(a2)))
a2r <- sweep(a2, 1, rowSums(a2), "/")
a2c <- sweep(a2, 2, colSums(a2), "/")
delta <- seq(bx, 1-bx, length=100)

library(drc)
par(mar=c(1, 2, 1, 2))
plot(0:1, 0:1, type="n", frame=FALSE, axes=FALSE, xlab="", ylab="")
mtext("誕生月", 2, cex=lcex)
mtext("誕生日", 4, cex=lcex)
for(i in 1:(length(y1)-1)){
	polygon(rep(c(0, bx), each=2), c(1-y1[i], 1-y1[i+1]+b, 1-y1[i+1]+b, 1-y1[i]), border=NA, col=col_r1[i])
}
for(i in 1:(length(y2)-1)){
	polygon(rep(c(1, 1-bx), each=2), c(1-y2[i], 1-y2[i+1]+b, 1-y2[i+1]+b, 1-y2[i]), border=NA, col=col_r2[i])
}
for(i in 1:(length(y1)-1)){
	tmp_x <- c(0, cumsum(a2c[,i]))
	for(j in 1:(length(y2)-1)){
	tmp_x1 <- c(0, cumsum(a2r[j,]))
		if(a2[j, i] != 0){
			ul <- 1-y1[i]-tmp_x[j]*(diff(y1[i:(i+1)])-b)    # 左上
			ll <- 1-y1[i]-tmp_x[j+1]*(diff(y1[i:(i+1)])-b)  # 左下
			ur <- 1-y2[j]-tmp_x1[i]*(diff(y2[j:(j+1)])-b)   # 右上
			lr <- 1-y2[j]-tmp_x1[i+1]*(diff(y2[j:(j+1)])-b) # 右下
			#points(bx, ul, pch=16, cex=0.7, col=col_r2[j])
			#points(bx, ll, pch=16, cex=0.7, col=col_r2[j])
			#points(1-bx, ur, pch=16, cex=0.7, col=col_r1[i])
			#points(1-bx, lr, pch=16, cex=0.7, col=col_r1[i])
			d_u <- try(drm(c(ul, ur) ~ c(bx, 1-bx), fct=L.4()), silent=TRUE)
			if(class(d_u) == "try-error"){
				d_u <- lm(c(ul, ur) ~ c(bx, 1-bx))
				p_u <- delta*d_u$coefficients[2] + d_u$coefficients[1]
			} else {
				p_u <- predict(d_u, data.frame(delta))
			}
			d_l <- try(drm(c(ll, lr) ~ c(bx, 1-bx), fct=L.4()), silent=TRUE)
			if(class(d_l) == "try-error"){
				d_l <- lm(c(ll, lr) ~ c(bx, 1-bx))
				p_l <- delta*d_l$coefficients[2] + d_l$coefficients[1]
			} else {
				p_l <- predict(d_l, data.frame(delta))
			}
			polygon(c(delta, rev(delta)), c(p_u, rev(p_l)), col=col_l, border=NA)
		}
	}
}

for(i in 1:(length(y1)-1)){
	text(bx, mean((1-y1)[i:(i+1)]), colnames(a2)[i], pos=4, cex=tcex)
}
for(j in 1:(length(y2)-1)){
	text(1-bx, mean((1-y2)[j:(j+1)]), rownames(a2)[j], pos=2, cex=tcex)
}

# 声優誕生日
idx1 <- 31*(0:11)+1
idx2 <- 31*(1:12)
par(mar=c(3, 4, 2, 0))
ba <- barplot(c(a2), col=rep(rainbow(12), each=31), ylab="誕生日の人数")
text((ba[idx1]+ba[idx2])/2, par()$usr[3]-1, colnames(a2), xpd=TRUE, cex=2)
abline(h=median(a2[a2>0]), lty=3)