R言語による電子カルテデータの二次利用

~R言語初心者がデータ処理を楽しめるように基本的内容中心のサイトです~

ggplot2で人口ピラミッドを描く

まずはサンプルデータを落としてきます。

library(dplyr)
library(RCurl)
# This is not actual patient data.
url <- getURL("https://raw.githubusercontent.com/Algo1970/EHR_data/master/pt_master_sample.csv")
df <- read.csv(text = url, header = TRUE)

よく見る患者マスターは、こんなデータ構成かと思います。(下記はコンピューターでランダムに作成されたデータです)

> df %>% tail()
     facility_code    id        name            kana gender
995        1234567 10995   金田 彩華   かねだ あやか     女
996        1234567 10996     藤森 兼   ふじもり けん     男
997        1234567 10997   木村 誠一 きむら せいいち     男
998        1234567 10998     足立 優     あだち ゆう     女
999        1234567 10999   小堀 一輝   こほり かずき     男
1000       1234567 11000 阿久津 一代   あくつ かずよ     女
         birth address
995  1948/6/19  茨城県
996  1956/1/24  岐阜県
997   1961/8/2  沖縄県
998  1984/7/13  愛知県
999  1978/4/30  岐阜県
1000  1968/9/5  岩手県
>

年齢が必要なので、誕生日から年齢を計算します。性別も文字列を変更しておきます(性別データも電子カルテメーカー毎にまちまちなので下記コードで…)。

# birthベクトルから年齢計算
make.age <- function(birth){
  df <- NULL
  for(i in 1:length(birth)){
    df[i] <- (length(seq(as.Date(birth[i]), Sys.Date(), "year"))-1)
  }
  df
}
df$age <- make.age(df$birth)

# genderを変更
df$gender <-as.character(df$gender) 
df$gender[df$gender %in% c("m","male","M","Male","MALE","男","男性")] <- c("male")
df$gender[df$gender %in% c("f","female","F","Female","FEMALE","女","女性")] <- c("female")

いちど確認。

> df %>% tail()
     facility_code    id        name            kana gender
995        1234567 10995   金田 彩華   かねだ あやか female
996        1234567 10996     藤森 兼   ふじもり けん   male
997        1234567 10997   木村 誠一 きむら せいいち   male
998        1234567 10998     足立 優     あだち ゆう female
999        1234567 10999   小堀 一輝   こほり かずき   male
1000       1234567 11000 阿久津 一代   あくつ かずよ female
         birth address age
995  1948/6/19  茨城県  68
996  1956/1/24  岐阜県  60
997   1961/8/2  沖縄県  55
998  1984/7/13  愛知県  32
999  1978/4/30  岐阜県  38
1000  1968/9/5  岩手県  48
>

年齢も、性別文字列変更もできてますね。
必要なデータのみtempにいれて、まず集計。

df[,c("gender","age")]->temp
result <- tapply(temp$age,temp$gender, hist, breaks=seq(0,100,5))

結果はこんな感じで返します。

> result
$female
$breaks
 [1]   0   5  10  15  20  25  30  35  40  45  50  55  60  65
[15]  70  75  80  85  90  95 100

$counts
 [1]  0  0  0  7 56 43 50 45 43 45 33 41 38 34 33 47  0  0  0
[20]  0

$density
 [1] 0.000000000 0.000000000 0.000000000 0.002718447
 [5] 0.021747573 0.016699029 0.019417476 0.017475728
 [9] 0.016699029 0.017475728 0.012815534 0.015922330
[13] 0.014757282 0.013203883 0.012815534 0.018252427
[17] 0.000000000 0.000000000 0.000000000 0.000000000

$mids
 [1]  2.5  7.5 12.5 17.5 22.5 27.5 32.5 37.5 42.5 47.5 52.5
[12] 57.5 62.5 67.5 72.5 77.5 82.5 87.5 92.5 97.5

$xname
[1] "X[[i]]"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$male
$breaks
 [1]   0   5  10  15  20  25  30  35  40  45  50  55  60  65
[15]  70  75  80  85  90  95 100

$counts
 [1]  0  0  0  8 42 43 35 43 32 38 40 44 47 39 35 39  0  0  0
[20]  0

$density
 [1] 0.000000000 0.000000000 0.000000000 0.003298969
 [5] 0.017319588 0.017731959 0.014432990 0.017731959
 [9] 0.013195876 0.015670103 0.016494845 0.018144330
[13] 0.019381443 0.016082474 0.014432990 0.016082474
[17] 0.000000000 0.000000000 0.000000000 0.000000000

$mids
 [1]  2.5  7.5 12.5 17.5 22.5 27.5 32.5 37.5 42.5 47.5 52.5
[12] 57.5 62.5 67.5 72.5 77.5 82.5 87.5 92.5 97.5

$xname
[1] "X[[i]]"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

> 

必要なのは、result$(性別)$breaksとresult$(性別)$countsなので、取り出してデータフレームを作成します。
breaksは最後のデータが不要なので削除。カテゴリーの順序はlevelsで指定しておかないと後で順番がおかしくなるかも…
男女別に作ってrbindで結合。

paste(result$male$breaks,"-")->category
category[1:(length(category)-1)]->category
category <- factor(category,levels = category)

temp_male <- data.frame(category=category,
                        counts=result$male$counts,
                        gender=rep("male",length(result$male$counts)))
temp_female <- data.frame(category=category,
                        counts=result$female$counts*(-1),
                        gender=rep("female",length(result$female$counts)))
temp_all <- rbind(temp_male,temp_female)

できたデータはこんな感じ。

> temp_all
   category counts gender
1       0 -      0   male
2       5 -      0   male
3      10 -      0   male
4      15 -      8   male
5      20 -     42   male
6      25 -     43   male
7      30 -     35   male
8      35 -     43   male
9      40 -     32   male
10     45 -     38   male
11     50 -     40   male
12     55 -     44   male
13     60 -     47   male
14     65 -     39   male
15     70 -     35   male
16     75 -     39   male
17     80 -      0   male
18     85 -      0   male
19     90 -      0   male
20     95 -      0   male
21      0 -      0 female
22      5 -      0 female
23     10 -      0 female
24     15 -     -7 female
25     20 -    -56 female
26     25 -    -43 female
27     30 -    -50 female
28     35 -    -45 female
29     40 -    -43 female
30     45 -    -45 female
31     50 -    -33 female
32     55 -    -41 female
33     60 -    -38 female
34     65 -    -34 female
35     70 -    -33 female
36     75 -    -47 female
37     80 -      0 female
38     85 -      0 female
39     90 -      0 female
40     95 -      0 female
> 

女性のカウント数に-1をかけているのは、あとでggplotで作図するためです。

brks <- seq(-60, 60, 20)
lbls = paste0(as.character(c(seq(60, 0, -20), seq(20, 60, 20))))

ggplot(temp_all, aes(x = category, y = counts, fill = gender)) +   
  geom_bar(stat = "identity", width = .8) +   
  scale_y_continuous(breaks = brks,labels = lbls) + 
  coord_flip() +  
  labs(title="Population pyramid") +
  theme_tufte() +  
  theme(text=element_text(size=12, family="Comic Sans MS"), 
        plot.title = element_text(hjust = .5), 
        axis.ticks = element_blank()) +   
  scale_fill_manual(values=c("#9999CC","#CC6666"))

f:id:r_beginner:20170115152512j:plain

サンプルデータは適当に作った生年月日データなので、実際の分布とは大きく異なりますが、作図のイメージはできたかと思います。
時間があったら、ちゃんとしたサンプルデータに変更します。

〈追記〉
総務省統計局に年齢各性別人口(平成26年)のデータがありましたので、それからサンプリングして、20000人分の性別、年齢データを作ってみます。
ホームページのエクセルファイルから年齢別の人口を取り出し、CSVファイルにしてありますので、読み込んでみてください。

url <- getURL("https://raw.githubusercontent.com/Algo1970/EHR_data/master/Population_in_Japan.csv")
df <- read.csv(text = url, header = TRUE)

データの途中を確認してみます。

> df[40:50,]
   age  male female
40  39   954    928
41  40 1,006    979
42  41 1,022    998
43  42 1,004    977
44  43   976    955
45  44   948    928
46  45   931    917
47  46   910    898
48  47   907    894
49  48   707    702
50  49   874    867
> 

先程作ったデータと似ていますが、数字にカンマが入っています。これはマズイです。

df$male <- df$male %>% as.character() %>% gsub(",","",.) %>% as.numeric()
df$female <- df$female %>% as.character() %>% gsub(",","",.) %>% as.numeric()

gsub関数でカンマを外してから数値に変換します。
つぎに年齢ごと、性別ごとの確率を求めて、新しいカラムを作成します。

df$male.ratio <- round(df$male/sum(df$male),6)
df$female.ratio <- round(df$female/sum(df$female),6)

新しい確率カラムを確認してみます。

> df %>% head()
  age male female male.ratio female.ratio
1   0  524    496   0.008480     0.007604
2   1  533    508   0.008626     0.007788
3   2  534    508   0.008642     0.007788
4   3  548    519   0.008868     0.007956
5   4  534    509   0.008642     0.007803
6   5  534    510   0.008642     0.007818
> 

ちゃんとできていますね。
きれいなピラミッドが見たいので、10000人ずつサンプルと取ってみます。

sample(df$age,size=10000,replace = T,df$male.ratio) -> male.age
sample(df$age,size=10000,replace = T,df$female.ratio) -> female.age
male.df <- data.frame(gender=rep("male",length(male.age)),age=male.age)
female.df <- data.frame(gender=rep("female",length(female.age)),age=female.age)
temp <- rbind(male.df,female.df)

先ほどと同じ形のデータフレームができているか確認しましょう。

> temp %>% tail()
      gender age
19995 female  41
19996 female   0
19997 female   4
19998 female  38
19999 female  45
20000 female   5
>

大丈夫です。あとは先程と同じながれで、データを整形してプロットします。

result <- tapply(temp$age,temp$gender, hist, breaks=seq(0,100,5))
paste(result$male$breaks,"-")->category
category[1:(length(category)-1)]->category
category <- factor(category,levels = category)
temp_male <- data.frame(category=category,
                        counts=result$male$counts,
                        gender=rep("male",length(result$male$counts)))
temp_female <- data.frame(category=category,
                          counts=result$female$counts*(-1),
                          gender=rep("female",length(result$female$counts)))
temp_all <- rbind(temp_male,temp_female)
brks <- seq(-600, 600, 200)
lbls = paste0(as.character(c(seq(600, 0, -200), seq(200, 600, 200))))
ggplot(temp_all, aes(x = category, y = counts, fill = gender)) +   
  geom_bar(stat = "identity", width = .8) +   
  scale_y_continuous(breaks = brks,labels = lbls) + 
  coord_flip() +  
  labs(title="Population pyramid") +
  theme_tufte() +  
  theme(text=element_text(size=12, family="Comic Sans MS"), 
        plot.title = element_text(hjust = .5), 
        axis.ticks = element_blank()) +   
  scale_fill_manual(values=c("#9999CC","#CC6666"))

f:id:r_beginner:20170115180915j:plain
こんどは美しい人口ピラミッドが作成されました。


参考