Python医療データ分析入門 1-1
やること
擬似健康保険データを作ってみる。 会員全体のデータを作るために、出生と死亡のシミュレーションを行う。
まずは出生。 年齢と性別を決めると、それぞれ何人生まれたか…は、公的統計の年次男女出生データを使えばいい。 そこから、乱数シミュレーションによりN人発生させる。
1-1では、その下準備のために、出生と死亡のデータを作成する。
IPSSが公開している、年次出生数のデータを利用する。
library(readr)
library(dplyr)
url = 'http://www.ipss.go.jp/p-toukei/JMD/00/STATS/Births.txt'
dat = read.table(url, skip=2, header = TRUE)
dat %>% head## Year Female Male Total
## 1 1947 1301806 1376986 2678792
## 2 1948 1303060 1378564 2681624
## 3 1949 1316630 1380008 2696638
## 4 1950 1134396 1203111 2337507
## 5 1951 1043048 1094641 2137689
## 6 1952 977101 1028061 2005162ある年に、FemaleとMaleが何人生まれたかのデータになっている。
性別列を作って縦に持っておく。
library(tidyr)
library(magrittr)
dat %>%
pivot_longer(cols = c("Male", "Female"), names_to = "Sex", values_to = "Life") %>%
mutate(Sex = if_else(Sex == "Female", "F", "M")) -> dat
dat## # A tibble: 140 x 4
## Year Total Sex Life
## <int> <int> <chr> <int>
## 1 1947 2678792 M 1376986
## 2 1947 2678792 F 1301806
## 3 1948 2681624 M 1378564
## 4 1948 2681624 F 1303060
## 5 1949 2696638 M 1380008
## 6 1949 2696638 F 1316630
## 7 1950 2337507 M 1203111
## 8 1950 2337507 F 1134396
## 9 1951 2137689 M 1094641
## 10 1951 2137689 F 1043048
## # … with 130 more rows乱数でシミュレーションをするために、累積比率を出しておく
dat %<>%
arrange(Sex, Year)
dat %<>%
mutate(ratio = Life / sum(Life)) %>%
mutate(cum_sum = cumsum(ratio))
dat %>% head## # A tibble: 6 x 6
## Year Total Sex Life ratio cum_sum
## <int> <int> <chr> <int> <dbl> <dbl>
## 1 1947 2678792 F 1301806 0.0121 0.0121
## 2 1948 2681624 F 1303060 0.0122 0.0243
## 3 1949 2696638 F 1316630 0.0123 0.0366
## 4 1950 2337507 F 1134396 0.0106 0.0472
## 5 1951 2137689 F 1043048 0.00973 0.0569
## 6 1952 2005162 F 977101 0.00912 0.0660dat %>%
select(-Total) %>%
write.csv("../data/ipss_birth.csv", row.names=FALSE, quote = FALSE)累積確率分布を作る。これにより、一様乱数(y軸:0-1)から、対応する行番号が得られて、 行番号から、(出生年,性別)を得ることができる。
欲しい人数分だけ一様乱数を取れば、出生のデータが得られる。
library(ggplot2)
dat = read_csv("../data/ipss_birth.csv")## Parsed with column specification:
## cols(
## Year = col_double(),
## Sex = col_character(),
## Life = col_double(),
## ratio = col_double(),
## cum_sum = col_double()
## )dat %>%
mutate(x_axis = paste(Year, Sex, sep="")) %>%
mutate(row_num = 1:nrow(.)) %>%
ggplot() +
geom_line(aes(x = row_num, y = cum_sum), stat = "identity") +
annotate("segment",x=1,xend=90,y=0.7,yend=0.7,colour="blue",
size=1, arrow=arrow()) +
annotate("segment",x=92,xend=92,y=0.65,yend=0.0,colour="blue",
size=1,arrow=arrow()) +
annotate("text", x=20, y=0.75, parse=TRUE, label="'Random Number 0.7'") +
annotate("text", x=120, y=0.2, parse=TRUE, label="'Rownum: 95 -> 1969-Female'") +
ggtitle("Cumulative ratio of Population Birth")
これを利用して、出生のシミュレーションができる。
次は、死亡率を算出して、死亡のシミュレーションを行う。
url_death_rate = "http://www.ipss.go.jp/p-toukei/JMD/00/STATS/Mx_1x1.txt"
dat = read.table(url_death_rate, skip=2, header = TRUE)
dat %>% str## 'data.frame': 7770 obs. of 5 variables:
## $ Year : int 1947 1947 1947 1947 1947 1947 1947 1947 1947 1947 ...
## $ Age : chr "0" "1" "2" "3" ...
## $ Female: chr "0.087401" "0.033723" "0.016994" "0.011412" ...
## $ Male : chr "0.099181" "0.034697" "0.016804" "0.011461" ...
## $ Total : chr "0.093432" "0.034220" "0.016897" "0.011437" ...死亡率のデータは、Year年にX才だった人が1年後に死んでいる確率を表している。(本当に？)
今回は簡単のために、2016年の死亡率でシミュレーションを行う。
dat %<>%
filter(Year == 2016) %>%
select(-Year) %>%
select(-Total)Ageがcharになっているので、確認する。
dat$Age %>% table## .
## 0 1 10 100 101 102 103 104 105 106 107 108 109 11 110+ 12
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 13 14 15 16 17 18 19 2 20 21 22 23 24 25 26 27
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 28 29 3 30 31 32 33 34 35 36 37 38 39 4 40 41
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 42 43 44 45 46 47 48 49 5 50 51 52 53 54 55 56
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 57 58 59 6 60 61 62 63 64 65 66 67 68 69 7 70
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 71 72 73 74 75 76 77 78 79 8 80 81 82 83 84 85
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 86 87 88 89 9 90 91 92 93 94 95 96 97 98 99
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1110+があるので、これは簡単のために111才にしておく。
dat %<>%
mutate(Age = if_else(Age == "110+", "111", Age)) %>%
mutate(Age = as.integer(Age))
dat %<>%
mutate(Anb = Age) %>%
select(-Age)
dat %>% head## Female Male Anb
## 1 0.002028 0.001995 0
## 2 0.000313 0.000340 1
## 3 0.000174 0.000178 2
## 4 0.000098 0.000133 3
## 5 0.000087 0.000095 4
## 6 0.000084 0.000101 5AnbとAlbについては、あとで追記する。Anbなので、Albに直す。
\(x\)才のAnb死亡率\(q_x\)とすると、Alb死亡率は、\(\frac{q_x+q_{x+1}}{2}\)になる。
dat %<>%
mutate(Female = as.numeric(Female)) %>%
mutate(Male = as.numeric(Male)) %>%
mutate(lead_Female = lead(Female)) %>%
mutate(lead_Male = lead(Male)) %>%
mutate(F = (Female + lead_Female)/2) %>%
mutate(M = (Male + lead_Male)/2) %>%
mutate(Alb = Anb) %>%
select(Alb,F,M)
dat %>% head## Alb F M
## 1 0 0.0011705 0.0011675
## 2 1 0.0002435 0.0002590
## 3 2 0.0001360 0.0001555
## 4 3 0.0000925 0.0001140
## 5 4 0.0000855 0.0000980
## 6 5 0.0000815 0.0001060年次死亡率を可視化してみる。 横軸が年齢、縦軸が1年以内に死亡する確率。
dat %>%
pivot_longer(cols=c("F","M"), names_to = "Sex", values_to = "Mortarity") %>%
ggplot(aes(x = Alb, y = Mortarity, group = Sex, color = Sex)) +
geom_line() +
geom_hline(yintercept = 1.0, linetype = "dashed") +
annotate("text", x = 10, y = 1.05, label = 'Mortarity: 1.0') +
ggtitle("Mortarity")## Warning: Removed 2 row(s) containing missing values (geom_path).
高齢部分が不自然だが、そういうものらしい。100才以上は使わない、として切り捨てる。
また、年次死亡率を月次死亡率に変換する。
\(y = 1 - (1-x)^{12}\)
\(x = 1 - (1-y)^{0.08...}\)
dat %<>%
filter(Alb < 100)
dat %<>%
mutate(F = 1 - (1-F)**(1/12)) %>%
mutate(M = 1 - (1-M)**(1/12))
dat %>% write.csv("../data/ipss_mortality.csv", quote=F, row.names = F)
dat %>% head## Alb F M
## 1 0 9.759403e-05 9.734377e-05
## 2 1 2.029393e-05 2.158590e-05
## 3 2 1.133404e-05 1.295926e-05
## 4 3 7.708660e-06 9.500496e-06
## 5 4 7.125279e-06 8.167034e-06
## 6 5 6.791920e-06 8.833763e-06ここまでやったこと
- 公的統計の取得
- 出生率データの作成
- 死亡率データの作成
つぎにやること
擬似レセプトデータを作るために… - 出生のシミュレーション - 死亡のシミュレーション