こんなツイートがあった。
この画像をみたときに、出生数の月日分布なんてデータとしてあるの? と思っていたら、estat の人口動態データとしてあった。
このうち、出生数，出生年月日時・出生の場所別 というデータが、月日時刻別に集計しているのでこれを整形する。なぜか1995年がPDF だったり、csv も全角数値が入っていたり行数が年度ごとに揃っていなかったりとネ申エクセルではないけどデータ管理の下手くそさを感じさせるが、せっかくのデータを公開してくれているので文句は言わない(言うけど。
これをPython で整形して、2015年から1996年まで、1月1日から12月31日までデータを取ってくる。
# ターミナルで文字コード変換 for i in `ls *csv`; do nkf -w --overwrite $i; done # Python import os import glob import re import codecs input_wd = "/birth/" files = glob.glob(input_wd + "*csv") digit = re.compile("\d+") month = re.compile(".月") month = "月" w0 = open("birth_estat.csv", "w") for f in files: g = open(f, "rU") res = [] hoge = "" while not (month in hoge and "総" in hoge): hoge = g.readline() print hoge while "病" not in hoge: hoge = g.readline() if "日" in hoge: res += [ hoge.rstrip().split(",")[1] ] year = digit.findall(f) dat = year + res text = ",".join(dat) + "\n" w0.write(text) w0.close()
20年の間に年間出生数は年々減ってきているため、比率にして全体を揃える。
365日間の、全出生数に対する各日の出生比率だが、下向きのスパイクがいくつかある。これはツイートでも言われているように、ハッピーマンデーに寄らない祝日である。各年度は1周間単位で振動しているが、これは週末に出生数が少ないことによるもので、これを20年分重ねれば週末分は消える。これでも残る成分は、例えば2月11日だったり、4月1日だったり、年末正月、GW、お盆である。
各年度の出生比率分布をカレンダープロットすると、国民休日や、特に週末日曜日に出生比率が低いことがわかる。出生は自然に生まれるパターンと、緊急で帝王切開で生まれるパターンと、予定手術の帝王切開で生まれるパターンがあるが、近年は帝王切開が増えていることと、帝王切開を予定手術でやるならば、普通は、人手の多い平日にやろうというのが普通なので、祝日週末は避ける。
4月1日も学年の変わり目になるので普通は避けるだろう。
さて、クリスマスになると、声優監視スレがあるように、クリスマスに何が行われているかというと、ナニである。ここで、出生日と出生までの妊娠週数データがあるので、これからナニの日を推定したい。いま、妊娠週数は厚生労働省のデータを取ってきて、満N週になっているので線形に補間して日数にする。
38-39週が最も多いが、妊娠最小限界の22週から、過産期の43週までこんな感じになる。
厳密には、妊娠週数の数え方は最終月経からということになっているが、生物学的な胎児の週数は受精してからになるので、2週間のずれが生じる。ここでは、受精とナニは同時にあったとして2週間引いている。
出生月日から単純に上の週数をランダムサンプリング(ただし-2週)して、生まれた児の受精日を逆算したら、それつまりナニの日ということをやってみる。20年間の出生数データがあるので、最初と最後以外の18年分をやってみる。
結果としては、年末年始が最も多いようだった。6-8月の夏は少なめ。意外と開放感がないが、5月のGW っぽいところでヒト山あるので、連休中はアレなのだろう。
同じようにカレンダープロットを作った。出生の生データでは週末や祝日の影響があったが、妊娠週数の分布データからサンプリングするとこれは消えて、連続する日はほとんど滑らかになる。
結局のところ、妊娠38週(胎児は36週)が最も頻度が高いので、9月に出生数が多いことから単純に36週引いたら年末年始になる、という簡単な結果になった。年末年始はもちろんクリスマスを含むので、まあアレだね!
本当ならば、出生以外にも、残念ながら死産や人工中絶というものがあるので、本来の受精数の365日ベクトル というものを念頭に置いて、死産や中絶の分布を考慮して真のナニ分布というものを推定したかった。特に中絶については、この20年で40-20万件あり、出生数に対して20% 程度を占めることから、真の分布に及ぼす影響は小さくないだろうと思われたが、件数と週分布はあっても月日分布はさすがに自治体、病院レベルでも見当たらなかったので、出産によるデータのみによって推定せざるを得なかった。
あとは、出生120万/年くらいとして、図では最大と最小の差でも0.0004 くらいしかないが、これは400-500くらいの差でしかない。365日均一だとしても1/365=0.0027 なので、絶対的な差としてはそんなにない。
結論：命は愛おしい
# 妊娠週数 bw <- 22:43 bwn <- c(38, 68, 106, 130, 142, 196, 208, 259, 265, 268, 351, 444, 620, 767, 1399, 4103, 10002, 16123, 14954, 7324, 798, 29) pbw <-bwn/sum(bwn) # 妊娠日数 pbd <- rep(pbw, each=7) + rep((c(pbw[-1], 0) - pbw)/7, each=7)*(0:6) pbd <- pbd/sum(pbd) bd <- c(t(outer(bw*7, 0:6, "+"))) plot(bd, pbd, type="o", pch=16, lwd=2, xlab="Birth day", ylab="Probability") abline(v=bw*7, lty=3) text((c(bw[-1], max(bw)+1)*7+bw*7)/2, par()$usr[4], bw, xpd=TRUE, pos=3) # estat を処理したデータ dat <- read.csv("birth_estat.csv", header=FALSE, row=1) dat <- t(dat) dat <- dat[!apply(dat == "-", 1, all),] dat[dat == "-"] <- NA dat <- apply(dat, 2, as.numeric) days <- c(31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) md <- unlist(mapply(rep, 1:12, days)) dd <- unlist(mapply(seq, days)) ymd <- t(outer(colnames(dat), paste(md, dd, sep="-"), paste, sep="-")) # うるう年の処理 ymdate <- mapply(function(z) as.Date(ymd[!is.na(dat[,z]), z]), 1:ncol(ymd)) datlist <- apply(dat, 2, na.omit) cols <- rainbow(ncol(dat)) datmedian <- sweep(dat, 2, colSums(dat, na.rm=TRUE), "/") matplot(datmedian, type="l", xaxt="n", lty=1, lwd=0.6, col=cols, xlab="Month", ylab="Median birth", xlim=c(1, nrow(datmedian)+50)) lines(apply(datmedian, 1, median, na.rm=TRUE), lwd=3) month <- c(1, cumsum(days)) abline(v=month, h=1/365, lty=3) text((head(month, -1)+tail(month, -1))/2, par()$usr[3], 1:12, xpd=TRUE, pos=1) legend("right", legend=colnames(dat), col=cols, lty=1, lwd=3, bg="white", bty="o", box.col=NA) for(i in seq(ymdate)){ plotdat <- data.frame(date=ymdate[[i]], Birth=datlist[[i]]/sum(datlist[[i]])) calendarPlot(plotdat, pollutant="Birth", year=names(datlist)[i], cols="jet") } # list を vector にする alldate <- as.Date(unlist(ymdate), origin="1970-1-1") allbirth <- unlist(datlist) years <- as.numeric(unique(format(alldate, "%Y"))) # 妊娠日数からランダムサンプリングしてナニの日を推定するシミュレーションを1回行う関数 sim <- function(){ res <- vector("list", length(alldate)) for(i in length(alldate):360){ diffday <- sample(bd-2*7, size=allbirth[i], prob=pbd, replace=TRUE) res[[i]] <- alldate[i] - diffday } res1 <- as.Date(unlist(res), origin="1970-1-1") tab <- table(res1) d2 <- as.Date(names(tab)) tab <- as.numeric(tab) names(tab) <- NULL dat2 <- data.frame(date=d2, birth=tab) years <- as.numeric(unique(format(alldate, "%Y"))) # カレンダープロット用に年度ごとにまとめる dat3 <- mapply(function(y){ idx <- d2 >= as.Date(paste0(y, "-1-1")) & d2 <= as.Date(paste0(y, "-12-31")) dat22 <- dat2[idx,] }, years[2:(length(years)-1)], SIMPLIFY=FALSE) return(dat3) } # 1回あたりにかなり時間がかかる iter <- 50 result <- vector("list", iter) for(i in seq(iter)){ result[[i]] <- sim() } result_date <- mapply(function(z) z$date, result[[1]]) result_birth <- mapply(function(i) mapply(function(z) z[[i]]$birth/sum(z[[i]]$birth), result), 1:(length(years)-2), SIMPLIFY=FALSE) # うるう年の処理 tmp <- lapply(result_birth, apply, 1, median) for(i in seq(tmp)){ if(length(tmp[[i]]) < 366){ tmp[[i]] <- c(tmp[[i]][1:59], NA, tmp[[i]][60:length(tmp[[i]])]) } } for(i in seq(tmp)){ if(length(tmp[[i]]) < 366){ tmp[[i]] <- c(tmp[[i]][1:59], NA, tmp[[i]][60:length(tmp[[i]])]) } } tmp <- do.call(cbind, tmp) cols <- rainbow(ncol(tmp)) matplot(tmp, type="l", xaxt="n", lty=1, col=cols, xlab="Month", ylab="No. estimated fertilization", xlim=c(1, nrow(tmp)+50)) month <- c(1, cumsum(days)) abline(v=month, h=1/365, lty=3) text((head(month, -1)+tail(month, -1))/2, par()$usr[3], 1:12, xpd=TRUE, pos=1) legend("right", legend=years[2:(length(years)-1)], col=cols, lty=1, lwd=3, bg="white", bty="o", box.col=NA) for(i in 1:ncol(tmp)){ plotdat <- data.frame(date=result_date[[i]], Fertilization=apply(result_birth[[i]], 1, median)) calendarPlot(plotdat, pollutant="Fertilization", year=years[i+1], cols="jet") }
