# 下载至临时文件
if (FALSE) {
tmpxlsxpath <- file.path(tempdir(), "airquality.xlsx")
download.file("https://drwater.rcees.ac.cn/git/course/RWEP/raw/branch/PUB/data/airquality.xlsx",
destfile = tmpxlsxpath)
airqualitydf <- readxl::read_xlsx(tmpxlsxpath, sheet = 2)
metadf <- readxl::read_xlsx(tmpxlsxpath, sheet = 1)
saveRDS(airqualitydf, "./airqualitydf.RDS")
saveRDS(metadf, "./metadf.RDS")
}
airqualitydf <- readRDS("./airqualitydf.RDS")
metadf <- readRDS("./metadf.RDS")
课后作业8
数据
下载airquality.xlsx,并读取数据。
描述统计
根据airqualitydf.xlsx
,按采样点统计白天(8:00-20:00)与夜晚(20:00-8:00)中空气质量指数(AQI)中位数,按城市统计低于所有采样点AQI30%分位值的采样点占比,列出上述占比最高的10个城市(不考虑采样点数低于5个的城市)。
require(tidyverse)
airqualitydf |>
select(datetime, site, AQI) |>
filter(!is.na(AQI)) |>
group_by(site) |>
summarize(AQI.median = median(AQI, na.rm = TRUE)) |>
left_join(metadf |> select(site, city = Area)) |>
group_by(city) |>
filter(n() > 5) |>
summarize(p = sum(AQI.median < quantile(airqualitydf$AQI, probs = 0.5, na.rm = TRUE)) / n()) |>
top_n(10, p)
# A tibble: 29 × 2
city p
<chr> <dbl>
1 东莞市 1
2 伊春市 1
3 佛山市 1
4 包头市 1
5 南充市 1
6 南宁市 1
7 吉林市 1
8 呼和浩特市 1
9 哈尔滨市 1
10 大庆市 1
# ℹ 19 more rows
airqualitydf |>
select(datetime, site, AQI) |>
filter(!is.na(AQI)) |>
group_by(site) |>
summarize(AQI.median = median(AQI, na.rm = TRUE))
# A tibble: 1,626 × 2
site AQI.median
<chr> <dbl>
1 1001A 69
2 1003A 63
3 1004A 69
4 1005A 57
5 1006A 71
6 1007A 66
7 1008A 55.5
8 1009A 53
9 1010A 65
10 1011A 63.5
# ℹ 1,616 more rows
airqualitydf |>
select(datetime, site, AQI) |>
filter(!is.na(AQI)) |>
left_join(metadf |> select(site, city = Area)) |>
group_by(city) |>
filter(length(unique(site)) >= 5) |>
summarize(p = sum(AQI < quantile(airqualitydf$AQI, probs = 0.2,
na.rm = TRUE)) / n()) |>
slice_max(p, n = 10) |>
knitr::kable()
city | p |
---|---|
阜新市 | 0.9245283 |
清远市 | 0.8833333 |
玉林市 | 0.8421053 |
肇庆市 | 0.8305085 |
韶关市 | 0.7962963 |
海口市 | 0.7857143 |
伊春市 | 0.7375000 |
贵港市 | 0.7288136 |
东莞市 | 0.7285714 |
锦州市 | 0.7192982 |
统计检验
按照不同城市分组,统计白天与夜晚AQI中位数是否具有显著差异。
if (FALSE) {
require(infer)
require(tidyverse)
testdf <- airqualitydf |>
select(datetime, site, AQI) |>
filter(!is.na(AQI)) |>
left_join(metadf |> select(site, city = Area)) |>
group_by(city) |>
filter(length(unique(site)) >= 5) |>
mutate(dayornight = factor(ifelse(between(hour(datetime), 8, 20), "day", "night"),
levels = c("day", "night"))
) |>
group_by(city) |>
nest(citydf = -city) |>
mutate(median_diff = purrr::map_dbl(citydf, ~
.x |>
specify(AQI ~ dayornight) |>
calculate(stat = "diff in medians", order = c("day", "night")) |>
pull(stat)
)) |>
ungroup() |>
# slice_sample(n = 12) |>
mutate(null_dist = purrr::map(citydf, ~
.x |>
specify(AQI ~ dayornight) |>
hypothesize(null = "independence") |>
generate(reps = 1000, type = "permute") |>
calculate(stat = "diff in medians", order = c("day", "night"))
)) |>
mutate(p_value = purrr::map2_dbl(null_dist, median_diff,
~ get_p_value(.x, obs_stat = .y, direction = "both") |>
pull(p_value)
)) |>
mutate(sigdiff = ifelse(p_value < 0.01, "显著差异", "无显著差异")) |>
mutate(fig = purrr::pmap(list(null_dist, median_diff, city, sigdiff),
~ visualize(..1) +
shade_p_value(obs_stat = ..2, direction = "both") +
ggtitle(paste0(..3, ":", ..4)) +
theme_sci(2, 2)
)) |>
arrange(p_value)
saveRDS(testdf, "./testdf.RDS")
}
lang <- "cn"
require(dwfun)
require(rmdify)
require(drwateR)
dwfun::init()
rmdify::rmd_init()
testdf <- readRDS("./testdf.RDS")
require(tidyverse)
testdf |>
select(city, median_diff, p_value, sigdiff) |>
knitr::kable()
city | median_diff | p_value | sigdiff |
---|---|---|---|
北京市 | -11.0 | 0.000 | 显著差异 |
石家庄市 | -23.5 | 0.000 | 显著差异 |
太原市 | -27.5 | 0.000 | 显著差异 |
阳泉市 | -28.5 | 0.000 | 显著差异 |
长治市 | -22.5 | 0.000 | 显著差异 |
晋城市 | -30.5 | 0.000 | 显著差异 |
临汾市 | -28.0 | 0.000 | 显著差异 |
赤峰市 | 12.0 | 0.000 | 显著差异 |
鄂尔多斯市 | -23.0 | 0.000 | 显著差异 |
大连市 | -29.0 | 0.000 | 显著差异 |
营口市 | -25.0 | 0.000 | 显著差异 |
长春市 | -22.5 | 0.000 | 显著差异 |
齐齐哈尔市 | -14.5 | 0.000 | 显著差异 |
牡丹江市 | -13.5 | 0.000 | 显著差异 |
上海市 | 10.0 | 0.000 | 显著差异 |
南京市 | 34.5 | 0.000 | 显著差异 |
徐州市 | 18.5 | 0.000 | 显著差异 |
常州市 | -9.0 | 0.000 | 显著差异 |
连云港市 | 20.5 | 0.000 | 显著差异 |
淮安市 | 61.5 | 0.000 | 显著差异 |
杭州市 | -12.0 | 0.000 | 显著差异 |
宁波市 | 17.0 | 0.000 | 显著差异 |
温州市 | -20.0 | 0.000 | 显著差异 |
合肥市 | 50.5 | 0.000 | 显著差异 |
蚌埠市 | 48.0 | 0.000 | 显著差异 |
淮南市 | 46.0 | 0.000 | 显著差异 |
阜阳市 | 39.5 | 0.000 | 显著差异 |
景德镇市 | 18.5 | 0.000 | 显著差异 |
济南市 | -32.0 | 0.000 | 显著差异 |
青岛市 | -16.5 | 0.000 | 显著差异 |
淄博市 | -52.0 | 0.000 | 显著差异 |
枣庄市 | 23.0 | 0.000 | 显著差异 |
东营市 | -48.5 | 0.000 | 显著差异 |
烟台市 | -44.0 | 0.000 | 显著差异 |
威海市 | -38.0 | 0.000 | 显著差异 |
临沂市 | 19.5 | 0.000 | 显著差异 |
德州市 | -32.0 | 0.000 | 显著差异 |
聊城市 | -36.0 | 0.000 | 显著差异 |
滨州市 | -93.0 | 0.000 | 显著差异 |
郑州市 | -57.5 | 0.000 | 显著差异 |
洛阳市 | -49.5 | 0.000 | 显著差异 |
新乡市 | -38.0 | 0.000 | 显著差异 |
焦作市 | -34.5 | 0.000 | 显著差异 |
漯河市 | 31.0 | 0.000 | 显著差异 |
南阳市 | 79.0 | 0.000 | 显著差异 |
荆门市 | 74.0 | 0.000 | 显著差异 |
衡阳市 | 15.0 | 0.000 | 显著差异 |
湛江市 | -11.5 | 0.000 | 显著差异 |
惠州市 | -12.0 | 0.000 | 显著差异 |
揭阳市 | -36.0 | 0.000 | 显著差异 |
南宁市 | 8.0 | 0.000 | 显著差异 |
凉山彝族自治州 | 8.0 | 0.000 | 显著差异 |
宝鸡市 | 29.5 | 0.000 | 显著差异 |
天津市 | -16.5 | 0.000 | 显著差异 |
绍兴市 | -9.5 | 0.000 | 显著差异 |
运城市 | -36.0 | 0.002 | 显著差异 |
吉林市 | -16.5 | 0.002 | 显著差异 |
盐城市 | 50.0 | 0.002 | 显著差异 |
潍坊市 | -38.5 | 0.002 | 显著差异 |
邯郸市 | -16.0 | 0.004 | 显著差异 |
永州市 | 13.5 | 0.004 | 显著差异 |
贵港市 | 14.0 | 0.004 | 显著差异 |
海口市 | 7.0 | 0.004 | 显著差异 |
泸州市 | 11.0 | 0.004 | 显著差异 |
达州市 | 28.0 | 0.004 | 显著差异 |
宜昌市 | 62.5 | 0.006 | 显著差异 |
茂名市 | -11.0 | 0.006 | 显著差异 |
兰州市 | 14.0 | 0.006 | 显著差异 |
朔州市 | -11.0 | 0.008 | 显著差异 |
南昌市 | 13.0 | 0.008 | 显著差异 |
泰安市 | -9.5 | 0.008 | 显著差异 |
清远市 | 11.5 | 0.008 | 显著差异 |
苏州市 | -7.0 | 0.012 | 无显著差异 |
扬州市 | 21.0 | 0.012 | 无显著差异 |
鹤壁市 | -34.5 | 0.012 | 无显著差异 |
赣州市 | 22.5 | 0.014 | 无显著差异 |
江门市 | -9.0 | 0.014 | 无显著差异 |
襄阳市 | 71.5 | 0.016 | 无显著差异 |
哈尔滨市 | -8.0 | 0.020 | 无显著差异 |
秦皇岛市 | -25.0 | 0.022 | 无显著差异 |
鞍山市 | -22.0 | 0.022 | 无显著差异 |
沈阳市 | -6.0 | 0.024 | 无显著差异 |
宜春市 | 15.5 | 0.024 | 无显著差异 |
深圳市 | -5.0 | 0.024 | 无显著差异 |
柳州市 | 23.5 | 0.026 | 无显著差异 |
贵阳市 | 37.0 | 0.028 | 无显著差异 |
珠海市 | -6.5 | 0.030 | 无显著差异 |
承德市 | 14.0 | 0.032 | 无显著差异 |
本溪市 | -10.0 | 0.032 | 无显著差异 |
娄底市 | 12.5 | 0.036 | 无显著差异 |
银川市 | -10.5 | 0.036 | 无显著差异 |
芜湖市 | 18.0 | 0.042 | 无显著差异 |
重庆市 | 3.5 | 0.044 | 无显著差异 |
平顶山市 | -14.5 | 0.048 | 无显著差异 |
南充市 | 11.0 | 0.048 | 无显著差异 |
马鞍山市 | 32.0 | 0.058 | 无显著差异 |
张家口市 | -8.0 | 0.058 | 无显著差异 |
抚顺市 | -10.0 | 0.064 | 无显著差异 |
邵阳市 | 22.5 | 0.076 | 无显著差异 |
韶关市 | 10.0 | 0.076 | 无显著差异 |
桂林市 | 20.5 | 0.078 | 无显著差异 |
绵阳市 | 3.5 | 0.080 | 无显著差异 |
湘潭市 | 6.0 | 0.088 | 无显著差异 |
许昌市 | -19.0 | 0.094 | 无显著差异 |
长沙市 | 4.5 | 0.096 | 无显著差异 |
遵义市 | 11.0 | 0.100 | 无显著差异 |
泰州市 | 10.0 | 0.104 | 无显著差异 |
开封市 | -28.0 | 0.106 | 无显著差异 |
常德市 | 23.5 | 0.108 | 无显著差异 |
乌鲁木齐市 | 6.0 | 0.152 | 无显著差异 |
武汉市 | 4.0 | 0.158 | 无显著差异 |
佛山市 | -10.0 | 0.158 | 无显著差异 |
怀化市 | 9.5 | 0.172 | 无显著差异 |
镇江市 | 5.0 | 0.220 | 无显著差异 |
阜新市 | 3.0 | 0.222 | 无显著差异 |
大同市 | 5.5 | 0.258 | 无显著差异 |
南通市 | 9.5 | 0.266 | 无显著差异 |
安阳市 | -9.5 | 0.280 | 无显著差异 |
昆明市 | 2.0 | 0.280 | 无显著差异 |
成都市 | 1.5 | 0.376 | 无显著差异 |
九江市 | 3.0 | 0.380 | 无显著差异 |
株洲市 | 2.0 | 0.388 | 无显著差异 |
大庆市 | -3.0 | 0.428 | 无显著差异 |
宜宾市 | 1.5 | 0.438 | 无显著差异 |
黄石市 | 4.0 | 0.450 | 无显著差异 |
厦门市 | 2.0 | 0.454 | 无显著差异 |
佳木斯市 | -2.5 | 0.472 | 无显著差异 |
石嘴山市 | -2.0 | 0.472 | 无显著差异 |
铜陵市 | 5.5 | 0.558 | 无显著差异 |
郴州市 | 3.0 | 0.564 | 无显著差异 |
岳阳市 | 3.0 | 0.566 | 无显著差异 |
包头市 | -8.5 | 0.580 | 无显著差异 |
玉林市 | 2.0 | 0.618 | 无显著差异 |
伊春市 | 1.0 | 0.688 | 无显著差异 |
锦州市 | -0.5 | 0.702 | 无显著差异 |
萍乡市 | -1.5 | 0.764 | 无显著差异 |
唐山市 | -3.0 | 0.778 | 无显著差异 |
益阳市 | 1.5 | 0.782 | 无显著差异 |
抚州市 | 1.0 | 0.806 | 无显著差异 |
广州市 | -1.0 | 0.814 | 无显著差异 |
乐山市 | 0.5 | 0.816 | 无显著差异 |
济宁市 | 1.5 | 0.826 | 无显著差异 |
无锡市 | 1.0 | 0.896 | 无显著差异 |
资阳市 | 0.5 | 0.934 | 无显著差异 |
西安市 | 2.0 | 0.954 | 无显著差异 |
呼和浩特市 | -0.5 | 0.980 | 无显著差异 |
福州市 | 0.5 | 0.982 | 无显著差异 |
保定市 | 0.0 | 1.000 | 无显著差异 |
日照市 | -0.5 | 1.000 | 无显著差异 |
汕头市 | 0.0 | 1.000 | 无显著差异 |
肇庆市 | 0.5 | 1.000 | 无显著差异 |
东莞市 | -1.0 | 1.000 | 无显著差异 |
攀枝花市 | 0.0 | 1.000 | 无显著差异 |
拉萨市 | 0.0 | 1.000 | 无显著差异 |
克拉玛依市 | 0.0 | 1.000 | 无显著差异 |
testdf |>
mutate(grp = (row_number() - 1)%/% 12) |>
group_by(grp) |>
nest(grpdf = -grp) |>
ungroup() |>
# slice(1) |>
mutate(gp = purrr::map(grpdf,
~(.x |>
pull(fig)) |>
patchwork::wrap_plots(ncol = 3) +
dwfun::theme_sci(5, 7))) |>
pull(gp)
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[9]]
[[10]]
[[11]]
[[12]]
[[13]]