在R Commander中导入数据
在经济学分析中不可避免的要和数据打交道，而目前数据分析中最主要的工具就是计量经济学。数据源于现实，而对待数据的态度方面，我更欣赏凯恩斯的观点：从数据中寻找直觉。既不是单纯的从计量的结果中寻求观点的佐证，也不是从归纳的角度来推理因果关系。这有些和“散点图是最好的统计图形”的观点有些不谋而合。但是数据本身的特性并不是简简单单的可以肉眼扫视原始数据就可以得出的，这个时候借助计量这个分析工具更有利于我们发现隐藏在原始数据背后的蛛丝马迹，进而寻求灵感。因此，玩转数据是做经济学研究必不可少的一个环节。有句话说得好：Let’s get our hands dirty with data first!
当然做计量的时候很依赖计算机软件，常用的有Eviews、Stata、SPSS、SAS等。可以看出，这和统计学中常用的软件惊人的一致。追根溯源，计量经济学本来就是从数理统计学中的回归分析等渐渐延伸出来的，所以其方法在统计软件中可以很容易的实现。近几年R的快速蓬勃发展使之成为了最前沿的统计软件，由于其良好的拓展性，大量的免费的包(package) 的出现使得R足以胜任最潮流的统计分析工作。因此，R也足以作为一个计量分析软件来处理计量经济学的问题。
我作为一个经济学专业的学生，机缘巧合接触到了R，并为之深深沉迷。2009年冬天在第二届中国R语言会议做了一个简单的“在计量和经济学中使用R”的报告后，感到有必要写一个简单的小册子，介绍各种计量经济学方法在R中的实现，也希望借此从丰富的实例数据中找寻更多的直觉。
这个小册子主要希望能对下列使用者有所裨益：
因水平所限，这本小册子将会比较简单，着重于介绍各种方法对应的R包和实现，帮助从未使用过R的朋友们尽快的熟悉、了解和应用这款软件。学习一个软件最好的方法无非是多多使用，因此除了囊括大量的实例，我想不出更好的办法。这些例子有些来源于现有出版的计量经济学书籍（例如伍德里奇的计量经济学导论2），也有些摘取于公开发表的论文。当然，这对我来说是一项浩瀚而繁重的工作，因此诸位朋友的帮助显得格外的珍贵。
从现有的关于计量经济学和R的书籍来说，从网上能找到几本英文的，大都是免费发行并具有非常高的质量。只是国内中文的资料还颇为零散。在撰写这个小册子的过程中，我参考了大量已有的成果并受益匪浅，也建议英文较好的朋友们直接去阅读相关的英文材料尤其是R包自带的介绍，相信会更深入的了解R。在这里特别要说的是AER（全称：Applied Econometrics with R）这个包，是配合同名的书发行的包。不过通过demo可以详尽的看到各个例子的R源代码，也带有丰富的数据集（来自格林的《计量经济分析》等有名的著作），是非常好的练手的包。
最后需要说明的是，这本小册子是我在担任统计之都中文论坛（http://cos.name/cn/）“经济统计版”版主的时候所撰写的。承蒙站长谢益辉兄和诸位骨干成员的大力帮助，此册子凝结了COS诸多成员的心血，换言之我只是一个代笔者而已。我们通过GIT这个多人协作平台共同完善，也借助了knitr包来结合R与markdown。这样高效且免费的开源平台使得我在撰写过程中受益匪浅，也使得本册子避免潜在的问题得以实现在互联网上的免费发行。
关于R的中文入门有很多，在这里不再一一枚举（也不是本册子的职责所在）。但我也不能在这里说，诸位客官回去学完R的基本结构再来吧。一是有点枯燥，二则我个人窃以为学习一个软件最好的办法就是不断的拿例子来磨练，遇到问题去网上搜寻，否则看再多的入门引导也只似过眼烟云，那些命令一会儿就抛到九霄云外了。所以我把会用到的一些东西直接融入在下面这个简单的例子中，有一些必要的说明，希望可以快速的熟悉R。
下面是来自Papke(1995)的一个例子。他研究的是一个退休金计划和计划的慷慨度。
在401K.DTA这个数据集中，我们关心两个变量。prate是在合法的工人中拥有活跃帐户的比例。mrate是用来衡量这个计划的匹配程度（用来代表慷慨度），即如果mrate = 0.5，则表示工人付出了$10，其工作单位相应的付出了$5。接下来，我们需要面对这么几个问题：
做计量分析，离不开的就是数据。所以我们第一步先来导入需要的数据。
获取数据有很多办法，在R 里面通过foreign包可以读/写Minitab, S, SAS, SPSS, Stata, Systat等等格式的数据。当然，R本身是支持从文本文件（包括CSV格式）和剪贴板中直接读取数据的。此外，对于R包里面自带的数据集，我们可以直接用data("name")来加载数据集。这里我采取的是读取Stata的数据（DTA格式）。
当然，我们首先要加载foreign包，可以在R中直接点击“加载程序包”，也可以手动输入：
library(foreign)然后就可以使用read.dta()命令：
Papke_1995 = read.dta("data/401K.DTA")## Warning message: cannot read factor labels from Stata 5 filessummary(Papke_1995)## prate mrate totpart totelg
## Min. : 3.0 Min. :0.010 Min. : 50 Min. : 51
## 1st Qu.: 78.0 1st Qu.:0.300 1st Qu.: 156 1st Qu.: 176
## Median : 95.7 Median :0.460 Median : 276 Median : 330
## Mean : 87.4 Mean :0.732 Mean : 1354 Mean : 1629
## 3rd Qu.:100.0 3rd Qu.:0.830 3rd Qu.: 750 3rd Qu.: 890
## Max. :100.0 Max. :4.910 Max. :58811 Max. :70429
## age totemp sole ltotemp
## Min. : 4.0 Min. : 58 Min. :0.000 Min. : 4.06
## 1st Qu.: 7.0 1st Qu.: 261 1st Qu.:0.000 1st Qu.: 5.57
## Median : 9.0 Median : 588 Median :0.000 Median : 6.38
## Mean :13.2 Mean : 3568 Mean :0.488 Mean : 6.69
## 3rd Qu.:18.0 3rd Qu.: 1804 3rd Qu.:1.000 3rd Qu.: 7.50
## Max. :51.0 Max. :144387 Max. :1.000 Max. :11.88 Papke_1995是我们赋值后在R里使用的数据表的名字。因为R是基于对象(object)的，所以我们需要在读取数据的时候指定数据存储的对象。同样的，后面会不断的用到对象这一概念。
如果觉得这些东西记起来比较麻烦，一个个字母的打起来也挺麻烦的，怎么办？好在有个包叫做Rcmdr。加载这个包之后就会出现图形界面，可以通过点击的方式来操作。
在R Commander中导入数据
之后，R Commander会自动记录每一步对应的代码，可供下次重复使用。
在介绍关于平均值的函数前，先介绍另一个有用的函数names()。这个函数的作用是显示数据表中所有的变量名称。用法和效果见后面的代码例子。
我们可以使用summary()来获取该数据表的摘要信息，里面包含平均值、最大最小值 、中位数等。不过我们这里只关心两个变量prate和mrate ，所以也可以使用numSummary()（需加载abind包）。
load("data/Papke_1995.rdata")
names(Papke_1995)## [1] "prate" "mrate" "totpart" "totelg" "age" "totemp" "sole"
## [8] "ltotemp"summary(Papke_1995)## prate mrate totpart totelg
## Min. : 3.0 Min. :0.010 Min. : 50 Min. : 51
## 1st Qu.: 78.0 1st Qu.:0.300 1st Qu.: 156 1st Qu.: 176
## Median : 95.7 Median :0.460 Median : 276 Median : 330
## Mean : 87.4 Mean :0.732 Mean : 1354 Mean : 1629
## 3rd Qu.:100.0 3rd Qu.:0.830 3rd Qu.: 750 3rd Qu.: 890
## Max. :100.0 Max. :4.910 Max. :58811 Max. :70429
## age totemp sole ltotemp
## Min. : 4.0 Min. : 58 Min. :0.000 Min. : 4.06
## 1st Qu.: 7.0 1st Qu.: 261 1st Qu.:0.000 1st Qu.: 5.57
## Median : 9.0 Median : 588 Median :0.000 Median : 6.38
## Mean :13.2 Mean : 3568 Mean :0.488 Mean : 6.69
## 3rd Qu.:18.0 3rd Qu.: 1804 3rd Qu.:1.000 3rd Qu.: 7.50
## Max. :51.0 Max. :144387 Max. :1.000 Max. :11.88 可以从上表中读出prate和mrate的平均值。
sumSummary()也可以通过R Commander的图形界面实现。 
在R里面进行线性回归还是比较容易的，直接使用lm()就可以。值得注意的是，由于R的面向对象特性，我们需要不断的赋值。对于赋值，有三种基本方法，分别可以用->、<-和=实现，其中前两个是有方向的赋值，所以一般来说更为常用。比如我们可以对变量mrate和prate 求乘积，并将结果赋予一个新变量mp，则只需写成mp <- mrate*prate。
因此在做回归的时候写成：
RegModel <- lm(prate ~ mrate, data = Papke_1995)
summary(RegModel)##
## Call:
## lm(formula = prate ~ mrate, data = Papke_1995)
##
## Residuals:
## Min 1Q Median 3Q Max
## -82.30 -8.18 5.18 12.71 16.81
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 83.075 0.563 147.5 <2e-16 ***
## mrate 5.861 0.527 11.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.1 on 1532 degrees of freedom
## Multiple R-squared: 0.0747, Adjusted R-squared: 0.0741
## F-statistic: 124 on 1 and 1532 DF, p-value: <2e-16
## 这样RegModel里面就存储了这次回归所得的数据。
我们还可以采用attach()命令，这样就不用每次都指定回归向量所在的数据集了，直接写RegModel<- lm(prate~mrate)，然后就可以用summary(RegModel)来看回归的结果了。(注：通常情况下不建议使用attach()，可能会导致变量名的一定程度混乱，尤其是在函数封装的时候。)
可以看出估计后的回归方程应为：
\(\hat{prate}=83.0755+5.8611mrate\)
其中\(R^{2}\)为0.0747。呃，这个\(R^{2}\)为什么这么小？看看散点图就知道了。
我们可以直接用最简单的plot()命令作图(当然更好的一个选择可能是ggplot2)，用法如下：
plot(Papke_1995$mrate, Papke_1995$prate)
abline(RegModel, col = "red")第二行命令是添加了那条回归拟合线。 可见这个图本来就很散，也难怪线性拟合效果这么差了。
最后，就是依赖估计方程做预测了。这里需要的是做一个点预测。R里面需要依据另一个数据集来预测，而且这个数据集中必须含有mrate 这个变量。新建一个数据集并赋值的办法有许多，最简单的就是直接赋值，方法如下：
mrate_new <- data.frame(mrate = 3.5)另外一种不推荐的方式是，利用数据编辑框来手动输入:mrate_new <- edit(as.data.frame(NULL))。不推荐的原因是难以有效的追溯数据的来源和数值，尤其是违背“可重复的研究”精神。
之后再利用predict()就可得到所需的预测值了。
mrate_new <- data.frame(mrate = 3.5)
predict(RegModel, mrate_new)## 1
## 103.6 当然现实中我们很少做一元的线性回归，解释变量往往是两个或者更多。这可以依旧用上面的lm()。如下面这个例子，研究的是出勤率和ACT测试成绩、学习成绩之间的关系。
(@ATTEND)
在ATTEND.DTA这个数据集中，atndrte 是出勤率（采用百分比表示），ACT 为ACT测试的成绩，priGPA 是之前的学习平均分。我们需要估计如下的方程：
\(atndrte=\beta_{0}+\beta_{1}priGPA+\beta_{2}ACT+u\)
很显然，这里我们和上面的例子一样，代码和结果如下：
library(foreign)
Attend <- read.dta("data/attend.dta")
Reg2 <- lm(atndrte ~ priGPA + ACT, data = Attend)
summary(Reg2)##
## Call:
## lm(formula = atndrte ~ priGPA + ACT, data = Attend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -65.37 -6.77 2.12 9.63 29.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 75.700 3.884 19.5 <2e-16 ***
## priGPA 17.261 1.083 15.9 <2e-16 ***
## ACT -1.717 0.169 -10.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.4 on 677 degrees of freedom
## Multiple R-squared: 0.291, Adjusted R-squared: 0.288
## F-statistic: 139 on 2 and 677 DF, p-value: <2e-16
## 虽然我们有 RCommander 创造的图形界面，但是每次都指定参数也是件很烦的事儿。因此养成一个好习惯，保存好上次运行的代码，下次直接在R里面调用就可以了，有什么修改的也只需要稍作调整即可。RCommander里面本身就有File -> Save Script，可以把Script Window里面所有代码存储为***.R的格式，从而方便下次调用。Script Window里面也是可以直接编辑代码的，删掉一些自己不想要的，调整个别的参数都是很方便的。
需要说明的是，.R文件就是告诉R应该怎么运行的文件，所以可以直接用文本编辑器软件打开并编辑。现在NotePad++, UltraEdit等等文本编辑软件都有支持R的插件，可以方便的把代码传送到R里面调用。R的基本界面中也是可以直接打开.R的脚本文件运行的。此外，推荐一个新兴的R编辑器——RStuidio，集成了R的各种窗口（
有了上述的例子，相信大家已经基本熟悉R了。那么遇到问题怎么办呢？比如summary()这个函数，对于不同的模型会有不同的用法，那么我们就需要去查看原始的帮助。在R中，最简单的办法就算再想要查看的命令前加一个?号。例如?summary之后就会蹦出来帮助页面了。这是查看某一包作者撰写原始文档的最快捷方式。此外也可以用两个连续的问号??来搜索所有相关的资料。
但是如果根本不知道有哪些命令，则需要去找包内原始的资料。可以直接在Google等搜索引擎里面搜寻，也可以查看R包自带的说明，亦可以参照各种书籍。总之方法很多，多多利用互联网是最好的办法。国内最佳的地方自然是统计之都论坛的R版，里面有丰富的资料和资深的UseR为大家解惑。
下面，我们将介绍一些R的基础知识，包括软件的安装、配置和基本的数据清理工作，为后续的分析打下坚实的基础（曾经某位自身的数据分析专家说过，“在数据正式进入模型之前，最关键的就是数据整理”。这一步可能会耗费很多时间，但把数据整理为理想的格式会为后面的统计分析提供极大的便利）。
R可以从www.r-project.org下载，里面点击CRAN后可以选择CHINA的几个镜像，然后依据操作系统选择对应的版本即可。
安装完R之后，最好再选择一个顺手的编辑器，用于书写对应的R代码。这里推荐的是RStudio，可以从rstudio.org下载。安装完成后，双击打开，界面如图[fig:RStudio]所示。
[fig:RStudio]RStudio界面
image
RStudio界面，分为四块。左上角是脚本编辑框和数据浏览框；右上角是当前空间中数据情况和使用命令的历史记录；左下角是实际的R的代码执行界面和相应返回的结果；右下角是文件目录列表、画图展示区、R包目录区和帮助区。
目前我们的工作仅仅搭好了R的架子，还需要依据分析任务下载并安装对应的R包。这部分内容在后续中会详细依据案例介绍。
在R中，默认的工作目录依系统配置而变化，可以在直接启动R之后，通过getwd()命令来查看。
getwd()另外，如果是通过后缀名为.R的脚本文件来直接调用R，那么工作目录就为该脚本文件所在的目录。如果对于任何命令的参数等希望得到进一步的说明，那么可以在命令前加上“?”来直接调用帮助。比如，
`?`(getwd)这个时候R会弹出帮助文档窗口。
在我们撰写R脚本之前，往往希望事前指定一个工作目录，这个时候就可以利用setwd(file_path).
R里面大多数操作都是面向对象的，所以上述代码的含义就是给file_path这个object赋上了"E:/example/"这个文本串作为其值。然后利用setwd()函数来指定工作目录。显然，这个代码和直接调用setwd("E:/example/")是等同的。
R作为一个开源软件，最大的特性就是有很多人在不断的贡献各种分析包，从基础的数据整理到高级的统计模型、可视化实现都有着相应的支持。R里面现在可以支持分析包括但不仅限于：
Bayesian： Bayesian Inference（贝叶斯推断）
ChemPhys：Chemometrics and Computational Physics（化学计量和计算物理）
ClinicalTrials：Clinical Trial Design, Monitoring, and Analysis（临床试验设计、监测和分析）
Cluster：Cluster Analysis & Finite Mixture Models（聚类分析和有限混合模型）
Distributions： Probability Distributions（概率分布）
Econometrics： Computational Econometrics（计量经济学）
Environmetrics： Analysis of Ecological and Environmental Data（生态和环境学数据分析）
ExperimentalDesign： Design of Experiments (DoE) & Analysis of Experimental Data（实验设计和分析）
Finance：Empirical Finance（实证金融分析）
Genetics：Statistical Genetics（统计遗传学）
Graphics：Graphic Displays & Dynamic Graphics & Graphic Devices & Visualization（图形显示、动态图形、图形设备和可视化）
HighPerformanceComputing： High-Performance and Parallel Computing with R（高性能并行计算）
MachineLearning：Machine Learning & Statistical Learning（机器学习和统计学习）
MedicalImaging：Medical Image Analysis（医学图像分析）
Multivariate：Multivariate Statistics（多元统计）
NaturalLanguageProcessing：Natural Language Processing（自然语言处理）
OfficialStatistics：Official Statistics & Survey Methodology（官方统计和普查方法）
Optimization：Optimization and Mathematical Programming（最优化和数学规划）
Pharmacokinetics：Analysis of Pharmacokinetic Data（药代动力学数据分析）
Phylogenetics：Phylogenetics, Especially Comparative Methods（系统发育、特别是比较方法）
Psychometrics：Psychometric Models and Methods（心理测量模型和方法）
ReproducibleResearch：Reproducible Research（可重复的研究）
Robust：Robust Statistical Methods（稳健性统计模型）
SocialSciences：Statistics for the Social Sciences（社会科学统计）
Spatial：Analysis of Spatial Data（空间数据分析）
Survival：Survival Analysis（生存分析）
TimeSeries：Time Series Analysis（时间序列分析）
gR：gRaphical Models in R（图形模型）
每一项的具体说明请参见：http://cran.r-project.org/web/views/。
针对每一项任务，我们往往都需要加载不同的包。在第一次加载一个包以前，我们往往需要先安装。在Rstudio界面中，加载包的命令位于Tools -> Install Packages，然后可以输入包的名字来安装，如图[install packages]所示。
[install packages]安装包（在Rstudio中）
image
image
安装完包之后，可以直接通过library()命令加载。如我们常用的data.table这个包：
library(data.table)注意：R里面严格区分大小写，所以大小写不对的话包是无法加载成功的。
CRAN下的Task Views不仅仅罗列了各个主题下的R包情况，其本身亦作为一个R包ctv出现。在R中安装这个包之后，可以很方便的安装所有主题下的包。比如，我们希望安装Econometrics这个主题下所有提及的包，那么:
install.packages("ctv")
library("ctv")
install.views("Econometrics") #安装Econometrics主题下的包
update.views("Econometrics") #更新Econometrics主题下的包对于刚刚接触R的用户来说，这样一次性下载并安装一系列R包，可以省去调用某些包时候需要现下载、安装的等待时间，也方便测试某些代码和函数。
文本为最常见的数据存储格式，包括以.txt、.csv、.tsv等一系列扩展名结尾的文件。文本文件可以通过windows自带的记事本或者Notepad++、UltraEdit等文本编辑软件直接打开。
在R中，文本文件的读取依赖read.table()等一类命令，使用?read.table可以看到，里面可以指定很多参数，其中常用的有file，header， sep ，quote等。
file：需要读入的文件名
header: 第一列是否为变量名
sep: 变量之间的分隔符
quote: 文本被包裹的符号
比如在当前工作目录下，我们有一个制表符（Tab或\t）分割的文本文件sample.txt，第一行含有英文变量名（中文变量名可能会出错，依系统而异），然后文本没有被任何符号包裹，那么我们读入它的时候需要采用：
sample_tab_data <- read.table("data/sample_book_purschase.txt", header = TRUE,
sep = "\t")而read.csv()、read.csv2()、read.delim()、read.delim2()都是read.table()不同默认参数的变形：
library("formatR")
usage(read.csv)## read.csv(file, header = TRUE, sep = ",", quote = "\"",
## dec = ".", fill = TRUE, comment.char = "", ...)usage(read.csv2)## read.csv2(file, header = TRUE, sep = ";", quote = "\"",
## dec = ",", fill = TRUE, comment.char = "", ...)usage(read.delim)## read.delim(file, header = TRUE, sep = "\t", quote = "\"",
## dec = ".", fill = TRUE, comment.char = "", ...)usage(read.delim2)## read.delim2(file, header = TRUE, sep = "\t", quote = "\"",
## dec = ",", fill = TRUE, comment.char = "", ...)当数据不整齐的时候，R会在读入过程中报错，并给出出错的行数。当然我们也可以通过更改参数来强制读入。read.table()的常用参数定义如下：
header：是否把第一行读为变量名
sep：列与列之间的分隔符
quote：引号的格式
col.names：可以指定一个文本向量，作为变量的名字，其中所含名字的个数需与列数相一致
colClasses：可以为每一列指定格式，如果不需要读入的话置为NULL即可。
fill：是否强制每一行都有相同的列数
其他的参数如as.is、na.strings等，可参见?read.table的说明。
Excel格式除了可以采用excel里面导出文本文件或者csv文件的方式外，还可以采取ODBC方式读入(Windows和Mac下)。如果采用这种方式，需要加载RODBC这个包。
library(RODBC)
excel_channel <- odbcConnectExcel("data/sample.xls")
sample_excel <- sqlFetch(excel_channel, "sample") #参数为要导入的excel数据表的名字
odbcClose(channel)
sample_excel_2007 <- odbcConnectExcel2007("data/sample.xlsx") #对于07版excel文件当然，除了excel之外，所有基于ODBC接口的数据都可以读入，包括常见的MySQL、Access等。在Linux下，MySQL数据库建议使用另外的RMySQL包连接。
除了ODBC之外，另一种简单的方式则是调用xlsx包（在linux平台下也可运行）。该包不仅仅可以读入excel数据，还可以“将Excel读取为数据框，以及将数据框写入为Excel文件都不是问题，而更加强大的是它能处理Excel中的格式，比如合并单元格，设置列的宽度，设置字体和颜色等等”（具体请参见yixuan的介绍）。这样，就为通过R脚本重复生成excel格式的报表等铺平了道路。
如果只是读入数据，那么直接调用read.xlsx()就可以了。
library(xlsx)
file <- system.file("tests", "test_import.xlsx", package = "xlsx")
sample_sheet1 <- read.xlsx(file, 1) # 读取第一个工作表对于该包中write.xlsx()等写入excel文件函数感兴趣的读者，可以借助?write.xlsx来阅读详细的使用说明。注：该包只支持excel 07以上格式，对于.xls格式还请先用高版本的Excel等软件进行转换。
其他数据库可以通过对应的包连接 <待补充>
大多数常见的数据都可以通过foreign这个包读入：
SPSS： read.spss()
Minitab：read.mtp()
SAS：read.ssd()或者read.xport()
Stata：read.dta()
导入数据之后，可以在RStudio的右上方Workspace那里看到一个新的数据，名为sample，共有针对17个变量的99行记录。单击可以在左上角的数据浏览窗内看到对应的数据样本。
这个时候，可以用names()来查看变量名：
sample_tab_data <- read.table("data/sample_book_purschase.txt", header = TRUE,
sep = "\t")
names(sample_tab_data)## [1] "CUSTOMER" "DAY" "LOCATION" "LOCATION2" "FROM"
## [6] "NOTE" "BOOK_ID" "NOTE2" "BOOK_TYPE" "CLASS_ID"
## [11] "CLASS_NAME" "AMOUNT" "LENGTH" "LENGTH2" "PAGES"
## [16] "PAGES2" "BOOK_NAME" 而后，可以对变量进行重新命名：
比如，我们想把第二个字段重命名为RECORD_DAY：
names(sample_tab_data)[2] <- "RECORD_DAY"这里把names(sample_tab_data)返回的第二个元素重新定义为了RECORD_DAY，故而实现了变量的重命名。
或者，我们希望对导入的没有变量名的数据集进行重命名（一般这种情况下对应的默认变量名是V1、V2等），那么可以直接对整个数据集操作。
names(sample_tab_data) <- c("CUSTOMER", "RECORD_DAY", "LOCATION",
"LOCATION2", "FROM", "NOTE", "BOOK_ID", "NOTE2", "BOOK_TYPE", "CLASS_ID",
"CLASS_NAME", "AMOUNT", "LENGTH", "LENGTH2", "PAGES", "PAGES2", "BOOK_NAME")在一个data.frame中，我们可以直接用$来调用其中的一个变量，是最简单的调用列的格式。如果希望调用某些行列，则需要分别指定调用条件，比如：
sub_sample <- sample_tab_data["BOOK_ID" == 348918047, c("CUSTOMER",
"RECORD_DAY", "BOOK_ID")]那么sub_sample里面现在就得到了购买过编号为348918047这本书的所有顾客购买记录，包括顾客ID、购买日期和书籍ID。即，对于任何一个data.frame对象，都可以在中括号内，逗号之前指定行选择条件，逗号之后指定要选择的列（变量）。c()为向量生成函数。
CUSTOMER_ID <- sample_tab_data$CUSTOMER这个时候就会在workspace那里出现一个新向量CUSTOMER_ID。各个格式之间可以直接转换，比如CUSTOMER_ID <- as.data.frame(CUSTOMER_ID)那么这个时候customer_ID就成为了只有一个变量的data.frame格式了。logic, character, integer, numeric, factor
对于向量中的元素，记录的格式则可能是逻辑型、文本型、整数型、数值型、因素型等等。各个向量格式之间可以直接转换，比如对于CLASS_ID这个变量，虽然是数值记录的但数值本身没有任何意义，只是一个相互区别和识别的代码，因此可以考虑专为character或者factor格式：
sample_tab_data$CLASS_ID <- as.factor(sample_tab_data$CLASS_ID)之后R里面就会识别其为文本或者因素型数据了。值得注意的是，如果需要把一个因素型的数据重新转换为整数型，则需要经过文本型过渡：sample_tab_data$CLASS_ID <- as.character(sample_tab_data$CLASS_ID)
sample_tab_data$CLASS_ID <- as.integer(sample_tab_data$CLASS_ID)否则会被直接重编码，丢失原有的数据串信息。（注：花开两朵、各表一枝。实际上，我们也可以应用这种特性来进行重编码工作。）逻辑型
如果我们基于一些记录判断生成新的变量，比如基于如果用户购买量>0，则我们认为其在当日有购买行为，那么可以使用：
sample_tab_data$purchase <- sample_tab_data$AMOUNT > 0这样就生成了一个新的逻辑型变量purchase（取值为TRUE 或者FALSE）。逻辑型变量的一大用处就是可以直接通过相乘操作来进行多个行为之间的交集运算，比如除了是否购买之外，我们还关心购买的书籍是不是在标号为200的书店购买的，那么就可以：sample_tab_data$book_store_200 <- sample_tab_data$LOCATION == 200
sample_tab_data$purchase_bs200 <- sample_tab_data$book_store_200 *
sample_tab_data$purchase最后得到的变量purchase_bs200为TRUE则该用户是在200号书店购买的图书。类似的，我们可以统计是不是每个月都有购买行为等。
这里还一并介绍一个有用的`%in%`运算符，表示一个元素是否属于一个给定的集合，比如：sample_tab_data$book_store <- sample_tab_data$LOCATION %in% c(200,
300)表示用户在200或者300号书店进行了购买。`c()`为向量生成函数，故得到的向量含有200和300两个元素。数值型
数值型变量的操作一般就是基本的运算，R里面的话，
+：加法
-：减法
*：乘法
/：除法
%% ：取余
^n ：n次方
字符操作
字符操作最常见的就是字符串生成操作，比如我们有CUSTOMER、LOCATION、和BOOK_NAME三个变量，希望批量生成一个变量，然后发送给顾客作为反馈记录，希望的格式为“CUSTOMER顾客您好，您在编号为LOCATION的书店购买了书籍BOOK_NAME，仅供确认。”，那么我们可以使用paste()这个函数：
sample_tab_data$message <- paste(sample_tab_data$CUSTOMER, "顾客您好，您在编号为",
sample_tab_data$LOCATION, "的书店购买了书籍", sample_tab_data$BOOK_NAME,
"，仅供确认。", sep = "")这个时候就得到的相应的新变量。`paste()`函数有个参数是`sep`，用来指定各个部分之间的连接符，默认为空格，如果不需要任何额外的符号用一对双引号设置为空即可。
字符的其他操作亦包括查找、截取`substr()`等。列合并：对于两个含有完全一样变量的数据集，可以采用rbind()函数来直接将一个data.frame附加在另外一个后面。
行合并：如果两个数据集有相同的记录数（行数），那么可以采用cbind()这个函数来直接把变量附加在后面。
数据集合并：merge()函数提供了更强大的数据集合并操作命令，可以按照一个主键（即用来识别个体的变量）来合并，比如我们另有一个文件 BOOK_MAP.txt，里面记录的是重编码后的书籍ID和原编码对照表，则可以读入之后利用来merge()合并：
book_map <- read.delim("data/BOOK_MAP.txt", header = T)
sample_merged <- merge(sample_tab_data, book_map, by.x = "BOOK_ID",
by.y = "BOOK_ID", all.x = T, all.y = F)这样就有了一列新的变量，记录的是重编码之后的书籍ID。删除/提取变量
如果要删除某个变量，可以直接使用NULL值置空，即：
sample_merged$CLASS_NAME <- NULL会删除掉` CLASS_NAME`这个变量。在需要删除多个变量的时候，不如只保留几个变量，如：sample_merged <- sample_merged[, c("BOOK_ID", "CLASS_ID")]会只保留`BOOK_ID`, `CLASS_ID` 两个变量。除了这些基本的变量操作之外，还有一类可能的需求就是对于整个数据集做一个形状的转换，比如把“长数据集”转换为“宽数据集”。这样的过程类似于“揉面”，而帮我们玩转面团的利器便是reshape2这个包。
在正式介绍强大的reshape2包之前，需要先提到一个轻量级武器——reshape()函数。这个函数可以帮我们在数据的长、宽形状之间自由玩转，比如我们现在有一个用户逐月购买记录，为长格式，想把它变为宽格式：
load("data/reshape_sample.rdata") #载入样本数据集
summary(reshape_sample) #基本统计量，有CUSTOMER_ID（顾客ID）、MONTH（月份）、PURCHASE（消费额）三个变量## CUSTOMER_ID MONTH PURCHASE
## 5 : 8 201103 : 5927 Min. : 1
## 49 : 8 201104 : 5927 1st Qu.: 30
## 64 : 8 201105 : 5927 Median : 481
## 79 : 8 201106 : 5927 Mean : 1694
## 122 : 8 201107 : 5927 3rd Qu.: 1951
## 135 : 8 201108 : 5927 Max. :193547
## (Other):47368 (Other):11854 head(reshape_sample) #调用数据的前几行，显示为长格式## CUSTOMER_ID MONTH PURCHASE
## 8 5 201103 130
## 9 5 201104 79
## 10 5 201105 73
## 11 5 201106 139
## 12 5 201107 24
## 13 5 201108 4reshape_sample_wide <- reshape(reshape_sample, idvar = "CUSTOMER_ID",
timevar = "MONTH", direction = "wide") #变为宽格式
head(reshape_sample_wide) #宽格式展现## CUSTOMER_ID PURCHASE.201103 PURCHASE.201104 PURCHASE.201105
## 8 5 130 79 73
## 113 49 317 297 441
## 149 64 3211 3290 3523
## 178 79 17630 13564 2086
## 278 122 2140 4506 277
## 319 135 20097 13988 31453
## PURCHASE.201106 PURCHASE.201107 PURCHASE.201108 PURCHASE.201109
## 8 139 24 4 76
## 113 167 139 73 15
## 149 3760 1610 548 806
## 178 2923 1374 906 1359
## 278 351 1 1442 1663
## 319 84411 114838 107842 35610
## PURCHASE.201110
## 8 7
## 113 404
## 149 925
## 178 2119
## 278 328
## 319 110981summary(reshape_sample_wide) #宽格式下基本统计量## CUSTOMER_ID PURCHASE.201103 PURCHASE.201104 PURCHASE.201105
## 5 : 1 Min. : 1 Min. : 1 Min. : 1
## 49 : 1 1st Qu.: 18 1st Qu.: 24 1st Qu.: 26
## 64 : 1 Median : 406 Median : 488 Median : 544
## 79 : 1 Mean : 1573 Mean : 1646 Mean : 1769
## 122 : 1 3rd Qu.: 1792 3rd Qu.: 1922 3rd Qu.: 2114
## 135 : 1 Max. :193547 Max. :99872 Max. :100394
## (Other):5921
## PURCHASE.201106 PURCHASE.201107 PURCHASE.201108 PURCHASE.201109
## Min. : 1 Min. : 1 Min. : 1 Min. : 1
## 1st Qu.: 33 1st Qu.: 44 1st Qu.: 38 1st Qu.: 40
## Median : 579 Median : 510 Median : 492 Median : 477
## Mean : 1823 Mean : 1771 Mean : 1760 Mean : 1653
## 3rd Qu.: 2063 3rd Qu.: 2034 3rd Qu.: 1982 3rd Qu.: 1888
## Max. :112609 Max. :114838 Max. :107842 Max. :96798
##
## PURCHASE.201110
## Min. : 1
## 1st Qu.: 28
## Median : 374
## Mean : 1561
## 3rd Qu.: 1737
## Max. :110981
## 在reshape2包中，melt()函数进一步简化了这个过程。比如我们现在希望把宽数据转回长数据：
library(reshape2)
reshape_sample_long <- melt(reshape_sample_wide, id = c("CUSTOMER_ID")) #转回长数据格式
head(reshape_sample_long[order(reshape_sample_long$CUSTOMER_ID),
]) #长格式展示## CUSTOMER_ID variable value
## 1 5 PURCHASE.201103 130
## 5928 5 PURCHASE.201104 79
## 11855 5 PURCHASE.201105 73
## 17782 5 PURCHASE.201106 139
## 23709 5 PURCHASE.201107 24
## 29636 5 PURCHASE.201108 4此外，该包提供的acast()/dcast()函数可以进一步帮我们分类展现数据及其统计量，具体使用请参见函数包内帮助。
数据导出最常用的应该就是write.table函数。比如我们要输出book_map这个数据集为文本格式，那么使用：
write.table(book_map, file = "book_map_new.txt", row.names = F, col.names = T,
sep = "\t", quote = F)这个时候就会得到book_map_new.txt这个文本文件，以制表符分隔。write.table()的主要参数有：
file：指定文件名
row.names：是否输出记录行编号
col.names：是否输出变量名
sep：分隔符
quote：引号形式
append：是否附加在现有文件后面（如为FALSE则新文件覆盖原有文件）
日常分析工作中，最常用到的就是分类统计了。简单的说，就是按一个字段归类之后，统计其他字段的量。
比如，我们现在希望知道每个顾客购买的图书的总数。那么可以使用如下的代码：
library(data.table)
sample_tab_data_frame <- data.table(sample_tab_data)
sample_stats <- sample_tab_data_frame[, list(Amount = sum(AMOUNT),
Book = length(unique(BOOK_ID))), by = "CUSTOMER"]
sample_stats## CUSTOMER Amount Book
## [1,] 297017f0d704830c05774f455c5919e3 5876 76
## [2,] c6970366fd1f9f6de96e80ca77151c58 1 1
## [3,] fcd51818074819c811c9dd3f3b9eecbc 4 1
## [4,] ace38a66dcc1133cfec26d6adb870091 12 3
## [5,] 4ae62768ad286695948069b95a4ed36e 1 1
## [6,] fa57634cdf8f3a99ff894cab54084452 1 1
## [7,] f6461a358469ffa52dd9cd2e295dbc76 12 2
## [8,] 47fead3dfdd9b483b70f7ff925abc8a4 29 1
## [9,] ba6a3dde6168c6207a63b7f73559b0a7 491 5这样在新的对象sample_stats里面，就有了每个顾客购买书籍的总本书、以及不同图书的数量。简单的说，需要使用这样的分类统计的时候，就是先加载data.table这个包，然后利用data.table()函数转换一个data.frame为data.table格式，然后在原有的data.frame行、列操作基础上，增加了一个by参数，可以用来指定分类的依据。当然这里我们可以同时针对多个变量及其组合分类统计，比如
library(data.table)
sample_stat_by_day <- sample_tab_data_frame[, list(Amount = sum(AMOUNT),
Book = length(unique(BOOK_ID))), by = c("CUSTOMER", "RECORD_DAY")]
sample_stat_by_day## CUSTOMER RECORD_DAY Amount Book
## [1,] 297017f0d704830c05774f455c5919e3 201106 5876 76
## [2,] c6970366fd1f9f6de96e80ca77151c58 201106 1 1
## [3,] fcd51818074819c811c9dd3f3b9eecbc 201106 4 1
## [4,] ace38a66dcc1133cfec26d6adb870091 201106 12 3
## [5,] 4ae62768ad286695948069b95a4ed36e 201106 1 1
## [6,] fa57634cdf8f3a99ff894cab54084452 201106 1 1
## [7,] f6461a358469ffa52dd9cd2e295dbc76 201106 12 2
## [8,] 47fead3dfdd9b483b70f7ff925abc8a4 201106 29 1
## [9,] ba6a3dde6168c6207a63b7f73559b0a7 201106 491 5这样返回的就是逐日统计的每位顾客的购买数量了。值得多说的是，sum()函数代表求和，length()函数代表计数，而unique()函数则是去掉重复值。
我们常用的循环和判断有三种：
for
while
if
循环和判断是基本的逻辑操作语句。在R中，大部分常用功能都已经有现成的函数，所以极少用到循环，我们也非常不提倡在R里面写循环（效率一般很低，因为一个循环背后往往是现有函数内部的许多层循环）。有的时候，为了一些特殊的情况，知道怎么写循环还是有用的。比如，我们希望把统计好的顾客购买记录按照顾客ID写入不同的文本文件，这里就需要用到for循环或者while循环。
for (i in 1:nrow(sample_stats)) {
file_name <- paste("data/result_", sample_stats[i, ]$CUSTOMER, "_record.txt",
sep = "")
write.table(sample_stats[i, ], file = file_name, sep = "\t", row.names = F,
col.names = T, quote = F)
}这里用到了for循环。另，nrow()是用来统计一个data.frame行数的函数，同样的ncol()会返回其列数。符号:会生成一个向量，从1到nrow(sample_stat)。常用的循环还有while，同时可以使用if来进行一些判断。这里不再赘述。
上面简单的说了一下多元回归，下面则是一些我们在回归分析中常用分析的实现。
得到一个回归方程后，关心的第一件事就是系数和方程整体的显著性，分别由t检验和F检验实现。来看下面这个有关法学院的例子。
(@LAWSCH85) 在LAWSCH85.DTA这个数据集中，法学院应届生薪水的中位数由下面的方程决定：
\(log(salary)=\beta_{0}+\beta_{1}LAST+\beta_{2}GPA+\beta_{3}log(libvol)+\beta_{4}log(cost)+\beta_{5}rank+u\)
其中LSAT 是班级里LSAT成绩中位数，GPA 是班级学习成绩的中位数，libvol 是法学院图书馆藏书卷数，cost 是在法学院的年消费额，rank 是法学院的排名（正序）。 接下来我们估计这个方程:
LAW_Result <- lm(log(salary) ~ LSAT + GPA + log(libvol) + log(cost) + rank, data=LAW)
很容易得到回归结果如下：
LAW <- read.dta("data/LAWSCH85.DTA")## Warning message: cannot read factor labels from Stata 5 filesLAW_Result <- lm(log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
rank, data = LAW)
summary(LAW_Result)##
## Call:
## lm(formula = log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
## rank, data = LAW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30136 -0.08498 -0.00436 0.07794 0.28861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.343226 0.532519 15.67 <2e-16 ***
## LSAT 0.004696 0.004010 1.17 0.2437
## GPA 0.247524 0.090037 2.75 0.0068 **
## log(libvol) 0.094993 0.033254 2.86 0.0050 **
## log(cost) 0.037554 0.032106 1.17 0.2443
## rank -0.003325 0.000348 -9.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.112 on 130 degrees of freedom
## (20 observations deleted due to missingness)
## Multiple R-squared: 0.842, Adjusted R-squared: 0.836
## F-statistic: 138 on 5 and 130 DF, p-value: <2e-16
## 在回归结果中已经报告了各变量的t统计值，从而可知：rank 的估计值很显著（通过0.1%显著性水平检验）。而GPA 和(log)libvol 则通过了1%显著性水平检验。 而变量LSAT 系数估计值不显著。
考虑到新生入学的时候只有GPA 和LSAT 两个变量可以观测，所以接下来我们进行变量GPA 和LSAT 的联合检验即F检验。该检验属于线性假设检验，在RCommander下可以在Models -> Hypothesis Tests ->Linear hypothesis里面通过图形化界面完成。
其中设置为“2”行，而后分别在LAST和GPA输入1，保持右侧为0。相当于检验假设
\(H_{0}:\;\beta_{1}=\beta_{2}=0\)
可得到输出结果如下：
library("car")
.Hypothesis <- matrix(c(0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0), 2, 6,
byrow = TRUE)
.RHS <- c(0, 0)
linearHypothesis(LAW_Result, .Hypothesis, rhs = .RHS)## Linear hypothesis test
##
## Hypothesis:
## LSAT = 0
## GPA = 0
##
## Model 1: restricted model
## Model 2: log(salary) ~ LSAT + GPA + log(libvol) + log(cost) + rank
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 132 1.89
## 2 130 1.64 2 0.252 9.95 9.5e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 remove(.Hypothesis, .RHS)从上面结果可知，F=9.9517 且通过了0.1%显著性水平检验，即变量LSAT 和GPA 联合显著。 此外还可以分别加入clsize （班级容量）和faculty （教师规模）来进行回归，代码如下： [未完成]
由回归结果，这两个变量的系数均不显著，且方程F统计量有所下降（调整后的\(R^{2}\) 变动不大），故不用加入到方程中关于模型变量选择的问题，将在后面另行论述，这里只做简单的分析。
回归分析里还常用的一项分析就是得到某一显著性水平下的置信区间。[未完成]
比较简单的例子就是已经含有分组变量的数据，比如变量gender 有两个值male, female, 那么我们只需把它们变成factor形式就可以了。如我们可以把法学院例子中的变量north进行变化 [footnote-factor][]： [footnote-factor]: 其实这里不用这么麻烦也可以，因为north变量本身的赋值只有0和1，可以直接进行回归。在这里只是用这个例子来说明factor()函数的调用形式而已。 “footnote-factor”
LAW$north_true <- factor(LAW$north, labels = c("others", "north"))可以这样做的原因是：R在调用lm()函数的时候会自动把factor类型的变量作为虚拟变量进行回归。
依旧采用法学院的例子。很明显在上面的分析中我们把各个学院排名直接当作一个可以“测距”的变量其实它只是一个排序而已，并不代表实质的差距。来使用了，这可能会引起一些争议。因此，我们可以采取另一种模式，即引入虚拟变量，把学校分为六组：top 10，11-25名，26-40名，41-60名，61-100名，100名开外。引入五个虚拟变量\(top10\) , \(r11\_25\) , \(r26\_40\) , \(r41\_60\) , \(r61\_100\) 。 我想你不会手动的把所有的学校都赋一个虚拟变量值吧？在R里面，我们需要先通过recode()来依照分组创建一个新的factor形式的变量\(rank\_f\)，然后再进行回归。这样我们就不需要在原来的数据库里面新增加五个变量并赋逻辑值了。
LAW$rank_f <- recode(LAW$rank, "1:10=\"top10\"; 11:25=\"r11_25\"; 26:40=\"r26_40\"; 41:60=\"r41_60\"; 61:100=\"r61_100\"; else=\"r101_\"; ",
as.factor.result = TRUE)
LAW_Result2 <- lm(log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
rank_f, data = LAW)
summary(LAW_Result2)##
## Call:
## lm(formula = log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
## rank_f, data = LAW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29489 -0.03969 -0.00168 0.04389 0.27750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.165295 0.411424 22.28 < 2e-16 ***
## LSAT 0.005691 0.003063 1.86 0.066 .
## GPA 0.013726 0.074192 0.19 0.854
## log(libvol) 0.036362 0.026017 1.40 0.165
## log(cost) 0.000841 0.025136 0.03 0.973
## rank_fr11_25 0.593543 0.039440 15.05 < 2e-16 ***
## rank_fr26_40 0.375076 0.034081 11.01 < 2e-16 ***
## rank_fr41_60 0.262819 0.027962 9.40 3.2e-16 ***
## rank_fr61_100 0.131595 0.021042 6.25 5.7e-09 ***
## rank_ftop10 0.699566 0.053492 13.08 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0856 on 126 degrees of freedom
## (20 observations deleted due to missingness)
## Multiple R-squared: 0.911, Adjusted R-squared: 0.905
## F-statistic: 143 on 9 and 126 DF, p-value: <2e-16
## 此外，我们可以通过RCommander里面，在Data -> Recode Variables的方框里逐行输入。
Recode Varibles图形化界面
我不知道这样叫是不是足够确切，在微观计量里面我们会经常用到两个变量相乘的回归项，比如\(female\cdot single\)，即单身女士。相比而言这样的虚拟变量并不需要特别的处理，在回归方程里面直接写成相乘的形式即可。注意此时不需要再写female 和single 变量，lm()会默认加入这两个变量。例如，在法学院的例子中，我们可以对top10 和west 两个变量进行相乘回归（这里top10 变量来源于数据库本身自带的）。
LAW_Result3 <- lm(log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
top10 * west, data = LAW)
summary(LAW_Result3)##
## Call:
## lm(formula = log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
## top10 * west, data = LAW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3723 -0.0997 -0.0219 0.0943 0.4159
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.22585 0.55409 9.43 2.3e-16 ***
## LSAT 0.00778 0.00521 1.49 0.1383
## GPA 0.50411 0.11417 4.42 2.1e-05 ***
## log(libvol) 0.21985 0.04086 5.38 3.4e-07 ***
## log(cost) 0.12096 0.04036 3.00 0.0033 **
## top10 0.10506 0.06688 1.57 0.1187
## west 0.01635 0.03052 0.54 0.5931
## top10:west 0.01128 0.11961 0.09 0.9250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.146 on 128 degrees of freedom
## (20 observations deleted due to missingness)
## Multiple R-squared: 0.737, Adjusted R-squared: 0.723
## F-statistic: 51.4 on 7 and 128 DF, p-value: <2e-16
## 当然也可以使用factor变量和其他虚拟变量相乘进行回归，反馈的结果中包含所有的相乘项。
LAW_Result4 <- lm(log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
rank_f * west, data = LAW)
LAW_Result4##
## Call:
## lm(formula = log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
## rank_f * west, data = LAW)
##
## Coefficients:
## (Intercept) LSAT GPA
## 9.183816 0.005732 0.011030
## log(libvol) log(cost) rank_fr11_25
## 0.035169 -0.000485 0.591775
## rank_fr26_40 rank_fr41_60 rank_fr61_100
## 0.378440 0.260765 0.137090
## rank_ftop10 west rank_fr11_25:west
## 0.713412 0.012817 0.008311
## rank_fr26_40:west rank_fr41_60:west rank_fr61_100:west
## -0.011688 0.005488 -0.020113
## rank_ftop10:west
## -0.055886
## 当R处理factor形式的数据的时候，默认以数据中的第一个层次 (level) 作为参照组。比如上例中，我们想把top10 这一组作为参照组，那么则需要使用relevel()命令。
attach(LAW)
rank_f2 <- relevel(rank_f, ref = "top10")
LAW_Result5 <- lm(log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
rank_f2, data = LAW)
summary(LAW_Result5)##
## Call:
## lm(formula = log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
## rank_f2, data = LAW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29489 -0.03969 -0.00168 0.04389 0.27750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.864861 0.449840 21.93 < 2e-16 ***
## LSAT 0.005691 0.003063 1.86 0.0655 .
## GPA 0.013726 0.074192 0.19 0.8535
## log(libvol) 0.036362 0.026017 1.40 0.1647
## log(cost) 0.000841 0.025136 0.03 0.9734
## rank_f2r101_ -0.699566 0.053492 -13.08 < 2e-16 ***
## rank_f2r11_25 -0.106023 0.038716 -2.74 0.0071 **
## rank_f2r26_40 -0.324490 0.044339 -7.32 2.6e-11 ***
## rank_f2r41_60 -0.436747 0.045913 -9.51 < 2e-16 ***
## rank_f2r61_100 -0.567971 0.047180 -12.04 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0856 on 126 degrees of freedom
## (20 observations deleted due to missingness)
## Multiple R-squared: 0.911, Adjusted R-squared: 0.905
## F-statistic: 143 on 9 and 126 DF, p-value: <2e-16
## 有件很没办法的事儿，那就是要想让OLS回归出来的结果最佳，必须要符合那五条经典假设（尤其是小样本下）。但是事实中的数据那里会那么完美呢？首当其冲的问题就是异方差。 异方差检验方法很多，这里给出两种常用的：BP检验 (Breusch-Pagan Test)、怀特检验 (White test for heteroskedasticity)。
下面给出一个关于房价的例子。
(@Hprice1) 在Hprice1.dta这个数据集中，有price （房价，按套计算）、lotsize （地皮面积要知道人家住的都是小别墅啊，自然要先有地、后建房，卖房子也都是一套一套的卖。）、sqrft （房屋面积）、bdrms （卧室数量）。接下来我们估计这个方程：
\(\hat{price}=\hat{\beta_{0}}+\hat{\beta_{1}}lotsize+\hat{\beta_{2}}sqrft+\hat{\beta_{3}}bdrms\)
使用OLS估计结果如下：
HPrice <- read.dta("data/hprice1.dta")
Hprice_Result <- lm(price ~ bdrms + lotsize + sqrft, data = HPrice)
Hprice_Result##
## Call:
## lm(formula = price ~ bdrms + lotsize + sqrft, data = HPrice)
##
## Coefficients:
## (Intercept) bdrms lotsize sqrft
## -21.77031 13.85252 0.00207 0.12278
## 下面我们进行BP检验来测定是否有异方差。进行BP检验需要加载包lmtest，而后者需要加载包zoo。调用BP检验最简单的方法就是直接写回归结果变量。更多参数关于该命令各种参数设置请使用?bptest来查看帮助文档。可以通过RCommander里面的图形化界面设定，位于Models > Numberical diagnostics > Breusch-Pagan Test for heteroskedasticity。
library("lmtest")
bptest(Hprice_Result)##
## studentized Breusch-Pagan test
##
## data: Hprice_Result
## BP = 14.09, df = 3, p-value = 0.002782
## 由结果来看，存在异方差。由于对数形式是消除异方差（尤其针对价格数据）的常用方法，因此我们再对对数形式进行回归:
\(\hat{log(price)}=\hat{\beta_{0}}+\hat{\beta_{1}}log(lotsize)+\hat{\beta_{2}}log(sqrft)+\hat{\beta_{3}}bdrms\)
得到回归结果如下：
Hprice_Result2 <- lm(log(price) ~ bdrms + log(lotsize) + log(sqrft),
data = HPrice)
Hprice_Result2##
## Call:
## lm(formula = log(price) ~ bdrms + log(lotsize) + log(sqrft),
## data = HPrice)
##
## Coefficients:
## (Intercept) bdrms log(lotsize) log(sqrft)
## -1.297 0.037 0.168 0.700
## 此时再进行异方差检验。
bptest(Hprice_Result2)##
## studentized Breusch-Pagan test
##
## data: Hprice_Result2
## BP = 4.223, df = 3, p-value = 0.2383
## 因此可以接受原假设，已经不存在异方差。 此外还有一个bptest()非常接近的ncv.test()，需加载car*包。我们进行一下对比。
library(car)
ncvTest(Hprice_Result2)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 3.525 Df = 1 p = 0.06044 bptest(Hprice_Result2, studentize = FALSE, varformula = ~fitted.values(Hprice_Result2),
data = HPrice)##
## Breusch-Pagan test
##
## data: Hprice_Result2
## BP = 3.525, df = 1, p-value = 0.06044
## 可以看到ncvTest()相当于把bptest()里面参数studentize设为FLASE，且进行固定值的检验。事实上，bptest()默认该项为TRUE，在TRUE时会采用Koenker (1981)的学生氏检验算法，也是目前最广为接受的算法。
还是上面这个例子，我们改用怀特检验。其实，我们可以把怀特检验看作广义上的BP检验的一种特殊形式，因此可以通过在bptest()里面赋予更多的参数来实现，即加入各个变量平方和交叉相乘的项。
比如回归为fm <- lm(y ~ x + z, data = foo)，那么则应写成bptest(fm, ~ x * z + I(x^2) + I(z^2), data = foo) 。其余以此类推，不再举例。
另，也可以通过sandwich这个专门对付“三明治”的包来实现怀特检验。
当存在异方差的时候，我们可以采用稳健的标准差来替代原有的异方差矩阵。在R中，可以利用sandwich包的vcovHC()函数实现。依旧采用上面房价的例子。
library("sandwich")
vcovHC(Hprice_Result)## (Intercept) bdrms lotsize sqrft
## (Intercept) 1683.68200 -221.25362 -0.0616765 -0.2399245
## bdrms -221.25362 133.67499 0.0475633 -0.3056587
## lotsize -0.06168 0.04756 0.0000511 -0.0002524
## sqrft -0.23992 -0.30566 -0.0002524 0.0016591vcovHC()函数可以选择各种形式，在后面附加type = 即可。可选形式有“HC3”, “const”, “HC”, “HC0”, “HC1”, “HC2”, “HC4”，分别对应不同的异方差假设形式。如该函数的文档中所述，const对应 \(\sigma^{2}(X'X)^{-1}\) ，HC（或HC0）对应 \((X'X)^{-1}X'\Omega X(X'X)^{-1}\) ，更多解释请参照?vcovHC。 之后我们就可以采用该异方差矩阵进行参数检验，如t检验。这里调用lmtest包内另一个非常有用的函数coeftest()。
coeftest(Hprice_Result, vcov = vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -21.77031 41.03269 -0.53 0.5971
## bdrms 13.85252 11.56179 1.20 0.2342
## lotsize 0.00207 0.00715 0.29 0.7731
## sqrft 0.12278 0.04073 3.01 0.0034 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 此外，对于含有自回归异方差的情形，可以采用vcovHAC()函数。更多相关用法后面章节中详述。这两个函数均可以附加到参数检验coeftest()或者waldtest()一起使用。
有些情况下，我们可以写出加权的形式，比如扰动项服从\(Var(u_{i}|inc)=\sigma^{2}inc\) ，那么可以直接在lm()函数里附加一项weight来实现。
(@smoke) 下面是一个烟草需求的例子。在SMOKE.rda中有如下几个变量：每天吸烟的数量 (cigs)、年收入 (income)、该州烟的价格 (cigpric)、受访者年龄 (age)、受教育程度 (educ)、该州有无饭店内吸烟禁令 (restaurn)。而后我们需要研究决定烟草需求的因素，即cigs为被解释变量，其他为解释变量。
load("data/SMOKE.rda")
SMOKE_OLS <- lm(cigs ~ log(income) + log(cigpric) + educ + age +
I(age^2) + restaurn, data = SMOKE)lm()所得结果，residuals()存储的是残差项。SMOKE_auxreg <- lm(log(residuals(SMOKE_OLS)^2) ~ log(income) + log(cigpric) +
educ + age + I(age^2) + restaurn, data = SMOKE)SMOKE_FGLS <- lm(cigs ~ log(income) + log(cigpric) + educ + age +
I(age^2) + restaurn, data = SMOKE, weights = 1/exp(fitted(SMOKE_auxreg)))
summary(SMOKE_FGLS)##
## Call:
## lm(formula = cigs ~ log(income) + log(cigpric) + educ + age +
## I(age^2) + restaurn, data = SMOKE, weights = 1/exp(fitted(SMOKE_auxreg)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.904 -0.953 -0.810 0.842 9.856
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.635463 17.803139 0.32 0.75167
## log(income) 1.295239 0.437012 2.96 0.00313 **
## log(cigpric) -2.940312 4.460145 -0.66 0.50993
## educ -0.463446 0.120159 -3.86 0.00012 ***
## age 0.481948 0.096808 4.98 7.9e-07 ***
## I(age^2) -0.005627 0.000939 -5.99 3.2e-09 ***
## restaurn -3.461064 0.795505 -4.35 1.5e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.58 on 800 degrees of freedom
## Multiple R-squared: 0.113, Adjusted R-squared: 0.107
## F-statistic: 17.1 on 6 and 800 DF, p-value: <2e-16
## gamma2i <- coef(SMOKE_auxreg)[2]
gamma2 <- 0
while (abs((gamma2i - gamma2)/gamma2) > 1e-07) {
gamma2 <- gamma2i
SMOKE_FGLSi <- lm(cigs ~ log(income) + log(cigpric) + educ + age + I(age^2) +
restaurn, data = SMOKE, weights = 1/exp(fitted(SMOKE_auxreg)))
SMOKE_auxreg <- lm(log(residuals(SMOKE_FGLSi)^2) ~ log(income) + log(cigpric) +
educ + age + I(age^2) + restaurn, data = SMOKE)
gamma2i <- coef(SMOKE_auxreg)[2]
}
SMOKE_FGLS2 <- lm(cigs ~ log(income) + log(cigpric) + educ + age +
I(age^2) + restaurn, data = SMOKE, weights = 1/exp(fitted(SMOKE_auxreg)))
summary(SMOKE_FGLS2)##
## Call:
## lm(formula = cigs ~ log(income) + log(cigpric) + educ + age +
## I(age^2) + restaurn, data = SMOKE, weights = 1/exp(fitted(SMOKE_auxreg)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.487 -0.977 -0.864 0.835 7.972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.40314 19.54430 0.07 0.94278
## log(income) 1.12124 0.42827 2.62 0.00901 **
## log(cigpric) -1.92067 4.87588 -0.39 0.69375
## educ -0.46133 0.12895 -3.58 0.00037 ***
## age 0.57558 0.10313 5.58 3.3e-08 ***
## I(age^2) -0.00661 0.00102 -6.46 1.8e-10 ***
## restaurn -3.31660 0.79246 -4.19 3.2e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.57 on 800 degrees of freedom
## Multiple R-squared: 0.107, Adjusted R-squared: 0.0999
## F-statistic: 15.9 on 6 and 800 DF, p-value: <2e-16
## 在其中我们使用当log(income)的系数估计值收敛\((<10^{-7})\)当作循环的条件。
通常被解释变量并不一定服从正态分布，比如为0-1的离散情况，因而需要进一步借助Probit, Logit等模型来进行估计。在R中，在采取广义线性估计法（Generalized Linear Models, GLM）来估计的时候，我们可以调用glm包。
(@MORZ) 这里我们看一个关于已婚妇女劳动参与率的例子（MROZ.dra）。当然，一个人工不工作是一个二值变量 (inlf)，我们设1为工作，0表示不工作。这里我们不妨认为其劳动参与行为主要取决于其他的收入来源——丈夫的工资 (nwifeinc)，受教育年限 (educ)，工作经验 (exper)，年龄 (age)，小于六岁的孩子数 (kidslt6)，六到十八岁的孩子数 (kidsge6)。 首先我们使用OLS来估计线性概率模型（Linear Probability Model，LPM）。
load("data/MROZ.rda")
MROZ_LPM <- lm(inlf ~ nwifeinc + educ + exper + I(exper^2) + age +
kidslt6 + kidsge6, data = MROZ)
MROZ_LPM##
## Call:
## lm(formula = inlf ~ nwifeinc + educ + exper + I(exper^2) + age +
## kidslt6 + kidsge6, data = MROZ)
##
## Coefficients:
## (Intercept) nwifeinc educ exper I(exper^2)
## 0.585519 -0.003405 0.037995 0.039492 -0.000596
## age kidslt6 kidsge6
## -0.016091 -0.261810 0.013012
## 对于非线性的估计最常用的方法就是最大似然估计法（MLE）。在stats4包中有对应的函数mle()可进行相应的估计。但是由于在计量中单独用到最大似然估计的时候很少，大多数情况下都是用于估计特定模型（如GLM），有着特定的函数，所以在这里不再作特别介绍。
当然使用OLS来估计概率模型并不理想（线性模型需假设解释变量具有不变的边际效应），因此我们再分别使用非线性的Probit和Logit模型来估计上例。在使用glm()的时候，只需要指定回归的类型(family)，其他的用法和lm()类似。此外可以使用logLik()来获取Log-Likelihood统计量。
MROZ_Probit <- glm(inlf ~ nwifeinc + educ + exper + I(exper^2) +
age + kidsge6 + kidslt6, family = binomial(probit), data = MROZ)
MROZ_Probit##
## Call: glm(formula = inlf ~ nwifeinc + educ + exper + I(exper^2) + age +
## kidsge6 + kidslt6, family = binomial(probit), data = MROZ)
##
## Coefficients:
## (Intercept) nwifeinc educ exper I(exper^2)
## 0.27007 -0.01202 0.13090 0.12335 -0.00189
## age kidsge6 kidslt6
## -0.05285 0.03601 -0.86832
##
## Degrees of Freedom: 752 Total (i.e. Null); 745 Residual
## Null Deviance: 1030
## Residual Deviance: 803 AIC: 819 logLik(MROZ_Probit)## 'log Lik.' -401.3 (df=8)MROZ_Logit <- glm(inlf ~ nwifeinc + educ + exper + I(exper^2) + age +
kidsge6 + kidslt6, family = binomial(probit), data = MROZ)
MROZ_Logit##
## Call: glm(formula = inlf ~ nwifeinc + educ + exper + I(exper^2) + age +
## kidsge6 + kidslt6, family = binomial(probit), data = MROZ)
##
## Coefficients:
## (Intercept) nwifeinc educ exper I(exper^2)
## 0.27007 -0.01202 0.13090 0.12335 -0.00189
## age kidsge6 kidslt6
## -0.05285 0.03601 -0.86832
##
## Degrees of Freedom: 752 Total (i.e. Null); 745 Residual
## Null Deviance: 1030
## Residual Deviance: 803 AIC: 819 logLik(MROZ_Logit)## 'log Lik.' -401.3 (df=8)依旧是上例，我们观察到每年的工作时间变化很大：对于不工作的来说，该值为0。此时如果画散点图那么必有近一半的点（325人）聚集在数轴上。因此，年工作时间 (hours)呈现很强烈的“边角解 (coner solution)”特性，对于这种情况我们可以采取Tobit模型。Tobit模型也是被审查的回归（Censored Regression）的一个特例。 其实Tobit很类似于生存分析里面的情况，因此可以调用survival包的survreg()函数。这里我们有个更简单的办法，AER包的作者进行了一个简单的转换，在包内自带了一个tobit()函数，调用更方便。这里我们将结果与OLS回归的进行对比。
library("AER")
MROZ_Tobit <- tobit(hours ~ nwifeinc + educ + exper + I(exper^2) +
age + kidslt6 + kidsge6, data = MROZ)
MROZ_Tobit##
## Call:
## tobit(formula = hours ~ nwifeinc + educ + exper + I(exper^2) +
## age + kidslt6 + kidsge6, data = MROZ)
##
## Coefficients:
## (Intercept) nwifeinc educ exper I(exper^2)
## 965.31 -8.81 80.65 131.56 -1.86
## age kidslt6 kidsge6
## -54.41 -894.02 -16.22
##
## Scale: 1122
## MROZ_OLS <- lm(hours ~ nwifeinc + educ + exper + I(exper^2) + age +
kidslt6 + kidsge6, data = MROZ)
MROZ_OLS##
## Call:
## lm(formula = hours ~ nwifeinc + educ + exper + I(exper^2) + age +
## kidslt6 + kidsge6, data = MROZ)
##
## Coefficients:
## (Intercept) nwifeinc educ exper I(exper^2)
## 1330.48 -3.45 28.76 65.67 -0.70
## age kidslt6 kidsge6
## -30.51 -442.09 -32.78
## 对于有序的Logit/Probit模型，可以调用MASS包中的polr()函数。比如研究数据集BankWages里面男性的有无工作（job）和教育程度（education）、少数民族（minority）之间的关系。
library("MASS")
data(BankWages)
bank_polr <- polr(job ~ education + minority, data = BankWages, subset = gender ==
"male", Hess = TRUE)
coeftest(bank_polr)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## education 0.8700 0.0931 9.35 < 2e-16 ***
## minorityyes -1.0564 0.4120 -2.56 0.01 *
## custodial|admin 7.9514 1.0769 7.38 1.5e-13 ***
## admin|manage 14.1721 1.4744 9.61 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 对于呈现\(\{0,1,2,...\}\)这样特性的数据，我们采取计数模型。而计数模型的实质是利用泊松回归。泊松回归也是上面所说的GLM的一种。下面我们给出一个关于犯罪率的例子。
(@CRIME1) 在CRIME1.rda这个数据集中，我们把每个人被逮捕的次数 (narr86)作为被解释变量。对于其中大多数人来说，该变量值为0；剩下的也大都在\(1\sim5\) 之间，因此泊松回归会比较合适。 一般说来我们直接用泊松回归就可以，在不假设整体都服从泊松分布的时候，我们采用的实质是准最大似然估计（Quasi-maximum likelihood estimation, QMLE）。
load("data/CRIME1.rda")
CRIME1_poisson <- glm(narr86 ~ pcnv + avgsen + tottime + ptime86 +
qemp86 + inc86 + black + hispan + born60, family = poisson(log), data = CRIME1)
CRIME1_poisson##
## Call: glm(formula = narr86 ~ pcnv + avgsen + tottime + ptime86 + qemp86 +
## inc86 + black + hispan + born60, family = poisson(log), data = CRIME1)
##
## Coefficients:
## (Intercept) pcnv avgsen tottime ptime86
## -0.59959 -0.40157 -0.02377 0.02449 -0.09856
## qemp86 inc86 black hispan born60
## -0.03802 -0.00808 0.66084 0.49981 -0.05103
##
## Degrees of Freedom: 2724 Total (i.e. Null); 2715 Residual
## Null Deviance: 3210
## Residual Deviance: 2820 AIC: 4520 CRIME1_quasipoisson <- glm(narr86 ~ pcnv + avgsen + tottime + ptime86 +
qemp86 + inc86 + black + hispan + born60, family = quasipoisson(log), data = CRIME1)
CRIME1_quasipoisson##
## Call: glm(formula = narr86 ~ pcnv + avgsen + tottime + ptime86 + qemp86 +
## inc86 + black + hispan + born60, family = quasipoisson(log),
## data = CRIME1)
##
## Coefficients:
## (Intercept) pcnv avgsen tottime ptime86
## -0.59959 -0.40157 -0.02377 0.02449 -0.09856
## qemp86 inc86 black hispan born60
## -0.03802 -0.00808 0.66084 0.49981 -0.05103
##
## Degrees of Freedom: 2724 Total (i.e. Null); 2715 Residual
## Null Deviance: 3210
## Residual Deviance: 2820 AIC: NA 当数据呈现过度离散特性的时候，就违背了泊松回归的假设（方差等于期望），所以在进行泊松回归前需要进行一下离散度检验。AER包中提供了一个函数dispersiontest()可以用来进行该检验。
(@RecreationDemand) 在RecreationDemand这个数据集中，涵盖了1980年德克萨斯州划船度假村的相关数据。我们希望用调查所得的（主观排名）设施质量（变量quality）、被调查者是否参加了滑水项目（变量ski）、家庭收入（income）、被调查者是否为游览该湖付费（userfee）、和三个表示机会成本的变量（costC , costS , costH）来解释trips 变量。
首先我们需要进行泊松回归。
library(AER)
data("RecreationDemand")
rd_pois <- glm(trips ~ ., data = RecreationDemand, family = poisson)而后对回归模型rd_pois进行离散度检验。在调用dispersiontest()函数之时，可以指定参数trafo 的值，从而只对\(\alpha\)进行估计。具体解释请参见函数说明。
dispersiontest(rd_pois)##
## Overdispersion test
##
## data: rd_pois
## z = 2.412, p-value = 0.007941
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 6.566
## dispersiontest(rd_pois, trafo = 2)##
## Overdispersion test
##
## data: rd_pois
## z = 2.938, p-value = 0.001651
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 1.316
## 在数据呈现过度离散的情况下，我们可以采取负二项回归模型。MASS包中提供了相关函数negative.binomial()和glm.nb()。前者用于参数\(\theta\)已知，后者用于该参数未知。依旧采用上面的例子。这里参数\(\theta\) 未知，所以调用glm.nb()，而后进行参数检验。
library("MASS")
rd_nb <- glm.nb(trips ~ ., data = RecreationDemand)
coeftest(rd_nb)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.12194 0.21430 -5.24 1.6e-07 ***
## quality 0.72200 0.04012 18.00 < 2e-16 ***
## skiyes 0.61214 0.15030 4.07 4.6e-05 ***
## income -0.02606 0.04245 -0.61 0.539
## userfeeyes 0.66917 0.35302 1.90 0.058 .
## costC 0.04801 0.00918 5.23 1.7e-07 ***
## costS -0.09269 0.00665 -13.93 < 2e-16 ***
## costH 0.03884 0.00775 5.01 5.4e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 当然，我们还可以再调用logLik()来获取Log-Likelihood统计量。
logLik(rd_nb)## 'log Lik.' -825.6 (df=9)对于计数模型来说，当数据集中的0观测值众多，所占比例超出泊松回归允许的范围时，回归结果就不是那么令人满意。此时需采用ZIP（Zero-inflated Poisson Model）模型。在R中，pscl包提供了一个函数zeroinfl()，用来实现ZIP泊松回归。比如例2.6中，大多数人(\(>70\%\))并未被逮捕过。所以我们不妨采用ZIP模型（负二项回归形式，ZINB）。
library("pscl")
CRIME1_zip <- zeroinfl(narr86 ~ pcnv + avgsen + tottime + ptime86 +
qemp86 + inc86 + black + hispan + born60 | inc86, data = CRIME1, dist = "negbin")
summary(CRIME1_zip)##
## Call:
## zeroinfl(formula = narr86 ~ pcnv + avgsen + tottime + ptime86 +
## qemp86 + inc86 + black + hispan + born60 | inc86, data = CRIME1,
## dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.815 -0.543 -0.428 0.165 13.307
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.28061 0.09839 -2.85 0.00434 **
## pcnv -0.47205 0.10293 -4.59 4.5e-06 ***
## avgsen -0.01989 0.02647 -0.75 0.45240
## tottime 0.02255 0.01968 1.15 0.25188
## ptime86 -0.09613 0.02710 -3.55 0.00039 ***
## qemp86 -0.15531 0.03995 -3.89 0.00010 ***
## inc86 -0.00685 0.00112 -6.11 9.8e-10 ***
## black 0.63444 0.09173 6.92 4.6e-12 ***
## hispan 0.49618 0.08845 5.61 2.0e-08 ***
## born60 -0.05239 0.07702 -0.68 0.49636
## Log(theta) 0.40448 0.15207 2.66 0.00782 **
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.769 0.254 -3.03 0.0025 **
## inc86 -30.038 96.621 -0.31 0.7559
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 1.499
## Number of iterations in BFGS optimization: 87
## Log-likelihood: -2.15e+03 on 13 Df当样本不能保证随机性的时候，即面对样本自选择问题（最常见的为偶然断尾）之时，我们一般采用Heckit模型。回到那个已婚妇女参加工作的例子（2.5）
library(sampleSelection)
MROZ_Heckit <- heckit(inlf ~ nwifeinc + educ + exper + I(exper^2) +
age + kidslt6 + kidsge6, log(wage) ~ educ + exper + I(exper^2), data = MROZ,
method = "2step")
summary(MROZ_Heckit)## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 753 observations (325 censored and 428 observed)
## 15 free parameters (df = 739)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.27008 0.50859 0.53 0.5956
## nwifeinc -0.01202 0.00484 -2.48 0.0132 *
## educ 0.13090 0.02525 5.18 2.8e-07 ***
## exper 0.12335 0.01872 6.59 8.3e-11 ***
## I(exper^2) -0.00189 0.00060 -3.15 0.0017 **
## age -0.05285 0.00848 -6.23 7.6e-10 ***
## kidslt6 -0.86833 0.11852 -7.33 6.2e-13 ***
## kidsge6 0.03601 0.04348 0.83 0.4079
## Outcome equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.578103 0.305006 -1.90 0.0584 .
## educ 0.109066 0.015523 7.03 4.8e-12 ***
## exper 0.043887 0.016261 2.70 0.0071 **
## I(exper^2) -0.000859 0.000439 -1.96 0.0507 .
## Multiple R-Squared:0.157, Adjusted R-Squared:0.149
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## invMillsRatio 0.0323 0.1336 0.24 0.81
## sigma 0.6636 NA NA NA
## rho 0.0486 NA NA NA
## --------------------------------------------Heckit模型实际上就是先进行一个probit回归，而后再利用前面的得到的估计值加上解释变量对被解释变量回归。
对于应用在结构单一的方程上的两阶段最小二乘法，我们只需要调用sem包中的tsls()函数。
回到例2.5，我们想研究教育的回报率。这个时候需要估计工资log(wage)和受教育程度（educ）、经验（exper）之间的关系。但是我们怀疑受教育程度educ 是内生变量，但他父母的受教育程度\(motheduc,\;\; fatheduc\) 则为外生变量，所以可以作为educ 的工具变量。
library(sem)
MROZ_2SLS <- tsls(lwage ~ educ + exper + I(exper^2), ~exper + I(exper^2) +
motheduc + fatheduc, data = MROZ)
summary(MROZ_2SLS)##
## 2SLS Estimates
##
## Model Formula: lwage ~ educ + exper + I(exper^2)
## <environment: 0x10458a60>
##
## Instruments: ~exper + I(exper^2) + motheduc + fatheduc
## <environment: 0x10458a60>
##
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.1000 -0.3200 0.0551 0.0000 0.3690 2.3500
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.048100 0.400328 0.12 0.9044
## educ 0.061397 0.031437 1.95 0.0515 .
## exper 0.044170 0.013432 3.29 0.0011 **
## I(exper^2) -0.000899 0.000402 -2.24 0.0257 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6747 on 424 degrees of freedom
## 由于经济系统内的各变量往往是相互联系的，所以在出现联立方程的情况下OLS估计会因为内生性的问题而产生偏差，此时需要借助两阶段最小二乘法（2SLS）来估计方程（组）。
对于联立方程模型，终极的解决方案莫过于systemfit包，它提供了似不相关回归及其变形、两阶段最小二乘法和工具变量等多种估计方法。这里我们看一个似不相关回归法来估计联立方程模型的例子。
(@Grunfeld) 我们现在来看投资和企业价值及资本存量之间的关系。在数据集Grunfeld之中，invest 代表总投资值，value 代表企业的价值，capital 代表企业的资本存量，以上数据都是剔除了通胀因素之后的真实值。 为了简单起见，我们只考虑两个厂商Chrysler和IBM，因此首先要使用subset()函数来从原数据集中找出一个子集。厂商之间可能是有联系的，所以需要分别设定变量，这里我们就可以用到factor类型的数据。这里最方便的就是利用面板数据的类型，调用plm包里面的plm.data()函数。
data(Grunfeld, package = "AER")
library("systemfit")
library("plm")
gr2 <- subset(Grunfeld, firm %in% c("Chrysler", "IBM"))
pgr2 <- plm.data(gr2, c("firm", "year"))
gr_sur <- systemfit(invest ~ value + capital, method = "SUR", data = pgr2)
summary(gr_sur, residCov = FALSE, equations = FALSE)##
## systemfit results
## method: SUR
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 40 34 4114 11022 0.929 0.927
##
## N DF SSR MSE RMSE R2 Adj R2
## Chrysler 20 17 3002 176.6 13.29 0.913 0.903
## IBM 20 17 1112 65.4 8.09 0.952 0.946
##
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Chrysler_(Intercept) -5.7031 13.2774 -0.43 0.67293
## Chrysler_value 0.0780 0.0196 3.98 0.00096 ***
## Chrysler_capital 0.3115 0.0287 10.85 4.6e-09 ***
## IBM_(Intercept) -8.0908 4.5216 -1.79 0.09139 .
## IBM_value 0.1272 0.0306 4.16 0.00066 ***
## IBM_capital 0.0966 0.0983 0.98 0.33951
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 从结果中我们可以得到似不相关回归的估计值。这里也可以调整一下summary()的参数来看更详细的回归结果。 此外，也可以加参数2SLS来调用systemfit()来进行两阶段最小二乘估计，具体方法不再赘述，请参见函数说明。
其实代理变量并没有多少要特别说明的地方，只需要选择合适的变量代替不能观测的变量进行相应的回归分析就好。但有一种常用的而且很特别的代理变量，就是被解释变量的滞后项。这里简单的以一个例子说明。
(@CRIME) 这个例子是有关城市中犯罪率的。在CRIME.rda这个数据集中，含有1987年46个城市的犯罪率数据，可见是个标准的截面数据分析。只是除了犯罪率 (crmrte)，可用的变量只有失业率 (unem) 和人均用于维护法律权益的支出 (lawexpc) 。这样就有一些难以观测的变量导致回归的结果不甚理想。因此，我们可以采取1982年各个城市的犯罪率做一个代理变量，来但因那些难以观测的变量。方程如下：
\(log(crmrte_{87})=\beta_{0}+\beta_{1}unem+\beta_{2}lawexpc+\beta_{3}log(crmrte_{82})\)
但是在这个数据集中，每个观测值均有犯罪率和年份两个属性，也就是为面板数据的形式。好在对于87年的数据，已经给出对应82年的犯罪率lcrmrte.1 ，故选取87年数据直接回归即可。
## need to correct data source
load("data/CRIME2.rda")
OLS87 <- lm(lcrmrte ~ lcrmrt.1 + llawexpc + unem, data = CRIME2,
subset = (year == 87))
summary(OLS87)对于更多时间序列问题和面板数据问题将在后面几章中详细阐述。
这一章开始我们着重关注面板数据。面板数据在当代计量经济学尤其是应用计量经济学领域的重要性不必多说，然而正是由于它在数据结构上拥有了时间和个体二重维度，所以其估计方法也较截面数据更为复杂。但从另外一个层面来说，由于其时间维度的引入使得很多截面数据下难以估计的量可以实现一致估计，所以其应用的现实意义颇大。此章主要围绕plm包展开，该包提供了大量的成熟的面板数据分析工具。
在对付面板数据的时候，一阶差分法有时是一种简便易行且有效的方法。如 第14章所指出的，
在只有两期数据的情况下，一阶差分法（FD）和固定效应模型（FE）等价；当多于两期的时候，两者在满足假定的条件下都是小样本下无偏且大样本下一致的。然而在误差项\(u_{it}\)序列不相关时，FE比FD更为有效；当\(u_{it}\)为随机游走（Random Walk）时，FD显然更有效；当\(\Delta u_{it}\)呈现序列负相关特性时，FE更有效。 然而在长面板（即时间维度长而观测个体少）的情况下，FE更为敏感，故FD更有优势。更需值得注意的是，FD和FE都对解释变量是否服从经典假设很敏感，但是在解释变量和残差项不相关的情况下，即使其他假设被违反，FE估计量比FD的偏差会小些（除非时间\(T=2\)）。”
在R中，使用pml包中的plm()函数就可以完成该估计，由于该函数同时可同于估计多种面板数据模型（包括固定效应模型、混合模型pooled model、随机效应模型、一阶差分法、组间估计法between model等），所以用法在下面的例子中一并给出。
由于固定效应模型允许不可观测效应与解释变量之间存在相关性，而随机效应模型则不可以，所以FE被广泛认为在其他条件相同下更有说服力。不过随机效应模型在特定情况下也需得到应有，比如当关键的解释变量不随时间变化的时候，我们就无法用FE去估计它的效应。实际上在更多的时候，我们会两种估计都进行，然后进行进一步的检验（如Hausman检验，见后）。
(@WAGEPAN) 在劳动经济学中，我们经常关注各种因素对于工资的影响。WAGEPAN.rda这个数据集中，有对545个劳动者1980到1987年间从事工作的情况调查数据。其中一些变量会随着时间而改变，比如工作经验(\(exper\))、婚姻状况(\(married\))和工会状况(\(union\))。还有一些变量不会改变，比如种族(\(hispan\)和\(black\))、教育(\(educ\))。
在下面我们依次进行混合回归、一阶差分法、固定效应、随机效应来估计该例。然而在进行分析之前，我们需要先调用plm.data()命令来将现有数据集中的数据转换为面板数据模式。
load("data/WAGEPAN.rdata")
library("plm")
WAGE_data <- plm.data(WAGEPAN, index = c("nr", "year"))
# 混合OLS
WAGE_PLM_Pooled <- plm(lwage ~ educ + black + hisp + exper + I(exper^2) +
married + union, data = WAGE_data, model = "pooling")
summary(WAGE_PLM_Pooled)## Oneway (individual) effect Pooling Model
##
## Call:
## plm(formula = lwage ~ educ + black + hisp + exper + I(exper^2) +
## married + union, data = WAGE_data, model = "pooling")
##
## Balanced Panel: n=545, T=8, N=4360
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -5.2700 -0.2490 0.0332 0.2960 2.5600
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -0.034706 0.064569 -0.54 0.59
## educ 0.099388 0.004678 21.25 < 2e-16 ***
## black -0.143842 0.023560 -6.11 1.1e-09 ***
## hisp 0.015698 0.020811 0.75 0.45
## exper 0.089179 0.010111 8.82 < 2e-16 ***
## I(exper^2) -0.002849 0.000707 -4.03 5.7e-05 ***
## married 0.107666 0.015696 6.86 7.9e-12 ***
## union 0.180073 0.017121 10.52 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1240
## Residual Sum of Squares: 1010
## R-Squared : 0.187
## Adj. R-Squared : 0.186
## F-statistic: 142.613 on 7 and 4352 DF, p-value: <2e-16# 一阶差分
WAGE_PLM_FD <- plm(lwage ~ educ + black + hisp + exper + I(exper^2) +
married + union, data = WAGE_data, model = "fd")
summary(WAGE_PLM_FD)## Oneway (individual) effect First-Difference Model
##
## Call:
## plm(formula = lwage ~ educ + black + hisp + exper + I(exper^2) +
## married + union, data = WAGE_data, model = "fd")
##
## Balanced Panel: n=545, T=8, N=4360
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -4.5800 -0.1460 -0.0127 0.1330 4.8400
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (intercept) 0.11575 0.01959 5.91 3.7e-09 ***
## I(exper^2) -0.00388 0.00139 -2.80 0.0051 **
## married 0.03814 0.02293 1.66 0.0963 .
## union 0.04279 0.01966 2.18 0.0296 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 751
## Residual Sum of Squares: 748
## R-Squared : 0.0042
## Adj. R-Squared : 0.0042
## F-statistic: 5.36242 on 3 and 3811 DF, p-value: 0.0011# 固定效应模型
WAGE_PLM_fixed <- plm(lwage ~ educ + black + hisp + exper + I(exper^2) +
married + union, data = WAGE_data, model = "within")
summary(WAGE_PLM_fixed)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = lwage ~ educ + black + hisp + exper + I(exper^2) +
## married + union, data = WAGE_data, model = "within")
##
## Balanced Panel: n=545, T=8, N=4360
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -4.17000 -0.12600 0.00925 0.16000 1.47000
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## exper 0.116847 0.008420 13.88 < 2e-16 ***
## I(exper^2) -0.004301 0.000605 -7.11 1.4e-12 ***
## married 0.045303 0.018310 2.47 0.013 *
## union 0.082087 0.019291 4.26 2.1e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 572
## Residual Sum of Squares: 470
## R-Squared : 0.178
## Adj. R-Squared : 0.156
## F-statistic: 206.375 on 4 and 3811 DF, p-value: <2e-16# 随机效应模型
WAGE_PLM_random <- plm(lwage ~ educ + black + hisp + exper + I(exper^2) +
married + union, data = WAGE_data, model = "random")
summary(WAGE_PLM_random)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = lwage ~ educ + black + hisp + exper + I(exper^2) +
## married + union, data = WAGE_data, model = "random")
##
## Balanced Panel: n=545, T=8, N=4360
##
## Effects:
## var std.dev share
## idiosyncratic 0.123 0.351 0.54
## individual 0.105 0.325 0.46
## theta: 0.643
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -4.5800 -0.1450 0.0235 0.1860 1.5400
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -0.107464 0.110706 -0.97 0.33174
## educ 0.101225 0.008913 11.36 < 2e-16 ***
## black -0.144131 0.047615 -3.03 0.00248 **
## hisp 0.020151 0.042601 0.47 0.63622
## exper 0.112119 0.008261 13.57 < 2e-16 ***
## I(exper^2) -0.004069 0.000592 -6.88 7.1e-12 ***
## married 0.062795 0.016773 3.74 0.00018 ***
## union 0.107379 0.017830 6.02 1.9e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 657
## Residual Sum of Squares: 540
## R-Squared : 0.178
## Adj. R-Squared : 0.178
## F-statistic: 134.85 on 7 and 4352 DF, p-value: <2e-16# 组内模型
WAGE_PLM_between <- plm(lwage ~ educ + black + hisp + exper + I(exper^2) +
married + union, data = WAGE_data, model = "between")
summary(WAGE_PLM_between)## Oneway (individual) effect Between Model
##
## Call:
## plm(formula = lwage ~ educ + black + hisp + exper + I(exper^2) +
## married + union, data = WAGE_data, model = "between")
##
## Balanced Panel: n=545, T=8, N=4360
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.1200 -0.2390 0.0246 0.2300 1.7400
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 0.49231 0.22101 2.23 0.02632 *
## educ 0.09460 0.01090 8.68 < 2e-16 ***
## black -0.13881 0.04887 -2.84 0.00468 **
## hisp 0.00478 0.04269 0.11 0.91097
## exper -0.05044 0.05033 -1.00 0.31676
## I(exper^2) 0.00512 0.00321 1.60 0.11119
## married 0.14366 0.04120 3.49 0.00053 ***
## union 0.27068 0.04656 5.81 1.1e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 83.1
## Residual Sum of Squares: 64.9
## R-Squared : 0.219
## Adj. R-Squared : 0.216
## F-statistic: 21.5386 on 7 and 537 DF, p-value: <2e-16一般说来，OLS估计（混合回归）不如其他方法有效。而在上例中我们可以看出，为了估计不随时间变化的量（种族、受教育程度），我们必须用到随机效应模型。然而对于随时间变化的量（经验等），固定效应模型或一阶差分法显然更为有效。此外，对于随机效应模型，还可以使用fixef()函数来提取其中的固定效应。
实际上，plm()函数不仅仅支持这些基本的模型，而且对于每种模型还提供了不同参数以调整具体方法，如双向模型（two-way model）和时间效应等。例如，在使用双向模型的时候，只需在调用函数plm()中加一个参数effect = "twoways"即可，但此时如果使用fixef()函数来提取其中的固定效应则需写成fixef(YOUR_MODLE, effect = "time")。具体说来，随机效应模型支持增加参数random.method，可选择 amemiya、swar、walhus、nerlove，而对应的误差分量的变化可以直接调用ercomp()函数指定method和effect来估计。
如果在不同时期的个体是不同的（增加或减少），那么该面板数据就成为了非平衡面板。目前plm包中对于非平衡面板还不能进行双向效应估计，能提供的误差分量估计法也只有一种。Baltagi(2001)重新估计了房地产市场的特徽价格方程（hedonic housing prices function），代码如下：
data("Hedonic", package = "plm")
Hed <- plm(mv ~ crim + zn + indus + chas + nox + rm + age + dis +
rad + tax + ptratio + blacks + lstat, Hedonic, model = "random", index = "townid")
summary(Hed)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = mv ~ crim + zn + indus + chas + nox + rm + age +
## dis + rad + tax + ptratio + blacks + lstat, data = Hedonic,
## model = "random", index = "townid")
##
## Unbalanced Panel: n=92, T=1-30, N=506
##
## Effects:
## var std.dev share
## idiosyncratic 0.0170 0.1302 0.5
## individual 0.0168 0.1297 0.5
## theta :
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.292 0.590 0.666 0.650 0.745 0.820
##
## Residuals :
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.641 -0.066 -0.001 -0.002 0.070 0.527
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 9.68e+00 2.07e-01 46.72 < 2e-16 ***
## crim -7.23e-03 1.03e-03 -6.99 8.9e-12 ***
## zn 3.96e-05 6.88e-04 0.06 0.95414
## indus 2.08e-03 4.34e-03 0.48 0.63208
## chasyes -1.06e-02 2.90e-02 -0.37 0.71473
## nox -5.86e-03 1.25e-03 -4.71 3.3e-06 ***
## rm 9.18e-03 1.18e-03 7.78 4.2e-14 ***
## age -9.27e-04 4.65e-04 -2.00 0.04657 *
## dis -1.33e-01 4.57e-02 -2.91 0.00379 **
## rad 9.69e-02 2.83e-02 3.42 0.00069 ***
## tax -3.75e-04 1.89e-04 -1.98 0.04799 *
## ptratio -2.97e-02 9.75e-03 -3.05 0.00243 **
## blacks 5.75e-01 1.01e-01 5.69 2.2e-08 ***
## lstat -2.85e-01 2.39e-02 -11.95 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 893
## Residual Sum of Squares: 8.68
## R-Squared : 0.99
## Adj. R-Squared : 0.963
## F-statistic: 3854.18 on 13 and 492 DF, p-value: <2e-16在进行面板数据分析前，一般会考虑到一个问题：对于每个个体，解释变量系数都是一样的吗？换言之，个体效应是否存在？当然，一般说来个体效应都会存在的，否则后面的分析就没有意义了，直接采用混合模型就好。可混合性检验就是用来检验个体效应的。我们这里需要调用pooltest()函数进行检验。回到例2.8，继续看公司的投资与真实价值和资本存量的关系。最简单的调用方法如下：
data(Grunfeld, package = "plm")
pooltest(inv ~ value + capital, data = Grunfeld, model = "within")##
## F statistic
##
## data: inv ~ value + capital
## F = 5.78, df1 = 18, df2 = 170, p-value = 1.219e-10
## alternative hypothesis: unstability
## 其实，该函数需要的由两部分组成，一为通过plm()得到的模型，二为通过pvcm()得到的固定效应模型。因此，上面的代码实则等价于： znp <- pvcm(inv~value + capital, data = Grunfeld, model = "within") zplm <- plm(inv~value + capital, data = Grunfeld) pooltest(zplm, znp)
实际上，这里的pvcm()函数进行的正是变系数回归，而可混合性检验的实质也就是对这两种模型进行比较的F检验。
用拉格朗格乘子法可以检验相对于混合模型的个体或时间效应。这里可以调用plmtest()函数。调用方式有两种，其一可以先做一个混合回归再进行检验（如双向效应）：
g_pooled <- plm(inv ~ value + capital, data = Grunfeld, model = "pooling")
plmtest(g_pooled, effect = "twoways", type = "ghm")##
## Lagrange Multiplier Test - two-ways effects (Gourieroux, Holly
## and Monfort)
##
## data: inv ~ value + capital
## chisq = 798.2, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effects
## 此外也可以直接在调用检验函数的同时指定回归模型，如：
plmtest(inv ~ value + capital, data = Grunfeld, effect = "twoways",
type = "ghm")##
## Lagrange Multiplier Test - two-ways effects (Gourieroux, Holly
## and Monfort)
##
## data: inv ~ value + capital
## chisq = 798.2, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effects
## 对于该函数来说，还可以附加参数type，可选bp, honda, kw, ghm之一。至于必须指明的effect参数，可选individual, time, twoways之一。
此外，还可以做各种效应的F检验（基于混合模型与固定效应模型之间的比较）。依旧用例2.8：
g_fixed_twoways <- plm(inv ~ value + capital, data = Grunfeld, effect = "twoways",
model = "within")## Error: cannot open file 'cache-upload/sur_b4f13e8b693a5cd09902dc33ebfcfc23.rdb': No such file or directorypFtest(g_fixed_twoways, g_pooled)## Error: object 'g_fixed_twoways' not found当然，pFtest()函数也提供了另外一种直接的调用形式，如pFtest(inv value + capital, data = Grunfeld, effect = "twoways")。
为了判断到底是采用固定效应还是随机效应模型，我们一般采取经典的Hausman检验，这可以通过调用phtest()函数实现。注：这里我们只展示Hausman检验的用法，至于实际应用中该检验结果是否足以支撑模型选择，则需另作学术讨论。
gw <- plm(inv ~ value + capital, data = Grunfeld, model = "within")
gr <- plm(inv ~ value + capital, data = Grunfeld, model = "random")
phtest(gw, gr)##
## Hausman Test
##
## data: inv ~ value + capital
## chisq = 2.33, df = 2, p-value = 0.3119
## alternative hypothesis: one model is inconsistent
## 同理，对于例4.1，我们也可以采用该检验来判断应该使用哪个模型。因为我们上面已经分别作了固定效应和随机效应的回归，所以此时直接调用两个回归模型结果就可以。
phtest(WAGE_PLM_fixed, WAGE_PLM_random)##
## Hausman Test
##
## data: lwage ~ educ + black + hisp + exper + I(exper^2) + married + union
## chisq = 31.45, df = 4, p-value = 2.476e-06
## alternative hypothesis: one model is inconsistent
## 注意Hasuman检验的原假设是“随机效应模型优于固定效应模型”，所以对于例2.8，我们可以接受原假设，采用随机效应模型；对于例4.1，则应考虑固定效应模型。
一般说来在固定或随机效应模型中，系数都被假设为不变的而截距项是变化的。下面我们可以考虑系数也变化的情形。依旧以2.8（Grunfeld）为例，欲使用变系数模型我们只需要调用pvcm()函数即可。
g_varw <- pvcm(inv ~ value + capital, data = Grunfeld, model = "within")
g_varr <- pvcm(inv ~ value + capital, data = Grunfeld, model = "random")## [1] 3.340e-02 1.633e-03 -1.120e+03
## attentionsummary(g_varr)## Oneway (individual) effect Random coefficients model
##
## Call:
## pvcm(formula = inv ~ value + capital, data = Grunfeld, model = "random")
##
## Balanced Panel: n=10, T=20, N=200
##
## Residuals:
## total sum of squares : 2178000
## id time
## 0.67678 0.02974
##
## Estimated mean of the coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -9.6293 17.0350 -0.57 0.57189
## value 0.0846 0.0200 4.24 2.2e-05 ***
## capital 0.1994 0.0527 3.79 0.00015 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Estimated variance of the coefficients:
## (Intercept) value capital
## (Intercept) 2344.244 -0.68523 -4.02766
## value -0.685 0.00312 -0.00118
## capital -4.028 -0.00118 0.02448
##
## Total Sum of Squares: 4.74e+08
## Residual Sum of Squares: 2190000
## Multiple R-Squared: 0.995前面说到，使用固定效应模型的时候无法估计不随时间变化的量的影响。在没有理想的外部工具变量的情况下，我们可以采用Hausman-Taylor估计法来进行工具变量的运用。但是该方法有一个前提假设，就是有部分变量为外生，和残差项不相关；另一部分变量为内生，和残差项相关，在具体研究中需要证明这一点才能保证估计的一致性。
调用该方法的函数为pht()。我们下面再来看一个工资的问题，这个时候使用Wages数据集，其包含了1976到1982年美国595个个体的数据。
data("Wages", package = "plm")
ht <- plm(lwage ~ wks + south + smsa + married + exp + I(exp^2) +
bluecol + ind + union + sex + black + ed | sex + black + bluecol + south +
smsa + ind, data = Wages, model = "ht", index = 595)
summary(ht)## Oneway (individual) effect Hausman-Taylor Model
## Call:
## pht(formula = lwage ~ wks + south + smsa + married + exp + I(exp^2) +
## bluecol + ind + union + sex + black + ed | sex + black +
## bluecol + south + smsa + ind, data = Wages, index = 595)
##
## T.V. exo : bluecol,south,smsa,ind
## T.V. endo : wks,married,exp,I(exp^2),union
## T.I. exo : sex,black
## T.I. endo : ed
##
## Balanced Panel: n=595, T=7, N=4165
##
## Effects:
## var std.dev share
## idiosyncratic 0.023 0.152 0.03
## individual 0.887 0.942 0.97
## theta: 0.939
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.92000 -0.07070 0.00657 0.07970 2.03000
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 2.78e+00 3.08e-01 9.04 < 2e-16 ***
## wks 8.37e-04 6.00e-04 1.40 0.163
## southyes 7.44e-03 3.20e-02 0.23 0.816
## smsayes -4.18e-02 1.90e-02 -2.21 0.027 *
## marriedyes -2.99e-02 1.90e-02 -1.57 0.116
## exp 1.13e-01 2.47e-03 45.79 < 2e-16 ***
## I(exp^2) -4.19e-04 5.46e-05 -7.67 1.7e-14 ***
## bluecolyes -2.07e-02 1.38e-02 -1.50 0.133
## ind 1.36e-02 1.52e-02 0.89 0.372
## unionyes 3.28e-02 1.49e-02 2.20 0.028 *
## sexmale 1.31e-01 1.27e-01 1.03 0.301
## blackyes -2.86e-01 1.56e-01 -1.84 0.066 .
## ed 1.38e-01 2.12e-02 6.49 8.5e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 887
## Residual Sum of Squares: 95.9
## F-statistic: 2852.33 on 12 and 4152 DF, p-value: <2e-16当存在一些外部的工具变量满足工具变量的要求的时候，我们可以直接使用它们代替原有变量进行回归，只是需要在调用plm()函数的时候用 符号指定就好。
我们依旧来看那个犯罪率的问题。这里我们扩展到整个面板数据，调用plm包中的Crime数据集，里面包含了1981-1987年间美国90个郡的数据。我们这里依旧想研究影响犯罪率的因素。
因为变量比较多，所以不一一解释其含义，可以直接参照包中的说明。下面我们希望以\(log(taxpc)\)和\(log(mix)\)作为工具变量，分别代替\(log(prbarr)\)和\(log(polpc)\)，故只需要在原始回归方程中加入. - log(prbarr) - log(polpc) + log(taxpc) + log(mix)即可。
data("Crime", package = "plm")
cr <- plm(log(crmrte) ~ log(prbarr) + log(polpc) + log(prbconv) +
log(prbpris) + log(avgsen) + log(density) + log(wcon) + log(wtuc) + log(wtrd) +
log(wfir) + log(wser) + log(wmfg) + log(wfed) + log(wsta) + log(wloc) +
log(pctymle) + log(pctmin) + region + smsa + factor(year) | . - log(prbarr) -
log(polpc) + log(taxpc) + log(mix), data = Crime, model = "random")
summary(cr)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
## Instrumental variable estimation
## (Balestra-Varadharajan-Krishnakumar's transformation)
##
## Call:
## plm(formula = log(crmrte) ~ log(prbarr) + log(polpc) + log(prbconv) +
## log(prbpris) + log(avgsen) + log(density) + log(wcon) + log(wtuc) +
## log(wtrd) + log(wfir) + log(wser) + log(wmfg) + log(wfed) +
## log(wsta) + log(wloc) + log(pctymle) + log(pctmin) + region +
## smsa + factor(year) | . - log(prbarr) - log(polpc) + log(taxpc) +
## log(mix), data = Crime, model = "random")
##
## Balanced Panel: n=90, T=7, N=630
##
## Effects:
## var std.dev share
## idiosyncratic 0.0223 0.1492 0.33
## individual 0.0460 0.2146 0.67
## theta: 0.746
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -5.0200 -0.4760 0.0273 0.5260 3.1900
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -0.45382 1.70298 -0.27 0.7900
## log(prbarr) -0.41412 0.22105 -1.87 0.0615 .
## log(polpc) 0.50493 0.22778 2.22 0.0270 *
## log(prbconv) -0.34324 0.13247 -2.59 0.0098 **
## log(prbpris) -0.19004 0.07334 -2.59 0.0098 **
## log(avgsen) -0.00644 0.02894 -0.22 0.8241
## log(density) 0.43435 0.07115 6.10 1.8e-09 ***
## log(wcon) -0.00430 0.04142 -0.10 0.9174
## log(wtuc) 0.04446 0.02154 2.06 0.0395 *
## log(wtrd) -0.00856 0.04198 -0.20 0.8385
## log(wfir) -0.00403 0.02946 -0.14 0.8912
## log(wser) 0.01056 0.02158 0.49 0.6248
## log(wmfg) -0.20179 0.08394 -2.40 0.0165 *
## log(wfed) -0.21346 0.21511 -0.99 0.3214
## log(wsta) -0.06011 0.12031 -0.50 0.6175
## log(wloc) 0.18351 0.13967 1.31 0.1894
## log(pctymle) -0.14584 0.22681 -0.64 0.5205
## log(pctmin) 0.19488 0.04594 4.24 2.6e-05 ***
## regionwest -0.22818 0.10103 -2.26 0.0243 *
## regioncentral -0.19877 0.06075 -3.27 0.0011 **
## smsayes -0.25954 0.14998 -1.73 0.0840 .
## factor(year)82 0.01321 0.02999 0.44 0.6597
## factor(year)83 -0.08477 0.03200 -2.65 0.0083 **
## factor(year)84 -0.10620 0.03879 -2.74 0.0064 **
## factor(year)85 -0.09774 0.05117 -1.91 0.0566 .
## factor(year)86 -0.07194 0.06058 -1.19 0.2355
## factor(year)87 -0.03965 0.07585 -0.52 0.6013
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 30.2
## Residual Sum of Squares: 558
## R-Squared : 0.592
## Adj. R-Squared : 0.567
## F-statistic: -21.9376 on 26 and 603 DF, p-value: 1当回归方程中有被解释变量的滞后值作为解释变量的时候，便形成了动态面板数据。由于滞后项的存在，原来的估计成为了不一致的估计，即“动态面板偏差”。对于此种情况的出现，分别有差分GMM（difference GMM）和系统GMM（system GMM）可以用来解决。R中对应的函数为pgmm()。
在EmplUK数据集中，涵盖了1976到1984年间英国工厂的就业和工资情况的数据，且该面板为非平衡面板。我们这里关系就业情况。由于上一期的就业会很直接的影响当期就业水平，所以回归方程中包含了就业作为解释变量，该问题成为了动态面板分析。
在这里我们首先用差分GMM来估计，选择双向效应和两步法。
data("EmplUK", package = "plm")
emp.gmm <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) +
log(capital) + lag(log(output), 0:1) | lag(log(emp), 2:99), data = EmplUK,
effect = "twoways", model = "twosteps")
summary(emp.gmm)## Twoways effects Two steps model
##
## Call:
## pgmm(formula = log(emp) ~ lag(log(emp), 1:2) + lag(log(wage),
## 0:1) + log(capital) + lag(log(output), 0:1) | lag(log(emp),
## 2:99), data = EmplUK, effect = "twoways", model = "twosteps")
##
## Unbalanced Panel: n=140, T=7-9, N=1031
##
## Number of Observations Used: 611
##
## Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.6190 -0.0256 0.0000 -0.0001 0.0332 0.6410
##
## Coefficients
## Estimate Std. Error z-value Pr(>|z|)
## lag(log(emp), 1:2)1 0.4742 0.0853 5.56 2.7e-08 ***
## lag(log(emp), 1:2)2 -0.0530 0.0273 -1.94 0.05222 .
## lag(log(wage), 0:1)0 -0.5132 0.0493 -10.40 < 2e-16 ***
## lag(log(wage), 0:1)1 0.2246 0.0801 2.81 0.00502 **
## log(capital) 0.2927 0.0395 7.42 1.2e-13 ***
## lag(log(output), 0:1)0 0.6098 0.1085 5.62 1.9e-08 ***
## lag(log(output), 0:1)1 -0.4464 0.1248 -3.58 0.00035 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sargan Test: chisq(25) = 30.11 (p.value=0.22)
## Autocorrelation test (1): normal = -2.428 (p.value=0.00759)
## Autocorrelation test (2): normal = -0.3325 (p.value=0.37)
## Wald test for coefficients: chisq(7) = 372 (p.value=<2e-16)
## Wald test for time dummies: chisq(6) = 26.9 (p.value=0.000151)当然我们也可以用系统GMM来估计，只需要加上一个参数transformation = "ld"即可。最后可以在summary()中加入robust = TRUE来查看模型的稳健性。
z2 <- pgmm(log(emp) ~ lag(log(emp), 1) + lag(log(wage), 0:1) + lag(log(capital),
0:1) | lag(log(emp), 2:99) + lag(log(wage), 2:99) + lag(log(capital), 2:99),
data = EmplUK, effect = "twoways", model = "onestep", transformation = "ld")
summary(z2, robust = TRUE)## Twoways effects One step model
##
## Call:
## pgmm(formula = log(emp) ~ lag(log(emp), 1) + lag(log(wage), 0:1) +
## lag(log(capital), 0:1) | lag(log(emp), 2:99) + lag(log(wage),
## 2:99) + lag(log(capital), 2:99), data = EmplUK, effect = "twoways",
## model = "onestep", transformation = "ld")
##
## Unbalanced Panel: n=140, T=7-9, N=1031
##
## Number of Observations Used: 1642
##
## Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.7530 -0.0369 0.0000 0.0003 0.0466 0.6000
##
## Coefficients
## Estimate Std. Error z-value Pr(>|z|)
## lag(log(emp), 1) 0.9356 0.0263 35.58 < 2e-16 ***
## lag(log(wage), 0:1)0 -0.6310 0.1181 -5.34 9.1e-08 ***
## lag(log(wage), 0:1)1 0.4826 0.1369 3.53 0.00042 ***
## lag(log(capital), 0:1)0 0.4839 0.0539 8.98 < 2e-16 ***
## lag(log(capital), 0:1)1 -0.4244 0.0585 -7.26 4.0e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sargan Test: chisq(100) = 236 (p.value=5.21e-13)
## Autocorrelation test (1): normal = -4.808 (p.value=7.61e-07)
## Autocorrelation test (2): normal = -0.28 (p.value=0.39)
## Wald test for coefficients: chisq(5) = 11175 (p.value=<2e-16)
## Wald test for time dummies: chisq(7) = 14.71 (p.value=0.0399)outreg()