1 分析前的环境设置

1.1 下载安装套件包

本章学习者須先下载安装相關套件包，方能正常運行。

#下载安装套件包
#如果之前有下载过，会跳过不下载；如果之前没有下载过，会下载安装套件包
if (!require("showtext")) install.packages("showtext")
if (!require("anesrake")) install.packages("anesrake")
if (!require("sjlabelled")) install.packages("sjlabelled")
if (!require("sjmisc")) install.packages("sjmisc")
if (!require("sjPlot")) install.packages("sjPlot")
if (!require("sjstats")) install.packages("sjstats")

1.2 环境设置

我们通常会将R的工作环境先进行一些设置，这些设置包括了 “设定当前工作目录”、“设定系統中文文字编码”、以及“设定绘图物件中的中文文字”等项。

#设置工作目录(此处要改成你自己的工作目录)
setwd("/Users/simonfair/Desktop/量化研究数据分析/R Ch5 样本代表性与加权")

getwd() #查看目前的工作目录

[1] "/Users/simonfair/Desktop/量化研究数据分析/R Ch5 样本代表性与加权"

#设置系統中文文字编码(简体中文)
Sys.setlocale(category = "LC_ALL", locale = "zh_CN.UTF-8") #

[1] "zh_CN.UTF-8/zh_CN.UTF-8/zh_CN.UTF-8/C/zh_CN.UTF-8/zh_CN.UTF-8"

#使绘图物件中的中文文字能正确呈现
#调用加载"showtext"套件包
library(showtext) 
showtext_auto(enable = TRUE)

2 读取SPSS数据档，并存成R数据档

#加载"sjlabelled" 套件包
library(sjlabelled)

#将SPSS数据汇入到R(数据档要放在目前的工作目录下)
mydata <- read_spss("stfood2019w.sav")

Converting atomic to factors. Please wait...

#看一下数据档中的前6笔数据
head(mydata)

  sex class dep Q14        w
1   2     1   2   4 1.229556
2   1     1   2   4 0.570904
3   2     2   2   4 1.239425
4   2     1   2   4 1.229556
5   1     1   2   4 0.570904
6   2     1   2   2 1.229556

#存成R数据檔
save(mydata, file = "mydata.rda")

#把R数据档叫进来
load("mydata.rda")

3 检视样本数据档

View(mydata)

样本数据档中的 w 变量是使用SPSS加权产生的 加权变量

4 呈现原始数据分析变量的频数分布表(未加权)

4.1 “sex”变量(性别)

可以使用 sjmisc 套件包的 frq( ) 函数，进行分析变量的频数分布。得知男性39.22%, 女性60.78%。

#呈现数据档(stfood2019w)中的变量名称
names(mydata)

[1] "sex"   "class" "dep"   "Q14"   "w"

#加载"sjmisc"套件包
library(sjmisc)

#sex变量的频数分布
frq(mydata$sex) #男39.22%, 女60.78%。

性别 (x) <categorical> 
# total N=1058 valid N=1058 mean=1.61 sd=0.49

Value | Label |   N | Raw % | Valid % | Cum. %
----------------------------------------------
    1 |    男 | 415 | 39.22 |   39.22 |  39.22
    2 |    女 | 643 | 60.78 |   60.78 | 100.00
 <NA> |  <NA> |   0 |  0.00 |    <NA> |   <NA>

4.2 “class”变量(年级)

从频数分布表得知，大一＋大二51.23%, 大三＋大四48.77%。

#加载"sjmisc"套件包
library(sjmisc)

#class变量的频数分布
frq(mydata$class) #大一＋大二51.23%, 大三＋大四48.77。

年级 (x) <categorical> 
# total N=1058 valid N=1058 mean=1.49 sd=0.50

Value |     Label |   N | Raw % | Valid % | Cum. %
--------------------------------------------------
    1 | 大一+大二 | 542 | 51.23 |   51.23 |  51.23
    2 | 大三+大四 | 516 | 48.77 |   48.77 | 100.00
 <NA> |      <NA> |   0 |  0.00 |    <NA> |   <NA>

4.3 “dep”变量(学系)

从频数分布表得知，新闻系36.77%,传播系39.88%, 设计系20.32%,智能媒体系2.93%。

#加载"sjmisc"套件包
library(sjmisc)

#dep变量的频数分布
frq(mydata$dep)

学系 (x) <categorical> 
# total N=1058 valid N=1058 mean=1.89 sd=0.82

Value |    Label |   N | Raw % | Valid % | Cum. %
-------------------------------------------------
    1 |     新闻 | 389 | 36.77 |   36.77 |  36.77
    2 |     传播 | 423 | 39.98 |   39.98 |  76.75
    3 |     设计 | 215 | 20.32 |   20.32 |  97.07
    4 | 智能媒体 |  31 |  2.93 |    2.93 | 100.00
 <NA> |     <NA> |   0 |  0.00 |    <NA> |   <NA>

#新闻系36.77%,传播系39.88%, 设计系20.32%,智能媒体系2.93%。

5 制造一个list档，内含分析变量应有的总体比例，并赋予标签

分析变量包括：

sex(性别) 1:“男性” 0.2423, 2:“女性” 0.7577
class(年级) 1:“大一＋大二” 0.4937, 2:“大三＋大四” 0.5063
dep(年级)变量 1:“新闻系” 0.3039, 2:“传媒系” 0.3915, 3:“设计系” 0.2793, 4:“智能媒体系” 0.0253

#新增观察值"ID"变量
library(sjmisc)
mydata <- add_id(mydata, var = "ID") 

#输入分析变量应有的总体比例
#sex(性别)变量 1:男性0.2423, 2:女性0.7577
sex <- c(.2423, .7577) 

#class(年级)变量 1:大一＋大二0.4937, 2:大三＋大四0.5062
class <- c(.4937, .5063)

#dep(年级)变量 1:新闻系0.3039, 2:传媒系0.3915, 3:设计系0.2793, 4:智能媒体系0.0252
dep <- c(.3039, .3915, .2793, .0253)

#将变量的总体比例存成一个list档，并加以命名(targets)
targets <- list(sex,class,dep)

#在targets档中，将各个变量的总体比例，个别给予变量标签。
#此一变量标签要和样本中分析变量的标签一致。
names(targets) <- c("sex", "class", "dep")

#查看数据档的属性以及内含的变量
attributes(targets)

$names
[1] "sex"   "class" "dep"

  #属性：list ; 名称："sex"  "class" "dep"  

#此一总体比例list档(档名：targets)在下面的加权分析中，将会使用到。

检视总体比例数据list档(无变量标签vs.有变量标签)

View(targets)

6 进行多变量反覆加权 anesrake::anesrake( )

所谓多变量反覆加权(raking)，就是将某一个样本分析变量比例(例如：样本性别比例)和总体比例（例如：总体性别比例）做一比较对照后，研究者根据总体比例，将样本中分析变量的比例调整和总体比例一致。

调整的方法是给予一个权数（总体比例除以样本比例）。调整完第一个变量之后(例如：性别变量），接续调整第二个变量比例(例如：年级变量）；调整完第二个变量之后，继续调整第三个变量比例(例如：学系变量）。

如此经过多个变量不断的反覆加权，直到研究者所设定的变量(例如：性别、年级、学系)，在经过加权后，都呈现出和总体比例一致的面貌。

此时，我们就认为该样本能代表总体，是一个具有样本代表性的样本。

欲检验的分析变量通常是可以获得总体统计数据的变量。

其步骤大略如下：

使用 sjlabelled 套件包的 to_numeric( )函数，将分析变量sex(性别), class(年级), dep(学系)设置为 numeric 数字变量
使用 anesrake 套件包的 anesrake( ) 函数，进行sex,class,dep 三个变量的多变量反覆加权，产生加权变量。
将前述步骤所产生的加权变量与个案ID，组合成一个数据框档案(data.frame）。
将前述步骤产生的加权数据档与原始数据档，进行纵向合并，成为一个包含有加权变量在内的新数据框档案。

具体的分析步骤如以下程式代码。

6.1 将欲检验的分析变量，转化为数字(numeric)变量

#加载sjlabelled 套件包
library(sjlabelled)

##sex(性别)
mydata$sex <- to_numeric(mydata$sex)
class(mydata$sex) #numeric

[1] "numeric"

##class(年级)
mydata$class <- to_numeric(mydata$class)
class(mydata$class) #numeric

[1] "numeric"

##dep(学系)
mydata$dep <- to_numeric(mydata$dep)
class(mydata$dep) #numeric

[1] "numeric"

6.2 使用anesrake套件包的anesrake( )函数，进行多变量反覆加权，产生加权变量档

#加载anesrake套件
library(anesrake)

#使用anesrake( )函数，将总体比例list档(例如：targets)，与样本数据档(例如：mydata)，根据样本数据档(mydata)的编号变量(例如：mydata$ID)，进行多变量反覆加权，并将结果另存为新的数据档（例如：outsave）

outsave <- anesrake(targets, #target:总体比例list档
                    mydata, #mydata:样本数据档
                    caseid = mydata$ID,
     #caseid:指明样本数据档(mydata)观察值编号变量名称(ID)
                    verbose = TRUE #verbose=TRUE:显示加权过程 
)

[1] "Raking...Iteration 1"
[1] "Current iteration changed total weights by 297.106358138516"
[1] "Raking...Iteration 2"
[1] "Current iteration changed total weights by 24.4630570814103"
[1] "Raking...Iteration 3"
[1] "Current iteration changed total weights by 0.374783430163922"
[1] "Raking...Iteration 4"
[1] "Current iteration changed total weights by 0.00574630827282396"
[1] "Raking...Iteration 5"
[1] "Current iteration changed total weights by 8.81054261507863e-05"
[1] "Raking...Iteration 6"
[1] "Current iteration changed total weights by 1.35087852592974e-06"
[1] "Raking...Iteration 7"
[1] "Current iteration changed total weights by 2.07161960918967e-08"
[1] "Raking...Iteration 8"
[1] "Current iteration changed total weights by 3.24115734251507e-10"
[1] "Raking...Iteration 9"
[1] "Current iteration changed total weights by 5.5543347698972e-12"
[1] "Raking...Iteration 10"
[1] "Current iteration changed total weights by 2.17215134767912e-12"
[1] "Raking...Iteration 11"
[1] "Current iteration changed total weights by 2.93687296704093e-12"
[1] "Raking converged in 11 iterations"

上述过程显示。总共反覆加权了11次，才达到收敛。

6.3 将 “样本编号” (例如：ID) 和 “加权权数”(例如：Raking_W) 组合成一个 “数据框”（data.frame），并给予名称（例如：caseweights）

新的加权结果数据档（outsave）包含有两个重要变量（caseid 与 weightvec）。

我们将这两个变量重新命名，将 caseid* 变量重新命名为 ID，将 weightvec 变量重新命名为 Raking_W。

这个 Raking_W变量就是我们经过多变量反覆加权后所得到的权数。

然后将这两个变量组合成一个数据框(data.frame)，并给予名称（例如：caseweights）。

#将变量重新命名并组合成一个数据框(data.frame)
caseweights <- data.frame(ID = outsave$caseid, 
                          Raking_W = outsave$weightvec)

#检视加权结果数据框
View(caseweights)

6.4 将原始样本数据档(例如:mydata)，与加权变量数据档(例如：caseweights) 进行纵向合并

接下来，我们可以使用 sjmisc 套件包的 replace_columns( )函数，将 原始样本数据档(例如mydata)，与 加权变量数据档(例如：caseweights)，根据共同关键变量(例如：ID)进行纵向合并，合并后的物件另外命名(例如：mydata_merge）

#加载sjmisc套件包
library(sjmisc)

#将两个数据框(mydata, caseweights)进行纵向合并
mydata_merge <- replace_columns(mydata,caseweights) 
                               #mydata:原始样本数据档
                               #caseweights:加权变量数据档

#检视合并后的数据档
View(mydata_merge)

可以发现，透过SPSS产生的加权变量 w，与透过R的anesrake套件包所产生的加权变量 Raking_W，两者的数值几乎是一样的。

6.5 将分析变量转换为因子变量(factor)

library(sjlabelled)
mydata_merge$sex <- to_factor(mydata_merge$sex)
mydata_merge$class <- to_factor(mydata_merge$class)
mydata_merge$dep <- to_factor(mydata_merge$dep)

6.6 呈现加权后分析变量的频数分布表

#加载sjmisc套件包
library(sjmisc)

#呈现加权后分析变量的频数分布表
frq(mydata_merge$sex, weights = mydata_merge$Raking_W)

性别 (xw) <categorical> 
# total N=1058 valid N=1058 mean=1.76 sd=0.43

Value | Label |   N | Raw % | Valid % | Cum. %
----------------------------------------------
    1 |    男 | 256 | 24.20 |   24.20 |  24.20
    2 |    女 | 802 | 75.80 |   75.80 | 100.00
 <NA> |  <NA> |   0 |  0.00 |    <NA> |   <NA>

frq(mydata_merge$class, weights = mydata_merge$Raking_W)

年级 (xw) <categorical> 
# total N=1058 valid N=1058 mean=1.50 sd=0.50

Value |     Label |   N | Raw % | Valid % | Cum. %
--------------------------------------------------
    1 | 大一+大二 | 524 | 49.53 |   49.53 |  49.53
    2 | 大三+大四 | 534 | 50.47 |   50.47 | 100.00
 <NA> |      <NA> |   0 |  0.00 |    <NA> |   <NA>

frq(mydata_merge$dep, weights = mydata_merge$Raking_W)

学系 (xw) <categorical> 
# total N=1058 valid N=1058 mean=2.03 sd=0.83

Value |    Label |   N | Raw % | Valid % | Cum. %
-------------------------------------------------
    1 |     新闻 | 322 | 30.43 |   30.43 |  30.43
    2 |     传播 | 414 | 39.13 |   39.13 |  69.57
    3 |     设计 | 295 | 27.88 |   27.88 |  97.45
    4 | 智能媒体 |  27 |  2.55 |    2.55 | 100.00
 <NA> |     <NA> |   0 |  0.00 |    <NA> |   <NA>

7 样本代表性检验

7.1 原始样本(未加权)样本代表性检验

可以使用 sjstats 套件包的 chi_squared_test( ) 函数，加上 probabilities( )参数，针对原始样本(未加权)的 sex(性别)、class(年级) 与dep(学系)变量，进行样本代表性检验。

检验结果发现，原始样本的sex(性别)与dep(学系)变量的分布与总体不一致，无法代表总体。而原始样本的class(年级)变量，虽通过统计检验，但比例方向不一致，所以总的来说，样本的class(年级)结构亦无法代表总体。

具体代码如下。

#original data原始样本代表检验
#加载sjstats套件包
library(sjstats)

#检验sex(性别)变量
chi_squared_test(mydata_merge, "sex", 
                 probabilities = c(.2423, .7577)
          #性别应有的总体比例       #男    #女  
#χ² = 129.576, df = 1, p < 0.001 #不具代表性
)

# Chi-squared test for given probabilities

  Data: sex against probabilities 24% and 76% (n = 1058)

  χ² = 129.576, פ‎ = 0.198 (small effect), df = 1, p < .001

#检验class(年级)变量
chi_squared_test(mydata_merge, "class", 
                 probabilities = c(.4937, .5063)
        #年级应有的总体比例 #大一＋大二 #大三＋大四 
#χ² = 1.462,  df = 1, p = .227 >0.05 
#样本的class变量，虽通过统计检验，但比例方向不一致，所以总的来说，样本的class结构无法代表总体。
)

# Chi-squared test for given probabilities

  Data: class against probabilities 49% and 51% (n = 1058)

  χ² = 1.462, פ‎ = 0.037 (tiny effect), df = 1, p = 0.227

#检验dep(学系)变量
chi_squared_test(mydata_merge, "dep", 
         probabilities = c(.3039, .3915, .2793, .0253) 
       #学系应有的总体比例 #新闻  #传媒  #设计  #智能媒体
#χ² = 36.945, df = 3, p < 0.001 #不具代表性
)

# Chi-squared test for given probabilities

  Data: dep against probabilities 30%, 39%, 28% and 3% (n = 1058)

  χ² = 36.945, פ‎ = 0.030 (tiny effect), df = 3, p < .001

7.2 加权后样本代表性检验

同样的，我们可以使用 sjstats套件包的 chi_squared_test( ) 函数，与 probabilities( ) 参数，再加上 weights( ) 这一个参数，就可得到加权后样本代表性检验结果。

检验结果发现，加权后样本的sex(性别)、class(年级)与dep(学系)变量的分布皆通过统计检验，与总体一致。所以总的来说，加权后的样本结构，是一个具有代表性的样本，可以代表总体。