mikuhatsune.hatenadiary.com
叡王戦で破れたため八冠の独占状態ではなくなった。
253日だったがそもそも八冠になることがすごい。
過去記事のスクリプトは一部修正した。
library(survival) library(prodlim) txt <- " class,title,start,end 出生,ama,2002-07-19,2012-09-21 六級,ama,2012-09-22,2012-11-09 五級,ama,2012-11-10,2013-05-04 四級,ama,2013-05-05,2013-06-01 三級,ama,2013-06-02,2013-09-14 二級,ama,2013-09-15,2014-04-04 一級,ama,2014-04-05,2014-06-20 一段,dan,2014-06-21,2015-02-27 二段,dan,2015-02-28,2015-10-17 三段,dan,2015-10-18,2016-09-30 四段,dan,2016-10-01,2018-01-31 五段,dan,2018-02-01,2018-02-16 六段,dan,2018-02-17,2018-05-17 七段,dan,2018-05-18,2020-08-19 八段,dan,2020-08-20,2021-07-02 九段,dan,2021-07-03, 棋聖,title,2020-07-16, 王位,title,2020-08-20, 叡王,title,2021-09-13,2024-06-20 竜王,title,2021-11-13, 王将,title,2022-02-12, 棋王,title,2023-03-19, 名人,title,2023-06-01, 竜王・名人,title,2023-06-01, 王座,title,2023-10-11, 二冠,kan,2020-08-20,2021-09-12 三冠,kan,2021-09-13,2021-11-12 四冠,kan,2021-11-13,2022-02-11 五冠,kan,2022-02-12,2023-03-18 六冠,kan,2023-03-19,2023-05-31 七冠,kan,2023-06-01,2023-10-10 八冠,kan,2023-10-11,2024-06-20 " sysd <- Sys.Date() # sysd <- as.Date("2023-12-07") dat <- read.csv(text=txt, stringsAsFactors=FALSE) dat$end <- ifelse(dat$end == "", as.character(sysd), dat$end) dat$days <- as.numeric(as.Date(dat$end) - as.Date(dat$start)) dat$event <- (dat$end != sysd) + 0 fit <- survfit(Surv(days, event=event, type="right") ~ 1, data=dat) sfit <- prodlim(Hist(days, event=event) ~ 1, data=dat) ss <- summary(sfit) s1 <- as.data.frame(ss)[ss$n.event == 0, ] s2 <- merge(subset(dat[dat$days%in%s1[, "time"],], title=="title"), s1, by.x="days", by.y="time") # MST の信頼区間の計算 pai <- 0.5 g <- switch((is.null(fit$strata))+1, unlist(mapply(rep, names(fit$strata), fit$strata)), rep(1, length(fit$n.event))) er <- tapply(fit$surv*fit$std.err, g, c) p <- tapply(fit$surv, g, c) t0 <- tapply(fit$time, g, c) i0 <- lapply(p, function(z) tail(which(z > pai), 1)) i1 <- lapply(p, function(z) head(which(z < pai), 1)) j0 <- mapply(function(z0, z1) z0[z1], t0, i0, SIMPLIFY=FALSE) j1 <- mapply(function(z0, z1) z0[z1], t0, i1, SIMPLIFY=FALSE) mst <- mapply(function(z0, z1, z2, z3){z0[z1]+(z0[z2]-z0[z1])*(z3[z1]-pai)/(z3[z1]-z3[z2])}, t0, i0, i1, p) cis <- do.call(rbind, list( mapply(function(z0, z1, z2, z3){z0[tail(which(z1+2*((z0>z2)-0.5)*qnorm(0.975)*z3 > pai & z0 < z2), 1)+1]}, t0, p, mst, er), mst, mapply(function(z0, z1, z2, z3){z0[head(which(z1+2*((z0>z2)-0.5)*qnorm(0.975)*z3 < pai & z0 > z2), 1)-1]}, t0, p, mst, er) ) ) half_time <- sfit$time[which(sfit$surv < 0.5)[1]] xl <- c(0, 365*5) axis1.at <- at.t <- floor(seq(0, ceiling(max(xl/365)), by=0.5)*365) axis1.labels <- sprintf(ifelse(at.t%%365 == 0, "%.f", "%.1f"), at.t/365) par(mar=c(4, 4, 1, 1), las=1, cex.lab=1.5) plot(sfit, lwd=3, xlim=xl, xlab="タイトル保持期間(年)", ylab="", axis1.at=axis1.at, axis1.labels=axis1.labels, atrisk.at=at.t, marktime=TRUE, background=TRUE, background.horizontal=NA, atrisk.labels="藤井棋士", confint.col="yellow") mtext("生存確率", 2, line=2, las=3, cex=1.5) p <- par()$usr # segments(half_time, 0.5, y1=p[3], lty=3) # segments(half_time, 0.5, x1=-100, lty=3) # text(half_time, p[3], sprintf("%d日", half_time), cex=2, adj=c(-0.1, -0.5)) arrows(cis[2, 1], 0.5, cis[1, 1], angle=90, length=0.1, col=4, lwd=3) arrows(cis[2, 1], 0.5, cis[3, 1], angle=90, length=0.1, col=4, lwd=3) segments(cis[2, 1], 0.5, y1=p[3], lty=3) txt <- sprintf("%.f日 [95%s CI: %.f ー %.f]", cis[2, 1], "%", cis[1, 1], cis[3, 1]) text(half_time, p[3], txt, cex=2, adj=c(0, -0.5)) txt <- sprintf("%s年%d月%d日時点", format(sysd, "%Y"), as.numeric(format(sysd, "%m")), as.numeric(format(sysd, "%d"))) legend("topright", legend=txt, bty="n", cex=2) # 生存確率が同じタイトルを同一x軸上にtext する xd <- 100 r <- rle(s2[, "surv"]) xx <- unlist(mapply(rep, s2[cumsum(r$lengths), "days"], r$lengths)) dup <- duplicated(s2[, "days"]) | duplicated(s2[, "days"], fromLast=TRUE) dif <- c(xd, diff(s2[, "days"])) tcex <- 1.2 for(i in seq(nrow(s2))){ x0 <- s2[i, "days"] y0 <- s2[i, "surv"] x1 <- xx[i] + xd if(i != 0){ # 王位と棋聖のように近いタイトルを上下離す flag <- 2*((!(dif < xd/2 & dif > 0))-0.5)[i] y1 <- y0 + (p[4]-p[3])*(x1-x0)/(p[2]-p[1])*flag } else { y1 <- y0 } if(!dup[i]){ txt <- sprintf("%s (%d日)", s2[i, "class"], s2[i, "days"]) segments(x0, y0, x1, y1, lty=3, lwd=2) text(x1, y1, txt, pos=4, cex=tcex) } else { if(duplicated(s2[, "days"])[i]){ # 竜王・名人がかぶる txt <- sprintf("%s (%d日)", paste(s2[dup, "class"], collapse="/"), s2[i, "days"]) segments(x0, y0, x1, y1, lty=3, lwd=2) text(x1, y1, txt, pos=4, cex=tcex) } } } d0 <- list(as.POSIXct(dat$start), as.POSIXct(dat$end)) interval <- list(seq(as.POSIXlt(format(min(d0[[1]]), "%Y-01-01")), max(d0[[2]]), by="year"), seq(as.POSIXlt(format(min(d0[[1]]), "%Y-01-01")), max(d0[[2]]), by="month") ) interval[[2]] <- setdiff(interval[[2]], interval[[1]]) y <- 1:nrow(dat) xl <- as.numeric(c(as.POSIXlt("2012-04-01"), max(d0[[2]]))) yl <- c(1, nrow(dat)) cols <- c(ama="red", dan="yellow", title="cyan", kan="orange") par(mar=c(2.5, 6, 2.5, 1.5)) plot(d0[[1]], y, type="n", xlim=xl, ylim=yl, xaxt="n", yaxt="n", xlab="", ylab="", bty="[") axis(2, at=y, labels=dat$class, las=1) axis(1, at=interval[[1]], labels=format(interval[[1]], "%Y")) axis(3, at=interval[[1]], labels=format(interval[[1]], "%Y")) abline(v=interval[[1]]) # abline(v=interval[[2]], lty=3) # 毎月 for(i in seq(nrow(dat))){ #segments(x0=d0[[1]][i], y0=i, x1=d0[[2]][i], lwd=20) px <- c(d0[[1]][i], d0[[2]][i])[c(1,2,2,1)] py <- i + c(-1, -1, 1, 1)*0.5 polygon(px, py, col=cols[dat$title[i]], border=NA) text(px[2], i, dat$day[i], pos=4, xpd=TRUE) } abline(h=c(0, y)+0.5, lty=3)