1 😄加载经常用的R包😄

library(pacman)
# 读数据
p_load(readxl,writexl,data.table,openxlsx,haven,rvest)
# 数据探索
p_load(tidyverse,DT,skimr,DataExplorer,explore,vtable,stringr,kableExtra,lubridate)
# 模型
p_load(grf,glmnet,caret,tidytext,fpp2,forecast,car,tseries,hdm,tidymodels,broom)
# 可视化
p_load(patchwork,ggrepel,ggcorrplot,gghighlight,ggthemes,shiny)
# 其它常用包
p_load(magrittr,listviewer,devtools,here,janitor,reticulate,jsonlite)

2 😄读入数据😄

第一步是读入数据，对数据进行初步了解。下面以小说《三生三世十里桃花》中的人物信息为背景，具体的变量解释表如表1所示，其中因变量Y为“决定”这个变量。

knitr::include_graphics(here::here("Machine_Learning_and_Causal_Inference/fig/640.png"))

data <- read.csv(here::here("Machine_Learning_and_Causal_Inference/data/相亲数据2.csv"),fileEncoding = "UTF-8")
data
data %>% str()
## 'data.frame':    15 obs. of  12 variables:
##  $ 名字              : chr  "夜华" "白浅" "东华" "素锦" ...
##  $ 决定              : int  1 0 0 0 1 0 0 0 1 1 ...
##  $ 性别              : int  1 0 1 0 1 0 1 1 1 0 ...
##  $ 是否喜欢矮矬穷    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ 年龄              : int  50000 140000 360000 70000 140000 300 370000 23 31 23 ...
##  $ 智力              : int  9 8 10 5 6 6 10 9 9 7 ...
##  $ 对方名字          : chr  "白浅" "离镜" "凤九" "离镜" ...
##  $ 对方决定          : int  1 1 1 0 1 1 0 1 1 1 ...
##  $ 对方性别          : int  0 1 0 1 0 1 0 0 1 1 ...
##  $ 对方是否喜欢矮矬穷: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ 对方年龄          : int  140000 140000 30000 140000 100000 90000 70000 20 31 30 ...
##  $ 对方智力          : int  8 6 5 6 4 7 5 7 7 8 ...
data %>% glimpse()
## Rows: 15
## Columns: 12
## $ 名字               <chr> "夜华", "白浅", "东华", "素锦", "离镜", "成玉", "折颜", "宇文玥", "...
## $ 决定               <int> 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0
## $ 性别               <int> 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1
## $ 是否喜欢矮矬穷     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## $ 年龄               <int> 50000, 140000, 360000, 70000, 140000, 300, 370000,...
## $ 智力               <int> 9, 8, 10, 5, 6, 6, 10, 9, 9, 7, 8, 5, 8, NA, 9
## $ 对方名字           <chr> "白浅", "离镜", "凤九", "离镜", "玄女", "连宋", "素锦", "大梁公主", "靖...
## $ 对方决定           <int> 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1
## $ 对方性别           <int> 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0
## $ 对方是否喜欢矮矬穷 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## $ 对方年龄           <int> 140000, 140000, 30000, 140000, 100000, 90000, 70000,...
## $ 对方智力           <int> 8, 6, 5, 6, 4, 7, 5, 7, 7, 8, 6, 9, 8, 8, 6
data %>% skimr::skim()
Table 2.1: Data summary
Name Piped data
Number of rows 15
Number of columns 12
_______________________
Column type frequency:
character 2
numeric 10
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
名字 0 1 2 3 0 15 0
对方名字 0 1 2 4 0 14 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
决定 0 1.00 0.47 0.52 0 0.00 0 1 1 ▇▁▁▁▇
性别 0 1.00 0.53 0.52 0 0.00 1 1 1 ▇▁▁▁▇
是否喜欢矮矬穷 0 1.00 0.00 0.00 0 0.00 0 0 0 ▁▁▇▁▁
年龄 0 1.00 75447.73 127536.03 21 33.00 300 105000 370000 ▇▂▁▁▂
智力 1 0.93 7.79 1.72 5 6.25 8 9 10 ▇▂▆▇▃
对方决定 0 1.00 0.67 0.49 0 0.00 1 1 1 ▃▁▁▁▇
对方性别 0 1.00 0.53 0.52 0 0.00 1 1 1 ▇▁▁▁▇
对方是否喜欢矮矬穷 0 1.00 0.00 0.00 0 0.00 0 0 0 ▁▁▇▁▁
对方年龄 0 1.00 47368.13 59187.35 20 31.00 300 95000 140000 ▇▁▁▂▃
对方智力 0 1.00 6.67 1.40 4 6.00 7 8 9 ▆▇▆▇▂

3 😄分割训练集和测试集😄

数据拿到之后，要先划分训练集和测试集。测试集是黑盒子，是不能触碰的炸弹，所做的任何处理，包括标准化、缺失值填补都只能基于训练集。下面介绍几个典型划分训练集和测试集的方法。

3.1 👏留出法👏

留出法分割是将样本分为两个互斥的子集，通常情况下，划分数据的80%为训练集，剩下的20%为测试集

之前提过，caret包中的createDataPartition()函数不仅可以实现这样的划分，而且可以保证训练集和测试集中Y的比例是一致的，简而言之就是按照Y进行分层抽样。

set.seed(1234)
data_id <- createDataPartition(data$决定,p = 0.7,list = FALSE,times = 1)
data_id
##       Resample1
##  [1,]         2
##  [2,]         4
##  [3,]         5
##  [4,]         6
##  [5,]         8
##  [6,]         9
##  [7,]        10
##  [8,]        11
##  [9,]        12
## [10,]        13
## [11,]        15
data_training <- data %>% slice(data_id)
data_testing <- data %>% slice(-data_id)
data_training
data_testing
data %$% table(决定) %>% prop.table() 
## 决定
##         0         1 
## 0.5333333 0.4666667
data_training %$% table(决定) %>% prop.table() 
## 决定
##         0         1 
## 0.5454545 0.4545455
data_testing %$% table(决定) %>% prop.table() 
## 决定
##   0   1 
## 0.5 0.5

3.2 👏交叉验证法👏

交叉验证法将原始数据分成K组（一般是均分），每次训练将其中一组作为测试集，另外K-1组作为训练集。

实际应用中一般十折交叉验证用得最多，但是这里由于数据量太少，就以3折交叉验证为例（见图2），展示代码如下：

set.seed(1234)
data_id1 <- createFolds(data$决定,k = 3,list = FALSE,returnTrain = TRUE)
data_id1
##  [1] 3 1 3 2 3 3 1 2 1 2 2 2 3 1 1
data_id1 %>% table()
## .
## 1 2 3 
## 5 5 5
data_training <- data[-which(data_id == 1),]
data_testing <- data[which(data_id == 1),]
data_training
data_testing

3.3 👏Bootstrap法👏

当数据量比较少时，Bootstrap抽样会成为“救命稻草”，它是一种从给定训练集中有放回的均匀抽样。也就是说，每当选中一个样本，它依然会被再次选中并被再次添加到训练集中。

createResample()函数中times参数用于设定生成几份随机样本，当times为3，意味着生成3份样本。不仅不同sample之间会有交叉，就连同一份sample中也会有重复的样本。

set.seed(1234)
data_id1 <- createResample(data$决定,times = 3,list = FALSE)
data_id1
##       Resample1 Resample2 Resample3
##  [1,]         2         3         2
##  [2,]         4         4         3
##  [3,]         5         4         4
##  [4,]         5         4         5
##  [5,]         6         4         7
##  [6,]         6         4         8
##  [7,]         6         5         9
##  [8,]         6         8        10
##  [9,]         7         8        11
## [10,]         9         8        12
## [11,]        10        14        13
## [12,]        10        14        14
## [13,]        12        14        15
## [14,]        12        14        15
## [15,]        15        15        15

上面这些划分训练集和测试集的方法都是针对横截面数据而言的，那么对于时间序列又该如何进行数据分割呢？

3.4 👏分割时间序列👏

data_growth <- data.frame(时间 = ymd(20160101) + months(1:6),
                            数值 = 1:6)
data_growth

接下来，利用caret包来分割时间序列。createTimeSlices()函数需要输入以下参数：initialWindow表示第一个训练集中的样本数；horizon参数表示每个测试集中的样本数；fixedWindow参数表示是否每个训练集中的样本数都相同。

从结果可以看出来，一共有2组训练集和测试集。第一组的训练集为1、2、3、4、5行观测，测试集为6、7行观测。那么第二组呢？从下面的数据就可以看出。

data_id <- createTimeSlices(data_growth$数值,initialWindow = 3,horizon = 2,fixedWindow = TRUE)
data_id
## $train
## $train$Training3
## [1] 1 2 3
## 
## $train$Training4
## [1] 2 3 4
## 
## 
## $test
## $test$Testing3
## [1] 4 5
## 
## $test$Testing4
## [1] 5 6
data_id$train
## $Training3
## [1] 1 2 3
## 
## $Training4
## [1] 2 3 4
data_id$test
## $Testing3
## [1] 4 5
## 
## $Testing4
## [1] 5 6
data_id1 <- createTimeSlices(data_growth$数值,initialWindow = 3,horizon = 2,fixedWindow = FALSE)
data_id1
## $train
## $train$Training3
## [1] 1 2 3
## 
## $train$Training4
## [1] 1 2 3 4
## 
## 
## $test
## $test$Testing3
## [1] 4 5
## 
## $test$Testing4
## [1] 5 6

4 😄处理缺失值😄

caret包中preProcess()函数实现了两种常用的缺失值处理方法：中位数填补法、K近邻方法。

4.1 👏中位数填补法👏

该方法直接用训练集的中位数代替缺失值，所以对于每个变量而言，填补的缺失值都相同，为训练集的中位数。该方法的优点是速度非常快，但填补的准确率有待验证。

data <- read.csv(here::here("Machine_Learning_and_Causal_Inference/data/相亲数据2.csv"),fileEncoding = "UTF-8")
data
set.seed(1234)
data_id <- createDataPartition(data$决定,p = 0.7,times = 1,list = FALSE)
data_training <- data[data_id,]
data_testing <- data[-data_id,]
data_training
data_testing

训练集中所有人的智力中位数值是不是7.5?

data_training$智力 %>% median()
## [1] 8
data_training %>% datatable()
data_training %>% is.na() %>% sum()
## [1] 0
data_testing %>% datatable()
data_testing %>% is.na() %>% sum()
## [1] 1
data_imputation_median <- preProcess(data_training,method = "medianImpute")
data_imputation_median
## Created from 11 samples and 12 variables
## 
## Pre-processing:
##   - ignored (2)
##   - median imputation (10)
data_training_imputation <- predict(data_imputation_median,data_training)


all.equal(data_training,data_training_imputation)  # 知道为什么是true吗？
## [1] TRUE
data_testing_imputation <- predict(data_imputation_median,data_testing)
data_testing_imputation

4.2 👏K近邻法👏

该方法的思想是“近朱者赤近墨者黑”。K近邻法对于需要插值的记录，基于欧氏距离计算k个和它最近的观测，然后接着利用k个近邻的数据来填补缺失值。

K近邻法会自动利用训练集的均值标准差信息对数据进行标准化，所以最后得到的数据是标准化之后的。如果想看原始值，那么还需要将其去标准化倒推回来。

data_imputation_knn <- preProcess(data_training,method = "knnImpute")
## Warning in preProcess.default(data_training, method = "knnImpute"): These
## variables have zero variances: 是否喜欢矮矬穷, 对方是否喜欢矮矬穷
data_imputation_knn
## Created from 11 samples and 12 variables
## 
## Pre-processing:
##   - centered (10)
##   - ignored (2)
##   - 5 nearest neighbor imputation (10)
##   - scaled (10)
data_training_imputation <- predict(data_imputation_knn,data_training)


all.equal(data_training,data_training_imputation)  # 知道为什么不是true吗？
## [1] "Component \"决定\": Mean relative difference: 1.088932"
## [2] "Component \"性别\": Mean relative difference: 1.088932"
## [3] "Component \"年龄\": Mean relative difference: 1"       
## [4] "Component \"智力\": Mean relative difference: 1"       
## [5] "Component \"对方决定\": Mean relative difference: 1"   
## [6] "Component \"对方性别\": Mean relative difference: 1"   
## [7] "Component \"对方年龄\": Mean relative difference: 1"   
## [8] "Component \"对方智力\": Mean relative difference: 1"
data_training$智力 <- data_training_imputation$智力 * sd(data_training$智力,na.rm = TRUE) +
  mean(data_training$智力,na.rm = TRUE)
data_testing_imputation <- predict(data_imputation_knn,data_testing)
data_testing_imputation
data_testing$智力 <- data_testing_imputation$智力 * sd(data_training$智力,na.rm = TRUE) + mean(data_training$智力,na.rm = TRUE)

data_testing

5 😄删除近零方差😄

零方差或者近零方差的变量传递不了什么信息，因为几乎所有人的取值都一样。可以利用caret包中的nearZeroVar()函数，一行代码就能找出近零方差的变量，操作过程非常简单。

nearZeroVar(data_training) # 4 和 10
## [1]  4 10
data_training_dropvariable <- data_training[,-nearZeroVar(data_training)]
data_testing_dropvariable <- data_testing[,-nearZeroVar(data_training)]
data_training_dropvariable
data_testing_dropvariable %>% str()
## 'data.frame':    4 obs. of  10 variables:
##  $ 名字    : chr  "夜华" "东华" "折颜" "罗子君"
##  $ 决定    : int  1 0 0 1
##  $ 性别    : int  1 1 1 0
##  $ 年龄    : int  50000 360000 370000 35
##  $ 智力    : num  9 10 10 7
##  $ 对方名字: chr  "白浅" "凤九" "素锦" "贺函"
##  $ 对方决定: int  1 1 0 1
##  $ 对方性别: int  0 0 0 1
##  $ 对方年龄: int  140000 30000 70000 40
##  $ 对方智力: int  8 5 5 8

6 😄删除共线性变量😄

caret包中的findCorrelation()函数会自动找到高度共线性的变量，并给出建议删除的变量。

但需要注意，这个函数对输入的数据要求比较高：

  • 首先，数据中不能有缺失值，所以在此之前需要先处理缺失值；
  • 其次，只能包含数值型变量。
data_training %<>% drop_na()
data_training %>% 
  select(-nearZeroVar(data_training)) %>% 
  select(where(is.numeric)) %>% 
  cor() -> data_cor
data_cor
##                 决定        性别        年龄        智力    对方决定
## 决定      1.00000000  0.26666667 -0.06595876 -0.16791512 -0.06900656
## 性别      0.26666667  1.00000000 -0.06198019  0.57091142  0.31052950
## 年龄     -0.06595876 -0.06198019  1.00000000 -0.26212272  0.19983184
## 智力     -0.16791512  0.57091142 -0.26212272  1.00000000  0.39396631
## 对方决定 -0.06900656  0.31052950  0.19983184  0.39396631  1.00000000
## 对方性别 -0.06900656 -0.82807867 -0.04628580 -0.37079182 -0.17857143
## 对方年龄 -0.35717624 -0.35636285  0.81215697 -0.48262196  0.10071829
## 对方智力  0.33565856 -0.23237900 -0.69247720 -0.05636215 -0.30735043
##             对方性别   对方年龄    对方智力
## 决定     -0.06900656 -0.3571762  0.33565856
## 性别     -0.82807867 -0.3563629 -0.23237900
## 年龄     -0.04628580  0.8121570 -0.69247720
## 智力     -0.37079182 -0.4826220 -0.05636215
## 对方决定 -0.17857143  0.1007183 -0.30735043
## 对方性别  1.00000000  0.2299345  0.28062430
## 对方年龄  0.22993447  1.0000000 -0.54965318
## 对方智力  0.28062430 -0.5496532  1.00000000
data_high_cor <- findCorrelation(data_cor,cutoff = 0.75,verbose = TRUE,names = TRUE)
## Compare row 7  and column  3 with corr  0.812 
##   Means:  0.413 vs 0.31 so flagging column 7 
## Compare row 2  and column  6 with corr  0.828 
##   Means:  0.378 vs 0.271 so flagging column 2 
## All correlations <= 0.75
data_high_cor
## [1] "对方年龄" "性别"
data_training %<>% select(-data_high_cor)
## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(data_high_cor)` instead of `data_high_cor` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
data_testing %<>% select(-data_high_cor)
data_training
data_testing

7 😄标准化😄

为什么要标准化？很简单，看看年龄，几十万岁，但是智力这个变量最高也才10分，这两列变量的量纲不同，为了防止年龄的权重过高，就需要将这些特征进行标准化才能学习各个变量真实的权重。需要注意的是：只能拿训练集的均值和标准差来对测试集进行标准化。

data_proprocess_value <- preProcess(data_training,method = c("scale","center"))
## Warning in preProcess.default(data_training, method = c("scale", "center")):
## These variables have zero variances: 是否喜欢矮矬穷, 对方是否喜欢矮矬穷
data_training_fin <- predict(data_proprocess_value,data_training)
data_testing_fin <- predict(data_proprocess_value,data_testing)
data_training
data_testing