R analysis of The US-born population in 2006
需要的R包
nutshell:此包出自书籍《R in a Nutshell, 2nd edition》
# install.packages("ggplot2") # install.packages("nutshell") |
library("ggplot2") library("nutshell") |
载入数据
这里面包含了2006年在美国出生的婴儿的记录，一次出生对应一个记录。注意后缀名.smpl(sample),births2006.smpl是随机抽取百分之十的样本。
数据来源：https://www.cdc.gov/nchs/data_access/Vitalstatsonline.htm
data(births2006.smpl) |
查看数据
head(births2006.smpl) |
## DOB_MM DOB_WK MAGER TBO_REC WTGAIN SEX APGAR5 ## 591430 9 1 25 2 NA F NA ## 1827276 2 6 28 2 26 M 9 ## 1705673 2 2 18 2 25 F 9 ## 3368269 10 5 21 2 6 M 9 ## 2990253 7 7 25 1 36 M 10 ## 966967 3 3 28 3 35 M 8 ## DMEDUC UPREVIS ESTGEST DMETH_REC DPLURAL DBWT ## 591430 NULL 10 99 Vaginal 1 Single 3800 ## 1827276 2 years of college 10 37 Vaginal 1 Single 3625 ## 1705673 NULL 14 38 Vaginal 1 Single 3650 ## 3368269 NULL 22 38 Vaginal 1 Single 3045 ## 2990253 2 years of high school 15 40 Vaginal 1 Single 3827 ## 966967 NULL 18 39 Vaginal 1 Single 3090 |
数据解释
可知样本一共有427323条记录，那么美国在2006年总共出生记录有4273230条，也就是2006年大约出生四百多万人口,其实现在也差不多.
我们先来看一下列名，总共有13个，分别是：
DOB_MM
Month of date of birth : 出生的月份 范围：1~12
DOB_WK
Day of week of birth : 出生的星期 注意的是，看一下日历，1为星期日，7为星期六 范围：1~7
MAGER
Mother’s age : 母亲的年龄 范围：12~50
TBO_REC
Total birth order : 出生顺序，其实就是第几胎 范围：1~8
WTGAIN
Weight gain by mother : 母亲体重增加值 范围:0~98 lb
SEX
factor with levels F M, representing the sex of the child : 性别
APGAR5
APGAR score :这是1953年一位叫Apgar的美国麻醉科医生提出的一种快速评价新生宝宝健康情况的指标, 满10分者为正常新生儿，评分7分以下的新生儿考虑患有轻度窒息，评分在4分以下考虑患有重度窒息 , 大部分新生儿的评分多在7 - 10分之间。
DMEDUC
Mother’s education level :母亲的教育水平 elementary school :小学 初中 high school: 高中 college :大学
UPREVIS
Number of prenatal visits : 出生以前去医院作了几次检查 范围:0~49
ESTGEST
Estimated weeks of gestation :估计的妊娠期 范围 : 12~51周
DMETH_REC
Delivery Method : C-Section/delivery(剖腹产) Vaginal(顺产)
Tips: Cesarean section :Caesar 凯撒大帝, 又称为君王剖腹术,据说他是第一个剖腹而生的,而凯撒在拉丁文中有“切开”的意思,所以叫做C-Section.
DPLURAL
Plural Births : 1 Single(就一个), 2 Twin(双胞胎), 3 Triplet(三胞胎), 4 Quadruplet(四胞胎), and 5 Quintuplet or higher(五胞胎或更多)
DBWT
Birth weight, in grams 出生体重 范围: 227~8165g
表中不存在的数据
数值数据用NA(Not Available/Not Applicable) ; 字符数据用NULL ; 或者其他可能的处理, 比如99等
数据预处理
# 变量太长了,简洁一点 births<-births2006.smpl # 判断婴儿出生那天是工作日(2~6)还是周末(1,7) # 6 %in% c(1, 7) + 1 - 2 births$WEEKEND <- c("Weekday","Weekend")[( births$DOB_WK %in% c(1, 7)) + 1] |
判断阿普加评分
CritLow(<=3),Low(4~6),Normal(>=7)
# findInterval(1:10,c(3.5,6.5)) # ---- 0 0 0 1 1 1 2 2 2 2 births$HEALTH <- c("CritLow","Low","Normal")[( findInterval(births$APGAR5, c(3.5, 6.5))) + 1] |
处理缺失值(NA)
因为R偏爱NA代表缺失值,所以染指一下数据
将妊娠期为99 weeks的替换为NA; 将产前检查为99 次的替换为NA
# replace(c(1:10),c(1:10)==5,NA) # 1 2 3 4 NA 6 7 8 9 10 births$ESTGEST <- replace(births$ESTGEST, births$ESTGEST==99, NA) births$UPREVIS <- replace(births$UPREVIS, births$UPREVIS==99, NA) |
最终结果
head(births) |
## DOB_MM DOB_WK MAGER TBO_REC WTGAIN SEX APGAR5 ## 591430 9 1 25 2 NA F NA ## 1827276 2 6 28 2 26 M 9 ## 1705673 2 2 18 2 25 F 9 ## 3368269 10 5 21 2 6 M 9 ## 2990253 7 7 25 1 36 M 10 ## 966967 3 3 28 3 35 M 8 ## DMEDUC UPREVIS ESTGEST DMETH_REC DPLURAL DBWT ## 591430 NULL 10 NA Vaginal 1 Single 3800 ## 1827276 2 years of college 10 37 Vaginal 1 Single 3625 ## 1705673 NULL 14 38 Vaginal 1 Single 3650 ## 3368269 NULL 22 38 Vaginal 1 Single 3045 ## 2990253 2 years of high school 15 40 Vaginal 1 Single 3827 ## 966967 NULL 18 39 Vaginal 1 Single 3090 ## WEEKEND HEALTH ## 591430 Weekend <NA> ## 1827276 Weekday Normal ## 1705673 Weekday Normal ## 3368269 Weekday Normal ## 2990253 Weekend Normal ## 966967 Weekday Normal |
使用 ggplot
先来看一下数据(births)
可知有427323个observations(行观测值)和15个variables(列变量)
嗯，数据太多了,减少一点,更快一些
births.small <- births[1:10000, ] bir <- ggplot(births.small) #产生一张灰色的白纸 bir |
画张柱状图
我们来看一下每个星期的周一到周日分别有多少birth,1 (星期天) ~ 7 (星期六).
Tips： factor(因子)能将数据进行classification(分类) DOB_WK：Day of week of birth
周末比平时出生的少嘛，因为医生可以延期分娩，比如Anti-contraction pills or Tocolytic agents.
bir + geom_bar(aes(x=factor(DOB_WK))) |
你还可以加些佐料，比如分娩方式,fill可以加一个颜色分类
bir + geom_bar(aes(x=factor(DOB_WK), fill=DMETH_REC),position="dodge") |
分面(Facet)
首先，我们画一个年龄的直方图，binwidth决定组距的大小，那么我们是否可以猜想，老年妇女不太可能接受宫缩抑制剂治疗(tocolytic)，因为他们可能已经有其他并发症了(complications)。
# range(births.small$MAGER) # 12 47 # 47 - 12 + 1 = 36 groups bir + geom_histogram(aes(x=MAGER), binwidth=1) |
可以验证一下分组 geom_vline(geometry verticle line)画线的好工具
bir + geom_histogram(aes(x=MAGER), binwidth=1) + geom_vline(xintercept=seq(12.5,46.5,by=1), alpha=0.2, color="white") |
还想到了正态分布，对吧 先来看一下总体密度曲线 geom_histogram中fill与col看图像理解即可 stat_density密度曲线，col(颜色)，fill(填充，NA表示不要填充)，stat_density对应y=..density..
p <- bir + geom_histogram(aes(x=MAGER,y=..density..), binwidth=1,fill = "purple",col = "yellow") + stat_density(aes(x = MAGER), col="red", fill=NA) p |
再加一个正态分布
p + stat_function(fun = dnorm, col = "blue", args=list(mean=mean(births.small$MAGER), sd=sd(births.small$MAGER))) |
Day of the week versus women’s age
先不要哪些 unknown or NA 的数据
bir <- ggplot(subset(births, DMETH_REC %in% c("C-section","Vaginal")& TBO_REC %in% c(1:8))) |
# only RHS(right hand side) bir + geom_histogram(aes(x=MAGER, fill=DMETH_REC), binwidth=1) + facet_wrap( ~ DOB_WK, scales = "free_y") |
labs(labels)可以设置标题、副标题什么的
# facet_grid:LHS + RHS bir + geom_histogram(aes(x=MAGER, fill=factor(TBO_REC)), binwidth=1) + facet_grid(WEEKEND ~ DMETH_REC, scale="free_y") + labs(title="Births in USA 2006", fill="Birth\nOrder") + geom_vline(xintercept=seq(15,45,by=5), alpha=0.2, color="white") |
线性回归
最后, 我们看一下完整数据集中的四胞胎与产前检查, 次数UPREVIS，母亲的年龄MAGER，预计妊娠期ESTGEST，分娩方式DMETH_REC和母亲的文化水平DMEDUC的相关性
ggplot(subset(births[!is.na(births$DMEDUC),], DPLURAL=="4 Quadruplet"), aes(x=UPREVIS, y=MAGER)) + geom_point(aes(size=ESTGEST, shape=DMETH_REC, col=DMEDUC)) + stat_smooth(aes(x=UPREVIS, y=MAGER, fill=DMETH_REC), method="lm") + scale_size(range=c(3, 10)) + scale_color_brewer(palette="Set1") + scale_shape(solid=TRUE) |
