给学生做R培训，我总会先讲数据码农编码三板斧：

• 代码文件必须有文件说明；
• 像写论文那样写注释；
• 设好工作目录；

说到第三点总是有点感觉对不起听众，一丝丝把人带到坑里去了的隐忧。因为在 Windows 系统下，R的文件路径总是这种画风：

setwd('c:\\users\\xxx\\Documents\\R\\projects\\yyy')
# 或者
setwd('c:/users/xxx/Documents/R/projects/yyy')

对于 Windows 的反斜杠(\)路径符号来说简直是逆行开车的节奏。Windows下码农兵器库中有一大堆武器可以很方便的提取路径名、文件名，最不济也可以从资源管理器直接复制路径名过来，但对于R的这种画风最后总不得不在 RStudio 中把 \ 换成 / 或 \\。查找替换用的再熟练，总归是项体力劳动。

那么我为什么还要坚持这项无聊的枯燥活动呢？因为我不需要啊！看动画！

再如图设个快捷键，Ctrl+Shift+V 搞定！

那么这么好用的码农居家神器哪里找呢？出门开车上 github 咯！

install_github("huashan/RSExpert")

这两年来 markdown 越发流行于各类网络客户端系统，各大博客系统、wiki系统、团队协作系统都纷纷支持 markdown 语法，支持 markdown 的编辑器也越来越多。但真正将 markdown 用于学术论文或出版物写作的实践却并不多。在此介绍一下本人自2012年来就开始打造的一套工作流程。

先介绍流程的两大基石：

knitr 背后的功臣——pandoc

markdown 起源于email写作中的简单格式，其语法非常简单，可能对 markdown 略有耳闻的人会觉得它不过是个玩具。其实不然，markdown 本身并非一个统一的语法，事实上有很多版本的实现。比较有经验的R用户可能经常用 knitr 来做一些数据报告，甚至用 knitr book 来写数据类长篇大作，或许你们还不知道，knitr 背后有个默默无名的英雄—— pandoc。这是一款由一名伯克利大学哲学系教授领衔开发的 markdown 超级语法及文档格式转换工具。与普通网页版的简易 markdown 语法不同，pandoc 可以说是学术 markdown，除了常规的 markdown 元素以外，它还支持学术论文中常见的公式、表格、参考文献、脚注、尾注的格式化输出。 不仅如此，它还是一款超级格式机，能与大量的文档格式进行转换，它支持的格式是如此之多我就不贴大图了，仅手绘一张 pandoc 流程图做一示意：

所以，有了 pandoc 这个超级工具，社科量化小白用户再也无需羡慕理工科学生手里能输出漂亮公式的 latex 了，用 markdown 同样能输出漂亮公式，beamer 幻灯片也不在话下，再也不用啃三个月 latex 手册了。而且，印刷论文的时候可以用漂亮的pdf输出，要求论文电子版存档就交个 word 交差，何其方便！

思维脑图写作——mindmanager

前面说了简单数据报告可以采用 knitr ，但完整的社科类论文写作还是很难完全依赖它的。除了社科类论文中，数据分析在文章中占比不大以外，更重要的是思维差异。写作思维与编程思维是完全不同的！数据分析过程，强调数据结构和处理逻辑的一致性，强调代码的重用性，尽量避免重复代码，而对计算结果进行分析和写作则是基于分类和演绎的逻辑来进行，很多时候，我们需要用多种口径检查计算结果，但并不会将其都放在论文中。个人经验，最有效的方式还是把这二者分开，各行其是才好。

脑图类软件无疑是整理构思、研究素材、论文结构的极好工具。脑图类软件也有很多，收费的、开源的、在线的，那么是不是随便搞个看着顺眼的就可以了？当然不是了。整理构思、管理素材固然是用脑图软件的一大目的，但是还得将它与 makrdown 文字相结合起来。该软件比如符合两个要求：

• 能够添加笔记：这条就杀掉了大多数在线脑图；
• 提供二次开发 API：添加笔记的目的是用来书写内容，而 API 则是为这些笔记的输出提供各种可能性。

所以我选择了 mindmanager：丰富的 OLE 接口，可以用来任意自定义脑图结构和笔记的输出方式。例如如何划分章、节和小节，第四级以下的脑图节点是否还作为小小节；不同输出版本的生成，例如下图中，标记有 advanced 的单独输出为一个版本。

有 API 的好处是，还可以对笔记内容在变成 markdown 正文前进行处理。例如可以将笔记中的图片直接转为 [img](path_to_image) 这样的语法，充分利用笔记本身所见即所得的功能。例如分枝标题上的 - 表示不显示该小节标题：

更进一步，在笔记内部进行再分页、分章节也完全不是问题。同一个脑图，即可以输出为文章，也可以输出为幻灯片，一切取决于你的相像！最后，mindmanager 的 OLE 接口语法是如此简单，只需要简单的 visual basic 知识就足以搞定。

流程整合：Paper Composer

有了上面两大基石的加持，只需小小地在 mindmanager OLE 接口上下点功夫，论文编译器就出台了！

主界面来了，嗒哒！

功能也不多，随想随加：

• 输出的 markdown 文件存放路径；
• 导出脑图范围：全部 vs. 分枝；
• 章节层级；
• 是否作为 Rmd 输出，以供 knitr 运算；
• 参考文献格式化；
• 根据 tag 标识选择不同分枝输出；
• 最终格式选择：Word, HTML, pdf；

还剩了一些问题，例如：

• 参考文献格式化；
• R 的输出如何整合；
• 导出到 Word 时的问题：模板设置、表格图片自动编号；表格边框设置；

请持续关注。

目的

• 在导入 SPSS Stata 等格式时，提供统一的处理标签的接口；
• 在输出表格时提供对标签的自动化处理；

什么是标签

SPSS 和 Stata 的用户最熟悉标签：

• 变量标签

gender `性别`
age    `年龄`

• 数值标签

gender:

    1 = 男性
    2 = 女性

R 中如何处理标签

首先，data.frame 没有变量标签的概念（attr 另外再说），变量名行使标签的功能：

table(dt$`性别`)
# 如果标签有空格
table(dt$`a variable`)

其次，可用 factor 提供数值标签功能（factor 其实就等价于字符型）：

levels(df$`性别`) <- c('男', '女')

R 包对标签的处理

通过对象属性设置来保存标签：

• foreign

attr(df, 'variable.labels') # 变量标签
attr(var, 'value.labels')   # 数值标签

• haven

和 foreign 道理一样，只不过命名方式不同。此外，haven 把变量标签和数值标签都作为变量的属性，而在 foreign 中，变量标签是 data.frame 的属性，数值标签才是变量的属性。

attr(var, 'label')
attr(var, 'labels')

显示在 RStudio 中是这样的：

ezdf 提供统一的数据导入接口

• 对数据导入提供统一接口，封装 foreign、haven 等包的导入函数；
• 修正一些 bug，例如：

`haven`: “Error: `x` and `labels` must be same type”

导入 Stata 数据

导入 Stata 数据使用 readStata() 函数：

library(ezdf)
dat <- readStata('CGSS2013（居民问卷）发布版_2014.dta', encoding = 'GB2312')
# View(dat)

# 参数 `encoding` 设置 Stata 标签的编码，该参数默认值为 UTF-8。
# 有的 Stata 数据对变量名以及字符变量（string）的值都采用不同编码，对于这种情况，
# 需分别设置 `varNameEncoding` 和 `charEncoding`。
dat <- readStata('CGSS2013（居民问卷）发布版_2014.dta', encoding = 'GB2312', 
                 varNameEncoding = 'UTF-8', charEncodin = 'UTF-8')

导入 SPSS 数据

导入 SPSS 数据使用 readSPSS() 函数

# 参数 `lib` 设置导入所使用的 R 包，目前支持 `foreign` 和 `haven`。
readSPSS(file, lib = "foreign", ...)

将 data.frame 转换为 ezdf

用 as.ez(dt, meta) 创建一个新的 ez.data.frame 对象

data(iris)
library(ezdf)
d1 = as.ez(iris)
class(d1)
## [1] "ez.data.frame" "data.table"    "data.frame"

ezdf 对标签的设置

变量标签

• 变量标签存储在 meta 属性当中；
• meta 可为 data.frame 或 matrix 类型对象：至少包括两列：第一列为变量名，第二列为变量标签。

两个辅助函数：

- `setmeta()`
- `getmeta()`

d1$test = sample(5, size = nrow(iris), replace = T)

# 对新变量 test 设置变量标签
setmeta(d1, data.frame(var= 'test', lbl = '这是新变量标签'))

# 显示数据 d1 的全部变量标签
attr(d1, 'meta')
##     var            lbl
## 1: test 这是新变量标签

# 或者
getmeta(d1)
##     var            lbl
## 1: test 这是新变量标签

varLabels(d1, c('Species', 'test'))
## [1] ""               "这是新变量标签"

# 用 default = "var" 只输出带有标签的变量
varLabels(d1, c('Species', 'test'), default = "var")
## [1] "Species"        "这是新变量标签"

# 返回所有已定义的变量标签
varLabels(d1)
##     var            lbl
## 1: test 这是新变量标签

单独设置部分变量标签

# 设置变量标签
varLabels(d1, "test") <- "新标签"
varLabels(d1, "test")
## [1] "新标签"

数值标签

数值标签的存储采用命名整数向量作为变量的 labels 属性

# 定义一个数值标签
c(C1 = 1, C2 = 2, C3 = 3, MI = 9)

• valueLabels()

vl1 = valueLabels(d1, 'test')
vl1
## list()
## attr(,"class")
## [1] "value.labels"
## attr(,"ez")
## [1] "d1"
## attr(,"col")
## [1] "test"

# 数值标签可以“加减”。
# 注意： MI=9 设了一个不存在的值标签
vl2 = vl1 + c("Class1"=1, "Class2"=2, "Class3"=3, 'Class4'=4, 'Class5' = 5, MI = 9, MM = 8)
valueLabels(d1, 'test') = vl2

制表函数

tbl()

tbl(d1, ~test)
##    test\t新标签  N
## 1:    1++Class1 25
## 2:    2++Class2 37
## 3:    3++Class3 23
## 4:    4++Class4 25
## 5:    5++Class5 40

# 分组求均值，添加样本数
tbl(d1, Sepal.Length ~ Species + test, 'mean', N = T)
##        Species test\t新标签 Sepal.Length  N
##  1:     setosa    1++Class1     5.141667 12
##  2:     setosa    2++Class2     4.923077 13
##  3:     setosa    3++Class3     5.200000  5
##  4:     setosa    4++Class4     5.044444  9
##  5:     setosa    5++Class5     4.836364 11
##  6: versicolor    1++Class1     5.900000  4
##  7: versicolor    2++Class2     5.941667 12
##  8: versicolor    3++Class3     5.908333 12
##  9: versicolor    4++Class4     6.142857  7
## 10: versicolor    5++Class5     5.866667 15
## 11:  virginica    1++Class1     6.533333  9
## 12:  virginica    2++Class2     6.583333 12
## 13:  virginica    3++Class3     6.883333  6
## 14:  virginica    4++Class4     6.344444  9
## 15:  virginica    5++Class5     6.657143 14

# tbl() 默认按照公式右端 x 的值排序，如果取消排序
tbl(d1, Sepal.Length ~ Species + test, 'mean', N = T, sort = F)
##        Species test\t新标签 Sepal.Length  N
##  1:     setosa    4++Class4     5.044444  9
##  2:     setosa    5++Class5     4.836364 11
##  3:     setosa    1++Class1     5.141667 12
##  4:     setosa    2++Class2     4.923077 13
##  5:     setosa    3++Class3     5.200000  5
##  6: versicolor    4++Class4     6.142857  7
##  7: versicolor    2++Class2     5.941667 12
##  8: versicolor    5++Class5     5.866667 15
##  9: versicolor    3++Class3     5.908333 12
## 10: versicolor    1++Class1     5.900000  4
## 11:  virginica    4++Class4     6.344444  9
## 12:  virginica    5++Class5     6.657143 14
## 13:  virginica    2++Class2     6.583333 12
## 14:  virginica    1++Class1     6.533333  9
## 15:  virginica    3++Class3     6.883333  6

ctbl()

ctbl() 是对 table() 的封装，采用 ctbl(ez, expr) 的调用方式。

ctbl(d1, Sepal.Length ~ Species + test)
# 等价于
table(d1$Sepal.Length, d1$Species, d1$test)

ftable()

ftable.ez.data.frame() 方法是对 ftable() 的封装

ftable(ez, formula, style = 1, prop_margin = 1, ...)

• prop_margin： 行百分比 / 列百分比；
• style = 1：输出频次；
• style = 2：输出百分比；
• style = 3：输出百分比和行加总频次。

ftable(d1, Species~test)
##            
##             setosa versicolor virginica
##   1++Class1     12          4         9
##   2++Class2     13         12        12
##   3++Class3      5         12         6
##   4++Class4      9          7         9
##   5++Class5     11         15        14

ftable(d1, Species~test, style = 2)
##            
##                setosa versicolor virginica
##   1++Class1 0.4800000  0.1600000 0.3600000
##   2++Class2 0.3513514  0.3243243 0.3243243
##   3++Class3 0.2173913  0.5217391 0.2608696
##   4++Class4 0.3600000  0.2800000 0.3600000
##   5++Class5 0.2750000  0.3750000 0.3500000

(t1 = ftable(d1, Species~test, style = 3))
##              setosa versicolor virginica  N
## 1++Class1 0.4800000  0.1600000 0.3600000 25
## 2++Class2 0.3513514  0.3243243 0.3243243 37
## 3++Class3 0.2173913  0.5217391 0.2608696 23
## 4++Class4 0.3600000  0.2800000 0.3600000 25
## 5++Class5 0.2750000  0.3750000 0.3500000 40

与 markdown 流程整合

• pander 是用于 markdown 格式输出的 R 包，提供了非常丰富的表格输出功能
• 在加载 ezdf 包之后，会自动与 pander 包结合，实现自动标签输出

# pander 输出
library(pander)

pander(t1, ez = d1)

这是输出的 markdown 结果：

-----------------------------------------------------
            setosa   versicolor   virginica   N 
--------------- -------- ------------ ----------- ---
 **1++Class1**    0.48       0.16        0.36     25 

 **2++Class2**   0.3514     0.3243      0.3243    37 

 **3++Class3**   0.2174     0.5217      0.2609    23 

 **4++Class4**    0.36       0.28        0.36     25 

 **5++Class5**   0.275      0.375        0.35     40 
-----------------------------------------------------

最终输出效果：

pander 与回归结果输出：

# 加上数值标签
options('ezdfKeepVal' = T)
pander(tbl(dat, a66 ~ s5a, 'mean'))

---------------------------------------------------------
 s5a	受访者居住的地区类型是   a66	您家是否拥有家用小汽车 
---------------------------- ----------------------------
    1++市/县城的中心地区                1.761            

    2++市/县城的边缘地区                1.775            

   3++市/县城的城乡结合部               1.809            

    4++市/县城区以外的镇                1.859            

          5++农村                       1.919            
---------------------------------------------------------

数值与标签之间分隔符

options('ezdfValueLabelSep' = '=')
pander(tbl(dat, a66 ~ s5a, 'mean'))

---------------------------------------------------------
 s5a	受访者居住的地区类型是   a66	您家是否拥有家用小汽车 
---------------------------- ----------------------------
    1=市/县城的中心地区                 1.761            

    2=市/县城的边缘地区                 1.775            

   3=市/县城的城乡结合部                1.809            

    4=市/县城区以外的镇                 1.859            

           5=农村                       1.919            
---------------------------------------------------------

回归模型的输出：

m1 = lm(a6 ~ a2 + a10, dat)
pander(m1)

这是通过 markdown 输出转成 pdf 后的效果，没有任何手工干预（pandoc 在输出小数点时，还有一点瑕疵，比如小数点位数不统一，不过本人 github 上的版本已经修正了这个问题）。

表格输出选项

目前提供三个选项：

• options('ezdfKeepVal' = T)
• options('ezdfValueLabelSep' = '=')
• options('ezdfKeepVarName' = T)

options('ezdfKeepVal' = T)
options('ezdfValueLabelSep' = '=')
options('ezdfKeepVarName' =  F)

tbl(d1, ~test)
##      新标签  N
## 1: 1=Class1 25
## 2: 2=Class2 37
## 3: 3=Class3 23
## 4: 4=Class4 25
## 5: 5=Class5 40
pander(tbl(d1, ~test))

options('ezdfKeepVarName' =  T)
options('ezdfValueLabelSep' = '++')
pander(tbl(d1, ~test))

下载与安装

github: https://github.com/huashan/ezdf

devtools::install_github('huashan/pander')

前面讲过，R 中字符向量可以有多种编码。一般情况下，对于混合编码，R 都能很好的自动处理。例如：

x1 <- '中国'
x2 <- iconv(x1, 'GB2312', 'UTF-8')
## == 判断也没问题
x2 == '中国'
## [1] TRUE

x <- c(x1, x2)
grep('中', x)
## [1] 1 2
grep(iconv('中', 'GB2312', 'UTF-8'), x)
## [1] 1 2

在写命令时，无需特意为了目标字符将匹配条件转为相应的编码，R 自动就在后台做了转换。然而，总有些让人意想不到的的角落，例如：

a1 <- iconv('中文 aa', 'gb2312', 'UTF-8')
Encoding(a1)
## [1] "UTF-8"

# 当待替换字符非 UTF-8 编码，替换字符为 UTF-8 编码，结果乱码
(r1 <- sub('aa', a1, 'aa ba'))
## [1] "涓枃 aa ba"
Encoding(r1)
## [1] "unknown"
# 需要手工纠正
Encoding(r1) <- 'UTF-8'
r1
## [1] "中文 aa ba"

(r2 <- sub('aa', a1, 'aa ba', useBytes = F, fixed = T))
## [1] "涓枃 aa ba"
Encoding(r2)
## [1] "unknown"
Encoding(r2) <- 'UTF-8'
r2
## [1] "中文 aa ba"

在线帮助里是这么写的：

For sub and gsub...,..., If useBytes = FALSE a non-ASCII substituted result will often be in UTF-8 with a marked encoding (e.g., if there is a UTF-8 input, and in a multibyte locale unless fixed = TRUE). Such strings can be re-encoded by enc2native.

嗯，这是少数几处帮助文本与实际情况脱离的个例。

待替换字符为 UTF-8 编码，倒是没有问题：

sub('aa', 'bb', a1, fixed = F)
## [1] "中文 bb"

如果是我们自己写的代码，可以小心谨慎，避免跳坑，很容易就能找到替代方案。博主随便列举一下，就有两种：

方式一：坚持极简主义不用第三方包，使用正则表达式之动力引擎 regexpr() 进行替换

s <- 'aa ba'
m <- regexpr('aa', s)
regmatches(s, m) <- a1
s
## [1] "中文 aa ba"

方式二：用包也不羞耻之 stringr 包

library(stringr)
str_replace('aa', 'aa', a1)
## [1] "中文 aa"
str_replace(a1, 'aa', '+')
## [1] "中文 +"

但是 sub() 和 gsub() 作为R的基础函数，在大量的包中广泛存在，怎么避免让这些包跳坑呢？

本来还想介绍一个与混合编码有关的奇葩案例（来自性能无敌的 data.table 包），但是实在太奇葩了，所以今天就提前收尾吧。

替第四篇做个预告。用过 tm 包的童鞋想必经常碰到为什么 tm 包里词都连起来了？ 这样的问题。问题就出在 scan() 这个基础函数（怎么又是基础函数！）。scan() 是R中用来拆分文本的函数，读取文本数据（read.table()）都要用。我记得以前是正常的，但不记得是哪个R版本升级之后，就变成下面这样了，“的”和“分词”被当作了同一个词。scan()的这个奇葩结果就是中文文本分了词让 tm 建矩阵结果还是分的乱七八糟的原因。

x = '测试 中文 的 分词'
scan(text = x, what='', sep=' ')
## [1] "测试"    "中文"    "的 分词"

第四篇将介绍如何修理它。

罚似然图模型扩展

基于罚似然回归方法的社会网络关系测度不仅适用于小群体的网络数据，更适用于大规模的社会网络数据。罚似然回归其本质上是回归估计和模型变量选择，统计学家们通过模拟分析，已经证明其具有非常好的稳健性，对于几千甚至上万的自变量选择具有一致性(Tibshirani, 1996)。另外，使用罚似然图模型进行社会关系网络测度，可以根据无向无权的二分双模网络数据估计得到无向有权的关系网络矩阵，从而不仅可以对关系进行有无的判定，还可以进行强度的比较，大大丰富了分析内容。除此之外，罚似然图模型还具有很强的扩展性，本节将对此展开介绍。

针对不同数据类型的图模型

经典的高斯图模型对变量的假设为多元正态分布。但在社会科学研究中，往往会碰到多种类型的数据，甚至是混合类型的数据，包括二分数据、定类数据、定序数据、计数数据、有偏分布的连续数据等等。例如，前述的美国南方妇女数据即为二分变量，网民在论坛的发帖回帖为计数型变量，在某个场所停留的时间可算是计数型变量或有偏的连续变量。再举个例子，对于健康领域的社会学大数据研究则可能要考虑性别（二分）、年龄（连续）、行为模式（计数）、事件发生的场所（类别）、用药的计量（连续）等各类数据类型之间的关联模式。基于此，统计学家们发展了多种特殊模型予以解决。例如，针对二项分布数据的估计问题可进一步参考(Banerjee & El Ghaoui et al., 2008; Ravikumar & Raskutti et al., 2008; van Borkulo & Borsboom et al., 2014)，针对泊松分布数据可参考(Allen & Liu, 2012; Allen & Liu, 2013)，针对多分类分布可参考(Dai & Ding et al., 2013)，针对混合数据类型的估计问题可参考。略有遗憾的是，目前为止，尚未有一个软件包将所有数据类型综合到一个框架下进行处理。

带协变量的罚似然图模型

在罚似然模型中，除了对所有变量加罚，还可以仅对部分变量加罚。将公式三的罚则项

$$\lambda\sum_{p = 1}^{p}\left| b_{p} \right|$$

改写为：

$$\lambda\sum_{p = 1}^{m}\left| b_{p} \right|$$

其中$m < p$，即表示仅对部分自变量加罚。因此，很容易引入其它协变量进入模型。以DGG数据为例，由于所记录的事件来自多种聚会类型，尽管人类学家们没有记录事件的具体类型从而缺失了相关信息，但是可以假设，不同的活动类型和参与规模相关，并进而影响不同人的参与程度。因此，在本示例中，将参与活动人数作为协变量引入图模型，得到的结果如图2所示。相比起图一，在控制了参与人数规模这个因素之后，图二仍然保留了基本相同的网络结构，编号17和18通过编号16与其它成员相关联。但与图一不同的是，编号1、编号8和编号16处在了网络桥的位置，而编号9不再成为网络桥，而是成了第一网络子群的成员。在弗里曼的元分析中，编号8的分组其实存有争议，21个方法中有4个将其判定为第二分组，另有7个方法无法处理编号8只能将其剔除。从原始数据来看，编号9所参与的4次活动均是这个群体参与人数最多的活动。因此，针对编号9的网络地位，除了技术问题以外（由于观测数太少而导致估计不稳定），可能的推论有两个：从众的边缘成员或者比较重大事件才出席的核心人物。选择何种推论就取决于对活动信息的了解。遗憾的是，原始数据缺乏相关信息。

图2：控制聚会规模以后的美国南方妇女网络

多组罚似然图模型

若协变量为类别变量，也就是说，观测样本可能来自不同的子总体，那么有两种策略：1、用不加罚的方式将协变量引入模型，此时估计得到的是一个总体网络，消除了不同类别之间的异质性。2、对子总体分别进行建模，从而得到多个网络。该方式的缺点在于无法进一步分析网络之间的共性。

除了不同子总体的样本之外，在时点观测数据中往往需要假设一定的异质性：在一个随时点变化的观测中，存在一个公共的网络结构，在不同时间段网络结构发生缓慢变化或者突变。例如，对于学术引文网络来说，在上个世纪六七十年代，由于布劳邓肯地位获得模型的成功，社会流动研究领域的引文可能会更多涉及路径模型和结构方程模型方面的文献，而在八九十年代之后引文中可能更多出现对数线性模型方面的文献。在这个统计技术变迁的过程当中，核心的关注主题并没有发生变化，不同时段的引文仍然具有一定的共性。

对于这种异质性数据，有不同的分析策略：一为在考虑异质性的条件下，估计一个平滑的共同网络结构(Zhou & Lafferty et al., 2010; Kolar & Xing, 2011)。二为假设不同子总体之间存在一个公共网络结构，但每个子总体由于其自身的结构特殊性而具有独特的网络结构，需同时估计多个子网络结构。对于后一种策略来说，针对观测的独立同分布（i.i.d）假设，将glasso模型进一步扩展为联合glasso模型，可以在一个统一的分析框架下，考察在同一个群体中多个性质的网络关系如何叠加和扩展。对于该问题，需要使用两个惩罚因子，一个惩罚因子用来控制所有子样本中的公共因子$\theta_{j,j^{'}}$的稀疏度，另一个惩罚因子用来控制子样本内部的稀疏度。朱等人(Zhu & Shen et al., 2014)所提的方案是为每一个子总体估计一个稀疏图结构，也同时估计跨子图的网络凝聚点。郭等人(Guo & Levina et al., 2011)使用分层罚模型估计来保留公共的图结构，同时允许组间差异。当$p\log\left( p \right)/n$趋向于0时，其中$p$为变量个数，$n$为样本规模，该方法可实现弗罗贝尼乌斯范数（Frobenius norm）收敛。但这也意味着当$p>n$时，并不能获得稳定的估计。

达纳赫等人(Danaher & Wang et al., 2014)提出两种联合估计算法，FGL（fused Lasso）和GGL（grouped Lasso）来使罚似然对数最大化，但并未给出其估计量统计收敛的理论验证。蔡天文等人(Cai & Li et al., 2015) 提出一个改进型模型（MPE）联合估计$K$个稀疏精度矩阵，并对其统计收敛属性给予了理论验证。

上述几种方法的视角是将网络边作为分析的核心，即假设网络中某些边是公共的或是特殊的。莫汉等人则提出一个以网络节点为核心的视角(Mohan & London et al., 2014)，即某些节点的连结性在子图中具有共性，而另一些节点的连结性在不同子图中则具有特殊性。

本节的示例采用学术文献关键词关联网络，数据来自《社会学研究》和《社会》这两个杂志从2006年至2015年所发表的所有文章的关键词。在社会学研究文章中，关键词是体现一篇文章现实关注点、理论切入口和方法流派的重要指标。通过构建关键词关联网络，除了可以发现社会学研究中比较重要的核心关注要素之外，还可以比较这两个杂志在文章题材选择方面的偏好差异。在这个数据，共有1348篇文章，492个关键词。由于本示例主要为了展示模型特点，在此仅提供初步的模拟拟合结果。使用GGL法拟合的结果表明《社会学研究》的关键词关联有79对，《社会》83对，其中两个杂志共同的关键词关联41对。如图3-5所示。

## 多分组罚似然图模型
## 本示例数据采用《社会学研究》和《社会》杂志2006年至2015年全部发表论文的关键词。

kw <- readRDS('soc_keywords.RDS')
lapply(kw, dim)

library(JGL)
jgl.grp<- JGL(Y = kw, penalty='group', screening= 'fast', penalize.diagonal = F,
              lambda1=.001, lambda2=.001)
print.jgl(jgl.grp)

theta <- jgl.grp$theta

library(igraph)

K=length(theta)

adj = JGL:::make.adj.matrix(theta)
diag(adj)=0
gadj = graph.adjacency(adj, mode="upper", weighted=TRUE)
lay <- NULL
if (is.null(lay)) lay <- layout.fruchterman.reingold(gadj)

# class 1 class2 shared
# 38     42    41
table(E(gadj)$weight)

# 公共子图
gcomm <- igraph::delete.edges(gadj, which(E(gadj)$weight != 3))
plot(gcomm, layout=lay, vertex.size=0, vertex.color = 'white', 
     edge.color = 'black', main = "Common Network")

# 两个子图
nets <- JGL:::make.adj.matrix(theta, separate = TRUE)
for (k in 1:K) {
  g1 <- graph.adjacency(nets[[k]], mode = "undirected", diag=F)
  #png(filename = sprintf('x_keywords_%d.png', k), width = 800, height = 800)
  plot.igraph(g1, layout = lay, vertex.color = 'white', vertex.size = 0, 
              edge.color = 'black', main = paste("Sub Network", k))
  #dev.off()
}

图3：学术论文关键词关联网络（公共部分）

图4：学术论文关键词关联网络（社会学研究）

图5：学术论文关键词关联网络（社会）

更多的扩展

请见原文

按：原文见 《社会》官网，2017年第2期。代码及数据在此处。

题图致敬 DGG、Linton Freeman、R. Breiger。

导言

随着互联网和智能设备愈来愈多的介入人们的日常生活以及大数据概念的提出，在社会科学研究领域，研究者们面临一个新的非常巨大的数据源。不同于传统的问卷调查数据，这个新的数据源来自于各类智能设备所记录的数据，例如手机信号塔所记录的在某个范围内的人群聚集、摄像头所捕捉到的人们在各个场所的出席、人们在互联网使用过程中所留下的轨迹或积累的信息，例如在微博上面的评论或转发，在网络论坛的发帖和回帖。也有一些数据早已进入研究者的视野，但是由于信息化手段的丰富，研究者们无需再大费周章专门进行数据录入或转换，例如人们的日常消费记录、学术文献的作者信息和引文信息、公司间的联动交易行为等等。对于上述这类数据，研究者们往往关注其中的共现关系，并探讨其潜在的社会机制。例如共同在一个论坛帖子里进行讨论的用户可能具有相同的兴趣话题，在科学文献作品中科学家之间如何形成合作关系，图书购买记录背后所蕴含的共同的政治态度和价值观等等。

一般来说，对上述数据的分析大多采用社会网络分析方法进行。从数据分析的角度，这类互动数据可以采用一个发生矩阵（incidence matrix）来表示，例如一个$n \times m$的二进制矩阵$P$，其中矩阵的行表示某个场所或事件，例如微博的博文，学术文章或者购物清单等等，矩阵的列则表示参与该事件的基本单位或成员，例如转发微博的用户、文章作者或者购物商品名称。若 $p_{\text{ij}} = 1$ 则表示第 $j$ 个成员参与了第 $i$ 次事件，反之则表示没有参与。矩阵 $P$ 也可以采用权值的方式表示，即矩阵元素的取值表示参与的权值，例如回帖或转发的次数、所购买的商品数量、停留的时长等等。

在社会网络分析方法中，将上述社会网络结构称为双模网络，又称双重网络 （bipartite）或"隶属网"。

对双模数据的分析方法已经有很多，既有直接对双模网络的关联模式进行分析的方法，也有采用降模法将双模网络变为单模网络的分析方法。对于行为观测数据而言，研究者通常只关心参与者这个模所潜在的社会网络模式，事件模通常作为协变量来考虑。将二进制双模数据做一个简单的矩阵映射，就可得到一个单模矩阵，又称共生矩阵（co-occurrence matrix）。这个单模矩阵既可以是表示列模的数据（$P^{T} \times P$，$m \times m$），也可以是表示行模的数据（$P \times P^{T}$，$n \times n$），取决于研究者的需要。单模矩阵的数值则表示参与者两两之间共同出现的频次，因此该矩阵可视为有权网。从计算的角度来说，使用矩阵映射进行降模可以得到参与者的行为频次信息，因此非常简单高效。

然而对于带权值的发生矩阵，降模映射则不太适用。并且，其最大的缺点在于降模之后所得到的是一个密集矩阵。在原双模网络中，若节点的中心度为d的话，则降模之后变为$(d \times (d - 1)/2$，从而放大了网络密度，并在某种程度上扭曲网络结构(Latapy & Magnien et al., 2008)。因此，通常来说还需进一步处理，将其转换为二进制以作为社会网络关系的量度。然而，这种二值化量度，其麻烦在于如何确定阈值。若采用一个单一阈值对矩阵权值进行转换，则其潜在假设是参与者有相同的分布。举个例子，在网络社区中，网络成员之间的发帖数和回复数存在非常大的变差。假设成员 $x$ 与 $z$ 共有3次回复，成员 $y$ 与 $z$ 共有5次回复，但$x$的总回帖量是3次，而$y$的总回帖量是100次。我们该如何判定或比较x-z和y-z这两对关系呢？显而易见，基于单一阈值的二值化网络测度难以处理此类情形。在此基础上，可以考虑相对比例，例如3/3与5/100，巴拉特等人就分别提出改进型的加权方案(Barrat & Barthélemy et al., 2004; Newman, 2004; Barthélemy & Barrat et al., 2005)，但其重点往往在于探测网络的社区结构，而非节点间的关系测度。

除了降模映射法，典型的方式包括直接计算列模的相关系数矩阵并作为社会网络测量。相关系数法的优点在于能够控制不同参与者的活跃程度，但无法识别虚假相关，同时相关系数矩阵作为一个稠密矩阵也不太适合作为大规模网络的测量。也有学者将双模数据当作"购物篮"问题，采用数据挖掘手段来发现列模之间的关联模式(Raeder & Chawla, 2011; Zweig & Kaufmann, 2011)，然而其可信度和解释力不太能够得到保证。

在过去十几年中，在许多学科中，特别是在生物学(Friedman & Linial et al., 2000)、基因学(Ghazalpour & Doss et al., 2006)、神经科学(Huang & Li et al., 2010)等领域，图模型已经成为非常流行的对复杂系统进行抽象，并获得关于大规模观测变量的关联模式的一种处理手段。相比起前述的所介绍的降模映射法、相关系数法等处理方法而言，图模型的计算结果除了避免了前述几种处理方法的缺点之外，能够较好地探测出真实网络结构特征，还具有可解释性强、扩展性高的特点，在面对不同问题时能够展现出强大的解决能力。然而在社会科学领域，相关的研究尚不多见，仅有个别学者用图模型研究美国参议院投票网络、在线论坛发帖网络(陈华珊, 2015)等。相比起图模型在自然科学领域应用的流行性，社会科学领域对它的认识和使用还非常粗浅。在此，本文尝试对图模型进行一番引介，以期引起社会科学界同仁重视，并推动相关研究与应用。

高斯图模型

高斯图模型的基本形式

将观测数据的发生矩阵用一个$n \times p$的矩阵$X$来表示，$n$为观测数，$p$为变量数，假设观测之间相互独立，且$X$为多元正态分布随机变量，

$$X = (X_{1},...,X_{p}) \sim N(\mu,\Sigma)\backslash n$$

在此，对矩阵$X$中节点的两两关系的估计，也称为"邻域选择"（neighborhood selection），其实质是协方差选择问题。邻域选择的目的是对于给定的$n$个i.i.d观测$X$，分别估计每个变量（节点）的相邻变量。即对于集合$\Gamma$中的一个节点$a$，$a \in \Gamma$，它的邻域变量集合用 $X_{ne_{a}}$ 表示，邻域选择的目标是让 $X_{ne_{a}}$ 成为 $\Gamma\backslash{ a}$ 的一个最小子集，使得对于给定的 $X_{ne_{a}}$，$X_{a}$ 条件独立于所有其它变量。从而邻域选择可以被转化为标准的回归问题并求解。

但在数学求解上，一般不直接计算协方差矩阵，而是估计其逆协方差矩阵。这是因为逆协方差矩阵具有独特的性质。假设一个从多元正态分布中独立抽取的$n$个样本，其协方差矩阵为$\Sigma$，则表征样本变量之间条件依赖关系的高斯图模型可由逆协方差矩阵$\Theta = \Sigma^{- 1}$来表示。首先，逆协方差矩阵$\Theta$与协方差矩阵$\Sigma$具有对偶性，由于协方差矩阵为正定矩阵，那么逆协方差矩阵也为正定矩阵，因此它们互为对偶范数（dual norm）。其次，逆协方差矩阵具有稀疏的特质(Mardia & Kent et al., 1980; Lauritzen, 1996)，也就是说，当且仅当$\Sigma_{\text{ij}}^{- 1} = 0$时，变量$i$与变量$j$条件独立，反之，变量$i$与变量$j$存在条件依赖关系。逆协方差矩阵在图模型中又称"精度矩阵"（Precision Matrix）或"聚集矩阵"（Concentration Matrix）。

图：逆协方差矩阵与协方差矩阵示例

$\Theta$矩阵与偏相关系数有如下关系：

$$\rho_{ij \mid { i,j}} = - \frac{\omega_{\text{ij}}}{\sqrt{\omega_{\text{ii}}\omega_{\text{jj}}}}$$

对于社会关系网络测量来说，当该偏相关系数矩阵的元素大于0，即表示所对应的两个网络节点之间存在联带关系，且该数值可表示联带关系的强弱；反之则不存在联带关系。因此，根据观测数据$X$，计算的步骤为：估计其样本逆协方差矩阵，再转换为偏相关系数矩阵，就可得到该网络的关系测度。

一般采用最大似然法来估计"精度矩阵"$\Sigma^{- 1}$。用$S$表示$X$的经验协方差矩阵，高斯对数最大似然的公式表达如下：

$$\text{log\ }\text{det\ Θ} - trace(\text{SΘ})$$ （1）

其中 $\Theta$ 表示逆协方差矩阵，即 $\Theta = \Sigma^{- 1}$。使公式一最大化可得最大似然估计$\widehat{\Theta} = S^{- 1}$。但是对于大规模观测数据来说，存在两个基本特征：1、高维性。社会网络数据通常包含大量的节点（变量），用矩阵表示的话就是变量数$p$大于观测数$n$，在此情况下，经验协方差矩阵S为奇异矩阵，并不可逆，从而无法估计$\Theta$矩阵。即使 $p \approx n$ ，并且$S$不为奇异矩阵，$\Theta$ 的最大似然估计也会由于过高的方差从而失去效力；2、稀疏性。用图模型所表示的社会网络数据，存在大量的两两条件独立变量，即 $\Theta$ 中存在很多0元素。而根据使公式一最大化所估计得到的 $\text{Θ\ }$一般来说不存在值为0的元素。基于这两个性质，样本协方差矩阵不可逆，估计逆协方差矩阵时存在不稳定、计算成本高和不精确等问题。

罚似然估计法

罚似然估计法

近几十年来，统计学家针对高维稀疏数据提出了很多解决方案，其中蒂施莱尼所提出的罚似然回归法(Tibshirani, 1996)取得了主流地位，并被其它研究者进一步扩展和引进到图模型中(Meinshausen & Bühlmann, 2006; Yuan & Lin, 2007; Peng & Wang et al., 2009)。罚似然法是在线性回归公式中引入一个约束项 （regularizer）或惩罚项（penalty term）$\Theta$，并由一个非负的优化参数(tuning parameter) $\lambda$来控制，当 $\lambda$ 足够大时，$\Theta$的一些元素的值将等于零，也就是说 $\lambda$ 值越大，所估计的逆协方差矩阵越稀疏。并且，即使在 $p \gg n$ 的情形下，公式仍然能够求解，其表达式如下：

$$\text{maximiz}e_{\Theta}{ log\text{detΘ} - trace(\text{SΘ}) - \lambda \mid \mid \Theta \mid \mid_{1}}$$

其中，$\mid \mid \Theta \mid \mid_{1}$ 为 $l_{1}$ 罚则，表示对矩阵$\Theta$的所有元素的绝对值求和。将公式二用社会统计学教材常用的残差最小化拟合公式写法来表示，就是将：

$$RSS=\sum_{i=1}^{n}(y_i - \beta_0 - \sum_{j=1}^{p}\beta_jx_{ij})^2$$

改写为：

$$RSS + \lambda \sum_{j=1}^{p} |(\beta_j)|$$

在上式中，当$\lambda = 0$时，即为常规的OLS回归残差项。由于$\lambda$非负，因此当整个回归模型保留的变量越多，残差惩罚越大，反之则残差惩罚越小，从而$\lambda$作为模型超参数能够控制模型中所保留变量的稀疏程度。

梅豪森和布尔曼最早将罚似然回归应用到图模型中(Meinshausen & Bühlmann, 2006)，他们实际上是将网络的每一个节点作为因变量，其它所有节点作为自变量来构建一系列（$p$个）回归方程，从而得到一个近似解。其后，许多研究者提出了不同的求解法。有的学者借用万德伯格等人所提出的"内点搜索法"（interior-point）(Vandenberghe & Boyd et al., 1998)进行求解(Yuan & Lin, 2007)。贝纳杰等人则提出用"分块坐标递降法"（blockwise coordinate descent approach）来求解(Banerjee & El Ghaoui et al., 2008)，弗里德曼和哈斯蒂等人在此基础上进一步提出用坐标递降法（coordinate descent procedure）来求解(Friedman & Hastie et al., 2008)，并证明，当$p>n$的时候，坐标递降法具有很高的计算效率。

所有采用罚则对图模型进行稀疏求解的算法都可被称为图形罚极大似然法或罚似然图模型（以下简称glasso或图模型）。glasso模型近年来在基因研究、流行病学等领域有了很多应用研究，并且模型进一步从单一高斯图模型扩展成动态图模型(Ahmed & Xing, 2009; Song & Kolar et al., 2009)、多组图模型(Guo & Levina et al., 2011; Danaher & Wang et al., 2014)以及多层次图模型和潜变量图模型(Ambroise & Chiquet et al., 2009; Chandrasekaran & Parrilo et al., 2012)等。本文将在第二节展开介绍扩展模型。

最优参数选择与模型评估

在公式二中，参数$\lambda$未知，且无法通过样本数据进行推断，因此$\lambda$也称为超参数（hyper parameter），一般采用穷举法（对于多个超参数则可使用网格搜索（Grid Search）等方法）进行搜索。为了更好地评估模型以及避免模型的过度拟合，在机器学习理论中，一般采用交叉验证（Cross Validation）的方式来进行，即将样本数据集分为训练集和测试集，前者用以建立模型，后者则用来评估模型对未知样本进行预测时的精确度。也有学者采用贝叶斯信息准则（BIC）来评估模型，并针对稀疏约束的特点提出扩展贝叶斯信息准则（eBIC）(Chen & Chen, 2008; Foygel & Drton, 2010)。

应用与示例

在社会科学领域，最为著名的数据集恐怕要算美国南方黑人妇女数据集以下简称DGG），该数据集被很多研究者所使用(Freeman, 2003; Neal, 2013)。该数据由人类学家戴维斯和加纳等人使用访谈、观察记录、访客名单以及报纸记载来收集社区妇女参与社区活动的信息。他们的数据包括18名参与者，14次社会事件。研究者们用他们的人类学观察直觉及经验洞察力对这些妇女的社会网络进行了归纳，把她们分成两个子群体，并且在每个组中区分出核心成员、主要成员和边缘成员三个层次。在他们汇报的结果中，编号1至编号8的妇女被分到第一组，其中编号1、2、3、4被称为核心成员，编号5、6、7为主要成员，编号8为边缘成员。编号10到18被归为第二组，其中编号13、14、15是核心成员，11、12为主要成员，编号10、16、17、18为边缘成员。编号9被标识为同时属于两个组，且都作为边缘成员。

表1 美国南方妇女社会活动日常参与记录数据（原始版）：

表1.1 美国南方妇女社会活动日常参与记录数据（整理版）：

数据来源： (Davis & Gardner et al., 1941)

根据罚似然图模型计算结果，可有两种方式来构建社会关系网络矩阵。方式一为根据所估计的样本逆协方差矩阵，将非零元素转换为1，就可得到常规的社会关系网络表示矩阵，用这种测量方式所得到的网络为无向网络。方式二为根据样本逆协方差矩阵进一步计算偏相关系数矩阵，作为社会关系网络的测量，其中偏相关系数可作为关系的权重，由此，可得到无向有权网络（undirected valued network）。在实际应用中，上述方式所得到的关系矩阵很可能不是对称矩阵，还需做对称化处理。对于样本逆协方差矩阵可采用"或法则"（OR rule，即矩阵中每一对对角元素若任一个值不等于零则视为存在条件依赖）或"且法则"（AND rule，每一对对角元素均不等于零才视为存在条件依赖）；对于偏相关系数矩阵则可对每一对对角元素求平均值、最大值或最小值等方式来处理。

使用glasso法对这18位妇女的社会关系网络判定的结果见图一，随着罚则系数 $\rho$ 值的增大，所估计的网络密度愈加稀疏。根据eBIC法则，选择$\rho=0.1$的模型为最优模型。可以看到，图1-4区分出三个群体：编号1至7为一组，编号8和编号10至16为一组，编号17和18为第三组。编号9被判定为同时属于两个组，也就是说它承担了网桥的作用，连接两个群体。弗里曼(Freeman, 2003)汇总了21中计算方法对DGG数据进行元分析，glasso法的判定结果与这21种分析方法的绝大多数判定结果是一致的。稍微有所不同的是，glasso法单独将编号17和编号18两人判定为第三个组别，从原始数据上可以看到，她们两人仅共同出席了两次活动。在弗里曼所进行的分析中，BGR74 和OSB00 这两个方法也都将他们判定为单独的组别，在戴维斯和加纳的人类学分析中虽然将他们与编号10至16合为一组，但是将她们判定为边缘成员。由此可见，glasso法对于小群体估计也是具有敏感性的。

图1：用glasso法计算的美国南方妇女网络关联

代码如下：

# DGG 数据及示例

## 数据载入
library(igraph)
dgg<-read.csv(file = 'DGG.csv', header=TRUE)
dgg[is.na(dgg)]<-0
dgg<-t(dgg[,-1])
colnames(dgg)<-c(1:18)
dim(dgg)
dgg

## 基本 glasso 模型

library(glasso)

result = glasso(cov(dgg), rho = 0.1, penalize.diagonal=F)

# 计算偏相关系数矩阵
wi_to_pcor <- function(wi) {
  p = -wi
  d = 1/sqrt(diag(wi))
  pcor = diag(d)%*%p%*%diag(d)
  diag(pcor) = 0
  colnames(pcor) <- seq_len(ncol(pcor))
  pcor
}

pcor = wi_to_pcor(result$wi)

# 偏相关系数矩阵
pcor

g = graph.adjacency(pcor > 0, mode="undirected", diag = FALSE)
plot(g, vertex.color = 'white', vertex.size = 20, edge.color = 'black')

rhos <- c(0.05, 0.08, 0.1, 0.2, 0.8)
ret <- lapply(rhos, glasso, s = cov(dgg), penalize.diagonal = F)

wi <- lapply(ret, getElement, 'wi')
ret.pcor <- lapply(wi, wi_to_pcor)
ret.mat <- lapply(ret.pcor, function(x) x > 0)
ret.g <- lapply(ret.mat, graph.adjacency, mode="undirected", diag = FALSE)

for (i in seq_along(ret.g)) {
    plot(ret.g[[i]], vertex.color = 'white', vertex.size = 20, edge.color = 'black', main = sprintf("rho = %.2f", rhos[i]))
}

#' @param x list of glasso result
get_ebic <- function(x, n, gamma = 0.01) {
  df <- (sum(abs(x$wi) > 1e-5) - ncol(x$wi)) / 2
  d <- ncol(x$wi)
  -2 * x$loglik + df * log(n) + 4 * df * gamma * log(d)
}
ebic <- sapply(ret, get_ebic, n = 14, gamma = 0.1)
plot(rhos[1:4], ebic[1:4], type='l', xlab='rho', ylab='eBIC')

保真的读写方式

这是什么意思呢？先看两个代码示例：

示例一

来自 sjPlot 的示例，青椒群主丁哥也提到了（青椒），不过遇到障碍就略过，不符合我们 rebuild the wheel 的风格。中文是我的第二母语，我不跟中文较劲谁来较劲。

library(sjPlot)
library(sjmisc)
data(efc)

# 我要跟中文较劲之1
attr(efc$c160age, 'label') <- '年龄'
fit1 <- lm(barthtot ~ c160age + c12hour + c161sex + c172code, data = efc)
sjt.lm(fit1)

# 我要跟中文较劲之2
attr(efc$c160age, 'label') <- iconv('年龄', to = 'UTF-8')
fit1 <- lm(barthtot ~ c160age + c12hour + c161sex + c172code, data = efc)
sjt.lm(fit1)

显示的表格是这样的：

统统可耻地失败了！！！夜深人静之时，熬到双眼通红，改代码改到崩溃也没有解决，对不对？

示例二

示例二需要一些工具准备工作。博主用的是 UltraEdit 这款文本处理专业助手。顺便说一句，做数据处理的人，必须得心应手的运用一款这样的专业工具，否则就是毛笔刷水桶，筷子搅大缸，根本就是数据处理的门外汉。

不同编码的真实面目

这里用 UltraEdit 来干什么呢？查看文本的字符编码。编码的原理不再重复介绍了，网上资源大把。就看几个典型的文本编码是什么样子：

图1：这是人畜无害的记事本显示的文本

图2：十六进制下 GB 码的样子

图3：十六进制下 UTF-8 没有 BOM 的样子

图4：十六进制下 UTF-8 有 BOM 的样子，多了 EFBBBF 的标记

其它的编码比如 UTF-16 这些就不再一一展示了。我们就是要用 UltraEdit 来帮助R检查输入输出有没有问题，确保 UltraEdit 看到的编码和 R 读写的一致。目标很简单，保真，也就是要求 R 读进来的文本编码与原始的一致，R 写出来的文本编码与 R 内存中的一致。对于有一点编程经验的人来说，这个要求似乎一点都不高。对于随机数都是“伪”的软件编程来说，消除不确定性不是天经地义的嘛？

示例二

事实真是如此嘛？让我们检验一下。

options('encoding')
## $encoding
## [1] "native.enc"

(x1 <- readLines('datacarpenter_ansi.txt'))
## [1] "数据匠"
(x2 <- readLines('datacarpenter_utf8.txt', encoding= 'UTF-8'))
## [1] "数据匠"
(x3 <- readLines('datacarpenter_utf8bom.txt', encoding = 'UTF-8'))
## [1] "数据匠"
Encoding(x1)
## [1] "unknown"
Encoding(x2)
## [1] "UTF-8"
Encoding(x3)
## [1] "UTF-8"

目前为止还都符合预期。那么在 UTF-8 系统下是什么样子呢？

Sys.setlocale('LC_CTYPE', locale = "English_United States.1252")
## [1] "English_United States.1252"
options(encoding='UTF-8')
 
(x1 <- readLines('datacarpenter_ansi.txt', encoding = 'unknown'))
## [1] ""
## Warning messages:
## 1: In readLines("datacarpenter_ansi.txt", encoding = "unknown") :
##   invalid input found on input connection 'datacarpenter_ansi.txt'
## 2: In readLines("datacarpenter_ansi.txt", encoding = "unknown") :
##   line 1 appears to contain an embedded nul
(x2 <- readLines('datacarpenter_utf8.txt', encoding= 'UTF-8'))
## [1] "?"
## Warning messages:
## 1: In readLines("datacarpenter_utf8.txt", encoding = "UTF-8") :
##   invalid input found on input connection 'datacarpenter_utf8.txt'
## 2: In readLines("datacarpenter_utf8.txt", encoding = "UTF-8") :
##  incomplete final line found on 'datacarpenter_utf8.txt'
Encoding(x1)
## [1] "unknown"
Encoding(x2)
## [1] "unknown"

惨不忍睹啊！还有一点，再次告诉我们，Windows 本身就是支持Unicode的系统，只要系统的语言和区域设置正确（在控制面板中），就用不着设为统一的 UTF-8 体系，否则就是摁下一个小葫芦起个超级大瓢。

那么写入呢：

x1 <- '数据匠'
writeLines(x1, 'chn_ansi.txt')

x2 <- iconv(x1, to = 'UTF-8')
writeLines(x2, 'chn_utf8.txt')

看上去还不错，无论是 GB 码还是 UTF-8 码都显示正确。

除了一点，为何这两个文件的大小是一模一样的的呢？用 UltraEdit 查看它们的十六进制内容，也是一模一样的，有没有！

喜欢读帮助文档的读者就会提醒了，加上 useByte 参数可破！确实，useByte = TRUE 参数可确保写入 UTF-8 文本，不过有个前提，就是当编码与 options('encoding') 不一致的时候。如果碰到改了 options() 参数的时候，就不一样了：

options('encoding' = 'UTF-8')
x = iconv('中华人民共和国', to = 'UTF-8')

# 写入乱码
writeLines(x, con = 'test.txt', useBytes = T)
## Warning message:
## In writeLines(x, con = "test.txt", useBytes = T) :
##   invalid char string in output conversion
readLines("test.txt", encoding="UTF-8")
## [1] "中华人民共和\xe5\u009b"

在上面这个例子中，文本是 UTF-8 编码，options() 设置是 UTF-8 编码，写入是字节码，一切可能的设置都已经做到perfect，完美统一，可是忒么的写到磁盘就是乱码！如果你是个较真的人，是不是会觉得很别扭，很糟心？

所以，尽管针对每个特定的情形，都能找到相应的解决方案，但是这些方案都取决于特定的系统设置。因此，同样的代码在不同人的设备不同的环境下面跑，出错就在所难免了。并且，像上面这样的代码书写起来不是非常的单调、重复、丑陋无比嘛？

那么不用 readLines() writeLines() 这些函数就好了嘛，还有更高效的 data.table 的 fread()，readr 的各种 read_*** 函数可用。是这样，仅仅只是数据导入还好办，不过大量的输出相关 R 包，比如示例一中的 sjplot、stargazer，幻灯输出的 slidify，动态网页和D3.js 输出的 clickme等等，大量的包都依赖于 readLines()、writeLines()以及cat() 这些基础函数。用 R 就是要享受这种全流程工作链的感觉。想象一下，把你的博士论文全程用 R 计算并生成，这真的不是幻想。

本系列博文的目标就是要找到一种可靠的不受各种外部系统设置影响的方式来实现保真的编码读写，同时，要有足够的兼容性，通过一点简单的操作就能一举消除大量的第三方包所存在的乱码 bug，无需一个个去提醒各个包的开发者去做修改。你能体会跟老外去沟通，让他们去理解这些乱码问题是多么费劲的一件事情嘛？更何况下个版本升级不小心又改回去了呢？以上都是博主遭遇过的血泪史。

这个泡泡是不是吹的有点大了？

在R教学中，首先要跨过去的一座大山就是乱码问题。很多学生在装好R和RStudio之后，刚刚运行RStudio，还未尝个鲜，写出R的第一段甚至人生第一段hello world代码，乱码就来立个下马威了（此处应该有乱码翻车现场截图，但对系统干净清爽的博主来说要制作这么一张图还是有难度的，哪位读者愿意贡献一张？）。

对于许多从 Stata 转过来的社科领域研究人员来说，在学习和使用R的过程中是如此频繁地遭遇乱码问题更是一件很不可思议的事情。乱码无处不在，令人抓狂，令人崩溃，茶饭不思，错过DDL，成为一块死肉。的确，相对于封闭的Stata来说，作为开放系统的R不仅要处理作为数据的文本，还需同各种系统进行数据交换，结果输出到各类格式中，必须要支持各种字符编码。然而文本编码作为计算机领域的专门知识，不仅普通使用者不懂，大多数R包的开发者也不甚了解，代码可能写的不规范，再加上 R 语言与 c 语言的交互调用，于是就出现了异常复杂的R生态乱码问题。不仅导入数据时乱码，输出结果时乱码，制图时乱码，甚至原本正确的编码到了中间某一步就乱码了。

博主作为资深文科生，今天却要自不量力地对R的乱码问题做一些解析，提供一些解决思路，以飨后来者。

先让RStudio运行起来

对于那些第一次运行RStudio就遭遇车祸的童鞋来说，大部分都是 Windows 用户，而且 Windows 的账户用的是中文名。

在此，先强调两个傻瓜原则，能给将来避免很多麻烦：

• R 以及 RStudio 的安装路径不要有中文和空格，必须全英文。
• Windows 账户要用英文！

那么已经设置了中文帐户名的怎么办？

• 一劳永逸的办法，当然是新建一个英文账户登录了；
• 若是不会上面的操作，或者觉得切换账户太麻烦，可以改一下 Windows 的临时目录：

右键“我的电脑”-“属性”-“高级系统设置”-“环境变量”，将 temp 和 tmp 都设为 c:\windows\temp\（下图红框所示）：

一路确定之后，重启 RStudio 就可以去敲写人生第一行R代码了。

那么，在 MacOS 下翻车的用户怎么办？博主不用 Mac，所以现在没法复原车祸现场，请快递一台 MacBook Pro 以便博主更好地服务大众。

R 的字符编码

在介绍R的字符编码之前，还是得有一些关于文本编码的基本知识，建议先读读这篇，包括下面的评论。

这个世界上的文本编码方案千百种，有些 linux 的支持者及R包开发者号称要用 UTF-8 一统天下，可惜这只能是一厢情愿，因为很多字符超出了 UTF-8所能涵盖的范围。不得不承认，大多数R包开发者用的是 linux 或 Mac 系统，得益于以 UTF-8 作为整个体系的基础，避免了很多乱码的情形。但也正是他们写的不规范，才导致了 Windows 下罄竹难书的灾难性后果啊。字符编码不仅有方案之争，还有系统之争，关于 BOM 的吵架就有一罗锅。作为实用主义者的文科僧，博主只关心如何保证文本编码的正确传递、显示和输出。

先提供两个基本的函数：

1、纠正乱码的一招鲜

碰到乱码，一般可以先尝试一下，给它设置一个正确的编码：

Encoding(someX)<-'UTF-8'
# 或者
Encoding(someX)<-'GB2312'

前一句告诉R，这是个 UTF-8编码的字符，后一句则是告诉 R 这是简体中文国标码编码的字符，你给我按照这个方式打印出来！如果这一招能搞定，则万事大吉，继续码农之旅吧。否则，则真正的麻烦开始了。

由此可见，R 之所以能正确处理字符编码，靠的是给不同编码方案的字符提供一个标记，例如上面代码中的 UTF-8 以及 GB2312。若是这个标记设置错误，不仅显示乱码，后续的字符查找、匹配、替换等操作也会跟着出问题。例如下面这段代码强行将 x 的编码设为国标码，就造成了显示输出的乱码。

Encoding(x)
## [1] "UTF-8"
Encoding(x) <- 'GB2312'
x
## [1] "涓枃"

当然，对于这种编码设置错误的情形，不必过于担心，在对字符进行任何__改动__操作之前都是安全的，只需将正确的编码标记告诉 R 即可。

2、编码转换引擎

已知一个字符对象的编码，需要将它转换为别的编码，用 iconv() 实现：

x <- '中文'
Encoding(x)
## [1] "unknown"
x2 <- iconv(x, 'GB2312', 'UTF-8')
Encoding(x2)
## [1] "UTF-8"

iconv() 有两个比较重要的参数：from 指原始的编码方案，to 指目标编码方案。iconv() 也不管 x 原来是什么编码，只会按照 from 的指示行事。

3、如何判断字符的正确编码？

那么，若是导入数据时遇到乱码，除了用不同的编码值去猜测之外，有没有什么办法去自动识别字符的编码呢？从编码规则上来说似乎并没有非常确切的识别方法（结合编码规则以及文本规律，可以考虑tellenc），除了那些带有 BOM 标识的 Unicode 文本。根据文件头的 BOM 标识，可以自动识别文本是 UTF-8 还是 UTF-16，是大头（big endian）还是小头（little endian）。听上去很完美对不对？可惜 linux 的拥趸不这么干。于是，我们要面对这个混杂的世界。

更重要的的是，R 对编码具有很好的包容能力，尽管目前力有不逮但仍努力在背后去自动应对不同的编码，甚至混合编码。在第二篇，再来解决混合字符编码的问题。我们首先要知道，在同一个 R 字符对象中，是可以有不同编码方案的。

x <- '中文'
x[2] <- iconv('中文', to = 'UTF-8')
x
## [1] "中文" "中文"
Encoding(x)
## [1] "unknown" "UTF-8" 

那么最坏的情况就是我们需要导入一个外部数据文本，可能是国标码，也可能是UTF-8编码，甚至是混合的。万幸的是我们可以对 UTF-8 做一个小小的自动判定（大约99%的情况下是对的）。来看看一个文本导入的示例。为了这个示例，博主特意准备了一个比较乱的原始数据，既包含国标码的中文也包括 UTF-8 编码的中文。同时，再随机地设置其中一些为 UTF-8。这样，测试数据中无论是国标码的文本还是 UTF-8 编码的文本同时都有编码声明正确和错误的两种情形。

xx <- readLines("e:/Huashan/Stat/R/Demo/encoding/data/mixed.txt")
xx
## [1] "这是一段中文gb码"             "杩欐槸涓€娈典腑鏂噓tf8鐮\x81"
# 复制10遍再乱序
xx <- rep(xx, times=10)
Encoding(xx)[sample(length(xx))] <- 'UTF-8'
xx <- xx[sample(length(xx))]

Encoding(xx) <- 'GB2312'
Encoding(xx)[huashan::isUTF8(xx)] <- 'UTF-8'
xx
##  [1] "这是一段中文gb码"   "这是一段中文gb码"   "这是一段中文gb码"  
##  [4] "这是一段中文utf8码" "这是一段中文gb码"   "这是一段中文utf8码"
##  [7] "这是一段中文utf8码" "这是一段中文gb码"   "这是一段中文utf8码"
## [10] "这是一段中文utf8码" "这是一段中文gb码"   "这是一段中文gb码"  
## [13] "这是一段中文gb码"   "这是一段中文gb码"   "这是一段中文utf8码"
## [16] "这是一段中文utf8码" "这是一段中文gb码"   "这是一段中文utf8码"
## [19] "这是一段中文utf8码" "这是一段中文utf8码"

对了，秘诀就是 huashan 包中的 isUTF8() 函数。

4、默认字符编码

让我们再看这段最简单的示例：

x <- '中文'
Encoding(x)
## [1] "unknown"

上述命令的结果可能各有不同，一般情况下，Windows 用户得到 unknown 而 Mac 用户得到 UTF-8。如果 Mac 用户得到的结果不是 UTF-8，那么恭喜你，你可能会有麻烦了。

上述代码传达了两个信息：

• 同样的代码，在不同的系统中跑，可能是两个结果。对于不喜欢不确定性的程序员来说，这是一件可怕的事情。类似的情形在上机课或者交作业中也经常见到，往往是同学B找同学A要解决方案，结果发现在自己的机上跑不过，连抄个代码都忒么这么为难我！
• 万幸的是，R 的默认字符编码取决于某种系统设置。

这第二点在于操作系统的区域设置，和 Sys.getlocale()、Sys.setlocale() 这两个函数有关。

Sys.getlocale()
## [1] "LC_COLLATE=Chinese (Simplified)_China.936;LC_CTYPE=Chinese (Simplified)_China.936;LC_MONETARY=Chinese (Simplified)_China.936;LC_NUMERIC=C;LC_TIME=Chinese (Simplified)_China.936"
x <- '中文'
Encoding(x)
## [1] "unknown"

Sys.setlocale('LC_CTYPE', locale = "English_United States.1252")
## [1] "English_United States.1252"
x2 <- '中文'
Encoding(x2)
## [1] "UTF-8"

尽管通过更改区域设置实现了创建一个默认字符编码为 UTF-8的字符对象x2，输出也没有问题，但是作为本地编码的 x 却不能正常显示了。由此可见，在默认 UTF-8 编码的系统下，只能正确显示 UTF-8 一种编码，其它编码都不能正常显示，必须通通转成 UTF-8 才行。Mac 系统就是如此。但在 Windows 下作为国标编码的中文和作为 UTF-8 编码的中文均能同时正常显示。这是一种优点，也是缺点，后续将有例子展示。

另外，在RStudio中，x2 的显示也不那么让人舒心，如下图。除非将 Sys.setlocale() 设在 RStudio 启动之前。所以在 Windows 下，一般不建议更改系统的 locale。