-
Notifications
You must be signed in to change notification settings - Fork 61
/
knitr-book.html
2005 lines (2003 loc) · 174 KB
/
knitr-book.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta http-equiv="Content-Style-Type" content="text/css" />
<meta name="generator" content="pandoc" />
<meta name="author" content="陈丽云" />
<title>在R中玩转计量</title>
<style type="text/css">
table.sourceCode, tr.sourceCode, td.lineNumbers, td.sourceCode {
margin: 0; padding: 0; vertical-align: baseline; border: none; }
table.sourceCode { width: 100%; }
td.lineNumbers { text-align: right; padding-right: 4px; padding-left: 4px; color: #aaaaaa; border-right: 1px solid #aaaaaa; }
td.sourceCode { padding-left: 5px; }
code > span.kw { color: #007020; font-weight: bold; }
code > span.dt { color: #902000; }
code > span.dv { color: #40a070; }
code > span.bn { color: #40a070; }
code > span.fl { color: #40a070; }
code > span.ch { color: #4070a0; }
code > span.st { color: #4070a0; }
code > span.co { color: #60a0b0; font-style: italic; }
code > span.ot { color: #007020; }
code > span.al { color: #ff0000; font-weight: bold; }
code > span.fu { color: #06287e; }
code > span.er { color: #ff0000; font-weight: bold; }
</style>
<script src="https://d3eoax9i5htok0.cloudfront.net/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"></script>
</head>
<body>
<div id="header">
<h1 class="title">在R中玩转计量</h1>
<h2 class="author">陈丽云</h2>
</div>
<div id="TOC">
<ul>
<li><a href="#序言"><span class="toc-section-number">1</span> 序言</a></li>
<li><a href="#熟悉r"><span class="toc-section-number">2</span> 熟悉R</a><ul>
<li><a href="#数据的导入"><span class="toc-section-number">2.1</span> 数据的导入</a></li>
<li><a href="#数据分析"><span class="toc-section-number">2.2</span> 数据分析</a><ul>
<li><a href="#平均值"><span class="toc-section-number">2.2.1</span> 平均值</a></li>
<li><a href="#线性回归普通最小二乘法ols"><span class="toc-section-number">2.2.2</span> 线性回归（普通最小二乘法，OLS）</a></li>
<li><a href="#作回归图像"><span class="toc-section-number">2.2.3</span> 作回归图像</a></li>
<li><a href="#点预测"><span class="toc-section-number">2.2.4</span> 点预测</a></li>
<li><a href="#多元线性回归"><span class="toc-section-number">2.2.5</span> 多元线性回归</a></li>
<li><a href="#保存和编辑代码"><span class="toc-section-number">2.2.6</span> 保存和编辑代码</a></li>
</ul></li>
<li><a href="#寻求帮助"><span class="toc-section-number">2.3</span> 寻求帮助</a></li>
</ul></li>
<li><a href="#基础数据整理与分析"><span class="toc-section-number">3</span> 基础数据整理与分析</a><ul>
<li><a href="#软件安装"><span class="toc-section-number">3.1</span> 软件安装</a><ul>
<li><a href="#r核心"><span class="toc-section-number">3.1.1</span> R核心</a></li>
<li><a href="#r的ide编辑器"><span class="toc-section-number">3.1.2</span> R的IDE编辑器</a></li>
<li><a href="#安装r包"><span class="toc-section-number">3.1.3</span> 安装R包</a></li>
</ul></li>
<li><a href="#数据的导入-1"><span class="toc-section-number">3.2</span> 数据的导入</a><ul>
<li><a href="#指定工作目录"><span class="toc-section-number">3.2.1</span> 指定工作目录</a></li>
<li><a href="#加载分析包"><span class="toc-section-number">3.2.2</span> 加载分析包</a></li>
<li><a href="#订阅某个主题下的r包及其更新"><span class="toc-section-number">3.2.3</span> 订阅某个主题下的R包及其更新</a></li>
</ul></li>
<li><a href="#数据来源"><span class="toc-section-number">3.3</span> 数据来源</a><ul>
<li><a href="#文本格式"><span class="toc-section-number">3.3.1</span> 文本格式</a></li>
<li><a href="#excel-格式通过odbc接口"><span class="toc-section-number">3.3.2</span> Excel 格式：通过ODBC接口</a></li>
<li><a href="#excel-格式通过xlsx包"><span class="toc-section-number">3.3.3</span> Excel 格式：通过xlsx包</a></li>
<li><a href="#数据库sql"><span class="toc-section-number">3.3.4</span> 数据库（SQL）</a></li>
<li><a href="#其他数据格式"><span class="toc-section-number">3.3.5</span> 其他数据格式</a></li>
</ul></li>
<li><a href="#变量操作"><span class="toc-section-number">3.4</span> 变量操作</a><ul>
<li><a href="#变量重命名"><span class="toc-section-number">3.4.1</span> 变量重命名</a></li>
</ul></li>
<li><a href="#data.frame的行列操作"><span class="toc-section-number">3.5</span> data.frame的行、列操作</a><ul>
<li><a href="#数据格式的转换"><span class="toc-section-number">3.5.1</span> 数据格式的转换</a></li>
<li><a href="#新变量生成"><span class="toc-section-number">3.5.2</span> 新变量生成</a></li>
<li><a href="#数据集合并操作"><span class="toc-section-number">3.5.3</span> 数据集合并操作</a></li>
<li><a href="#数据集形状的转换reshape2包"><span class="toc-section-number">3.5.4</span> 数据集形状的转换（reshape2包）</a></li>
</ul></li>
<li><a href="#数据的导出"><span class="toc-section-number">3.6</span> 数据的导出</a></li>
<li><a href="#分类统计data.table包"><span class="toc-section-number">3.7</span> 分类统计（<strong>data.table</strong>包）</a></li>
<li><a href="#循环和判断"><span class="toc-section-number">3.8</span> 循环和判断</a></li>
</ul></li>
<li><a href="#从截面数据分析说起"><span class="toc-section-number">4</span> 从截面数据分析说起</a><ul>
<li><a href="#参数检验"><span class="toc-section-number">4.1</span> 参数检验</a><ul>
<li><a href="#单变量显著性t检验"><span class="toc-section-number">4.1.1</span> 单变量显著性：t检验</a></li>
<li><a href="#多个变量显著性f检验"><span class="toc-section-number">4.1.2</span> 多个变量显著性：F检验</a></li>
</ul></li>
<li><a href="#置信区间"><span class="toc-section-number">4.2</span> 置信区间</a></li>
<li><a href="#虚拟变量"><span class="toc-section-number">4.3</span> 虚拟变量</a><ul>
<li><a href="#按性质分组"><span class="toc-section-number">4.3.1</span> 按性质分组</a></li>
<li><a href="#按数量值分组"><span class="toc-section-number">4.3.2</span> 按数量值分组</a></li>
<li><a href="#交叉项intersections"><span class="toc-section-number">4.3.3</span> 交叉项（intersections）</a></li>
<li><a href="#指定参照组"><span class="toc-section-number">4.3.4</span> 指定参照组</a></li>
</ul></li>
<li><a href="#异方差检验"><span class="toc-section-number">4.4</span> 异方差检验</a><ul>
<li><a href="#bp检验-breusch-pagan-test"><span class="toc-section-number">4.4.1</span> BP检验 (Breusch-Pagan Test)</a></li>
<li><a href="#怀特检验-white-test-for-heteroskedasticity"><span class="toc-section-number">4.4.2</span> 怀特检验 (White test for heteroskedasticity)</a></li>
</ul></li>
<li><a href="#稳健标准差"><span class="toc-section-number">4.5</span> 稳健标准差</a></li>
<li><a href="#加权最小二乘估计-wls"><span class="toc-section-number">4.6</span> 加权最小二乘估计 (WLS)</a><ul>
<li><a href="#扰动项形式已知"><span class="toc-section-number">4.6.1</span> 扰动项形式已知<需添加例子></a></li>
<li><a href="#扰动项形式未知feasible-linear-regression"><span class="toc-section-number">4.6.2</span> 扰动项形式未知（Feasible Linear Regression）</a></li>
</ul></li>
<li><a href="#广义线性估计glm"><span class="toc-section-number">4.7</span> 广义线性估计(GLM)</a><ul>
<li><a href="#最大似然估计-maximum-likelihood-estimation-mle"><span class="toc-section-number">4.7.1</span> 最大似然估计 (Maximum Likelihood Estimation, MLE)</a></li>
<li><a href="#二元被解释变量probit和logit模型"><span class="toc-section-number">4.7.2</span> 二元被解释变量：Probit和Logit模型</a></li>
<li><a href="#离散被解释变量"><span class="toc-section-number">4.7.3</span> 离散被解释变量</a></li>
<li><a href="#边角解tobit模型"><span class="toc-section-number">4.7.4</span> 边角解：Tobit模型</a></li>
<li><a href="#有序的probitlogit模型ordered-logitprobit"><span class="toc-section-number">4.7.5</span> 有序的probit/logit模型（Ordered Logit/Probit）</a></li>
</ul></li>
<li><a href="#计数模型-count-model"><span class="toc-section-number">4.8</span> 计数模型 (Count Model)</a><ul>
<li><a href="#泊松回归-poisson-regression-model"><span class="toc-section-number">4.8.1</span> 泊松回归 (Poisson Regression Model)</a></li>
<li><a href="#过度离散数据检验"><span class="toc-section-number">4.8.2</span> 过度离散数据检验</a></li>
<li><a href="#负二项回归模型negative-binomial-regression"><span class="toc-section-number">4.8.3</span> 负二项回归模型（negative binomial regression）</a></li>
<li><a href="#零膨胀泊松模型-zero-inflated-poisson-model-zip"><span class="toc-section-number">4.8.4</span> 零膨胀泊松模型 (Zero-inflated Poisson Model, ZIP)</a></li>
</ul></li>
<li><a href="#选择性样本问题"><span class="toc-section-number">4.9</span> 选择性样本问题</a><ul>
<li><a href="#被审查的样本censored-regression-models"><span class="toc-section-number">4.9.1</span> 被审查的样本（Censored Regression Models）</a></li>
<li><a href="#被截断的样本truncated-regression-models"><span class="toc-section-number">4.9.2</span> 被截断的样本（Truncated Regression Models）</a></li>
<li><a href="#heckit模型"><span class="toc-section-number">4.9.3</span> Heckit模型</a></li>
</ul></li>
<li><a href="#工具变量法instrument-variable-estimation"><span class="toc-section-number">4.10</span> 工具变量法（Instrument Variable Estimation）</a><ul>
<li><a href="#内生性问题与两阶段最小二乘法2sls"><span class="toc-section-number">4.10.1</span> 内生性问题与两阶段最小二乘法(2SLS)</a></li>
<li><a href="#联立方程模型simultaneous-equations"><span class="toc-section-number">4.10.2</span> 联立方程模型（Simultaneous Equations）</a></li>
<li><a href="#内生性检验test-for-endogeneity"><span class="toc-section-number">4.10.3</span> 内生性检验（Test for Endogeneity）</a></li>
<li><a href="#过度识别约束检验test-for-overidentification-restrictions"><span class="toc-section-number">4.10.4</span> 过度识别约束检验（Test for Overidentification Restrictions）</a></li>
<li><a href="#联立方程模型估计似不相关回归法seemingly-unrelated-regression"><span class="toc-section-number">4.10.5</span> 联立方程模型估计：似不相关回归法（Seemingly Unrelated Regression）</a></li>
</ul></li>
<li><a href="#代理变量-proxy-variables"><span class="toc-section-number">4.11</span> 代理变量 (Proxy Variables)</a></li>
</ul></li>
<li><a href="#面板数据分析"><span class="toc-section-number">5</span> 面板数据分析</a><ul>
<li><a href="#一阶差分法first-differenced"><span class="toc-section-number">5.1</span> 一阶差分法（First-Differenced）</a></li>
<li><a href="#固定效应模型fixed-effects-model和随机效应模型random-effects-model"><span class="toc-section-number">5.2</span> 固定效应模型（Fixed Effects Model）和随机效应模型（Random Effects Model）</a></li>
<li><a href="#非平衡面板unbalanced-panels"><span class="toc-section-number">5.3</span> 非平衡面板（Unbalanced panels）</a></li>
<li><a href="#面板数据相关检验"><span class="toc-section-number">5.4</span> 面板数据相关检验</a><ul>
<li><a href="#可混合性检验tests-of-poolability"><span class="toc-section-number">5.4.1</span> 可混合性检验（Tests of poolability）</a></li>
<li><a href="#个体效应或时间效应检验tests-for-individual-and-time-effects"><span class="toc-section-number">5.4.2</span> 个体效应或时间效应检验（Tests for individual and time effects）</a></li>
<li><a href="#hausman检验"><span class="toc-section-number">5.4.3</span> Hausman检验</a></li>
</ul></li>
<li><a href="#变系数模型variable-coeffcients-model"><span class="toc-section-number">5.5</span> 变系数模型（Variable coeffcients model）</a></li>
<li><a href="#面板工具变量法"><span class="toc-section-number">5.6</span> 面板工具变量法</a><ul>
<li><a href="#hausman-taylor估计法"><span class="toc-section-number">5.6.1</span> Hausman-Taylor估计法</a></li>
<li><a href="#外部工具变量"><span class="toc-section-number">5.6.2</span> 外部工具变量</a></li>
</ul></li>
<li><a href="#序列相关性检验tests-of-serial-correlation"><span class="toc-section-number">5.7</span> 序列相关性检验（Tests of serial correlation）</a></li>
<li><a href="#横截面的依赖性检验tests-for-cross-sectional-dependence"><span class="toc-section-number">5.8</span> 横截面的依赖性检验（Tests for cross-sectional dependence）</a></li>
<li><a href="#动态面板分析广义矩估计gmm"><span class="toc-section-number">5.9</span> 动态面板分析：广义矩估计（GMM）</a></li>
</ul></li>
<li><a href="#生存分析survival-analysis"><span class="toc-section-number">6</span> 生存分析（Survival Analysis）</a></li>
<li><a href="#分位数回归quantile-regression"><span class="toc-section-number">7</span> 分位数回归（Quantile Regression）</a></li>
<li><a href="#断点回归regression-discontinuity-designsrdd"><span class="toc-section-number">8</span> 断点回归（Regression Discontinuity Designs，RDD）</a></li>
<li><a href="#回归结果输出"><span class="toc-section-number">9</span> 回归结果输出</a><ul>
<li><a href="#函数outreg"><span class="toc-section-number">9.1</span> 函数<code>outreg()</code></a></li>
</ul></li>
<li><a href="#可重复的研究"><span class="toc-section-number">10</span> 可重复的研究</a><ul>
<li><a href="#r与latex等文档的结合"><span class="toc-section-number">10.1</span> R与latex等文档的结合</a><ul>
<li><a href="#knitr"><span class="toc-section-number">10.1.1</span> Knitr</a></li>
<li><a href="#markdown"><span class="toc-section-number">10.1.2</span> MarkDown</a></li>
<li><a href="#latex"><span class="toc-section-number">10.1.3</span> LaTeX</a></li>
<li><a href="#pandoc"><span class="toc-section-number">10.1.4</span> Pandoc</a></li>
</ul></li>
</ul></li>
<li><a href="#多人协作的研究"><span class="toc-section-number">11</span> 多人协作的研究</a><ul>
<li><a href="#github"><span class="toc-section-number">11.1</span> github</a></li>
</ul></li>
</ul>
</div>
<h1 id="序言"><a href="#TOC"><span class="header-section-number">1</span> 序言</a></h1>
<p>在经济学分析中不可避免的要和数据打交道，而目前数据分析中最主要的工具就是计量经济学。数据源于现实，而对待数据的态度方面，我更欣赏凯恩斯的观点：从数据中寻找直觉。既不是单纯的从计量的结果中寻求观点的佐证，也不是从归纳的角度来推理因果关系。这有些和“散点图是最好的统计图形”的观点有些不谋而合。但是数据本身的特性并不是简简单单的可以肉眼扫视原始数据就可以得出的，这个时候借助计量这个分析工具更有利于我们发现隐藏在原始数据背后的蛛丝马迹，进而寻求灵感。因此，玩转数据是做经济学研究必不可少的一个环节。有句话说得好：Let’s get our hands dirty with data first!</p>
<p>当然做计量的时候很依赖计算机软件，常用的有Eviews、Stata、SPSS、SAS等。可以看出，这和统计学中常用的软件惊人的一致。追根溯源，计量经济学本来就是从数理统计学中的回归分析等渐渐延伸出来的，所以其方法在统计软件中可以很容易的实现。近几年R的快速蓬勃发展使之成为了最前沿的统计软件，由于其良好的拓展性，大量的免费的包(package) 的出现使得R足以胜任最潮流的统计分析工作。因此，R也足以作为一个计量分析软件来处理计量经济学的问题。</p>
<p>我作为一个经济学专业的学生，机缘巧合接触到了R，并为之深深沉迷。2009年冬天在第二届中国R语言会议做了一个简单的“在计量和经济学中使用R”的报告后，感到有必要写一个简单的小册子，介绍各种计量经济学方法在R中的实现，也希望借此从丰富的实例数据中找寻更多的直觉。</p>
<p>这个小册子主要希望能对下列使用者有所裨益：</p>
<ol style="list-style-type: decimal">
<li>想了解经济学和计量经济学分析方法的统计学学生，尤其是有至于转到经济学方向的。</li>
<li>想使用更先进的统计软件R来分析计量经济问题的用户，尤其是想从Stata等转到R的。</li>
</ol>
<p>因水平所限，这本小册子将会比较简单，着重于介绍各种方法对应的R包和实现，帮助从未使用过R的朋友们尽快的熟悉、了解和应用这款软件。学习一个软件最好的方法无非是多多使用，因此除了囊括大量的实例，我想不出更好的办法。这些例子有些来源于现有出版的计量经济学书籍（例如伍德里奇的计量经济学导论2），也有些摘取于公开发表的论文。当然，这对我来说是一项浩瀚而繁重的工作，因此诸位朋友的帮助显得格外的珍贵。</p>
<p>从现有的关于计量经济学和R的书籍来说，从网上能找到几本英文的，大都是免费发行并具有非常高的质量。只是国内中文的资料还颇为零散。在撰写这个小册子的过程中，我参考了大量已有的成果并受益匪浅，也建议英文较好的朋友们直接去阅读相关的英文材料尤其是R包自带的介绍，相信会更深入的了解R。在这里特别要说的是AER（全称：Applied Econometrics with R）这个包，是配合同名的书发行的包。不过通过demo可以详尽的看到各个例子的R源代码，也带有丰富的数据集（来自格林的《计量经济分析》等有名的著作），是非常好的练手的包。</p>
<p>最后需要说明的是，这本小册子是我在担任统计之都中文论坛（<a href="http://cos.name/cn/"><code class="url">http://cos.name/cn/</code></a>）“经济统计版”版主的时候所撰写的。承蒙站长谢益辉兄和诸位骨干成员的大力帮助，此册子凝结了COS诸多成员的心血，换言之我只是一个代笔者而已。我们通过GIT这个多人协作平台共同完善，也借助了<strong>knitr</strong>包来结合R与markdown。这样高效且免费的开源平台使得我在撰写过程中受益匪浅，也使得本册子避免潜在的问题得以实现在互联网上的免费发行。</p>
<h1 id="熟悉r"><a href="#TOC"><span class="header-section-number">2</span> 熟悉R</a></h1>
<p>关于R的中文入门有很多，在这里不再一一枚举（也不是本册子的职责所在）。但我也不能在这里说，诸位客官回去学完R的基本结构再来吧。一是有点枯燥，二则我个人窃以为学习一个软件最好的办法就是不断的拿例子来磨练，遇到问题去网上搜寻，否则看再多的入门引导也只似过眼烟云，那些命令一会儿就抛到九霄云外了。所以我把会用到的一些东西直接融入在下面这个简单的例子中，有一些必要的说明，希望可以快速的熟悉R。</p>
<p>下面是来自Papke(1995)的一个例子。他研究的是一个退休金计划和计划的慷慨度。</p>
<ol style="list-style-type: example">
<li><p>在401K.DTA这个数据集中，我们关心两个变量。prate是在合法的工人中拥有活跃帐户的比例。mrate是用来衡量这个计划的匹配程度（用来代表慷慨度），即如果mrate = 0.5，则表示工人付出了$10，其工作单位相应的付出了$5。接下来，我们需要面对这么几个问题：</p>
<ol style="list-style-type: decimal">
<li>找到prate和mrate这两个变量的平均值。</li>
<li>对下面这个方程进行最简单的OLS回归：\(\hat{prate} = \hat{\beta}_0 + \hat{\beta}_1mrate\)，并报告\(R^2\)。</li>
<li>找到当mrate = 3.5的时候，prate的预测值。</li>
</ol></li>
</ol>
<p>做计量分析，离不开的就是数据。所以我们第一步先来导入需要的数据。</p>
<h2 id="数据的导入"><a href="#TOC"><span class="header-section-number">2.1</span> 数据的导入</a></h2>
<p>获取数据有很多办法，在R 里面通过<strong>foreign</strong>包可以读/写Minitab, S, SAS, SPSS, Stata, Systat等等格式的数据。当然，R本身是支持从文本文件（包括CSV格式）和剪贴板中直接读取数据的。此外，对于R包里面自带的数据集，我们可以直接用<code>data("name")</code>来加载数据集。这里我采取的是读取Stata的数据（DTA格式）。</p>
<p>当然，我们首先要加载<strong>foreign</strong>包，可以在R中直接点击“加载程序包”，也可以手动输入：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(foreign)</code></pre>
<p>然后就可以使用<code>read.dta()</code>命令：</p>
<pre class="sourceCode r"><code class="sourceCode r">Papke_1995 = <span class="kw">read.dta</span>(<span class="st">"data/401K.DTA"</span>)</code></pre>
<pre><code>## Warning message: cannot read factor labels from Stata 5 files</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(Papke_1995)</code></pre>
<pre><code>## 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 </code></pre>
<p>Papke_1995是我们赋值后在R里使用的数据表的名字。因为R是基于对象(object)的，所以我们需要在读取数据的时候指定数据存储的对象。同样的，后面会不断的用到对象这一概念。</p>
<p>如果觉得这些东西记起来比较麻烦，一个个字母的打起来也挺麻烦的，怎么办？好在有个包叫做<strong>Rcmdr</strong>。加载这个包之后就会出现图形界面，可以通过点击的方式来操作。</p>
<div class="figure">
<img src="imgs/import_data_from_rcmdr.JPG" alt="在R Commander中导入数据" /><p class="caption">在R Commander中导入数据</p>
</div>
<p>之后，R Commander会自动记录每一步对应的代码，可供下次重复使用。</p>
<h2 id="数据分析"><a href="#TOC"><span class="header-section-number">2.2</span> 数据分析</a></h2>
<h3 id="平均值"><a href="#TOC"><span class="header-section-number">2.2.1</span> 平均值</a></h3>
<p>在介绍关于平均值的函数前，先介绍另一个有用的函数<code>names()</code>。这个函数的作用是显示数据表中所有的变量名称。用法和效果见后面的代码例子。</p>
<p>我们可以使用<code>summary()</code>来获取该数据表的摘要信息，里面包含平均值、最大最小值 、中位数等。不过我们这里只关心两个变量<code>prate</code>和<code>mrate</code> ，所以也可以使用<code>numSummary()</code>（需加载<strong>abind</strong>包）。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">load</span>(<span class="st">"data/Papke_1995.rdata"</span>)
<span class="kw">names</span>(Papke_1995)</code></pre>
<pre><code>## [1] "prate" "mrate" "totpart" "totelg" "age" "totemp" "sole"
## [8] "ltotemp"</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(Papke_1995)</code></pre>
<pre><code>## 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 </code></pre>
<p>可以从上表中读出<code>prate</code>和<code>mrate</code>的平均值。</p>
<p><code>sumSummary()</code>也可以通过R Commander的图形界面实现。 <img src="imgs/numerical_summary_rcmdr.JPG" alt="R Commander里调用sumSummary()分析数据" /></p>
<h3 id="线性回归普通最小二乘法ols"><a href="#TOC"><span class="header-section-number">2.2.2</span> 线性回归（普通最小二乘法，OLS）</a></h3>
<p>在R里面进行线性回归还是比较容易的，直接使用<code>lm()</code>就可以。值得注意的是，由于R的面向对象特性，我们需要不断的赋值。对于赋值，有三种基本方法，分别可以用<code>-></code>、<code><-</code>和<code>=</code>实现，其中前两个是有方向的赋值，所以一般来说更为常用。比如我们可以对变量<code>mrate</code>和<code>prate</code> 求乘积，并将结果赋予一个新变量<code>mp</code>，则只需写成<code>mp <- mrate*prate</code>。</p>
<p>因此在做回归的时候写成：</p>
<pre class="sourceCode r"><code class="sourceCode r">RegModel <- <span class="kw">lm</span>(prate ~ mrate, <span class="dt">data =</span> Papke_1995)
<span class="kw">summary</span>(RegModel)</code></pre>
<pre><code>##
## 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
## </code></pre>
<p>这样RegModel里面就存储了这次回归所得的数据。</p>
<p>我们还可以采用<code>attach()</code>命令，这样就不用每次都指定回归向量所在的数据集了，直接写<code>RegModel<- lm(prate~mrate)</code>，然后就可以用<code>summary(RegModel)</code>来看回归的结果了。(注：通常情况下不建议使用<code>attach()</code>，可能会导致变量名的一定程度混乱，尤其是在函数封装的时候。)</p>
<p>可以看出估计后的回归方程应为：</p>
<p>\(\hat{prate}=83.0755+5.8611mrate\)</p>
<p>其中\(R^{2}\)为0.0747。呃，这个\(R^{2}\)为什么这么小？看看散点图就知道了。</p>
<h3 id="作回归图像"><a href="#TOC"><span class="header-section-number">2.2.3</span> 作回归图像</a></h3>
<p>我们可以直接用最简单的<code>plot()</code>命令作图(当然更好的一个选择可能是<em><a href="ggplot2" title="ggplot2">ggplot2</a></em>)，用法如下：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(Papke_1995$mrate, Papke_1995$prate)
<span class="kw">abline</span>(RegModel, <span class="dt">col =</span> <span class="st">"red"</span>)</code></pre>
<p>第二行命令是添加了那条回归拟合线。 可见这个图本来就很散，也难怪线性拟合效果这么差了。</p>
<h3 id="点预测"><a href="#TOC"><span class="header-section-number">2.2.4</span> 点预测</a></h3>
<p>最后，就是依赖估计方程做预测了。这里需要的是做一个点预测。R里面需要依据另一个数据集来预测，而且这个数据集中必须含有mrate 这个变量。新建一个数据集并赋值的办法有许多，最简单的就是直接赋值，方法如下：</p>
<pre class="sourceCode r"><code class="sourceCode r">mrate_new <- <span class="kw">data.frame</span>(<span class="dt">mrate =</span> <span class="fl">3.5</span>)</code></pre>
<p>另外一种不推荐的方式是，利用数据编辑框来手动输入:<code>mrate_new <- edit(as.data.frame(NULL))</code>。不推荐的原因是难以有效的追溯数据的来源和数值，尤其是违背“可重复的研究”精神。</p>
<p>之后再利用<code>predict()</code>就可得到所需的预测值了。</p>
<pre class="sourceCode r"><code class="sourceCode r">mrate_new <- <span class="kw">data.frame</span>(<span class="dt">mrate =</span> <span class="fl">3.5</span>)
<span class="kw">predict</span>(RegModel, mrate_new)</code></pre>
<pre><code>## 1
## 103.6 </code></pre>
<h3 id="多元线性回归"><a href="#TOC"><span class="header-section-number">2.2.5</span> 多元线性回归</a></h3>
<p>当然现实中我们很少做一元的线性回归，解释变量往往是两个或者更多。这可以依旧用上面的<code>lm()</code>。如下面这个例子，研究的是出勤率和ACT测试成绩、学习成绩之间的关系。</p>
<p>(@ATTEND)<br />在ATTEND.DTA这个数据集中，atndrte 是出勤率（采用百分比表示），ACT 为ACT测试的成绩，priGPA 是之前的学习平均分。我们需要估计如下的方程：</p>
<p>\(atndrte=\beta_{0}+\beta_{1}priGPA+\beta_{2}ACT+u\)</p>
<p>很显然，这里我们和上面的例子一样，代码和结果如下：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(foreign)
Attend <- <span class="kw">read.dta</span>(<span class="st">"data/attend.dta"</span>)
Reg2 <- <span class="kw">lm</span>(atndrte ~ priGPA + ACT, <span class="dt">data =</span> Attend)
<span class="kw">summary</span>(Reg2)</code></pre>
<pre><code>##
## 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
## </code></pre>
<h3 id="保存和编辑代码"><a href="#TOC"><span class="header-section-number">2.2.6</span> 保存和编辑代码</a></h3>
<p>虽然我们有 <em>RCommander</em> 创造的图形界面，但是每次都指定参数也是件很烦的事儿。因此养成一个好习惯，保存好上次运行的代码，下次直接在R里面调用就可以了，有什么修改的也只需要稍作调整即可。<em>RCommander</em>里面本身就有<code>File -> Save Script</code>，可以把Script Window里面所有代码存储为***.R的格式，从而方便下次调用。Script Window里面也是可以直接编辑代码的，删掉一些自己不想要的，调整个别的参数都是很方便的。</p>
<p>需要说明的是，.R文件就是告诉R应该怎么运行的文件，所以可以直接用文本编辑器软件打开并编辑。现在NotePad++, UltraEdit等等文本编辑软件都有支持R的插件，可以方便的把代码传送到R里面调用。R的基本界面中也是可以直接打开.R的脚本文件运行的。此外，推荐一个新兴的R编辑器——RStuidio，集成了R的各种窗口（<red>将在后续章节详述</red>）。</p>
<h2 id="寻求帮助"><a href="#TOC"><span class="header-section-number">2.3</span> 寻求帮助</a></h2>
<p>有了上述的例子，相信大家已经基本熟悉R了。那么遇到问题怎么办呢？比如<code>summary()</code>这个函数，对于不同的模型会有不同的用法，那么我们就需要去查看原始的帮助。在R中，最简单的办法就算再想要查看的命令前加一个<code>?</code>号。例如<code>?summary</code>之后就会蹦出来帮助页面了。这是查看某一包作者撰写原始文档的最快捷方式。此外也可以用两个连续的问号<code>??</code>来搜索所有相关的资料。</p>
<p>但是如果根本不知道有哪些命令，则需要去找包内原始的资料。可以直接在Google等搜索引擎里面搜寻，也可以查看R包自带的说明，亦可以参照各种书籍。总之方法很多，多多利用互联网是最好的办法。国内最佳的地方自然是<a href="http://cos.name/cn/" title="统计之都论坛的R版">统计之都论坛的R版</a>，里面有丰富的资料和资深的UseR为大家解惑。</p>
<p>下面，我们将介绍一些R的基础知识，包括软件的安装、配置和基本的数据清理工作，为后续的分析打下坚实的基础（曾经某位自身的数据分析专家说过，“在数据正式进入模型之前，最关键的就是<strong>数据整理</strong>”。这一步可能会耗费很多时间，但把数据整理为理想的格式会为后面的统计分析提供极大的便利）。</p>
<h1 id="基础数据整理与分析"><a href="#TOC"><span class="header-section-number">3</span> 基础数据整理与分析</a></h1>
<h2 id="软件安装"><a href="#TOC"><span class="header-section-number">3.1</span> 软件安装</a></h2>
<h3 id="r核心"><a href="#TOC"><span class="header-section-number">3.1.1</span> R核心</a></h3>
<p>R可以从<a href="http://www.r-project.org">www.r-project.org</a>下载，里面点击CRAN后可以选择CHINA的几个镜像，然后依据操作系统选择对应的版本即可。</p>
<h3 id="r的ide编辑器"><a href="#TOC"><span class="header-section-number">3.1.2</span> R的IDE编辑器</a></h3>
<p>安装完R之后，最好再选择一个顺手的编辑器，用于书写对应的R代码。这里推荐的是RStudio，可以从<a href="http://rstudio.org/">rstudio.org</a>下载。安装完成后，双击打开，界面如图[fig:RStudio]所示。</p>
<p>[fig:RStudio]RStudio界面</p>
<div class="figure">
<img src="imgs/rstudio.jpg" alt="image" /><p class="caption">image</p>
</div>
<p>RStudio界面，分为四块。左上角是脚本编辑框和数据浏览框；右上角是当前空间中数据情况和使用命令的历史记录；左下角是实际的R的代码执行界面和相应返回的结果；右下角是文件目录列表、画图展示区、R包目录区和帮助区。</p>
<h3 id="安装r包"><a href="#TOC"><span class="header-section-number">3.1.3</span> 安装R包</a></h3>
<p>目前我们的工作仅仅搭好了R的架子，还需要依据分析任务下载并安装对应的R包。这部分内容在后续中会详细依据案例介绍。</p>
<h2 id="数据的导入-1"><a href="#TOC"><span class="header-section-number">3.2</span> 数据的导入</a></h2>
<h3 id="指定工作目录"><a href="#TOC"><span class="header-section-number">3.2.1</span> 指定工作目录</a></h3>
<p>在R中，默认的工作目录依系统配置而变化，可以在直接启动R之后，通过<code>getwd()</code>命令来查看。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">getwd</span>()</code></pre>
<p>另外，如果是通过后缀名为.R的脚本文件来直接调用R，那么工作目录就为该脚本文件所在的目录。如果对于任何命令的参数等希望得到进一步的说明，那么可以在命令前加上“?”来直接调用帮助。比如，</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="st">`</span><span class="dt">?</span><span class="st">`</span>(getwd)</code></pre>
<p>这个时候R会弹出帮助文档窗口。</p>
<p>在我们撰写R脚本之前，往往希望事前指定一个工作目录，这个时候就可以利用<code>setwd(file_path)</code>.</p>
<p>R里面大多数操作都是面向对象的，所以上述代码的含义就是给<code>file_path</code>这个object赋上了<code>"E:/example/"</code>这个文本串作为其值。然后利用<code>setwd()</code>函数来指定工作目录。显然，这个代码和直接调用<code>setwd("E:/example/")</code>是等同的。</p>
<h3 id="加载分析包"><a href="#TOC"><span class="header-section-number">3.2.2</span> 加载分析包</a></h3>
<p>R作为一个开源软件，最大的特性就是有很多人在不断的贡献各种分析包，从基础的数据整理到高级的统计模型、可视化实现都有着相应的支持。R里面现在可以支持分析包括但不仅限于：</p>
<ul>
<li><p>Bayesian： Bayesian Inference（贝叶斯推断）</p></li>
<li><p>ChemPhys：Chemometrics and Computational Physics（化学计量和计算物理）</p></li>
<li><p>ClinicalTrials：Clinical Trial Design, Monitoring, and Analysis（临床试验设计、监测和分析）</p></li>
<li><p>Cluster：Cluster Analysis & Finite Mixture Models（聚类分析和有限混合模型）</p></li>
<li><p>Distributions： Probability Distributions（概率分布）</p></li>
<li><p>Econometrics： Computational Econometrics（计量经济学）</p></li>
<li><p>Environmetrics： Analysis of Ecological and Environmental Data（生态和环境学数据分析）</p></li>
<li><p>ExperimentalDesign： Design of Experiments (DoE) & Analysis of Experimental Data（实验设计和分析）</p></li>
<li><p>Finance：Empirical Finance（实证金融分析）</p></li>
<li><p>Genetics：Statistical Genetics（统计遗传学）</p></li>
<li><p>Graphics：Graphic Displays & Dynamic Graphics & Graphic Devices & Visualization（图形显示、动态图形、图形设备和可视化）</p></li>
<li><p>HighPerformanceComputing： High-Performance and Parallel Computing with R（高性能并行计算）</p></li>
<li><p>MachineLearning：Machine Learning & Statistical Learning（机器学习和统计学习）</p></li>
<li><p>MedicalImaging：Medical Image Analysis（医学图像分析）</p></li>
<li><p>Multivariate：Multivariate Statistics（多元统计）</p></li>
<li><p>NaturalLanguageProcessing：Natural Language Processing（自然语言处理）</p></li>
<li><p>OfficialStatistics：Official Statistics & Survey Methodology（官方统计和普查方法）</p></li>
<li><p>Optimization：Optimization and Mathematical Programming（最优化和数学规划）</p></li>
<li><p>Pharmacokinetics：Analysis of Pharmacokinetic Data（药代动力学数据分析）</p></li>
<li><p>Phylogenetics：Phylogenetics, Especially Comparative Methods（系统发育、特别是比较方法）</p></li>
<li><p>Psychometrics：Psychometric Models and Methods（心理测量模型和方法）</p></li>
<li><p>ReproducibleResearch：Reproducible Research（可重复的研究）</p></li>
<li><p>Robust：Robust Statistical Methods（稳健性统计模型）</p></li>
<li><p>SocialSciences：Statistics for the Social Sciences（社会科学统计）</p></li>
<li><p>Spatial：Analysis of Spatial Data（空间数据分析）</p></li>
<li><p>Survival：Survival Analysis（生存分析）</p></li>
<li><p>TimeSeries：Time Series Analysis（时间序列分析）</p></li>
<li><p>gR：gRaphical Models in R（图形模型）</p></li>
</ul>
<p>每一项的具体说明请参见：<a href="http://cran.r-project.org/web/views/">http://cran.r-project.org/web/views/</a>。</p>
<p>针对每一项任务，我们往往都需要加载不同的包。在第一次加载一个包以前，我们往往需要先安装。在Rstudio界面中，加载包的命令位于<code>Tools -> Install Packages</code>，然后可以输入包的名字来安装，如图[install packages]所示。</p>
<p>[install packages]安装包（在Rstudio中）</p>
<div class="figure">
<img src="imgs/install-packages.jpg" alt="image" /><p class="caption">image</p>
</div>
<div class="figure">
<img src="imgs/install-packages-2.jpg" alt="image" /><p class="caption">image</p>
</div>
<p>安装完包之后，可以直接通过<code>library()</code>命令加载。如我们常用的<code>data.table</code>这个包：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(data.table)</code></pre>
<p>注意：<strong>R里面严格区分大小写</strong>，所以大小写不对的话包是无法加载成功的。</p>
<h3 id="订阅某个主题下的r包及其更新"><a href="#TOC"><span class="header-section-number">3.2.3</span> 订阅某个主题下的R包及其更新</a></h3>
<p><a href="http://cran.r-project.org/web/views/">CRAN下的Task Views</a>不仅仅罗列了各个主题下的R包情况，其本身亦作为一个R包<em>ctv</em>出现。在R中安装这个包之后，可以很方便的安装所有主题下的包。比如，我们希望安装Econometrics这个主题下所有提及的包，那么:</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">"ctv"</span>)
<span class="kw">library</span>(<span class="st">"ctv"</span>)
<span class="kw">install.views</span>(<span class="st">"Econometrics"</span>) <span class="co">#安装Econometrics主题下的包</span>
<span class="kw">update.views</span>(<span class="st">"Econometrics"</span>) <span class="co">#更新Econometrics主题下的包</span></code></pre>
<p>对于刚刚接触R的用户来说，这样一次性下载并安装一系列R包，可以省去调用某些包时候需要现下载、安装的等待时间，也方便测试某些代码和函数。</p>
<h2 id="数据来源"><a href="#TOC"><span class="header-section-number">3.3</span> 数据来源</a></h2>
<h3 id="文本格式"><a href="#TOC"><span class="header-section-number">3.3.1</span> 文本格式</a></h3>
<p>文本为最常见的数据存储格式，包括以<code>.txt</code>、<code>.csv</code>、<code>.tsv</code>等一系列扩展名结尾的文件。文本文件可以通过windows自带的记事本或者Notepad++、UltraEdit等文本编辑软件直接打开。</p>
<p>在R中，文本文件的读取依赖<code>read.table()</code>等一类命令，使用<code>?read.table</code>可以看到，里面可以指定很多参数，其中常用的有<code>file，header， sep ，quote</code>等。</p>
<ul>
<li><p>file：需要读入的文件名</p></li>
<li><p>header: 第一列是否为变量名</p></li>
<li><p>sep: 变量之间的分隔符</p></li>
<li><p>quote: 文本被包裹的符号</p></li>
</ul>
<p>比如在当前工作目录下，我们有一个制表符（Tab或<code>\t</code>）分割的文本文件<code>sample.txt</code>，第一行含有英文变量名（中文变量名可能会出错，依系统而异），然后文本没有被任何符号包裹，那么我们读入它的时候需要采用：</p>
<pre class="sourceCode r"><code class="sourceCode r">sample_tab_data <- <span class="kw">read.table</span>(<span class="st">"data/sample_book_purschase.txt"</span>, <span class="dt">header =</span> <span class="ot">TRUE</span>,
<span class="dt">sep =</span> <span class="st">"</span><span class="ch">\t</span><span class="st">"</span>)</code></pre>
<p>而<code>read.csv()</code>、<code>read.csv2()</code>、<code>read.delim()</code>、<code>read.delim2()</code>都是<code>read.table()</code>不同默认参数的变形：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(<span class="st">"formatR"</span>)
<span class="kw">usage</span>(read.csv)</code></pre>
<pre><code>## read.csv(file, header = TRUE, sep = ",", quote = "\"",
## dec = ".", fill = TRUE, comment.char = "", ...)</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">usage</span>(read.csv2)</code></pre>
<pre><code>## read.csv2(file, header = TRUE, sep = ";", quote = "\"",
## dec = ",", fill = TRUE, comment.char = "", ...)</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">usage</span>(read.delim)</code></pre>
<pre><code>## read.delim(file, header = TRUE, sep = "\t", quote = "\"",
## dec = ".", fill = TRUE, comment.char = "", ...)</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">usage</span>(read.delim2)</code></pre>
<pre><code>## read.delim2(file, header = TRUE, sep = "\t", quote = "\"",
## dec = ",", fill = TRUE, comment.char = "", ...)</code></pre>
<p>当数据不整齐的时候，R会在读入过程中报错，并给出出错的行数。当然我们也可以通过更改参数来强制读入。<code>read.table()</code>的常用参数定义如下：</p>
<ul>
<li><p><code>header</code>：是否把第一行读为变量名</p></li>
<li><p><code>sep</code>：列与列之间的分隔符</p></li>
<li><p><code>quote</code>：引号的格式</p></li>
<li><p><code>col.names</code>：可以指定一个文本向量，作为变量的名字，其中所含名字的个数需与列数相一致</p></li>
<li><p><code>colClasses</code>：可以为每一列指定格式，如果不需要读入的话置为NULL即可。</p></li>
<li><p><code>fill</code>：是否强制每一行都有相同的列数</p></li>
</ul>
<p>其他的参数如<code>as.is</code>、<code>na.strings</code>等，可参见<code>?read.table</code>的说明。</p>
<h3 id="excel-格式通过odbc接口"><a href="#TOC"><span class="header-section-number">3.3.2</span> Excel 格式：通过ODBC接口</a></h3>
<p>Excel格式除了可以采用excel里面导出文本文件或者csv文件的方式外，还可以采取<code>ODBC</code>方式读入(Windows和Mac下)。如果采用这种方式，需要加载<em>RODBC</em>这个包。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(RODBC)
excel_channel <- <span class="kw">odbcConnectExcel</span>(<span class="st">"data/sample.xls"</span>)
sample_excel <- <span class="kw">sqlFetch</span>(excel_channel, <span class="st">"sample"</span>) <span class="co">#参数为要导入的excel数据表的名字</span>
<span class="kw">odbcClose</span>(channel)
sample_excel_2007 <- <span class="kw">odbcConnectExcel2007</span>(<span class="st">"data/sample.xlsx"</span>) <span class="co">#对于07版excel文件</span></code></pre>
<p>当然，除了excel之外，所有基于ODBC接口的数据都可以读入，包括常见的MySQL、Access等。在Linux下，MySQL数据库建议使用另外的<code>RMySQL</code>包连接。</p>
<h3 id="excel-格式通过xlsx包"><a href="#TOC"><span class="header-section-number">3.3.3</span> Excel 格式：通过xlsx包</a></h3>
<p>除了ODBC之外，另一种简单的方式则是调用<strong>xlsx</strong>包（在linux平台下也可运行）。该包不仅仅可以读入excel数据，还可以“将Excel读取为数据框，以及将数据框写入为Excel文件都不是问题，而更加强大的是它能处理Excel中的格式，比如合并单元格，设置列的宽度，设置字体和颜色等等”（具体请参见<a href="http://yixuan.cos.name/cn/2012/01/new-method-to-read-excel-file-in-r/">yixuan的介绍</a>）。这样，就为通过R脚本重复生成excel格式的报表等铺平了道路。</p>
<p>如果只是读入数据，那么直接调用<code>read.xlsx()</code>就可以了。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(xlsx)
file <- <span class="kw">system.file</span>(<span class="st">"tests"</span>, <span class="st">"test_import.xlsx"</span>, <span class="dt">package =</span> <span class="st">"xlsx"</span>)
sample_sheet1 <- <span class="kw">read.xlsx</span>(file, <span class="dv">1</span>) <span class="co"># 读取第一个工作表</span></code></pre>
<p>对于该包中<code>write.xlsx()</code>等写入excel文件函数感兴趣的读者，可以借助<code>?write.xlsx</code>来阅读详细的使用说明。注：该包只支持excel 07以上格式，对于<code>.xls</code>格式还请先用高版本的Excel等软件进行转换。</p>
<h3 id="数据库sql"><a href="#TOC"><span class="header-section-number">3.3.4</span> 数据库（SQL）</a></h3>
<p>其他数据库可以通过对应的包连接 <待补充></p>
<h3 id="其他数据格式"><a href="#TOC"><span class="header-section-number">3.3.5</span> 其他数据格式</a></h3>
<p>大多数常见的数据都可以通过<strong>foreign</strong>这个包读入：</p>
<ul>
<li><p>SPSS： <code>read.spss()</code></p></li>
<li><p>Minitab：<code>read.mtp()</code></p></li>
<li><p>SAS：<code>read.ssd()</code>或者<code>read.xport()</code></p></li>
<li><p>Stata：<code>read.dta()</code></p></li>
</ul>
<h2 id="变量操作"><a href="#TOC"><span class="header-section-number">3.4</span> 变量操作</a></h2>
<h3 id="变量重命名"><a href="#TOC"><span class="header-section-number">3.4.1</span> 变量重命名</a></h3>
<p>导入数据之后，可以在RStudio的右上方Workspace那里看到一个新的数据，名为sample，共有针对17个变量的99行记录。单击可以在左上角的数据浏览窗内看到对应的数据样本。</p>
<p>这个时候，可以用<code>names()</code>来查看变量名：</p>
<pre class="sourceCode r"><code class="sourceCode r">sample_tab_data <- <span class="kw">read.table</span>(<span class="st">"data/sample_book_purschase.txt"</span>, <span class="dt">header =</span> <span class="ot">TRUE</span>,
<span class="dt">sep =</span> <span class="st">"</span><span class="ch">\t</span><span class="st">"</span>)
<span class="kw">names</span>(sample_tab_data)</code></pre>
<pre><code>## [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" </code></pre>
<p>而后，可以对变量进行重新命名：</p>
<p>比如，我们想把第二个字段重命名为RECORD_DAY：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">names</span>(sample_tab_data)[<span class="dv">2</span>] <- <span class="st">"RECORD_DAY"</span></code></pre>
<p>这里把<code>names(sample_tab_data)</code>返回的第二个元素重新定义为了RECORD_DAY，故而实现了变量的重命名。</p>
<p>或者，我们希望对导入的没有变量名的数据集进行重命名（一般这种情况下对应的默认变量名是V1、V2等），那么可以直接对整个数据集操作。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">names</span>(sample_tab_data) <- <span class="kw">c</span>(<span class="st">"CUSTOMER"</span>, <span class="st">"RECORD_DAY"</span>, <span class="st">"LOCATION"</span>,
<span class="st">"LOCATION2"</span>, <span class="st">"FROM"</span>, <span class="st">"NOTE"</span>, <span class="st">"BOOK_ID"</span>, <span class="st">"NOTE2"</span>, <span class="st">"BOOK_TYPE"</span>, <span class="st">"CLASS_ID"</span>,
<span class="st">"CLASS_NAME"</span>, <span class="st">"AMOUNT"</span>, <span class="st">"LENGTH"</span>, <span class="st">"LENGTH2"</span>, <span class="st">"PAGES"</span>, <span class="st">"PAGES2"</span>, <span class="st">"BOOK_NAME"</span>)</code></pre>
<h2 id="data.frame的行列操作"><a href="#TOC"><span class="header-section-number">3.5</span> data.frame的行、列操作</a></h2>
<p>在一个data.frame中，我们可以直接用$来调用其中的一个变量，是最简单的调用列的格式。如果希望调用某些行列，则需要分别指定调用条件，比如：</p>
<pre class="sourceCode r"><code class="sourceCode r">sub_sample <- sample_tab_data[<span class="st">"BOOK_ID"</span> == <span class="dv">348918047</span>, <span class="kw">c</span>(<span class="st">"CUSTOMER"</span>,
<span class="st">"RECORD_DAY"</span>, <span class="st">"BOOK_ID"</span>)]</code></pre>
<p>那么sub_sample里面现在就得到了购买过编号为348918047这本书的所有顾客购买记录，包括顾客ID、购买日期和书籍ID。即，对于任何一个data.frame对象，都可以在中括号内，逗号之前指定行选择条件，逗号之后指定要选择的列（变量）。<code>c()</code>为向量生成函数。</p>
<h3 id="数据格式的转换"><a href="#TOC"><span class="header-section-number">3.5.1</span> 数据格式的转换</a></h3>
<ul>
<li>vector, data.frame , matrix, list 对于各种数据集，一般读入之后默认的是data.frame格式。此外我们常用的还有向量格式vector、矩阵格式matrix、混合格式list等。简单的说，一个data.frame的一列就是一个vector，比如我们需要所有顾客的ID这个向量：</li>
</ul>
<pre class="sourceCode r"><code class="sourceCode r">CUSTOMER_ID <- sample_tab_data$CUSTOMER</code></pre>
<pre><code>这个时候就会在workspace那里出现一个新向量CUSTOMER_ID。各个格式之间可以直接转换，比如</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">CUSTOMER_ID <- <span class="kw">as.data.frame</span>(CUSTOMER_ID)</code></pre>
<pre><code>那么这个时候customer_ID就成为了只有一个变量的data.frame格式了。</code></pre>
<ul>
<li><p>logic, character, integer, numeric, factor</p>
<p>对于向量中的元素，记录的格式则可能是逻辑型、文本型、整数型、数值型、因素型等等。各个向量格式之间可以直接转换，比如对于CLASS_ID这个变量，虽然是数值记录的但数值本身没有任何意义，只是一个相互区别和识别的代码，因此可以考虑专为character或者factor格式：</p></li>
</ul>
<pre class="sourceCode r"><code class="sourceCode r">sample_tab_data$CLASS_ID <- <span class="kw">as.factor</span>(sample_tab_data$CLASS_ID)</code></pre>
<pre><code>之后R里面就会识别其为文本或者因素型数据了。值得注意的是，如果需要把一个因素型的数据重新转换为整数型，则需要经过文本型过渡：</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">sample_tab_data$CLASS_ID <- <span class="kw">as.character</span>(sample_tab_data$CLASS_ID)
sample_tab_data$CLASS_ID <- <span class="kw">as.integer</span>(sample_tab_data$CLASS_ID)</code></pre>
<pre><code>否则会被直接重编码，丢失原有的数据串信息。（注：花开两朵、各表一枝。实际上，我们也可以应用这种特性来进行重编码工作。）</code></pre>
<h3 id="新变量生成"><a href="#TOC"><span class="header-section-number">3.5.2</span> 新变量生成</a></h3>
<ul>
<li><p>逻辑型</p>
<p>如果我们基于一些记录判断生成新的变量，比如基于如果<code>用户购买量>0</code>，则我们认为其在当日有购买行为，那么可以使用：</p></li>
</ul>
<pre class="sourceCode r"><code class="sourceCode r">sample_tab_data$purchase <- sample_tab_data$AMOUNT > <span class="dv">0</span></code></pre>
<pre><code>这样就生成了一个新的逻辑型变量purchase（取值为TRUE 或者FALSE）。逻辑型变量的一大用处就是可以直接通过相乘操作来进行多个行为之间的交集运算，比如除了是否购买之外，我们还关心购买的书籍是不是在标号为200的书店购买的，那么就可以：</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">sample_tab_data$book_store_200 <- sample_tab_data$LOCATION == <span class="dv">200</span>
sample_tab_data$purchase_bs200 <- sample_tab_data$book_store_200 *
sample_tab_data$purchase</code></pre>
<pre><code>最后得到的变量purchase_bs200为TRUE则该用户是在200号书店购买的图书。类似的，我们可以统计是不是每个月都有购买行为等。
这里还一并介绍一个有用的`%in%`运算符，表示一个元素是否属于一个给定的集合，比如：</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">sample_tab_data$book_store <- sample_tab_data$LOCATION %in% <span class="kw">c</span>(<span class="dv">200</span>,
<span class="dv">300</span>)</code></pre>
<pre><code>表示用户在200或者300号书店进行了购买。`c()`为向量生成函数，故得到的向量含有200和300两个元素。</code></pre>
<ul>
<li><p>数值型</p>
<p>数值型变量的操作一般就是基本的运算，R里面的话，</p>
<ul>
<li><p>+：加法</p></li>
<li><p>-：减法</p></li>
<li><p>*：乘法</p></li>
<li><p>/：除法</p></li>
<li><p>%% ：取余</p></li>
<li><p>^n ：n次方</p></li>
</ul></li>
<li><p>字符操作</p>
<p>字符操作最常见的就是字符串生成操作，比如我们有CUSTOMER、LOCATION、和BOOK_NAME三个变量，希望批量生成一个变量，然后发送给顾客作为反馈记录，希望的格式为“CUSTOMER顾客您好，您在编号为LOCATION的书店购买了书籍BOOK_NAME，仅供确认。”，那么我们可以使用paste()这个函数：</p></li>
</ul>
<pre class="sourceCode r"><code class="sourceCode r">sample_tab_data$message <- <span class="kw">paste</span>(sample_tab_data$CUSTOMER, <span class="st">"顾客您好，您在编号为"</span>,
sample_tab_data$LOCATION, <span class="st">"的书店购买了书籍"</span>, sample_tab_data$BOOK_NAME,
<span class="st">"，仅供确认。"</span>, <span class="dt">sep =</span> <span class="st">""</span>)</code></pre>
<pre><code>这个时候就得到的相应的新变量。`paste()`函数有个参数是`sep`，用来指定各个部分之间的连接符，默认为空格，如果不需要任何额外的符号用一对双引号设置为空即可。
字符的其他操作亦包括查找、截取`substr()`等。</code></pre>
<h3 id="数据集合并操作"><a href="#TOC"><span class="header-section-number">3.5.3</span> 数据集合并操作</a></h3>
<ul>
<li><p>列合并：对于两个含有完全一样变量的数据集，可以采用<code>rbind()</code>函数来直接将一个data.frame附加在另外一个后面。</p></li>
<li><p>行合并：如果两个数据集有相同的记录数（行数），那么可以采用<code>cbind()</code>这个函数来直接把变量附加在后面。</p></li>
<li><p>数据集合并：<code>merge()</code>函数提供了更强大的数据集合并操作命令，可以按照一个主键（即用来识别个体的变量）来合并，比如我们另有一个文件 BOOK_MAP.txt，里面记录的是重编码后的书籍ID和原编码对照表，则可以读入之后利用来<code>merge()</code>合并：</p></li>
</ul>
<pre class="sourceCode r"><code class="sourceCode r">book_map <- <span class="kw">read.delim</span>(<span class="st">"data/BOOK_MAP.txt"</span>, <span class="dt">header =</span> T)
sample_merged <- <span class="kw">merge</span>(sample_tab_data, book_map, <span class="dt">by.x =</span> <span class="st">"BOOK_ID"</span>,
<span class="dt">by.y =</span> <span class="st">"BOOK_ID"</span>, <span class="dt">all.x =</span> T, <span class="dt">all.y =</span> F)</code></pre>
<pre><code>这样就有了一列新的变量，记录的是重编码之后的书籍ID。</code></pre>
<ul>
<li><p>删除/提取变量</p>
<p>如果要删除某个变量，可以直接使用NULL值置空，即：</p></li>
</ul>
<pre class="sourceCode r"><code class="sourceCode r">sample_merged$CLASS_NAME <- <span class="ot">NULL</span></code></pre>
<pre><code>会删除掉` CLASS_NAME`这个变量。在需要删除多个变量的时候，不如只保留几个变量，如：</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">sample_merged <- sample_merged[, <span class="kw">c</span>(<span class="st">"BOOK_ID"</span>, <span class="st">"CLASS_ID"</span>)]</code></pre>
<pre><code>会只保留`BOOK_ID`, `CLASS_ID` 两个变量。</code></pre>
<h3 id="数据集形状的转换reshape2包"><a href="#TOC"><span class="header-section-number">3.5.4</span> 数据集形状的转换（reshape2包）</a></h3>
<p>除了这些基本的变量操作之外，还有一类可能的需求就是对于整个数据集做一个形状的转换，比如把“长数据集”转换为“宽数据集”。这样的过程类似于“揉面”，而帮我们玩转面团的利器便是<strong>reshape2</strong>这个包。</p>
<p>在正式介绍强大的<strong>reshape2</strong>包之前，需要先提到一个轻量级武器——<code>reshape()</code>函数。这个函数可以帮我们在数据的长、宽形状之间自由玩转，比如我们现在有一个用户逐月购买记录，为长格式，想把它变为宽格式：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">load</span>(<span class="st">"data/reshape_sample.rdata"</span>) <span class="co">#载入样本数据集</span>
<span class="kw">summary</span>(reshape_sample) <span class="co">#基本统计量，有CUSTOMER_ID（顾客ID）、MONTH（月份）、PURCHASE（消费额）三个变量</span></code></pre>
<pre><code>## 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 </code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">head</span>(reshape_sample) <span class="co">#调用数据的前几行，显示为长格式</span></code></pre>
<pre><code>## 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 4</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">reshape_sample_wide <- <span class="kw">reshape</span>(reshape_sample, <span class="dt">idvar =</span> <span class="st">"CUSTOMER_ID"</span>,
<span class="dt">timevar =</span> <span class="st">"MONTH"</span>, <span class="dt">direction =</span> <span class="st">"wide"</span>) <span class="co">#变为宽格式</span>
<span class="kw">head</span>(reshape_sample_wide) <span class="co">#宽格式展现</span></code></pre>
<pre><code>## 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 110981</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">summary</span>(reshape_sample_wide) <span class="co">#宽格式下基本统计量</span></code></pre>
<pre><code>## 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
## </code></pre>
<p>在<strong>reshape2</strong>包中，<code>melt()</code>函数进一步简化了这个过程。比如我们现在希望把宽数据转回长数据：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(reshape2)
reshape_sample_long <- <span class="kw">melt</span>(reshape_sample_wide, <span class="dt">id =</span> <span class="kw">c</span>(<span class="st">"CUSTOMER_ID"</span>)) <span class="co">#转回长数据格式</span>
<span class="kw">head</span>(reshape_sample_long[<span class="kw">order</span>(reshape_sample_long$CUSTOMER_ID),
]) <span class="co">#长格式展示</span></code></pre>
<pre><code>## 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</code></pre>
<p>此外，该包提供的<code>acast()</code>/<code>dcast()</code>函数可以进一步帮我们分类展现数据及其统计量，具体使用请参见函数包内帮助。</p>
<h2 id="数据的导出"><a href="#TOC"><span class="header-section-number">3.6</span> 数据的导出</a></h2>
<p>数据导出最常用的应该就是<code>write.table</code>函数。比如我们要输出book_map这个数据集为文本格式，那么使用：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">write.table</span>(book_map, <span class="dt">file =</span> <span class="st">"book_map_new.txt"</span>, <span class="dt">row.names =</span> F, <span class="dt">col.names =</span> T,
<span class="dt">sep =</span> <span class="st">"</span><span class="ch">\t</span><span class="st">"</span>, <span class="dt">quote =</span> F)</code></pre>
<p>这个时候就会得到book_map_new.txt这个文本文件，以制表符分隔。<code>write.table()</code>的主要参数有：</p>
<ul>
<li><p>file：指定文件名</p></li>
<li><p>row.names：是否输出记录行编号</p></li>
<li><p>col.names：是否输出变量名</p></li>
<li><p>sep：分隔符</p></li>
<li><p>quote：引号形式</p></li>
<li><p>append：是否附加在现有文件后面（如为FALSE则新文件覆盖原有文件）</p></li>
</ul>
<h2 id="分类统计data.table包"><a href="#TOC"><span class="header-section-number">3.7</span> 分类统计（<strong>data.table</strong>包）</a></h2>
<p>日常分析工作中，最常用到的就是分类统计了。简单的说，就是按一个字段归类之后，统计其他字段的量。</p>
<p>比如，我们现在希望知道每个顾客购买的图书的总数。那么可以使用如下的代码：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(data.table)
sample_tab_data_frame <- <span class="kw">data.table</span>(sample_tab_data)
sample_stats <- sample_tab_data_frame[, <span class="kw">list</span>(<span class="dt">Amount =</span> <span class="kw">sum</span>(AMOUNT),
<span class="dt">Book =</span> <span class="kw">length</span>(<span class="kw">unique</span>(BOOK_ID))), by = <span class="st">"CUSTOMER"</span>]
sample_stats</code></pre>
<pre><code>## 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</code></pre>
<p>这样在新的对象sample_stats里面，就有了每个顾客购买书籍的总本书、以及不同图书的数量。简单的说，需要使用这样的分类统计的时候，就是先加载<code>data.table</code>这个包，然后利用<code>data.table()</code>函数转换一个data.frame为data.table格式，然后在原有的data.frame行、列操作基础上，增加了一个by参数，可以用来指定分类的依据。当然这里我们可以同时针对多个变量及其组合分类统计，比如</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(data.table)
sample_stat_by_day <- sample_tab_data_frame[, <span class="kw">list</span>(<span class="dt">Amount =</span> <span class="kw">sum</span>(AMOUNT),
<span class="dt">Book =</span> <span class="kw">length</span>(<span class="kw">unique</span>(BOOK_ID))), by = <span class="kw">c</span>(<span class="st">"CUSTOMER"</span>, <span class="st">"RECORD_DAY"</span>)]
sample_stat_by_day</code></pre>
<pre><code>## 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</code></pre>
<p>这样返回的就是逐日统计的每位顾客的购买数量了。值得多说的是，<code>sum()</code>函数代表求和，<code>length()</code>函数代表计数，而<code>unique()</code>函数则是去掉重复值。</p>
<h2 id="循环和判断"><a href="#TOC"><span class="header-section-number">3.8</span> 循环和判断</a></h2>
<p>我们常用的循环和判断有三种：</p>
<ul>
<li><p>for</p></li>
<li><p>while</p></li>
<li><p>if</p></li>
</ul>
<p>循环和判断是基本的逻辑操作语句。在R中，大部分常用功能都已经有现成的函数，所以极少用到循环，我们也非常不提倡在R里面写循环（效率一般很低，因为一个循环背后往往是现有函数内部的许多层循环）。有的时候，为了一些特殊的情况，知道怎么写循环还是有用的。比如，我们希望把统计好的顾客购买记录按照顾客ID写入不同的文本文件，这里就需要用到for循环或者while循环。</p>
<pre class="sourceCode r"><code class="sourceCode r">for (i in <span class="dv">1</span>:<span class="kw">nrow</span>(sample_stats)) {
file_name <- <span class="kw">paste</span>(<span class="st">"data/result_"</span>, sample_stats[i, ]$CUSTOMER, <span class="st">"_record.txt"</span>,
<span class="dt">sep =</span> <span class="st">""</span>)
<span class="kw">write.table</span>(sample_stats[i, ], <span class="dt">file =</span> file_name, <span class="dt">sep =</span> <span class="st">"</span><span class="ch">\t</span><span class="st">"</span>, <span class="dt">row.names =</span> F,
<span class="dt">col.names =</span> T, <span class="dt">quote =</span> F)
}</code></pre>
<p>这里用到了for循环。另，<code>nrow()</code>是用来统计一个data.frame行数的函数，同样的<code>ncol()</code>会返回其列数。符号:会生成一个向量，从1到<code>nrow(sample_stat)</code>。常用的循环还有while，同时可以使用if来进行一些判断。这里不再赘述。</p>
<h1 id="从截面数据分析说起"><a href="#TOC"><span class="header-section-number">4</span> 从截面数据分析说起</a></h1>
<p>上面简单的说了一下多元回归，下面则是一些我们在回归分析中常用分析的实现。</p>
<h2 id="参数检验"><a href="#TOC"><span class="header-section-number">4.1</span> 参数检验</a></h2>
<p>得到一个回归方程后，关心的第一件事就是系数和方程整体的显著性，分别由t检验和F检验实现。来看下面这个有关法学院的例子。</p>
<p>(@LAWSCH85) 在LAWSCH85.DTA这个数据集中，法学院应届生薪水的中位数由下面的方程决定：</p>
<p>\(log(salary)=\beta_{0}+\beta_{1}LAST+\beta_{2}GPA+\beta_{3}log(libvol)+\beta_{4}log(cost)+\beta_{5}rank+u\)</p>
<p>其中LSAT 是班级里LSAT成绩中位数，GPA 是班级学习成绩的中位数，libvol 是法学院图书馆藏书卷数，cost 是在法学院的年消费额，rank 是法学院的排名（正序）。 接下来我们估计这个方程:</p>
<p><code>LAW_Result <- lm(log(salary) ~ LSAT + GPA + log(libvol) + log(cost) + rank, data=LAW)</code></p>
<p>很容易得到回归结果如下：</p>
<pre class="sourceCode r"><code class="sourceCode r">LAW <- <span class="kw">read.dta</span>(<span class="st">"data/LAWSCH85.DTA"</span>)</code></pre>
<pre><code>## Warning message: cannot read factor labels from Stata 5 files</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">LAW_Result <- <span class="kw">lm</span>(<span class="kw">log</span>(salary) ~ LSAT + GPA + <span class="kw">log</span>(libvol) + <span class="kw">log</span>(cost) +
rank, <span class="dt">data =</span> LAW)
<span class="kw">summary</span>(LAW_Result)</code></pre>
<pre><code>##
## 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
## </code></pre>
<h3 id="单变量显著性t检验"><a href="#TOC"><span class="header-section-number">4.1.1</span> 单变量显著性：t检验</a></h3>
<p>在回归结果中已经报告了各变量的t统计值，从而可知：rank 的估计值很显著（通过0.1%显著性水平检验）。而GPA 和(log)libvol 则通过了1%显著性水平检验。 而变量LSAT 系数估计值不显著。</p>
<h3 id="多个变量显著性f检验"><a href="#TOC"><span class="header-section-number">4.1.2</span> 多个变量显著性：F检验</a></h3>
<p>考虑到新生入学的时候只有GPA 和LSAT 两个变量可以观测，所以接下来我们进行变量GPA 和LSAT 的联合检验即F检验。该检验属于线性假设检验，在<em>RCommander</em>下可以在<code>Models -> Hypothesis Tests ->Linear hypothesis</code>里面通过图形化界面完成。</p>
<p><img src="imgs/f_test.JPG" alt="F检验图形化界面" /> 其中设置为“2”行，而后分别在LAST和GPA输入1，保持右侧为0。相当于检验假设</p>
<p>\(H_{0}:\;\beta_{1}=\beta_{2}=0\)</p>
<p>可得到输出结果如下：</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(<span class="st">"car"</span>)
.Hypothesis <- <span class="kw">matrix</span>(<span class="kw">c</span>(<span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>), <span class="dv">2</span>, <span class="dv">6</span>,
<span class="dt">byrow =</span> <span class="ot">TRUE</span>)
.RHS <- <span class="kw">c</span>(<span class="dv">0</span>, <span class="dv">0</span>)
<span class="kw">linearHypothesis</span>(LAW_Result, .Hypothesis, <span class="dt">rhs =</span> .RHS)</code></pre>
<pre><code>## 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 </code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">remove</span>(.Hypothesis, .RHS)</code></pre>
<p>从上面结果可知，F=9.9517 且通过了0.1%显著性水平检验，即变量LSAT 和GPA 联合显著。 此外还可以分别加入clsize （班级容量）和faculty （教师规模）来进行回归，代码如下： [未完成]</p>
<p>由回归结果，这两个变量的系数均不显著，且方程F统计量有所下降（调整后的\(R^{2}\) 变动不大），故不用加入到方程中关于模型变量选择的问题，将在后面另行论述，这里只做简单的分析。</p>
<h2 id="置信区间"><a href="#TOC"><span class="header-section-number">4.2</span> 置信区间</a></h2>
<p>回归分析里还常用的一项分析就是得到某一显著性水平下的置信区间。[未完成]</p>
<h2 id="虚拟变量"><a href="#TOC"><span class="header-section-number">4.3</span> 虚拟变量</a></h2>
<h3 id="按性质分组"><a href="#TOC"><span class="header-section-number">4.3.1</span> 按性质分组</a></h3>
<p>比较简单的例子就是已经含有分组变量的数据，比如变量gender 有两个值male, female, 那么我们只需把它们变成factor形式就可以了。如我们可以把法学院例子中的变量north进行变化 [footnote-factor][]： [footnote-factor]: 其实这里不用这么麻烦也可以，因为north变量本身的赋值只有0和1，可以直接进行回归。在这里只是用这个例子来说明<code>factor()</code>函数的调用形式而已。 “footnote-factor”</p>
<pre class="sourceCode r"><code class="sourceCode r">LAW$north_true <- <span class="kw">factor</span>(LAW$north, <span class="dt">labels =</span> <span class="kw">c</span>(<span class="st">"others"</span>, <span class="st">"north"</span>))</code></pre>
<p>可以这样做的原因是：R在调用<code>lm()</code>函数的时候会自动把factor类型的变量作为虚拟变量进行回归。</p>
<h3 id="按数量值分组"><a href="#TOC"><span class="header-section-number">4.3.2</span> 按数量值分组</a></h3>
<p>依旧采用法学院的例子。很明显在上面的分析中我们把各个学院排名直接当作一个可以“测距”的变量其实它只是一个排序而已，并不代表实质的差距。来使用了，这可能会引起一些争议。因此，我们可以采取另一种模式，即引入虚拟变量，把学校分为六组：top 10，11-25名，26-40名，41-60名，61-100名，100名开外。引入五个虚拟变量\(top10\) , \(r11\_25\) , \(r26\_40\) , \(r41\_60\) , \(r61\_100\) 。 我想你不会手动的把所有的学校都赋一个虚拟变量值吧？在R里面，我们需要先通过<code>recode()</code>来依照分组创建一个新的factor形式的变量\(rank\_f\)，然后再进行回归。这样我们就不需要在原来的数据库里面新增加五个变量并赋逻辑值了。</p>
<pre class="sourceCode r"><code class="sourceCode r">LAW$rank_f <- <span class="kw">recode</span>(LAW$rank, <span class="st">"1:10=</span><span class="ch">\"</span><span class="st">top10</span><span class="ch">\"</span><span class="st">; 11:25=</span><span class="ch">\"</span><span class="st">r11_25</span><span class="ch">\"</span><span class="st">; 26:40=</span><span class="ch">\"</span><span class="st">r26_40</span><span class="ch">\"</span><span class="st">; 41:60=</span><span class="ch">\"</span><span class="st">r41_60</span><span class="ch">\"</span><span class="st">; 61:100=</span><span class="ch">\"</span><span class="st">r61_100</span><span class="ch">\"</span><span class="st">; else=</span><span class="ch">\"</span><span class="st">r101_</span><span class="ch">\"</span><span class="st">; "</span>,
<span class="dt">as.factor.result =</span> <span class="ot">TRUE</span>)
LAW_Result2 <- <span class="kw">lm</span>(<span class="kw">log</span>(salary) ~ LSAT + GPA + <span class="kw">log</span>(libvol) + <span class="kw">log</span>(cost) +
rank_f, <span class="dt">data =</span> LAW)
<span class="kw">summary</span>(LAW_Result2)</code></pre>
<pre><code>##
## 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
## </code></pre>
<p>此外，我们可以通过<em>RCommander</em>里面，在<code>Data -> Recode Variables</code>的方框里逐行输入。</p>
<div class="figure">
<img src="imgs/recode_variables.JPG" alt="Recode Varibles图形化界面" /><p class="caption">Recode Varibles图形化界面</p>
</div>
<h3 id="交叉项intersections"><a href="#TOC"><span class="header-section-number">4.3.3</span> 交叉项（intersections）</a></h3>
<p>我不知道这样叫是不是足够确切，在微观计量里面我们会经常用到两个变量相乘的回归项，比如\(female\cdot single\)，即单身女士。相比而言这样的虚拟变量并不需要特别的处理，在回归方程里面直接写成相乘的形式即可。注意此时不需要再写female 和single 变量，<code>lm()</code>会默认加入这两个变量。例如，在法学院的例子中，我们可以对top10 和west 两个变量进行相乘回归（这里top10 变量来源于数据库本身自带的）。</p>
<pre class="sourceCode r"><code class="sourceCode r">LAW_Result3 <- <span class="kw">lm</span>(<span class="kw">log</span>(salary) ~ LSAT + GPA + <span class="kw">log</span>(libvol) + <span class="kw">log</span>(cost) +
top10 * west, <span class="dt">data =</span> LAW)
<span class="kw">summary</span>(LAW_Result3)</code></pre>
<pre><code>##
## 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
## </code></pre>
<p>当然也可以使用factor变量和其他虚拟变量相乘进行回归，反馈的结果中包含所有的相乘项。</p>
<pre class="sourceCode r"><code class="sourceCode r">LAW_Result4 <- <span class="kw">lm</span>(<span class="kw">log</span>(salary) ~ LSAT + GPA + <span class="kw">log</span>(libvol) + <span class="kw">log</span>(cost) +
rank_f * west, <span class="dt">data =</span> LAW)
LAW_Result4</code></pre>
<pre><code>##
## 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
## </code></pre>
<h3 id="指定参照组"><a href="#TOC"><span class="header-section-number">4.3.4</span> 指定参照组</a></h3>
<p>当R处理factor形式的数据的时候，默认以数据中的第一个层次 (level) 作为参照组。比如上例中，我们想把top10 这一组作为参照组，那么则需要使用<code>relevel()</code>命令。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">attach</span>(LAW)
rank_f2 <- <span class="kw">relevel</span>(rank_f, <span class="dt">ref =</span> <span class="st">"top10"</span>)
LAW_Result5 <- <span class="kw">lm</span>(<span class="kw">log</span>(salary) ~ LSAT + GPA + <span class="kw">log</span>(libvol) + <span class="kw">log</span>(cost) +
rank_f2, <span class="dt">data =</span> LAW)
<span class="kw">summary</span>(LAW_Result5)</code></pre>
<pre><code>##
## 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
## </code></pre>
<h2 id="异方差检验"><a href="#TOC"><span class="header-section-number">4.4</span> 异方差检验</a></h2>
<p>有件很没办法的事儿，那就是要想让OLS回归出来的结果最佳，必须要符合那五条经典假设（尤其是小样本下）。但是事实中的数据那里会那么完美呢？首当其冲的问题就是异方差。 异方差检验方法很多，这里给出两种常用的：BP检验 (Breusch-Pagan Test)、怀特检验 (White test for heteroskedasticity)。</p>
<p>下面给出一个关于房价的例子。</p>
<p>(@Hprice1) 在Hprice1.dta这个数据集中，有price （房价，按套计算）、lotsize （地皮面积要知道人家住的都是小别墅啊，自然要先有地、后建房，卖房子也都是一套一套的卖。）、sqrft （房屋面积）、bdrms （卧室数量）。接下来我们估计这个方程：</p>
<p>\(\hat{price}=\hat{\beta_{0}}+\hat{\beta_{1}}lotsize+\hat{\beta_{2}}sqrft+\hat{\beta_{3}}bdrms\)</p>
<p>使用OLS估计结果如下：</p>
<pre class="sourceCode r"><code class="sourceCode r">HPrice <- <span class="kw">read.dta</span>(<span class="st">"data/hprice1.dta"</span>)
Hprice_Result <- <span class="kw">lm</span>(price ~ bdrms + lotsize + sqrft, <span class="dt">data =</span> HPrice)
Hprice_Result</code></pre>
<pre><code>##
## Call:
## lm(formula = price ~ bdrms + lotsize + sqrft, data = HPrice)
##
## Coefficients:
## (Intercept) bdrms lotsize sqrft
## -21.77031 13.85252 0.00207 0.12278
## </code></pre>
<h3 id="bp检验-breusch-pagan-test"><a href="#TOC"><span class="header-section-number">4.4.1</span> BP检验 (Breusch-Pagan Test)</a></h3>
<p>下面我们进行BP检验来测定是否有异方差。进行BP检验需要加载包<em>lmtest</em>，而后者需要加载包<em>zoo</em>。调用BP检验最简单的方法就是直接写回归结果变量。更多参数关于该命令各种参数设置请使用<code>?bptest</code>来查看帮助文档。可以通过<em>RCommander</em>里面的图形化界面设定，位于<code>Models > Numberical diagnostics > Breusch-Pagan Test for heteroskedasticity</code>。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(<span class="st">"lmtest"</span>)
<span class="kw">bptest</span>(Hprice_Result)</code></pre>
<pre><code>##
## studentized Breusch-Pagan test
##
## data: Hprice_Result
## BP = 14.09, df = 3, p-value = 0.002782
## </code></pre>
<p>由结果来看，存在异方差。由于对数形式是消除异方差（尤其针对价格数据）的常用方法，因此我们再对对数形式进行回归:</p>
<p>\(\hat{log(price)}=\hat{\beta_{0}}+\hat{\beta_{1}}log(lotsize)+\hat{\beta_{2}}log(sqrft)+\hat{\beta_{3}}bdrms\)</p>
<p>得到回归结果如下：</p>
<pre class="sourceCode r"><code class="sourceCode r">Hprice_Result2 <- <span class="kw">lm</span>(<span class="kw">log</span>(price) ~ bdrms + <span class="kw">log</span>(lotsize) + <span class="kw">log</span>(sqrft),
<span class="dt">data =</span> HPrice)
Hprice_Result2</code></pre>
<pre><code>##
## 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
## </code></pre>
<p>此时再进行异方差检验。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">bptest</span>(Hprice_Result2)</code></pre>
<pre><code>##
## studentized Breusch-Pagan test
##
## data: Hprice_Result2
## BP = 4.223, df = 3, p-value = 0.2383
## </code></pre>
<p>因此可以接受原假设，已经不存在异方差。 此外还有一个<code>bptest()</code>非常接近的<code>ncv.test()</code>，需加载<strong>car</strong>*包。我们进行一下对比。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(car)
<span class="kw">ncvTest</span>(Hprice_Result2)</code></pre>
<pre><code>## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 3.525 Df = 1 p = 0.06044 </code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">bptest</span>(Hprice_Result2, <span class="dt">studentize =</span> <span class="ot">FALSE</span>, <span class="dt">varformula =</span> ~<span class="kw">fitted.values</span>(Hprice_Result2),
<span class="dt">data =</span> HPrice)</code></pre>
<pre><code>##
## Breusch-Pagan test
##
## data: Hprice_Result2
## BP = 3.525, df = 1, p-value = 0.06044
## </code></pre>
<p>可以看到<code>ncvTest()</code>相当于把<code>bptest()</code>里面参数studentize设为FLASE，且进行固定值的检验。事实上，<code>bptest()</code>默认该项为TRUE，在TRUE时会采用Koenker (1981)的学生氏检验算法，也是目前最广为接受的算法。</p>
<h3 id="怀特检验-white-test-for-heteroskedasticity"><a href="#TOC"><span class="header-section-number">4.4.2</span> 怀特检验 (White test for heteroskedasticity)</a></h3>
<p>还是上面这个例子，我们改用怀特检验。其实，我们可以把怀特检验看作广义上的BP检验的一种特殊形式，因此可以通过在<code>bptest()</code>里面赋予更多的参数来实现，即加入各个变量平方和交叉相乘的项。</p>
<p>比如回归为<code>fm <- lm(y ~ x + z, data = foo)</code>，那么则应写成<code>bptest(fm, ~ x * z + I(x^2) + I(z^2), data = foo)</code> 。其余以此类推，不再举例。</p>
<p>另，也可以通过<strong>sandwich</strong>这个专门对付“三明治”的包来实现怀特检验。</p>
<h2 id="稳健标准差"><a href="#TOC"><span class="header-section-number">4.5</span> 稳健标准差</a></h2>
<p>当存在异方差的时候，我们可以采用稳健的标准差来替代原有的异方差矩阵。在R中，可以利用<strong>sandwich</strong>包的<code>vcovHC()</code>函数实现。依旧采用上面房价的例子。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(<span class="st">"sandwich"</span>)
<span class="kw">vcovHC</span>(Hprice_Result)</code></pre>
<pre><code>## (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.0016591</code></pre>
<p><code>vcovHC()</code>函数可以选择各种形式，在后面附加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}\) ，更多解释请参照<code>?vcovHC</code>。 之后我们就可以采用该异方差矩阵进行参数检验，如t检验。这里调用<em>lmtest</em>包内另一个非常有用的函数<code>coeftest()</code>。</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">coeftest</span>(Hprice_Result, <span class="dt">vcov =</span> vcovHC)</code></pre>
<pre><code>##
## 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
## </code></pre>
<p>此外，对于含有自回归异方差的情形，可以采用<code>vcovHAC()</code>函数。更多相关用法后面章节中详述。这两个函数均可以附加到参数检验<code>coeftest()</code>或者<code>waldtest()</code>一起使用。</p>
<h2 id="加权最小二乘估计-wls"><a href="#TOC"><span class="header-section-number">4.6</span> 加权最小二乘估计 (WLS)</a></h2>
<h3 id="扰动项形式已知"><a href="#TOC"><span class="header-section-number">4.6.1</span> 扰动项形式已知<需添加例子></a></h3>
<p>有些情况下，我们可以写出加权的形式，比如扰动项服从\(Var(u_{i}|inc)=\sigma^{2}inc\) ，那么可以直接在<code>lm()</code>函数里附加一项weight来实现。</p>
<需添加例子>
<h3 id="扰动项形式未知feasible-linear-regression"><a href="#TOC"><span class="header-section-number">4.6.2</span> 扰动项形式未知（Feasible Linear Regression）</a></h3>
<p>(@smoke) 下面是一个烟草需求的例子。在SMOKE.rda中有如下几个变量：每天吸烟的数量 (cigs)、年收入 (income)、该州烟的价格 (cigpric)、受访者年龄 (age)、受教育程度 (educ)、该州有无饭店内吸烟禁令 (restaurn)。而后我们需要研究决定烟草需求的因素，即cigs为被解释变量，其他为解释变量。</p>