灘高校の現役合格率

rstan Rを使いこなす R 医学

こんなツイートを観測した。


灘（兵庫） 主要大合格者数推移 - 進学校データ名鑑

2025, 2026年でなんとなく減少傾向に見えるので全体的になんとなく減少傾向に見えている。
が、実際にどうなのか。
東大、京大、医学部、東工大（東京科学大学）、一橋大学、のしかも現役、浪人混ぜのデータのようなので、灘高現役生の合格実績がどうかデータを取ってきた。
2008-2009年の東大理Ⅲ、京大医学部の現役合格者数が抜けていたのでこちらも参照した。
2008年度 東京大学合格者数 | 受験教育情報サイト：インターエデュ・ドットコム

東工大や一橋大の現役浪人まで確認していると面倒だったので、東大、京大、医学部、の現役浪人をなんとか調べた。

stan 的に考えると、各年度y に灘高3年生がN_y 人いるとして、合格者数Y_yは、合格率\theta_y から二項分布でサンプリングされるとする。
N\sim\rm{binom}(\theta_y, N_y)
合格率\theta_y は、モデルとしては難ありだが実装が簡単なので一次線形
\theta(t)=at+b
として、aが0より有意に離れているかが経時的な合格率の増減、b が切片としてベースの合格率、となる。
モデルとして難あり、というのは、合格率は以前の年度の影響を受けるだろうが説明変数を何も考えずに年度としているからであるが、許してほしい。

理Ⅲのみ、理Ⅲと京医学部、国公立医学部（理Ⅲと京医除く）、東大と京大（理Ⅲと京医除く）、東大と京大（理Ⅲと京医除く）と国公立医学部、で分けた結果が図となる。
結果として、経時的にはどの区分も合格率は95%信用区間的には増減していない。

理Ⅲの現役合格率は概ね6%、京医も合わせると13%が現役合格している。

n <- read.table(text="
入試年 卒業数 東大合格 東大現役 理3 理3現役 京大合格 京大現役 京医合格 京医現役 医学部合計 医学部現役
2003 217 88 66 6 4 55 32 15 5 48 23
2004 214 89 56 15 8 44 22 15 9 58 28
2005 215 101 72 21 18 54 34 22 19 82 53
2006 213 80 57 13 8 51 36 25 16 70 33
2007 216 100 73 11 9 33 18 15 9 65 31
2008 216 114 82 19 14 23 12 14 8 72 35
2009 217 103 66 15 9 37 24 22 15 68 39
2010 219 103 74 21 19 35 25 23 16 93 60
2011 219 99 75 17 14 45 28 25 15 88 51
2012 219 98 81 16 14 34 25 22 15 74 40
2013 224 105 80 27 21 41 27 24 19 91 55
2014 220 104 76 12 8 30 18 16 8 64 29
2015 219 94 72 15 11 36 22 21 14 66 37
2016 218 94 72 20 17 47 26 25 13 95 42
2017 220 95 75 19 17 39 30 21 18 83 52
2018 219 92 70 15 11 42 25 22 14 97 48
2019 219 74 59 21 17 48 33 26 19 90 58
2020 220 79 57 14 11 49 31 24 15 81 48
2021 216 97 75 12 10 34 21 14 9 50 38
2022 221 92 62 10 8 48 36 20 15 86 41
2023 220 86 66 15 14 42 31 17 14 76 46
2024 218 94 71 12 12 53 38 25 14 95 35
2025 214 77 59 9 7 50 37 19 14 84 46
2026 224 95 77 7 6 47 31 13 8 44 18
", header=TRUE)

library(rstan)
library(vioplot)
library(stringr)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

code <- "
data {
  int N;
  int Y[N]; // 合格者数
  int D[N]; // 一学年
}
parameters{
  real a;
  real b;
  real<lower=0> s;
}
transformed parameters{
  real<lower=0, upper=1> theta[N];
  for(i in 1:N){
    theta[i] = a*i + b;
  }
}
model{
  a ~ normal(0, 5);
  b ~ normal(0, 5);
  s ~ exponential(1);
  for(i in 1:N){
    Y[i] ~ binomial(D[i], theta[i]);
  }
}
"

m0 <- stan_model(model_code=code)
standata <- list(N=nrow(n), D=n$卒業数,
                 Y=list("東大・京大・医学部"=n$東大現役 + n$京大現役 - 2*n$理3現役 - 2*n$京医現役 + n$医学部現役,
                        "東大・京大"=n$東大現役 + n$京大現役 - n$理3現役 - n$京医現役,
                        "理Ⅲ・京医"=n$理3現役 + n$京医現役,
                        "理Ⅲ"=n$理3現役,
                        "国公立医学部"=n$医学部現役)
                )

# 上記パターンを一気にやる
fits <- mapply(function(z) sampling(m0, c(standata[c(1,2)], Y=list(z)), iter=2000, warmup=1000, chain=20), standata$Y)
ex <- mapply(function(z) extract(z, pars=head(z@model_pars, -1)), fits, SIMPLIFY=FALSE)
#p <- mapply(function(z) 100/(1+exp(-z$theta)), ex, SIMPLIFY=FALSE)
p <- mapply(function(z) 100*z$theta, ex, SIMPLIFY=FALSE)
cols <- c("red", "pink", "green", "skyblue", "orange")
alpha <- c(0.025, 0.5, 0.975)
x <- 1:nrow(n)

#png("nada.png", 840, 480)
par(mar=c(4, 5, 2, 10), cex.lab=1.5)
plot(range(x), c(0, 5)*10, type="n", xlab="", ylab="", frame=FALSE, xaxt="n", yaxt="n")
abline(h=seq(0, 5, by=1)*10, lty=3)
for(i in seq(standata$Y)){
vioplot(p[[i]], add=TRUE, ylim=100*c(0, 0.5), las=1, frame.plot=FALSE,
        rectCol=NA, lineCol=NA, colMed=NA, col=cols[i],
        xlab="", ylab="")
}
for(i in seq(standata$Y)){
  if(i == 1){
    mtext("20XX年度", side=1, line=2.5, cex=1.5)
    mtext("現役合格率 [%]", side=2, line=3, cex=1.5)
  }
  y <- 100*standata$Y[[i]]/standata$D
  lines(x, y, col=cols[i], lwd=3)
  points(x, y, col=cols[i], pch=15)
  lines(x, y, lwd=1)
  pa <- par()$usr
  y1 <- tail(colMeans(p[[i]]), 1)
  pp <- quantile(ex[[i]]$a, alpha)*100
  pb <- quantile(ex[[i]]$b, alpha)*100
  txt <- sprintf("%s %.1f%s [%.1f, %.1f]\n%.2f%s/年 [%.2f, %.2f]",
                 names(standata$Y)[i], pb[2], "%", pb[1], pb[3],
                 pp[2], "%", pp[1], pp[3])
  txt <- sprintf("%.1f%s [%.1f, %.1f]\n%.2f%s/年 [%.2f, %.2f]",
                 pb[2], "%", pb[1], pb[3],
                 pp[2], "%", pp[1], pp[3])
  #text(pa[2], y1, txt, xpd=TRUE, pos=4)
  #legend(pa[2], y1, legend=txt, xpd=TRUE, col=cols[i], pch=15, yjust=0.5, bty="n")
  legend(pa[2], y1, legend=txt, title=names(standata$Y)[i], xpd=TRUE, col=cols[i], pch=15, yjust=0.5, bty="n", pt.cex=2, title.adj=c(NA, 0.5))
}
text(pa[2], pa[4], "現役合格率 (%) [CI]\n年次増減 (%/年) [CI]", xpd=TRUE, pos=4)
axis(1, at=seq(nrow(n)), labels=str_extract(n$入試年, "\\d{2,2}$"), lwd=0, lwd.ticks=1)
axis(2, lwd=0, lwd.ticks=1, las=2)
box()
dev.off()