1 数据的读入与处理

代码提供 and 整理： 曾子尧 黄国政

1.1 数据读入：以.dta格式为例

数据读入包一般有haven、foreign和readstata13三种。包的详情可参见此处。

国内各大型数据库基本都有stata的.dta数据格式文件，我们这里以CGSS2021为案例，主要介绍haven包和readstata13包如何读入.dta数据格式的文件。其中，haven包的读入代码为read_dta()或read_stata()，readstata13包的读入代码为read.dta13()。

read_dta()和read_stata()的读入效果没有差别，我们选取read_dta()和readstata13()进行比较。

library(haven)
CGSS2021a <- read_dta("D://04_Rrender/Bookdown/2020sRBook/2020sRBook1/CGSS2021.dta")
library(readstata13)
CGSS2021c <- read.dta13("D://04_Rrender/Bookdown/2020sRBook/2020sRBook1/CGSS2021.dta")
## Warning in read.dta13("D://04_Rrender/Bookdown/2020sRBook/2020sRBook1/CGSS2021.dta"): 
##    Duplicated factor levels for variables
## 
##    XA1, XA6, XB1, XB6, XC1, XC6, isco08_a57e,
##    isco08_a59d, isco08_a60d, isco08_sp,
##    isco08_f, isco08_m
## 
##    Unique labels for these variables have been generated.
## Warning in read.dta13("D://04_Rrender/Bookdown/2020sRBook/2020sRBook1/CGSS2021.dta"): 
##    Factor codes of type double or float detected in variables
## 
##    type, A8a, A8b, L2, L4b, L5, L8a, L8b, L9,
##    A60m, A62, A75a, A75b
## 
##    No labels have been assigned.
##    Set option 'nonint.factors = TRUE' to assign labels anyway.
## Warning in read.dta13("D://04_Rrender/Bookdown/2020sRBook/2020sRBook1/CGSS2021.dta"): 
##    Missing factor labels for variables
## 
##    P18
## 
##    No labels have been assigned.
##    Set option 'generate.factors=TRUE' to generate labels.

typeof()函数和class()函数可用于检查数据集或列中的数据类型。比较观察两种读入代码对变量数据类型的影响有什么不同。

class(CGSS2021a$A2)
## [1] "haven_labelled" "vctrs_vctr"     "double"
class(CGSS2021c$A2)
## [1] "factor"

可以发现，stata的byte类型变量经由read_dta()读入后数据类型为double，同时由于haven包的影响，数据类型兼为haven_labelled和vctrs_vctr。前一种数据类型保留了数值型的性质，后两种数据类型则便于转换将数据为因子型，但有时也给数据分析带来一些问题，关于这一点，我们在后续分析时再具体进行讨论。

我们特别关心变量的类型，原因是变量的类型规定了可以对该变量执行的统计操作的类型。

概览数据库，可总结出两种语法存在以下两处明显的区别：

  1. read_dta读入了对应的变量标签内容和变量选项编码数值，而read.dta13没有对应的变量标签内容，且直接读出变量编码数值的内容，视觉上呈现了选项的含义。
  2. read_dta读入的数据是数值型的，而read.dta13读入的数据是因子型的。

我们认为，虽然read.dta13()会带来视觉上的方便，但是read_dta()和read_stata()在重新编码上更具有优势。

1.2 数据处理1

数据读入完毕。由于CGSS2021数据库庞大且复杂，为便于分析，我们选取关注的变量并创建小型的数据集，再进一步进行数据的处理和诊断，最后再进行清洗。根据CGSS2021年的调查问卷，我们选取以下几个变量创建一级数据集cgss2021（相对原始数据集）。

基本情况

  1. A27d 出生户口登记地
  2. A2 性别
  3. A3_1出生日期
  4. A13 身高
  5. A14 体重
  6. A7a 目前最高教育程度
  7. A10 政治面貌
  8. A8a 个人2020年总收入
  9. A62 2020年全年家庭总收入
  10. A15 身体健康感受
  11. A36 生活幸福感
cgss2021 <- CGSS2021a %>% select(
    stype = A27d,
    gender = A2,
    birth = A3_1,
    height = A13,
    weight = A14,
    edu = A7a,
    poli = A10,
    income = A8a,
    famincome = A62,
    health = A15,
    happy = A36
  )

使用mutate()函数，可以转换旧变量birth来增加新变量age，得出样本的具体年龄。

cgss2021 <- cgss2021 %>% mutate(age = 2021-birth)

使用head()函数，可以查看数据集的前6行，观察增加的新变量age。

library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
head(cgss2021) %>% kable("html") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "300px")
stype gender birth height weight edu poli income famincome health happy age
5 2 1967 160 118 7 1 96000 200000 5 3 54
5 2 1945 160 120 7 4 60000 120000 5 5 76
1 1 1949 172 156 11 4 96000 108000 4 NA 72
1 1 1986 172 130 11 1 160000 220000 4 4 35
5 2 1958 158 115 7 1 48000 100000 4 4 63
5 1 1996 180 150 9 1 600000 2000000 4 4 25

age属于连续型变量，对连续型变量的重新编码同样可以使用mutate()函数。

cgss2021 <- cgss2021 %>% mutate(agegroup = case_when(
  age <= 30 ~ "30岁以下",
  age > 30 & age <= 45 ~ "30~45岁",
  age > 45 & age <= 60 ~ "45~60岁",
  age > 60 ~ "60岁以上"
))

使用typeof()函数查看agegroup变量类型，会发现已转化为离散变量中的字符型变量。

typeof(cgss2021$agegroup)
## [1] "character"

使用factor()函数将agegroup转化为离散变量中的因子型变量。因子型变量属于类别，可以以数字的形式呈现，也可以以文本的形式呈现。因子型比字符型便于分析。

cgss2021$agegroup <- factor(cgss2021$agegroup,
                              levels = c("30岁以下","30~45岁","45~60岁","60岁以上"))

1.2.1 R中的变量类型

前面提及到了R中的连续变量和离散变量等概念，在对数据进行处理之前，需要对R中可用的变量类型进行介绍。其中，R的主要变量类型有两种：离散变量和连续变量。

1.2.1.1 离散变量

离散变量主要由可数的数字集合组成，例如从1到N的计数，或类别变量如女性和男性等，包括计数变量（整数）、顺序变量和类别变量，后两者又被称为分类变量。

1.2.1.2 连续变量

连续变量可以测量到任意精度，在刻度上的任意两个值之间加上额外的数。包括计数变量和数值变量。两者又被称为定量变量。

1.2.1.3 R的变量储存类型

导入R中的数据储存格式有如下几种：double、integer、factor、character和logical、complex和raw。haven包导入的数据中，若变量的数据为连续型，则存储格式不变；若变量的数据为离散型，则存储为haven_labelled格式。 double可以理解为连续变量中的数值变量，interger属于离散变量或连续变量中的计数变量，factor和character都属于离散变量中的类别变量。在R中，我们使用数字、整数、有序和因子分别表示数值变量、计数变量、顺序变量和类别变量。通常，我们可以使用class()或typeof()函数进行检验。

1.3 数据处理2

接下来我们对”最高教育程度”edu变量进行再编码，先查看变量类型。

typeof(cgss2021$edu)
## [1] "double"

虽然edu显示为double（双精度），但实际上其中的值可数，相当于整数型，我们还是采取对离散型变量重新编码的方式（这同样适用于字符型和因子型）。

cgss2021 <- cgss2021 %>% mutate(educ = case_when(
  edu == 1 ~ "未受教育",
  edu == 2 | edu == 3 ~ "初等教育",
  edu == 4 | edu == 5 | edu == 6 | edu == 7 | edu == 8 ~ "中等教育",
  edu == 9 | edu == 10 | edu == 11 | edu == 12 | edu == 13 ~ "高等教育",
  edu == 14 ~ "其他"
))
#显然，对字符型或文本形式的因子型进行重新编码比较麻烦。
#cgss2021<- cgss2021 %>% mutate(educ = case_when(
#    edu == "没有受过任何教育" ~ "未受教育",
#    edu == "私塾、扫盲班"| edu == "小学" ~ "初等教育",
#    edu == "初中"| edu == "职业高中"| edu == "普通高中"| edu == "中专"| edu == "技校" ~ "中等教育",
#    edu == "大学专科（成人高等教育）"| edu == "大学专科（正规高等教育）"| edu == "大学本科（成人高等教育）"| #edu == "大学本科（正规高等教育）"| edu == "研究生及以上" ~ "高等教育",
#    edu == "其他（请注明）" ~ "其他"
#  ))
cgss2021$educ <- factor(cgss2021$educ,
                        levels = c("未受教育","初等教育","中等教育","高等教育","其他"))

另一种编码方式。

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
cgss2021$stype <- as.numeric(cgss2021$stype)
typeof(cgss2021$stype)
## [1] "double"
cgss2021$stype <- 
  car::recode(cgss2021$stype, "1='1';2:5='2';6:99='99'")
cgss2021$stype <- 
  factor(cgss2021$stype,
         levels = c(1,2),
         labels = c("农村","城镇"))

1.4 数据处理3

建立好一级数据集，并对部分变量进行再编码后，接下来进行异常值、缺失值和重复值的诊断和处理。

1.4.1 诊断数据集

CGSS2021中往往存在很多异常值和缺失值。异常值是指CGSS问卷中数值特别大的选项，往往指”不知道”、“拒绝回答”，或者是填了不合理的数据，比如特别高或低的体重、身高、收入等等。缺失值往往指受访者没有回答的问题。诊断数据集的目的是为了发现数据集的异常状况和缺失状况，然后再进行下一步的处理。

我们先发现异常值，而后转化为缺失值，后续根据分析需要再将其与原有的缺失值一道删除。此处的数据集中变量较多，不适合立即删除所有的缺失值，以免影响其他变量。

ExpData()函数提供一份数据结构的分布情况，包含缺失值、唯一值等。

library(SmartEDA)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
res<- ExpData(data=cgss2021,type=2)

根据res，我们可以直接得到各变量的缺失值数量，但我们更关心的是潜在的异常值。根据变量的不同类型，我们对异常值的勘察方式也有所不同。

  1. 离散变量 对照离散变量的distincts值和问卷中的选项数量。distincts值是变量唯一取值的数量，包括NA。若distincts值大于问卷选项数量，说明存在异常值。
  2. 连续变量 连续变量一般查看其取值范围。

离散变量：

table(cgss2021$stype)
## 
## 农村 城镇 
## 5931 2163
table(cgss2021$gender)
## 
##    1    2 
## 3679 4469
table(cgss2021$edu)
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13 
##  858   47 1751 2311  127  960  374   28  261  361  194  738  117
table(cgss2021$poli)
## 
##    1    2    3    4   98   99 
## 6557  606    6  966   12    1
table(cgss2021$health)
## 
##    1    2    3    4    5   98 
##  429 1069 2288 2864 1492    6
table(cgss2021$happy)
## 
##    1    2    3    4    5   98   99 
##   75  243  698 3113 1319    7    3
table(cgss2021$agegroup)
## 
## 30岁以下  30~45岁  45~60岁 60岁以上 
##     1232     1684     2409     2823
table(cgss2021$educ)
## 
## 未受教育 初等教育 中等教育 高等教育     其他 
##      858     1798     3800     1671        0

连续变量：

library(sjmisc)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
## 
## Attaching package: 'sjmisc'
## The following object is masked from 'package:purrr':
## 
##     is_empty
## The following object is masked from 'package:tidyr':
## 
##     replace_na
## The following object is masked from 'package:tibble':
## 
##     add_case
descr(cgss2021$birth)
## 
## ## Basic descriptive statistics
## 
##  var    type                          label    n NA.prc    mean    sd   se   md
##   dd numeric [年]  A3. 您的出生日期是什么？ 8148      0 1969.36 17.57 0.19 1968
##  trimmed          range iqr skew
##  1969.06 81 (1922-2003)  29 0.14
descr(cgss2021$age)
## 
## ## Basic descriptive statistics
## 
##  var    type                          label    n NA.prc  mean    sd   se md
##   dd numeric [年]  A3. 您的出生日期是什么？ 8148      0 51.64 17.57 0.19 53
##  trimmed      range iqr  skew
##    51.94 81 (18-99)  29 -0.14
descr(cgss2021$height)
## 
## ## Basic descriptive statistics
## 
##  var    type                 label    n NA.prc   mean     sd   se  md trimmed
##   dd numeric A13. 您目前的身高是： 8148      0 178.46 110.25 1.22 163  178.46
##         range iqr skew
##  901 (98-999)  12 7.26
descr(cgss2021$weight)
## 
## ## Basic descriptive statistics
## 
##  var    type                 label    n NA.prc   mean    sd   se  md trimmed
##   dd numeric A14. 您目前的体重是： 8148      0 130.77 81.85 0.91 120  130.77
##         range iqr skew
##  959 (40-999)  31 9.64
descr(cgss2021$income)
## 
## ## Basic descriptive statistics
## 
##  var    type                              label    n NA.prc    mean      sd
##   dd numeric 您个人去年（2020年）全年的总收入是 8148      0 1046361 3001742
##        se    md trimmed               range   iqr skew
##  33254.31 30000 1046361 9999999 (0-9999999) 67000 2.64
descr(cgss2021$famincome)
## 
## ## Basic descriptive statistics
## 
##  var    type                          label    n NA.prc    mean      sd
##   dd numeric 您家2020年全年家庭总收入是多少 8148      0 2216411 4066856
##        se    md trimmed               range    iqr skew
##  45053.99 80000 2216411 9999999 (0-9999999) 270000 1.39

可以发现，变量height、weight、income和famincome都有异常的极大值。

1.4.2 处理异常值和缺失值

我们往往选择将异常值转化为缺失值，而后进行删除。

初始处理离散变量异常值的时候，我们的思路或许是将cgss2021中指定的值统一全部删掉。

cgss2021[cgss2021 == 98] <- NA
cgss2021[cgss2021 == 99] <- NA

这种方法用来处理离散变量的异常值很常用，但可能会损失连续变量中具有意义的值。例如98在“户口登记地”中是“不知道”，但在“体重”中是有意义的。

如果一一指定，似乎可以解决误删导致的信息损失问题。

cgss2021$stype[cgss2021$stype == 98] <- NA
cgss2021$stype[cgss2021$stype == 99] <- NA

但如果一个变量中存在多个异常值，一一指定会十分麻烦。

处理连续变量的异常值所要遭遇的问题与以上相近。

cgss2021[cgss2021 >= 999] <- NA

例如999在“身高”中指“拒绝回答”，但这样处理会使得收入中的大多数值也转换为NA。当然，我们可以具体指定只删除变量“身高”中的999以及999以上的数值，但当需要分析的连续变量较多时，也会十分麻烦。

综上，我们会发现通过指定数值来转换异常值通常有两个问题：

  1. 将其他变量中的相关有效数值转换为无效数值。
  2. 当特定变量的异常值比较多时，一一具体指定会十分麻烦。

沿着这个问题，我想到的第一个可能的解决方案是搭配使用sjlabelled和sjmisc两个包。

library(sjlabelled)
library(sjmisc)
cgss <- sjlabelled::drop_labels(cgss2021)
head(cgss)
cgss <- sjmisc::to_label(cgss)
head(cgss)

这两个包可以将haven_labelled类型变量转换为以文本形式显示的因子型变量，此时可以通过统一定义特定文本来转换异常值，例如：

cgss2021[cgss2021==“不知道”] <- NA
cgss2021[cgss2021==“拒绝回答”] <- NA

如此转换将更加精准。

但haven_labelled类型变量1中的纯数值会转换为NA，只有带有标签的数值会转换为文本显示（dbl类型的变量不会受到影响）

下面我们将数据集分为离散型和连续型子数据集来分别处理，最后再进行合并。

library(sjlabelled)
## 
## Attaching package: 'sjlabelled'
## The following objects are masked from 'package:haven':
## 
##     as_factor, read_sas, read_spss, read_stata, write_sas, zap_labels
## The following object is masked from 'package:forcats':
## 
##     as_factor
## The following object is masked from 'package:dplyr':
## 
##     as_label
## The following object is masked from 'package:ggplot2':
## 
##     as_label
library(sjmisc)

建立离散变量数据集cgss1

cgss1 <- cgss2021 %>% select(
  stype = stype,
  gender = gender,
  agegroup = agegroup,
  edu = edu,
  educ = educ,
  poli = poli,
  health = health,
  happy = happy
)

head(cgss1)
## # A tibble: 6 × 8
##   stype gender    agegroup edu                    educ  poli    health  happy   
##   <fct> <dbl+lbl> <fct>    <dbl+lbl>              <fct> <dbl+l> <dbl+l> <dbl+lb>
## 1 城镇  2 [女]    45~60岁   7 [中专]            … 中等… 1 [群… 5 [很…  3 [说…
## 2 城镇  2 [女]    60岁以上  7 [中专]            … 中等… 4 [共… 5 [很…  5 [非…
## 3 农村  1 [男]    60岁以上 11 [大学本科（成人高… 高等… 4 [共… 4 [比… NA     …
## 4 农村  1 [男]    30~45岁  11 [大学本科（成人高… 高等… 1 [群… 4 [比…  4 [比…
## 5 城镇  2 [女]    60岁以上  7 [中专]            … 中等… 1 [群… 4 [比…  4 [比…
## 6 城镇  1 [男]    30岁以下  9 [大学专科（成人高… 高等… 1 [群… 4 [比…  4 [比…

处理离散变量数据集cgss1的异常值

cgss1[cgss1 == 98] <- NA 
cgss1[cgss1 == 99] <- NA

cgss1 <- sjmisc::to_label(cgss1)

head(cgss1)
## # A tibble: 6 × 8
##   stype gender agegroup edu                      educ     poli     health  happy
##   <fct> <fct>  <fct>    <fct>                    <fct>    <fct>    <fct>   <fct>
## 1 城镇  女     45~60岁  中专                     中等教育 群众     很健康  说不…
## 2 城镇  女     60岁以上 中专                     中等教育 共产党员 很健康  非常…
## 3 农村  男     60岁以上 大学本科（成人高等教育） 高等教育 共产党员 比较健… <NA> 
## 4 农村  男     30~45岁  大学本科（成人高等教育） 高等教育 群众     比较健… 比较…
## 5 城镇  女     60岁以上 中专                     中等教育 群众     比较健… 比较…
## 6 城镇  男     30岁以下 大学专科（成人高等教育） 高等教育 群众     比较健… 比较…

!!!注意：本次分析的离散变量每个值都带有标签意义。但存在一种情况，例如询问满意度分值，分为1-5分，同时1分具有“非常不赞成”的标签含义，5分具有“非常赞成”的标签含义，而中间的2-4分不具有任何标签含义。此时若运行cgss1 <- sjmisc::to_label(cgss1)，会将2-4分值转换为NA。解决的方法是直接用命令factor()将这一类变量转换为因 子型变量，无需给“98”和“99”定义标签和排序，它们会自动被归为NA。例如2：

cgss1$x <- factor(cgss1$x,
                        labels = c(1,2,3,4,5),           
                        levels = c("1","2","3","4","5"))

处理以后，haven_lablled类型变量已全体转换因子型变量，并以文本的形式呈现。根据研究需要，我们后续也可以将文本形式的因子转换为数字形式的因子，或者将因子型变量转换为数值型变量。

建立连续变量数据集cgss2

cgss2 <- cgss2021 %>% select(
  age = age,
  height = height,
  weight = weight,
  income = income,
  famincome = famincome
)

head(cgss2)
## # A tibble: 6 × 5
##     age height    weight    income    famincome
##   <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl+lbl>
## 1    54 160       118        96000     200000  
## 2    76 160       120        60000     120000  
## 3    72 172       156        96000     108000  
## 4    35 172       130       160000     220000  
## 5    63 158       115        48000     100000  
## 6    25 180       150       600000    2000000

处理连续变量数据集cgss2的异常值

cgss2[cgss2 > 9999996] <- NA
cgss2$height[cgss2$height >= 998] <- NA
cgss2$weight[cgss2$weight >= 998] <- NA

cgss2 <- sjlabelled::drop_labels(cgss2)
cgss2 <- sjlabelled::zap_labels(cgss2)
## Warning: `na` is not allowed to be `NULL` or to contain `NA`-values.

## Warning: `na` is not allowed to be `NULL` or to contain `NA`-values.

## Warning: `na` is not allowed to be `NULL` or to contain `NA`-values.
head(cgss2)
## # A tibble: 6 × 5
##     age height weight income famincome
##   <dbl>  <dbl>  <dbl>  <dbl>     <dbl>
## 1    54    160    118  96000    200000
## 2    76    160    120  60000    120000
## 3    72    172    156  96000    108000
## 4    35    172    130 160000    220000
## 5    63    158    115  48000    100000
## 6    25    180    150 600000   2000000

可以看到，haven_labelled类型的变量全部转换为了dal类型的变量。

我们保留cgss1中变量的因子类型，合并cgss1和cgss2。

cbind()函数是按行合并。

cgss3 <- cbind(cgss1,cgss2)

使用select函数对数据集中的变量重新进行排序。everything()函数相当于一个占位符，用于标明数据集中之前未被选中的其他所有列。

cgss3 <- select(cgss3,
                stype,
                gender,
                age,
                agegroup,
                height,
                weight,
                edu,
                educ,
                poli,
                income,
                famincome,
                everything(),
                )
cgss3 <- as_tibble(cgss3)

tibble和data.frame各有长处，但tibble相对而言更稳定，例如从不改变输入的类型（不将字符串转换为因子），也从不改变变量的名称，也不创建行名。同时，tibble的打印只显示数据集的前10行（这使得处理大型数据更加容易），以及所有适合在屏幕上显示的列，除了列的名字，每一列还报告它的类型，相当于str()的功能。

如果遇上旧的函数不能与tibble工作，可以用as.data.frame()将tibble转换回data.frame的格式，转换的代码：as.data.frame(cgss3)

再进行一次数据诊断，查看缺失值状况。

library(SmartEDA)
res3<- ExpData(data=cgss3,type=2)

save()函数可以单独保存包含中间处理步骤的数据集。将新的数据集另存为新的.Rdata数据文件（由于.Rdata数据文件的内容会被加载并显示在Rstudio中，因此建议数据名与数据集的名称保持一致，避免数据混淆），这样在实际操作中，我们不必重复运行以上数据处理的代码便可导入我们已经处理好的数据集。

save(cgss3, file = "cgss3.Rdata")

至此，异常值清理完成，数据的处理初步完成。但由于缺失值较多，后续的分析仍然需要从cgss3中截取部分数据，建立新的数据集进行分析。

接下来，让我们基于由CGSS2021截取的小型数据集cgss3展开数据分析之旅吧！

  1. haven_labelled变量是指未明确类型的变量，往往同时含有没有标签的数值和具有标签的数值，即既体现出数值型变量的特性，也体现着因子型的变量，因而被认为是未明确类型的。↩︎

  2. 这里假设x是一个分值为1-5的变量，其中1为“非常不赞同”，5为“非常赞同”，2-4只具有数值意义，无标签意义。↩︎