课后作业8

作者

姓名

发布日期

2024年4月9日

数据

下载airquality.xlsx，并读取数据。

# 下载至临时文件
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")

描述统计

根据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]]