東京大学での式辞の統計問題を考える

こんな話を見かけた。ちなみに式辞は読んでない。
qiita.com

データ取得まではやってくれていて、解析もしていた。解析者では、東京大学鳥取大学島根大学の医学部での男女合格率には差がない、ということだった。
せっかく81大学の男女受験者数と合格者数を平成30年から25年まで(2018年から2013年まで)まとめてくれているので、もうちょっと解析したい。かつてこんなことをしたので考えてみよう。
mikuhatsune.hatenadiary.com

男女の受験者数と、合格者数が出ているので、単純には分割表検定のようなことができるが、これが81大学(と6年分)あるので、メタアナリシスのようなことをやってみよう。
受験者が男であることを対照として、受験者が女であることを介入というか治療群ということにみなしてデータを整形する。ここではメタアナリシス自体の説明はしない。

結論からいうと、女性であることは合格率が15%下がる(OR=0.85 [0.80, 0.90] でp値は10^{-7} くらいでる)。
meta パッケージの都合で字が非常に小さいが、Odds Ratio のみぎ側が女性有利、ひだり側が男性有利である。OR=1 の垂線は、男女で差がないことを示している。
Experimental とControl がそれぞれ女性、男性を示している。Event が合格者数である。男女がごちゃごちゃになってわかりにくい場合は、下3/4 くらいにある東京女子医科大学を見てもらうと、Control が0で男性であることがわかるし、その2つ上の東京医科大学を見てもらうと、forest plot が非常にひだり(つまり、男性有利)に偏っているのでそれが目印である。
f:id:MikuHatsune:20190415002531p:plain
各々の大学の受験者数は男性300、女性50くらいなので、単一の大学を見ても有意差が出るような合格率の差は生まれないが、81大学も集めると男性7.6万、女性4.8万人なので、流石に小さな差でも有意差は出る。各大学を見ると横線が長いが、Fixed effect model もしくはRandom effects model のところはほぼ丸で横線はほとんどないことがわかる。これがメタアナリシスのサンプルサイズの暴力である。

年度ごとにやっても0.1から0.15程度のORの低下があるので、女性が医学部に合格しにくいのはたぶん本当なんだろう(適当

te
     2018      2017      2016      2015      2014      2013 
0.8256398 0.9172780 0.8441605 0.8861341 0.8508194 0.8718301 


さて、東京大学での話なので、東京大学だけ取り出してみると、東京大学では男女で合否は差はないような感じである。いや、メタアナリシスで冒頭に全大学でやっているのにサブグループでしていいの?という疑問はある。
f:id:MikuHatsune:20190415002603p:plain

メタアナリシスは本当はfunnel plotとかデータの偏りとか検討しないといけないが、今回は省いた。

library(meta)
dat <- read.csv("data.csv")

gokei <- dat[,4] == "合計"
d0 <- dat[gokei, c(3,5,6,8,9)]

d1 <- cbind.data.frame(d0[,1], mapply(function(z) as.numeric(gsub(",", "", as.character(d0[,z]))), 2:5))
colnames(d1) <- c("Univ", "n.ctrl", "n.trt", "col.ctrl", "col.trt")
u <- as.character(d1$Univ)

b1 <- metabin(col.trt, n.trt, col.ctrl, n.ctrl, data=d1, sm="OR", studlab=Univ, comb.fixed=TRUE, comb.random=TRUE, label.e="女性", label.c="男性")

# png("result.png", 720, 1300)
forest.meta(b1, fontsize=12)
dev.off()
funnel(b1)
metabias(b1)

idx <- mapply(function(z) c(0,1,3,4) + z, grep("年度", colnames(dat)), SIMPLIFY=FALSE)
d2 <- mapply(function(w) mapply(function(z) as.numeric(gsub(",", "", as.character(dat[gokei, w][,z]))), 1:4), idx, SIMPLIFY=FALSE)

u0 <- c(mapply(function(z) sprintf("%d_%s", z, u), 2018:2013))
d3 <- cbind.data.frame(u0, do.call(rbind, d2))
colnames(d3) <- c("Univ", "n.ctrl", "n.trt", "col.ctrl", "col.trt")

b3 <- metabin(col.trt, n.trt, col.ctrl, n.ctrl, data=d3, sm="OR", studlab=Univ, comb.fixed=TRUE, comb.random=TRUE, label.e="女性", label.c="男性")

# 各年度
d4 <- mapply(function(z) cbind.data.frame(u, z), d2, SIMPLIFY=FALSE)
for(i in seq(d4)) colnames(d4[[i]]) <- c("Univ", "n.ctrl", "n.trt", "col.ctrl", "col.trt")
b4 <- mapply(function(z) metabin(col.trt, n.trt, col.ctrl, n.ctrl, data=z, sm="OR", studlab=Univ, comb.fixed=TRUE, comb.random=TRUE, label.e="女性", label.c="男性"), d4, SIMPLIFY=FALSE)
# TE
te <- 1 + mapply(function(z) z$TE.fixed, b4)

# 東京大学
d5 <- do.call(rbind.data.frame, mapply(function(z) z[u == "東京大学",], d4, SIMPLIFY=FALSE))
d5$Univ <- paste0(2018:2013, "_", as.character(d5$Univ))
b5 <- metabin(col.trt, n.trt, col.ctrl, n.ctrl, data=d5, sm="OR", studlab=Univ, comb.fixed=TRUE, comb.random=TRUE, label.e="女性", label.c="男性")

# png("result1.png", 720, 220)
forest.meta(b5)
dev.off()