4 数据的描述

4.1 引子

记不太清楚是从什么时候开始，在很多研究生的认识里描述性统计工作的重要性逐渐下降，越来越追求包含复杂模型的“高大上”研究，当然也包括我自己。简单分析之后，发现可能是这些年因果推断革命在社会科学领域的兴起，以及学术期刊对于“高大上”研究的偏好所引导的吧。

尽管这是现阶段很常见的一种现象，但是这并不代表数据描述在科研研究中就不重要了。通过回归模型获得因果推断证据的首要前提就是在现实世界中这个因果关系是客观存在的，而这就需要进行充分的数据描述和探索性分析去发现这种客观规律，如果数据描述过程中无法显示出较明显的因果关联，那么因果推断模型提示的显著性将值得怀疑。

数据的描述，或者是描述性统计分析，包含的内容其实比较简洁，就是数据的集中趋势（Central Tendency）、离散趋势（Variability）和分布（Frequency Distribution）。学习过统计学的人都知道均数标准差，但是想要做好一项描述性统计分析工作却是需要丰富的经验以及一定的耐心和精力，不过高质量的数据描述也会给整个研究带来非常大的收获。因此描述性统计分析也是探索性数据分析（Exploratory Data Analysis）的重要组成部分。

当然了，数据描述不单单只是为了弄清楚单个变量的集中、离散和分布情况，而是从杂乱的数据中，找到不同的变量之间的关联性，按照数据科学中的术语即是寻找个体不同属性或特征之间的关联性。尽管相关性不是因果关系，但是相关性往往是因果关系的前提，并且发现事物之间新的相关性也是十分有价值的，比如大家所熟知的“尿布与啤酒”的关系。

在本章中，依据数据的展现形式，我将数据描述分为表和图两大部分进行介绍。本来打算将数据可视化（Data Visualization）单独介绍，但是考虑到内容实在太多，精力有限无法面面俱到，又因这方面的数据和学习资料已经十分丰富了，故而放弃了单独介绍的想法。然而，一图胜千言，字不如表，表不如图，图形对于信息的展示能力是明显优于表格的，因此正确掌握可视化方法是描述统计一项重要的技能。

本章就主要内容：

数据的表格描述
数据的可视化
回归模型表格和图形的输出

4.2 数据的表格描述

在描述性统计分析中，表格常用于展示样本的基本信息、社会人口学特征等在不同组间等分布等，通常也是一篇医学或者社会科学方向学术论文的第一个表格，一般也称之为Table One或者Summary Table。此外，统统计表格也常用于比较研究结局变量在不同组别之间分布情况，比如展示参保病人和未参保病人的年卫生支出。

这一类表格具有共同的特点，一是展示的信息基本是固定的，对于连续型变量通常展示均数±标准差、中位数（四分位间距）、最大值、最小值等，对于分类变量通常展示频数（百分比）；二是展示的结果通常是需要通过聚合运算之后才能得出；三是重复工作量大。

这三个特点就使得完成这部分的描述工作简单但是费时。对于定量研究的初学者，面对这样的工作，最广泛采用的方法应该是：首先利用统计软件分别计算得出每个需要描述的变量结果之后，然后再在MS Word或者MS Excel中进行手动整理，最后根据研究报告或者学术论文的要求调整相应的格式。这样会带来一个十分明显的弊端，若数据库发生变化，那么所有的工作均需要重新来过。再就是对于数据清洗和探索性分析时极其不友好¹。分析数据集发生变化也是经常发生的状况，例如前期数据清洗不完整、计算方法错误等。以上所提及的这部分内容实际也属于可重复研究（Reproducible Research）的范畴。在本章中就将会介绍如何快速地完成这项工作。在本书中就将会介绍如何快速地完成这项工作。

4.2.1 R中表格描述的推荐工具

根据我自己的使用经验，整理出了几个在R中进行表格描述非常便利的工具，各个Package的特点可以在官方文档中进行查看，这里不赘述。近期我个人最推荐modelsummary、table1、gtsummary这三个包，以下内容也将利用它们分别进行演示。表 4.1 中统计描述是指支持自动完成变量的描述。

表 4.1: R中支持输出符合统计规范表格的第三方包情况一览

包名称	支持功能	推荐指数
modelsummary	统计描述，回归模型	***
gtsummary	统计描述，回归模型、自定义表格	***
table1	统计描述、自定义表格	**
stargazer	回归模型	*
gt	统计描述、自定义表格	**
kableExtra	自定义表格	**
flextable	自定义表格	**

4.2.2 快速生成Summary Table

这部分用卫生政策研究中常使用的公开微观数据库CHARLS进行演示，使用的是美国南加州大学社会经济研究中心（CESR）对原始数据进行整理后所提供的Harmonized CHARLS数据库。

首先，可以利用小节 3.3.3 中提到的datadictionary包生成Harmonized CHARLS数据库的数据字典，以便于对数据库的基本信息进行概览。当然了，在实际的分析过程中，数据字典应该是在分析之前就已经准备好了，而数据概览操作通常使用str()函数或者glimpse()函数来进行，在这里用str()函数进行演示。从结果可以看出有50658条记录，19296个受访者。

library(datadictionary)
library(readr)

harm <- 
  read_csv("dataset/harmonized_data_long_sub220522.csv")

str(harm[, 1:10])

tibble [50,658 × 10] (S3: tbl_df/tbl/data.frame)
 $ ...1       : num [1:50658] 1 2 3 4 5 6 7 8 9 10 ...
 $ ID         : num [1:50658] 1.01e+10 1.01e+10 1.01e+10 1.01e+10 1.01e+10 ...
 $ householdID: num [1:50658] 1.01e+08 1.01e+08 1.01e+08 1.01e+08 1.01e+08 ...
 $ communityID: num [1:50658] 101041 101041 101041 101041 101041 ...
 $ rdob       : Date[1:50658], format: "1965-05-01" "1965-05-01" ...
 $ ragender   : num [1:50658] 2 2 2 2 1 2 2 2 2 1 ...
 $ raeducl    : num [1:50658] 1 1 1 1 1 1 1 1 1 1 ...
 $ marriage   : num [1:50658] 1 1 1 1 1 1 1 1 1 1 ...
 $ wave       : num [1:50658] 1 2 3 4 1 1 2 3 4 4 ...
 $ ruralcounty: num [1:50658] 1 1 1 1 1 1 1 1 1 1 ...

这里简单选择4个典型的分类变量（ragender、raeducl、wave、ruralcounty），3个数值型变量（self_health、cesd10、mmse、hosfee）作为示例。由于选择的这四个分类变量在数据库中的类型为numeric²，这一点可以从以上的codebook结果中看出，因此为了方便进行后续的统计描述，建议先将分类变量的类型转换为character，代码如下。当然也可以无需提前进行变量类型转换，见后文中关于gtsummary包的例子。

library(dplyr)
harm <- harm %>%
    mutate_at(vars(c("ragender", "raeducl", 
                     "wave", "ruralcounty")), 
              as.character)

4.2.2.1 仅分类变量

只对分类变量进行描述的需求可能在论文撰写过程中不会经常出现，但是在对数据进行概览时，还是会经常需要。表 4.1 中提到的package基本都具备此功能，但更推荐使用table1和gtsummary包。

如果只用于数据概览，不考虑输出成规范的表格，推荐使用modelsummary包中的datasummary_skim()函数，结果见表 4.2 。

library(modelsummary)

datasummary_skim(harm[, c("ragender", "raeducl", 
                          "wave", "ruralcounty")], 
                 type="categorical")

表 4.2: 利用modelsummary包对分类变量进行描述

		N	%
ragender	1	12800	25.3
	2	37858	74.7
raeducl	1	43984	86.8
	2	5369	10.6
	3	1297	2.6
	NA	8	0.0
wave	1	11005	21.7
	2	12102	23.9
	3	13601	26.8
	4	13950	27.5
ruralcounty	0	22517	44.4
	1	28141	55.6

如果需要输出符合统计规范的表格，则推荐使用table1或者gtsummary包，结果见表 4.3 。gtsummary包tbl_summary()函数中的type参数即可用于指定在函数中应该予以数值型（Continuous）或是分类型（Categorical）进行统计，当需要改变统计表格中展示的数据类型，但不改变数据库中的变量类型时，可以使用type参数，即上文中提及的无需提前进行变量类型转换的处理办法。对于修改变量的标签，在table1中，可以通过其自带的label()函数，通过一个for循环向分析数据集添加标签，这些标签可一直存储在数据集中，但仅table1包中的函数可以调用。在gtsummary包tbl_summary()函数中，则可直接通过label参数以列表的形式指定，具体代码如下。

library(table1)
vars <- c("ragender", "raeducl", 
          "wave", "ruralcounty", 
          "self_health", "cesd10",
          "mmse", "hosfee")

lbs <- c("性别", "受教育程度", 
         "访次", "城乡", 
         "自评健康", "认知水平", 
         "自我行动能力评分", "住院费用")

for (i in 1:length(lbs)) {
    label(harm[[vars[i]]]) <- lbs[i]
}

table1(~ ragender + raeducl + 
         wave + ruralcounty, 
         data = harm)
#---#
library(gtsummary)
tbl_summary(harm[, c("ragender", "raeducl", 
                     "wave", "ruralcounty")],
           missing_text = "(Missing)",
           type = list(ragender ~ "categorical", raeducl ~ "categorical"),
           label = list(ragender ~ "Gender",
                        raeducl ~ "Education",
                        wave ~ "Wave",
                        ruralcounty ~ "Rural or Not"),
           ) %>% as_gt()

表 4.3: 利用table1和gtsummary包作仅分类变量描述表

(a) table1包

	Overall (N=50658)
性别
1	12800 (25.3%)
2	37858 (74.7%)
受教育程度
1	43984 (86.8%)
2	5369 (10.6%)
3	1297 (2.6%)
Missing	8 (0.0%)
访次
1	11005 (21.7%)
2	12102 (23.9%)
3	13601 (26.8%)
4	13950 (27.5%)
城乡
0	22517 (44.4%)
1	28141 (55.6%)

(b) gtsummary包

Characteristic	N = 50,658¹
Gender
1	12,800 (25%)
2	37,858 (75%)
Education
1	43,984 (87%)
2	5,369 (11%)
3	1,297 (2.6%)
(Missing)	8
Wave
1	11,005 (22%)
2	12,102 (24%)
3	13,601 (27%)
4	13,950 (28%)
Rural or Not
0	22,517 (44%)
1	28,141 (56%)
¹ n (%)

4.2.2.2 仅数值变量

modelsummary、table1、gtsummary均推荐。如果只用于数据概览，不考虑输出成规范的表格，推荐使用modelsummary包中的datasummary()或datasummary_skim()函数，其优点在于可以快速的了解数值型变量的分布情况，这两个函数之间其实是等价的，结果见表 4.4 。如果需要将输出结果用于报告或者论文写作，方法同分类变量，此处不再举例说明。

datasummary(self_health + cesd10 + mmse +
            hosfee ~ Mean + SD + Min + Max + Histogram,
            data = harm,
            output = "default" )

表 4.4: 利用modelsummary包作仅数值变量描述表

	Mean	SD	Min	Max	Histogram
self_health	2.97	0.97	1.00	5.00	▁▂▇▃▁
cesd10	8.37	6.38	0.00	30.00	▇▆▄▃▂▂▁▁
mmse	15.43	5.03	0.00	30.00	▁▂▄▆▇▅▂▁
hosfee	1724.53	11988.19	0.00	1400000.00	▇

4.2.2.3 分类和数值变量混合描述

将分类变量和数值变量同时进行描述是在论文写作中最常见的一种需求，推荐使用table1和gtsummary包。

table1，结果见表 4.5 。

table1(~ ragender + raeducl + 
         ruralcounty + 
         self_health + cesd10 + 
         mmse + hosfee, 
         data = harm)

表 4.5: 利用table1包生成分类和数值变量描述表（一）

	Overall (N=50658)
性别
1	12800 (25.3%)
2	37858 (74.7%)
受教育程度
1	43984 (86.8%)
2	5369 (10.6%)
3	1297 (2.6%)
Missing	8 (0.0%)
城乡
0	22517 (44.4%)
1	28141 (55.6%)
自评健康
Mean (SD)	2.97 (0.966)
Median [Min, Max]	3.00 [1.00, 5.00]
Missing	2440 (4.8%)
认知水平
Mean (SD)	8.37 (6.38)
Median [Min, Max]	7.00 [0, 30.0]
Missing	3890 (7.7%)
自我行动能力评分
Mean (SD)	15.4 (5.03)
Median [Min, Max]	16.0 [0, 30.0]
Missing	15173 (30.0%)
住院费用
Mean (SD)	1720 (12000)
Median [Min, Max]	0 [0, 1400000]
Missing	951 (1.9%)

gtsummary，结果见表 4.6 。

tbl_summary(harm[, c("ragender", "raeducl", 
                     "wave", "ruralcounty", 
                     "self_health", "cesd10",
                     "mmse", "hosfee")],
           missing_text = "(Missing)",
           label = list(ragender ~ "Gender",
                        raeducl ~ "Education",
                        wave ~ "Wave",
                        ruralcounty ~ "Rural or Not",
                        self_health ~ "Self report Health",
                        cesd10 ~ "CESD 10",
                        mmse ~ "MMSE",
                        hosfee ~ "Hospital expense"),
           statistic = list(all_continuous() ~ 
                             "{mean} ({sd})")) %>% as_gt()
# as_flex_table() %>%
#     save_as_docx(path = "tbl-1.docx") #导出word格式#

# as_gt() %>%
#     gtsave("tbl-1.html") #导出html格式#

表 4.6: 利用gtsummary包作分类和数值变量描述表（二）

Characteristic	N = 50,658¹
Gender
1	12,800 (25%)
2	37,858 (75%)
Education
1	43,984 (87%)
2	5,369 (11%)
3	1,297 (2.6%)
(Missing)	8
Wave
1	11,005 (22%)
2	12,102 (24%)
3	13,601 (27%)
4	13,950 (28%)
Rural or Not
0	22,517 (44%)
1	28,141 (56%)
Self report Health
1	4,626 (9.6%)
2	6,770 (14%)
3	24,569 (51%)
4	9,759 (20%)
5	2,494 (5.2%)
(Missing)	2,440
CESD 10	8 (6)
(Missing)	3,890
MMSE	15.4 (5.0)
(Missing)	15,173
Hospital expense	1,725 (11,988)
(Missing)	951
¹ n (%); Mean (SD)

4.2.2.4 按变量分组后分类和数值变量混合描述

除上述几种情况之外，还有一种使用更频繁场景，那就是分组描述并进行统计检验，推荐使用gtsummary和table1包。

gtsummary，默认是对数值变量进行Wilcoxon rank sum test，分类变量进行Pearson’s Chi-squared test。结果见表 4.7 。

tbl_mix <- 
    tbl_summary(harm[, c("ragender", "raeducl", 
                        "wave", "ruralcounty", 
                        "self_health", "cesd10",
                        "mmse", "hosfee")],
                by = ruralcounty,
                missing_text = "(Missing)",
                label = list(ragender ~ "Gender",
                            raeducl ~ "Education",
                            wave ~ "Wave",
                            ruralcounty ~ "Rural or Not",
                            self_health ~ "Self report Health",
                            cesd10 ~ "CESD 10",
                            mmse ~ "MMSE",
                            hosfee ~ "Hospital expense"),
                statistic = list(all_continuous() ~ 
                                "{mean} ({sd})"),  
            ) 
                     
tbl_mix %>% add_p() %>% as_gt()

表 4.7: 利用gtsummary包作分组分类和数值变量描述表（一）

Characteristic	0, N = 22,517¹	1, N = 28,141¹	p-value²
Gender			<0.001
1	6,996 (31%)	5,804 (21%)
2	15,521 (69%)	22,337 (79%)
Education			<0.001
1	17,409 (77%)	26,575 (94%)
2	3,905 (17%)	1,464 (5.2%)
3	1,197 (5.3%)	100 (0.4%)
(Missing)	6	2
Wave			0.8
1	4,933 (22%)	6,072 (22%)
2	5,386 (24%)	6,716 (24%)
3	6,026 (27%)	7,575 (27%)
4	6,172 (27%)	7,778 (28%)
Self report Health			<0.001
1	2,263 (11%)	2,363 (8.8%)
2	3,467 (16%)	3,303 (12%)
3	11,280 (53%)	13,289 (49%)
4	3,413 (16%)	6,346 (24%)
5	813 (3.8%)	1,681 (6.2%)
(Missing)	1,281	1,159
CESD 10	7 (6)	9 (7)	<0.001
(Missing)	1,970	1,920
MMSE	16.6 (4.8)	14.3 (5.0)	<0.001
(Missing)	5,216	9,957
Hospital expense	2,238 (15,889)	1,317 (7,548)	<0.001
(Missing)	514	437
¹ n (%); Mean (SD)
² Pearson’s Chi-squared test; Wilcoxon rank sum test

若需要将默认的Wilcoxon rank sum test改成Student’s t-test，可以先自编t-test函数，再利用add_stat函数，具体代码如下，结果见表 4.8 。

ttest_func <- function(df, var, by, ...) {

    if (is.numeric(df[[var]])) {
        t.test(df[[var]] ~ 
               as.factor(df[[by]])) %>%
        broom::tidy() %>%
        select(statistic, p.value)
    } else {
        chisq.test(table(df[[var]], 
                         as.factor(df[[by]]))) %>%
        broom::tidy() %>%
        select(statistic, p.value)
    }

}


tbl_mix %>%
  add_stat(fns = everything() ~ ttest_func) %>%
  as_gt()

表 4.8: 利用gtsummary包作分组分类和数值变量描述表（二）

Characteristic	0, N = 22,517¹	1, N = 28,141¹
Gender
1	6,996 (31%)	5,804 (21%)
2	15,521 (69%)	22,337 (79%)
Education
1	17,409 (77%)	26,575 (94%)
2	3,905 (17%)	1,464 (5.2%)
3	1,197 (5.3%)	100 (0.4%)
(Missing)	6	2
Wave
1	4,933 (22%)	6,072 (22%)
2	5,386 (24%)	6,716 (24%)
3	6,026 (27%)	7,575 (27%)
4	6,172 (27%)	7,778 (28%)
Self report Health
1	2,263 (11%)	2,363 (8.8%)
2	3,467 (16%)	3,303 (12%)
3	11,280 (53%)	13,289 (49%)
4	3,413 (16%)	6,346 (24%)
5	813 (3.8%)	1,681 (6.2%)
(Missing)	1,281	1,159
CESD 10	7 (6)	9 (7)
(Missing)	1,970	1,920
MMSE	16.6 (4.8)	14.3 (5.0)
(Missing)	5,216	9,957
Hospital expense	2,238 (15,889)	1,317 (7,548)
(Missing)	514	437
¹ n (%); Mean (SD)

table1，结果见表 4.9 。

pvalue <- function(x, ...) {

    y <- unlist(x)
    m <- factor(rep(1:length(x), 
                   times = sapply(x, length)))
    if (is.numeric(y)) {
        p <- t.test(y ~ m)$p.value
    } else {
        p <- chisq.test(table(y, m))$p.value
    }
    c("", sub("<", "&lt;", 
              format.pval(p, digits = 3, 
                          eps = 0.001)))
}
#---#
table1(~ ragender + raeducl + 
         self_health + cesd10 + 
         mmse + hosfee | ruralcounty, 
         overall = FALSE, 
         extra.col = list("P-value" = pvalue),
         data = harm)

表 4.9: 利用table1包生成分组后分类和数值变量描述表

	0 (N=22517)	1 (N=28141)	P-value
性别
1	6996 (31.1%)	5804 (20.6%)	<0.001
2	15521 (68.9%)	22337 (79.4%)
受教育程度
1	17409 (77.3%)	26575 (94.4%)	<0.001
2	3905 (17.3%)	1464 (5.2%)
3	1197 (5.3%)	100 (0.4%)
Missing	6 (0.0%)	2 (0.0%)
自评健康
Mean (SD)	2.86 (0.940)	3.06 (0.976)	<0.001
Median [Min, Max]	3.00 [1.00, 5.00]	3.00 [1.00, 5.00]
Missing	1281 (5.7%)	1159 (4.1%)
认知水平
Mean (SD)	7.16 (5.86)	9.31 (6.61)	<0.001
Median [Min, Max]	6.00 [0, 30.0]	8.00 [0, 30.0]
Missing	1970 (8.7%)	1920 (6.8%)
自我行动能力评分
Mean (SD)	16.6 (4.80)	14.3 (5.00)	<0.001
Median [Min, Max]	17.0 [0, 30.0]	15.0 [0, 30.0]
Missing	5216 (23.2%)	9957 (35.4%)
住院费用
Mean (SD)	2240 (15900)	1320 (7550)	<0.001
Median [Min, Max]	0 [0, 1400000]	0 [0, 250000]
Missing	514 (2.3%)	437 (1.6%)

如果进行数据探索，可以使用modelsummary中的datasummary_balance()的函数快速进行组间比较及描述，结果见表 4.10 。

datasummary_balance(~ ruralcounty,
                    dinm = FALSE,
                    dinm_statistic = "p.value",
                    data = harm[, 
                           c("ragender", "raeducl", 
                             "wave", "ruralcounty", 
                             "self_health", "cesd10",
                             "mmse", "hosfee")])

表 4.10: 利用modelsummary包作分组后分类和数值变量描述表

		0 (N=22517)		1 (N=28141)
		Mean	Std. Dev.	Mean	Std. Dev.
self_health		2.9	0.9	3.1	1.0
cesd10		7.2	5.9	9.3	6.6
mmse		16.6	4.8	14.3	5.0
hosfee		2238.3	15889.1	1316.5	7548.3
		N	Pct.	N	Pct.
ragender	1	6996	31.1	5804	20.6
	2	15521	68.9	22337	79.4
raeducl	1	17409	77.3	26575	94.4
	2	3905	17.3	1464	5.2
	3	1197	5.3	100	0.4
	NA	6	0.0	2	0.0
wave	1	4933	21.9	6072	21.6
	2	5386	23.9	6716	23.9
	3	6026	26.8	7575	26.9
	4	6172	27.4	7778	27.6

4.2.2.5 多变量嵌套分组

以上示例仅展示了单变量分组的描述性表格，然而多变量嵌套分组也是十分常见的分析需求。例如，当需要比较城市与农村地区不同性别居民的特征差异时，此时就需要进行地区与性别的嵌套分组。本节同样提供基于gtsummary和table1包，两者各有利弊，可根据实际研究需求灵活选择。具体代码如下，输出结果见表 4.11 和表 4.12。

#---gtsummary---#
harm %>%
    select(c("ragender", "raeducl", 
             "wave", "ruralcounty", 
             "self_health", "cesd10",
             "mmse", "hosfee")) %>%
    tbl_strata(
      strata = ruralcounty,
      .tbl_fun =
        ~ .x %>%
    tbl_summary(
                by = ragender,
                missing_text = "(Missing)",
                label = list(ragender ~ "Gender",
                             raeducl ~ "Education",
                             wave ~ "Wave",
                             self_health ~ "Self report Health",
                             cesd10 ~ "CESD 10",
                             mmse ~ "MMSE",
                             hosfee ~ "Hospital expense"),
                statistic = list(all_continuous() ~ "{mean} ({sd})")),
      .header = "**{strata}**, N = {n}"  ) %>% as_gt()
#---table1---#
table1(~ raeducl + 
         self_health + cesd10 + 
         mmse + hosfee | ruralcounty + ragender, 
         overall = FALSE, 
         data = harm)
# tbl_stack(list(tbl_9416_1, tbl_9416_2, tbl_9416_3), 
#           group_header = c("60~69", "70~79", "≥ 80"))  %>% as_gt() %>%
#     gtsave("Output/表9-4-16-分地区城乡.html")

表 4.11: 利用gtsummary包作分组后分类和数值变量描述表

Characteristic	0, N = 22517		1, N = 28141
Characteristic	1, N = 6,996¹	2, N = 15,521¹	1, N = 5,804¹	2, N = 22,337¹
Education
1	4,479 (64%)	12,930 (83%)	4,875 (84%)	21,700 (97%)
2	1,762 (25%)	2,143 (14%)	845 (15%)	619 (2.8%)
3	755 (11%)	442 (2.8%)	84 (1.4%)	16 (<0.1%)
(Missing)	0	6	0	2
Wave
1	1,372 (20%)	3,561 (23%)	1,035 (18%)	5,037 (23%)
2	1,632 (23%)	3,754 (24%)	1,307 (23%)	5,409 (24%)
3	1,922 (27%)	4,104 (26%)	1,649 (28%)	5,926 (27%)
4	2,070 (30%)	4,102 (26%)	1,813 (31%)	5,965 (27%)
Self report Health
1	875 (13%)	1,388 (9.5%)	727 (13%)	1,636 (7.7%)
2	1,315 (20%)	2,152 (15%)	824 (15%)	2,479 (12%)
3	3,468 (53%)	7,812 (53%)	3,043 (54%)	10,246 (48%)
4	749 (11%)	2,664 (18%)	847 (15%)	5,499 (26%)
5	190 (2.9%)	623 (4.3%)	201 (3.6%)	1,480 (6.9%)
(Missing)	399	882	162	997
CESD 10	5.5 (4.9)	7.9 (6.1)	7 (6)	10 (7)
(Missing)	590	1,380	244	1,676
MMSE	17.3 (4.2)	16.2 (5.1)	15.7 (4.4)	13.9 (5.1)
(Missing)	1,083	4,133	950	9,007
Hospital expense	2,972 (23,572)	1,902 (10,643)	1,458 (8,929)	1,279 (7,141)
(Missing)	84	430	35	402
¹ n (%); Mean (SD)

	0		1
	1 (N=6996)	2 (N=15521)	1 (N=5804)	2 (N=22337)
受教育程度
1	4479 (64.0%)	12930 (83.3%)	4875 (84.0%)	21700 (97.1%)
2	1762 (25.2%)	2143 (13.8%)	845 (14.6%)	619 (2.8%)
3	755 (10.8%)	442 (2.8%)	84 (1.4%)	16 (0.1%)
Missing	0 (0%)	6 (0.0%)	0 (0%)	2 (0.0%)
自评健康
Mean (SD)	2.71 (0.934)	2.93 (0.935)	2.82 (0.960)	3.13 (0.970)
Median [Min, Max]	3.00 [1.00, 5.00]	3.00 [1.00, 5.00]	3.00 [1.00, 5.00]	3.00 [1.00, 5.00]
Missing	399 (5.7%)	882 (5.7%)	162 (2.8%)	997 (4.5%)
认知水平
Mean (SD)	5.51 (4.85)	7.91 (6.11)	7.00 (5.51)	9.94 (6.74)
Median [Min, Max]	4.00 [0, 30.0]	7.00 [0, 30.0]	6.00 [0, 30.0]	9.00 [0, 30.0]
Missing	590 (8.4%)	1380 (8.9%)	244 (4.2%)	1676 (7.5%)
自我行动能力评分
Mean (SD)	17.3 (4.18)	16.2 (5.05)	15.7 (4.38)	13.9 (5.12)
Median [Min, Max]	18.0 [1.00, 30.0]	17.0 [0, 30.0]	16.0 [0, 30.0]	14.0 [0, 30.0]
Missing	1083 (15.5%)	4133 (26.6%)	950 (16.4%)	9007 (40.3%)
住院费用
Mean (SD)	2970 (23600)	1900 (10600)	1460 (8930)	1280 (7140)
Median [Min, Max]	0 [0, 1400000]	0 [0, 500000]	0 [0, 220000]	0 [0, 250000]
Missing	84 (1.2%)	430 (2.8%)	35 (0.6%)	402 (1.8%)

table1(~ raeducl + 
         self_health + cesd10 + 
         mmse + hosfee | ruralcounty + ragender, 
         overall = FALSE, 
         data = harm)

表 4.12: 利用table1包作分组后分类和数值变量描述表

	0		1
	1 (N=6996)	2 (N=15521)	1 (N=5804)	2 (N=22337)
受教育程度
1	4479 (64.0%)	12930 (83.3%)	4875 (84.0%)	21700 (97.1%)
2	1762 (25.2%)	2143 (13.8%)	845 (14.6%)	619 (2.8%)
3	755 (10.8%)	442 (2.8%)	84 (1.4%)	16 (0.1%)
Missing	0 (0%)	6 (0.0%)	0 (0%)	2 (0.0%)
自评健康
Mean (SD)	2.71 (0.934)	2.93 (0.935)	2.82 (0.960)	3.13 (0.970)
Median [Min, Max]	3.00 [1.00, 5.00]	3.00 [1.00, 5.00]	3.00 [1.00, 5.00]	3.00 [1.00, 5.00]
Missing	399 (5.7%)	882 (5.7%)	162 (2.8%)	997 (4.5%)
认知水平
Mean (SD)	5.51 (4.85)	7.91 (6.11)	7.00 (5.51)	9.94 (6.74)
Median [Min, Max]	4.00 [0, 30.0]	7.00 [0, 30.0]	6.00 [0, 30.0]	9.00 [0, 30.0]
Missing	590 (8.4%)	1380 (8.9%)	244 (4.2%)	1676 (7.5%)
自我行动能力评分
Mean (SD)	17.3 (4.18)	16.2 (5.05)	15.7 (4.38)	13.9 (5.12)
Median [Min, Max]	18.0 [1.00, 30.0]	17.0 [0, 30.0]	16.0 [0, 30.0]	14.0 [0, 30.0]
Missing	1083 (15.5%)	4133 (26.6%)	950 (16.4%)	9007 (40.3%)
住院费用
Mean (SD)	2970 (23600)	1900 (10600)	1460 (8930)	1280 (7140)
Median [Min, Max]	0 [0, 1400000]	0 [0, 500000]	0 [0, 220000]	0 [0, 250000]
Missing	84 (1.2%)	430 (2.8%)	35 (0.6%)	402 (1.8%)

4.2.2.6 将表格导出至文档中

在完成统计表格的分析后，通常需要将其整理到文档中，譬如MS Word或者Latex等。当然也可以利用quarto或者jupyter note等嵌入式编程工具，将数据分析、报告或者论文撰写的过程进行整合，但是这要求使用者熟悉各种文件，故其使用场景相对有限。对于使用table1包时，若使用的IDE为Rstudio，则在运行上述示例代码时，在Rstudio的输出框中会显示完整表格，此时只需点击导出按钮即可将表格存储为html格式文件。导出的html文件即是一个单独的表格，可以使用MS Word或者WPS软件打开，即可十分方便的将其转入MS Word文件中。对于使用gtsummary包时，可以将以下代码通过管道函数%>%与上述示例代码联接，即可将表格输出为docx或者html格式文件，但前提需要提前安装好gt和flextable包。根据可重复性研究的原则，在研究资料分析过程中，十分推荐将表格或者图形均以单独的文件进行存储，这样即方便研究资料的管理，也可以保证研究的可信度。

# as_flex_table() %>%
#     save_as_docx(path = "tbl-1.docx") #导出word格式#

# as_gt() %>%
#     gtsave("tbl-1.html") #导出html格式#

4.2.3 回归模型的输出

回归结果的表格输出虽然不属于描述性分析的内容，但是与本节内容关联性较大，因此在此一并介绍了。推荐modelsummary包，相比于另一个用于回归模型输出的stargazer包而言，其支持的模型函数更全面，并且输出的文件类型也更全面，比如html、latex、word等，示例结果如表 4.13。modelsummary包不仅仅只是将回归结果输出，也拥有较丰富的自定义功能，比如可以指定输出的回归变量、回归参数等。同时也可以直接通过vcov参数将值传递给sandwich包来实现稳健标准误的计算。关于其支持的模型类型和稳健标准误计算方法，可详见modelsummary 。

library(plm)
library(modelsummary)

data("Grunfeld", package="plm")

grun.fe <- plm(inv ~ value + capital, 
               index=c("firm","year"),
               data = Grunfeld, model = "within")
               
grun.re <- plm(inv ~ value + capital, 
               index=c("firm","year"),
               data = Grunfeld, model = "random")

modelsummary(list("Fixed Effect" = grun.fe,
                  "Random Effect" = grun.re), 
             vcov = "robust",
             stars = TRUE,
             output = "gt")

表 4.13: 利用modelsummary包输出回归结果

	Fixed Effect	Random Effect
value	0.110***	0.110***
	(0.016)	(0.015)
capital	0.310***	0.308***
	(0.062)	(0.068)
(Intercept)		-57.834*
		(27.793)
Num.Obs.	200	200
R2	0.767	0.770
R2 Adj.	0.753	0.767
AIC	2147.6	2159.0
BIC	2157.5	2172.2
RMSE	51.16	52.39
Std.Errors	HC3	HC3
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

4.2.4 小结

如果你和我一样有选择焦虑，可选项多了反而不知所措了，那么就直接看表 4.14 。

表 4.14: 支持输出符合统计规范表格的packages总结

需求	首推荐	次推荐	最次推荐
描述表格	gtsummary	table1	modelsummary
回归结果	modelsummary	stargazer	gtsummary

4.3 数据的可视化

统计图形的种类非常丰富，并且随着可视化工具的极大发展，目前数据可视化的展示方法越来越多样化。这样带来的好处是可以充分发挥出研究者的想象力，不好的地方是眼花缭乱的可视化方式容易造成研究人员选择困难，并且难以找到各类可视化方式的是用场景。在本节中，将进行总结并给出一定的示例，以下示例均基于R的ggplot2包实现。

统计图形可以归纳为四大类，这四类图形基本可以涵盖大部分的数据可视化需求。

分布型：Distribution of a single variable
关系型：Relationship between two variables
构成型：Composition of a single or multiple variables
对比型：Comparison between different categories/individuals

4.3.1 分布型图形

4.3.1.1 直方图

直方图（Histogram）是直观了解连续型变量分布情况的首选和最有效手段之一。一般情况下用于单个变量分布的观察，如图 4.1 (a) ，但是若需要进行同一个变量不同组间，或是不同变量之间的比较时，亦可在同一个图形中展示，如图 4.1 (b) 。

library(ggplot2)

set.seed(2022)
df <- data.frame("value" = c(rnorm(1000, 6, 2), 
                             rnorm(1000, 8, 4)),
                 "group" = c(rep("Rural", 1000), 
                             rep("Urban", 1000))
                )

# Univariate
ggplot(df) +
  geom_histogram(aes(x = value, 
                     y = ..density..),
                     bins = 50,
                     fill = "blue",
                     color = "black") +
  geom_density(aes(x = value), 
               alpha = 0.5, 
               size = 1.2,
               color = "red") +
  labs(x = "Indicator",
       y = "Density") +
  theme_bw()

# Multivariates
ggplot(df) +
  geom_histogram(aes(x = value, 
                     y = ..density.., 
                     fill = group),
                     bins = 50,
                     alpha = 0.5,
                     position = "identity",
                     color = "black") +
  geom_density(aes(x = value, color = group), 
               alpha = 0.5, 
               size = 1.2) +
  labs(x = "Indicator",
       y = "Density") +
  theme_bw() +
    theme(legend.position = c(0.1, 0.85))

4.3.1.2 箱图和小提琴图

箱图（Boxplot）是另一种观察数据分布情况的统计图形，箱子的上横线表示上四分位数（Q75），下横线表示下四分位数（Q25），中间的横线为中位数（Median），上下延长线的末端表示中位数加减1.5倍的四分位间距，如图 4.2 (a) ，箱图在识别异常值时可以提供很大的帮助。

箱图（Boxplot）的缺点也很明显，就是无法观测到数据的聚集和离散情况，如在哪一区间数据分布较为集中。为了弥补这一不足，可以将小提琴图（Violin plot）或是Jitter plot与箱图结合，也可以单独使用，如图 4.2 (b) ，图中较胖的区域说明数据分布较为集中。

p <- 
    ggplot(df, aes(x = group, y = value, 
                   fill = group)) + 
    geom_boxplot(width = 0.1) +
    theme_bw() +
      theme(legend.position="none")
p
#---#
p +
  geom_violin(trim = FALSE)+
  geom_jitter(shape = 16, 
              alpha = 0.2,
              position = position_jitter(0.2)) +
  scale_fill_brewer(palette = "Dark2") +
  labs(x = "Group", 
       y = "Value") +
  theme_bw() +
    theme(legend.position="none")

4.3.2 关系型图形

探索两个因素之间的关联性或是因果关系是社会科学研究中的关键问题，而这一类问题的分析起点就是通过图形探索两个因素之间是否存在相互依存关系，此时关系型图形的作用就非常大。常见的关系型图形包括散点图、线图以及热力图。

4.3.2.1 散点图

散点图（Scatter Plot）最长应用的场景就是描述两个连续型变量之间的关系，当然有时也可用于连续型变量与分类变量之间，示例结果如图 4.3 。

data(mtcars)
#---#
ggplot(mtcars, aes(x = hp, y = mpg)) +
  geom_point() +
  labs(x = "Horsepower",
       y = "Miles/(US) gallon") +
  theme_bw()

#---#
ggplot(mtcars, aes(x = hp, y = mpg)) +
  geom_point() +
  geom_smooth(method = "loess") +
  labs(x = "Horsepower",
       y = "Miles/(US) gallon") +
  theme_bw()

#---#
ggplot(mtcars, aes(x = hp, y = mpg)) +
  geom_point(shape = 15) +
  labs(x = "Horsepower",
       y = "Miles/(US) gallon") +
  theme_bw()

#---#
# change shape, color, fill, size
ggplot(mtcars, aes(x = hp, y = mpg)) +
  geom_point(shape=23, 
             fill = "blue", 
             color = "darkred", size=3) +
  labs(x = "Horsepower",
       y = "Miles/(US) gallon") +
  theme_bw()

4.3.2.2 线图

线图（Line Plot）也是最常见的用于观察两个变量之间关系的统计图形之一，例如股市的K线图，用于观察股价与时间的关系。示例结果如图 4.4 。

data(Orange)

ggplot(Orange) +
    geom_line(aes(x = age, 
                  y = circumference, 
                  linetype = Tree)) +
    theme_bw() +
      theme(legend.position = "bottom")

4.3.2.3 热力图

热力图（Heatmap）是通过颜色的深浅来区分某个连续型变量在两个分类变量之间的集中和离散情况，在生物信息学中应用较为广泛。目前，ggplot2包支持简单的热力图可视化，若需求更高，如需要同时了解不同分类在连续型变量上的聚类情况，可以使用pheatmap包。

lett <- LETTERS[1:20]
week <- paste0("Week", seq(1,20))
data <- expand.grid("lett" = lett, "week" = week)
data$value <- runif(400, 60, 100)
 
# Heatmap 
ggplot(data, aes(lett, week, fill = value)) + 
  geom_tile() +
  scale_fill_distiller(palette = "RdPu")

4.3.3 构成型图形

构成型图形相对比较简单，主要是饼图（Pie plot）。在ggplot2中，饼图的实现方法有两种思路，一是原始数据已经给出了某个分类的构成比，此时可以直接作图，二是需要经过计算之后方可获得构成比，如本书中的示例。一般而言，第二种情况更为常见。示例结果如图 4.6 。

好物推荐

推荐一个支持自动填充符合不同学术期刊配色要求的包ggsci。
其优点在于可以让你省去不少寻找RGB等配色的时间，其中预设了包括Lancet、JAMA、Nature等期刊常用的配色板。
尽管未能覆盖全部的应用场景，但是适用面已经相对较广了。

library(dplyr)
data(mtcars)

table(mtcars$cyl) %>% 
  prop.table(.) %>%
  as.data.frame(.) %>%
  ggplot(aes(x = "", y = Freq *  100, fill = Var1)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = paste0(Freq * 100, "%")),
            position = position_stack(vjust = 0.5)) +
  labs(x = "", y = "") +
  coord_polar(theta = "y") +
  ggsci::scale_fill_jama(name = "Cyl") +
  theme_bw() +
    theme(legend.position = "bottom")

4.3.4 对比型图形

对比型图形主要是条形图（Bar Plot），但是实际研究过程中线图和散点图等均可以实现对比的分析需求，这里所讲的对比型图形主要是指在进行数据对比时更为常用和直观的图形。条形图也称作柱状图，其变化形式较多样，根据柱子的不同摆放方式，可以分为堆积条图（Stack），并排条图（Dodge），百分条图（Fill）。

三种不同类型的偏好场景略有不同，堆积条图和百分条图常用于展示计数资料，而并排条图既常用于计数变量也常用于计量变量。在ggplot2中，可以通过geom_bar()或geom_col()函数实现。条形图的变化非常多，本书中仅对常见的情况进行演示，其他的希望读者自行探索。

ggplot(mtcars, aes(x = factor(cyl), y = ..count..)) +
  geom_bar(stat = "count", fill = "blue") +
  labs(x = "Cyl", y = "Count") +
  theme_bw()

#---#
mtcars %>% group_by(cyl, gear) %>%
  summarise(mpg_mean = mean(mpg)) %>%
  ggplot(aes(x = factor(cyl), 
             y = mpg_mean, 
             fill = factor(gear))) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = mpg_mean), 
            vjust = -0.8,
            size = 3,
            position = position_dodge(0.9)) +
  labs(x = "Cyl", y = "Mpg") +
  ggsci::scale_fill_aaas(name = "Gear") +
  theme_bw()  +
    theme(legend.position = c(0.85, 0.8))

#---#
ggplot(mtcars, aes(x = factor(cyl))) +
  geom_bar(aes(fill = factor(gear)), position = "stack") +
  labs(x = "Cyl", y = "Count") +
  ggsci::scale_fill_npg(name = "Gear") +
  theme_bw() +
    theme(legend.position = "bottom")

#---#
ggplot(mtcars, aes(x = factor(cyl))) +
  geom_bar(aes(fill = factor(gear)), position = "fill") +
  labs(x = "Cyl", y = "Percentage") +
  ggsci::scale_fill_lancet(name = "Gear") +
  theme_bw() +
    theme(legend.position = "bottom")

4.3.5 分面

一般而言，二维平面图形中最多能展示三个属性，假设有4个属性（Age、BMI、Sex、Region）想要一同展示在同一个图形中，常规方法是无法完成的，如横轴为Age，纵轴为BMI，那么只能利用颜色或者形状来区别Sex或者Region。当需要在一个图形中展现四种属性时，就需要利用到分面（Facet）功能，在ggplot2包中，可以利用 facet_wrap()和facet_grid()函数实现，后者多用于多个分面变量，使结果呈现为行列交叉的形式，结果如图 4.8 。

ggplot(mtcars) +
  geom_point(aes(x = hp, y = mpg, 
                 shape = factor(vs))) +
  facet_wrap(~ cyl) +
  labs(x = "Horsepower",
       y = "Miles/(US) gallon") +
  scale_shape_discrete(name = "Engine",
                       breaks = c(0, 1),
                       labels = c("V-shaped", 
                                  "Straight")) +
  theme_bw() +
    theme(legend.position = "bottom")

4.3.6 空间可视化

关于空间可视化的内容十分复杂，其本来也是空间数据分析或者空间统计学的主要内容，想要介绍清楚需要很多的篇幅和精力。在本书中，将仅仅介绍在数据描述过程中最常用的空间填色图（Choropleth）。

4.3.6.1 地图素材类型

空间可视化的基础素材就是想要研究的地区的地图图层，一般有三种常有类型：shp素材，geojson素材，地图包内置地图素材。

shp和geojson格式素材均可转换为以下两种地图数据格式：

sp:SpatialPolygonDataFrame
sf:Simple feature list column

这两种格式的数据集所描述的信息差不多是一致的。第一种格式（sp）是R语言绘图比较传统的数据格式，它将地理信息数据分割为两大块：描述层和映射层。在数据存放时，描述层记录各个地理区域的名称、ID、编号、简写、iOS编码，以及其他标识信息和度量变量，描述层是一个dataframe，我们可以用data@data来提取描述层的数据框。而对应的几何映射层，是每一个行政区域的多边形边界点，这些边界点按照order排序，按照group分组。多边形边界点信息是一个多层嵌套的list结构，但是我们仍然可以通过fortity()函数将其转化为数据框。即sp空间数据对象是一个dataframe（描述层）和polygons（几何映射层）两个对象的组合对象。

而sf对象将这种控件数据格式件进行了更加整齐的布局，使用st_read()导入的空间数据对象完全是一个整齐的数据框，拥有整齐的行列，这些行列中包含着数据描述和几何多边形的边界点信息。其中最大的特点是，它将每一个行政区划所对应的几何边界点封装成了一个list对象的记录，这条记录就像其他普通的文本记录、数值记录一样，被排列在对应行政区划描述的单元格中。

4.3.6.2 地图来源

世界地图基本在很多支持空间分析的软件包中内置，虽然方便，但是这些地图基本会存在更新不及时导致的过于老旧问题。获得正确的地图素材是进行空间研究的基本前提。

关于中国地图需要注意的

在这里重点要指出的是在适用中国地图时，不论是商业出版还是发表学术成果，必须适用具有官方审图号的地图素材，否在会带来潜在的不必要麻烦。
目前建议使用民政部网站提供的地图（http://xzqh.mca.gov.cn/map），其审图号为GS（2022）1873号。
该地图中，县级地图数据不包括香港和澳门特别行政区，市级地图数据不包括台湾省。

4.3.6.3 地图数据的读入方式

在R中目前有三种读取地图的方法，但是第一种方法目前不太推荐。

sp::readShapePoly()
rgdal::readOGR()
sf::st_read()

4.3.6.4 R中空间可视化常见三种方式

第一种，geom_polygon()函数。此方法必须通过@data及fortity()函数将描述层与几何层分别读取后，添加需要映射的变量，然后合并后进行画图。

ggplot(map_data) +
    geom_polygon(aes(x = long, y = lat, 
                   group = group, 
                   fill = fill_var))

第二种，geom_map()函数。fill_file为包含需要作图变量的数据框，map_file 为地图素材，即通过直接读取的地图文件，使用geom_map()函数不需要进行合并，因为合并的过程由map_id所指定的参数自动进行合并，merge_id为fill_file中的识别变量。示例如图 4.9 。

ids <- factor(c("1.1", "2.1", "1.2", 
                  "2.2", "1.3", "2.3"))
values <- data.frame(
      ids = ids,
      value = c(3, 3.1, 3.1, 
              3.2, 3.15, 3.5)
    )
    
positions <- data.frame(
      id = rep(ids, each = 4),
      x = c(2, 1, 1.1, 2.2, 1, 0, 0.3, 1.1, 
          2.2, 1.1, 1.2, 2.5, 1.1, 0.3,
            0.5, 1.2, 2.5, 1.2, 1.3, 2.7, 
          1.2, 0.5, 0.6, 1.3),
      y = c(-0.5, 0, 1, 0.5, 0, 0.5, 1.5, 1, 
          0.5, 1, 2.1, 1.7, 1, 1.5,
            2.2, 2.1, 1.7, 2.1, 3.2, 2.8, 
          2.1, 2.2, 3.3, 3.2)
    )

ggplot(values, aes(fill = value)) + 
    geom_map(aes(map_id = ids, fill = value),
           map = positions, color = "red") +
    expand_limits(positions)

注意

map_file中必须包含三个变量x或long、y或lat、region或id。
map_id指定的merge_id必须是能与region或id变量进行合并。

第三种，geom_sf()函数，此方法仍然需要一次合并过程，示例如图 4.10 。

library(ggplot2)

china_map <- 
  sf::st_read("dataset/MCA_China_province.geojson",
              quiet = TRUE)

nine_line <- 
  sf::st_read("dataset/MCA_Nine_line.geojson",
              quiet = TRUE)

china_map$pop <- 
  runif(length(china_map$NAME), 100, 1000)


ggplot() +
  geom_sf(aes(fill = pop), data = china_map) +
  geom_sf(data = nine_line)

但是，绘制中国地图有一个比较特殊的地方。由于南海九段线会占用较大的绘图空间，导致图片比例不太协调，因此常以小图框的形式展示。结合一些必要的美化后，中国地图填色图的示例如图 4.11 。

library(ggspatial)
library(cowplot)

fig_main <-
  ggplot() +
    geom_sf(aes(fill = pop), data = china_map) +
    geom_sf(data = nine_line) +
    labs(caption = "Map No.: GS(2022)1873") +
    scale_fill_gradient2(name = "Population",
                        low = "orange", 
                        mid = "yellow", 
                        high = "red",
                        n.breaks = 8,
      guide = guide_legend(
                        direction = "horizontal",
                        title.position = 'top',
                        title.hjust = 0.5,
                        label.hjust = 1,
                        nrow = 1,
                        byrow = T,
                        label.position = "bottom")) +
    annotation_scale(location = "bl") +
    annotation_north_arrow(location = "tl",
                          height = unit(0.8, "cm"),
                          width = unit(0.8, "cm")) +
    coord_sf(ylim = c(1869414.769862395, 
                      7187874.74616931), 
             crs = 3857) +
    theme_map() +
      theme(legend.position = c(0.5, 0.85))


#-南海九段线-#
fig_nine <- 
  ggplot() +
    geom_sf(aes(fill = pop), data = china_map) +
    geom_sf(data = nine_line) +
    scale_fill_gradient2(name = "Population",
                        low = "orange", 
                        mid = "yellow", 
                        high = "red",
                        n.breaks = 8) +
    coord_sf(ylim = c(278392.10080518876, 
                      2991845.069153875),
             xlim = c(11631734.185889415, 
                      13868701.579770062), 
             crs = 3857) +
    theme_map() +
      theme(legend.position = "none",
            plot.margin = unit(c(0, 0, 0, 0), "mm"),
            panel.border = element_rect(fill = NA,
                                        color = "grey10",
                                        linetype = 1,
                                        size = 0.5)
            )

# 使用cowplot包将主图和小图框合并
fig_china <- ggdraw(fig_main) +
             draw_plot(fig_nine, 
                       x = 0.84, 
                       y = 0.058, 
                       width = 0.15, 
                       height = 0.45)

# pdf("China map.pdf", width = 10, height = 14)
fig_china
# dev.off()

4.3.6.5 R中的tmap包

R中支持空间可视化的第三方包选项还是十分丰富的，比如上文演示的ggplot2包，还有latticeExtra包、tmap包、leaflet包等，这些包各有优劣，一般综合起来使用。这里再推荐tmap包，其具有几点优点：

与ggplot2包相比，在地理可视化方面其更加专业，可以提供比例尺、指北针等选项。当然结合ggspatial包，ggplot2也可以实现这些功能，如图 4.11 。
与leaflet包相比，其不仅同样支持reveal.js动态交互地图，还能输出符合学术要求的图片格式。

tm_shape(china_map) +
  tm_polygons(col = "pop", 
              pal = c("white", "skyblue")) +
  tm_scale_bar(position = c("left", "bottom"), 
               width = 0.1) +
  tm_compass(position = c("right", "bottom"),
             size = 1) +
  tm_layout(legend.width = 0.1)

4.3.6.6 关于坐标系

做空间可视化，坐标系是不可忽略的一个问题。投影的问题源于地球是球形，而地图是平面，因此就存在一个如何将球体转换为平面的问题，这里面的学问挺大的，我也讲不清楚就不展开了。投影可以理解为将地理坐标系转换到投影坐标系的过程，即将地球曲面转换为地图平面的过程。地理坐标系（Geographic Coordinate Systems）是球面坐标，由经度、纬度和高度组成，参考平面是椭球面。投影坐标系（Projected Coordinate Systems）是平面坐标系，参考平面是水平面。地理坐标系和投影坐标系统称为坐标参考系（Coordinate Reference System，CRS) 。只有当地理坐标系和投影坐标系采用同一个坐标系时，才能正确地将真实的地理坐标投影到地图中。

关于CRS只需要记住两个常用的就可以了：WGS84坐标系和伪墨卡托投影坐标系（Pseudo-Mercator）。

WGS84坐标系：EPSG编号为4326³，也就是以地心基准面，使用地球的质心作为原点，也称为地心大地坐标系，由卫星数据得到，美国GPS系统就是采用此坐标系。
伪墨卡托投影坐标系：EPSG编号为3857，是基于墨卡托投影，把WGS84坐标系投影到正方形上的一种坐标系。伪墨卡托坐标系是擅长显示数据，但不适合存储数据的，因此通常使用WGS84存储数据，使用伪墨卡托显示数据。伪墨卡托投影最早由Google提出，并在其地图产品中广泛使用，目前绝大多数的在线地图均采用伪墨卡托投影。但是，凡是都有例外，显示中国境内地理信息的在线地图使用的是GCJ02经纬度投影坐标系，包括Google Map，因为Google Map中国大陆地区的数据是由高德提供的。这样不影响日常使用，但是当通过利用高德等提供的坐标编码等功能获取某地的坐标时，如全国的医院，得到的经纬度不是基于WGS84而是GCJ-02。

GCJ-02是由中国国家测绘局（国测局）制订的地理信息系统的坐标系统，本质是对WGS84坐标系统进行了加偏处理，在WGS84坐标基础上随机加上偏倚量，并且这个偏移量为随机的，各地的偏移情况都会有所不同，广大网友戏称“火星坐标系统”。

4.3.6.7 关于在R中设置坐标系

一般可以使用sf包中的st_crs()函数，sf是R中开展空间分析的重要工具，对于空间对象的处理可非常丰富的函数。在读取完地图素材后，第一步就应该是查看其坐标参考系。例如，上文使用的china_map素材的自带CRS就是WGS84。

sf::st_crs(china_map)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

如果需要修改CRS，那么可以使用st_transform()函数，如下：

china_map_3875 = sf::st_transform(china_map, 3857)
sf::st_crs(china_map_3875)

Coordinate Reference System:
  User input: EPSG:3857 
  wkt:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]

而对于还未设置CRS的空间对象，可以使用st_crs()函数来设置，如下：

waterfalls <- 
  data.frame("name" = c("Iguazu Falls", "Niagara Falls",
                       "Victoria Falls"), 
             "lat" = c(-25.686785, 43.092461, -17.931805), 
             "lon" = c(-54.444981, -79.047150, 25.825558))

waterfalls_sf <- sf::st_as_sf(waterfalls, 
                              coords = c("lon", "lat"))

sf::st_crs(waterfalls_sf) <- 3875

此外，在geom_polygon和geom_map方法中，也可使用coord_map()函数指定坐标投影系，第三种可以使用coord_sf()函数，其中有参数crs用于指定投影方式，需要使用proj4string格式的字符串。如墨卡托投影用crs = “+proj=merc”。某些地图投影其字符串中包含经纬度等参数，如crs = “+proj=laea +lat_0=35 +lon_0=-100”。其它地图投影见Projection methods。crs本身包含坐标数据的，就不能再增加xlim, ylim等参数。

小结

以上示例展示的只是较为常见的统计图形可视化方法，所举示例也是相对比较简单，只是针对数据分析过程中较为常见的场景，并未能覆盖大部分乃至全部需求。数据可视化是一个充满了创新和创意的领域，读者可以根据自身需求和兴趣大胆尝试和探索。
如果你对基于ggplot2包的空间可视化感兴趣，Timo Grossenbacher的Blog也许对你有帮助。
没有画不出好看图的工具，只有画不出好看图的人，创意和探索非常重要。

数据清洗和探索性分析时通常需要快速地掌握数据的分布情况，如果逐个变量进行分析就会非常耗时。↩︎
这是很常见的一种做法，就如男和女在数据库中通常以整数型的1和0存储。↩︎
EPSG是欧洲石油调查组织(European Petroleum Survey Group)的缩写，其负责专门维护地球上所有的测量坐标系统，并且给每组坐标系统都赋予了一个编号和一组富文本描述（如”+proj=merc”）。↩︎