Categories
日常应用

七天搞定SAS(一):数据的导入、数据结构

标题有些噱头,不过这里的重点是: speak SAS in 7 days。也就是说,知识是现成的,我这里只是要学会如何讲这门语言,而不是如何边学SAS边学模型。顺便发现我最近喜欢写连载了,自从西藏回来后.....

之所以下定决定学SAS,是因为周围的人都在用SAS。为了和同事的沟通更有效率,还是多学一门语言吧。R再灵活,毕竟还是只有少数人能直接读懂。理论上语言是不应该成为障碍的~就像外语一样,多学一点总是好的,至少出门不发怵是不是?

最后一根稻草则是施老师传给我的一个link:http://blog.softwareadvice.com/articles/bi/3-career-secrets-for-data-scientists-1101712/,据说有数据分析师的职业秘笈...我就忍不住去看了看。其中一句话还是蛮有启发的:

如果有人问你要学什么工具,是SAS,R,EXCEL,SQL,SPSS还是?直接回答:所有。

这个答案一方面霸气,一方面也是,何必被工具束缚呢?

这东西宜突击不宜拖延,所以还是集中搞定吧。七天应该是个不错的时间段。

大致分配如下:
1. 熟悉SAS的数据结构,如基本的向量,数据集,数组;熟悉基本的数据类型,如文本,数字。
2. 熟悉基本的数据输入与输出。
3. 熟悉基本的逻辑语句:循环,判断
4. 熟悉基本的数据操作:筛选行列,筛选或计算变量,合并数据集,计算基本统计量,转置
5. 熟悉基本的文本操作函数
6. 熟悉基本的计量模型函数
7. 熟悉基本的macro编写,局部变量与全局变量

其实这大概也是按照我常用的R里面完成的任务来罗列的。基本计划是完成就可以大致了解SAS的语法了,其他的高级功能现用现学吧。

书籍方面,中文的抢了同事的一本《SAS编程与数据挖掘商业案例》,英文的找了一本「Applied Econometrics Using The SAS System」和「The Little SAS Book」,先这么看着吧。

后知后觉的补充:其实这一系列笔记都是先写再发布的,主要是方便我调整顺序什么的。事实证明绝大多数时间我在看(或者更直接的,抄)「The Little SAS Book」这本书,姚老师的《SAS编程与数据挖掘商业案例》简单看了一晚,作为对于SAS语法的预热。最后那本「Applied Econometrics Using The SAS System」更多是看具体模型的用法了,不是熟悉语法的问题了。例子都是第一本little book上的,很好用。

本系列连载文章:

-------笔记开始-------

SAS的数据类型

2013-12-09 16_41_10-The Little SAS Book(Fourth).PDF - Adobe Reader

首先,sas的编程大概就两块:Data和PROC,这个倒是蛮清晰的划分。然后目前关注data部分。

SAS的数据类型还真的只有两种:数字和文本。那么看来日期就要存成文本型了。变量名称后面加$代表文本型。

SAS的数据读入

手动输入这种就不考虑了,先是怎么从本地文件读入。比如我们有文本文件如下:

Lucky 2.3 1.9 . 3.0
Spot 4.6 2.5 3.1 .5
Tubs 7.1 . . 3.8
Hop 4.5 3.2 1.9 2.6
Noisy 3.8 1.3 1.8 1.5
Winner 5.7 . . .

然后SAS里面就可以用

* Create a SAS data set named toads;
* Read the data file ToadJump.dat using list input;
DATA toads;
INFILE ’c:\MyRawData\ToadJump.dat’;
INPUT ToadName $ Weight Jump1 Jump2 Jump3;
RUN;
* Print the data to make sure the file was read correctly;
PROC PRINT DATA = toads;
TITLE ’SAS Data Set Toads’;
RUN;

这样就建立了一个名为toads的临时数据集,然后读入外部文件ToadJump.dat,然后告诉SAS有四个变量,其中第一个是文本型。这样就OK了。缺失值用一个点.标记。

偶尔数据没那么规范,比如长成:

----+----1----+----2----+----3----+----4
Columbia Peaches 35 67 1 10 2 1
Plains Peanuts 210 2 5 0 2
Gilroy Garlics 151035 12 11 7 6
Sacramento Tomatoes 124 85 15 4 9 1

那么就要有点类似正则表达式的感觉,告诉SAS更多的参数:

* Create a SAS data set named sales;
* Read the data file OnionRing.dat using column input;
DATA sales;
INFILE ’c:\MyRawData\OnionRing.dat’;
INPUT VisitingTeam $ 1-20 ConcessionSales 21-24 BleacherSales 25-28
OurHits 29-31 TheirHits 32-34 OurRuns 35-37 TheirRuns 38-40;
RUN;
* Print the data to make sure the file was read correctly;
PROC PRINT DATA = sales;
TITLE ’SAS Data Set Sales’;
RUN;

这样SAS就可以正确的读数据了—类似于excel的导入文本-固定宽度分隔。

再不规则的话,比如有日期型的:

Alicia Grossman 13 c 10-28-2008 7.8 6.5 7.2 8.0 7.9
Matthew Lee 9 D 10-30-2008 6.5 5.9 6.8 6.0 8.1
Elizabeth Garcia 10 C 10-29-2008 8.9 7.9 8.5 9.0 8.8
Lori Newcombe 6 D 10-30-2008 6.7 5.6 4.9 5.2 6.1
Jose Martinez 7 d 10-31-2008 8.9 9.510.0 9.7 9.0
Brian Williams 11 C 10-29-2008 7.8 8.4 8.5 7.9 8.0

那么接下来就是:

* Create a SAS data set named contest;
* Read the file Pumpkin.dat using formatted input;
DATA contest;
INFILE ’c:\MyRawData\Pumpkin.dat’;
INPUT Name $16. Age 3. +1 Type $1. +1 Date MMDDYY10.
(Score1 Score2 Score3 Score4 Score5) (4.1);
RUN;
* Print the data set to make sure the file was read correctly;
PROC PRINT DATA = contest;
TITLE ’Pumpkin Carving Contest’;
RUN;

就是说,name是一个长度为16的字符;age是长度为3、无小数点的数字;+1跳过空列;type是长度为1的文本;date是MMDDYY长度为10的日期;score1-5是长度为4,小数部分为1位的数字。

还有若干更复杂的,可以遇到时侯回来查手册。此外还有@可用来直接指定开始读的列。鉴于我接触的数据一般比较规范,这些就不细看了。

此外SAS可以指定开始读的行数,读取的行数等。

DATA icecream;
INFILE ’c:\MyRawData\IceCreamSales.dat’ FIRSTOBS = 3;
INPUT Flavor $ 1-9 Location BoxesSold;
RUN;

SAS读取CSV数据

以我最关心的CSV文件为例,如下数据:

Lupine Lights,12/3/2007,45,63,70,
Awesome Octaves,12/15/2007,17,28,44,12
"Stop, Drop, and Rock-N-Roll",1/5/2008,34,62,77,91
The Silveyville Jazz Quartet,1/18/2008,38,30,42,43
Catalina Converts,1/31/2008,56,,65,34

只需要:

DATA music;
INFILE ’c:\MyRawData\Bands.csv’ DLM = ’,’ DSD MISSOVER;
INPUT BandName :$30. GigDate :MMDDYY10. EightPM NinePM TenPM ElevenPM;
RUN;
PROC PRINT DATA = music;
TITLE ’Customers at Each Gig’;
RUN;

其实,貌似更简单的办法是:

DATA music;
INFILE ’c:\MyRawData\Bands.csv’ DLM = ’,’ DSD MISSOVER;
INPUT BandName :$30. GigDate :MMDDYY10. EightPM NinePM TenPM ElevenPM;
RUN;
PROC PRINT DATA = music;
TITLE ’Customers at Each Gig’;
RUN;

好吧,import果然更直接一点...excel文件也可以如法炮制。

SAS读取excel数据

* Read an Excel spreadsheet using PROC IMPORT;
PROC IMPORT DATAFILE = 'c:\MyExcelFiles\OnionRing.xls' DBMS=XLS OUT = sales;
RUN;
PROC PRINT DATA = sales;
TITLE 'SAS Data Set Read From Excel File';
RUN;

如果需要SAS永久存着这些数据,则需要先指定libname:

LIBNAME plants ’c:\MySASLib’;
DATA plants.magnolia;
INFILE ’c:\MyRawData\Mag.dat’;
INPUT ScientificName $ 1-14 CommonName $ 16-32 MaximumHeight
AgeBloom Type $ Color $;
RUN;

后期就可以直接调用啦:

LIBNAME example ’c:\MySASLib’;
PROC PRINT DATA = example.magnolia;
TITLE ’Magnolias’;
RUN;

SAS 读取Teradata数据

最后就是从teradata里面读数据,可以利用teradata fastexport特性:

libname tra Teradata user=terauser pw=XXXXXX server=boom;
proc freq data=tra.big(dbsliceparm=all);
table x1-x3;
run;

等价于:

proc sql;
connect to teradata(user=terauser password=XXXXXX server=boom dbsliceparm=all);
select * from connection to teradata
(select * from big);
quit;

暂时没有fastload的需求,就先这样吧。可以参见SAS的TD手册:http://support.sas.com/resources/papers/teradata.pdf

本系列连载文章:

 

Categories
网络新发现

Coursera上的R语言课程

今天登上Coursera一看,随便点开几门课居然都是用R来辅助的...R是什么时候悄悄的渗透到这么多大学和行业的哇?孤陋寡闻了呢。

入门的,如专门的数据分析计算,有一门Computing for Data Analysis,是时长为4节的R语言课程。讲的貌似比较基础:

This course is about learning the fundamental computing skills necessary for effective data analysis. You will learn to program in R and to use R for reading data, writing functions, making informative graphs, and applying modern statistical methods.

还有一门类似的,Data Analysis(居然是Johns Hopkins的生物统计研究生院一年级的课程):

This course will focus on how to plan, carry out, and communicate analyses of real data sets. While we will cover the basics of how to use R to implement these analyses, the course will not cover specific programming skills. Computing for Data Analysis will cover some statistical programming topics that will be useful for this class, but it is not a prerequisite for the course.

当然,基础的统计课程也是R的天下了:Statistics One

Statistics One also provides an introduction to the R programming language. All the examples and assignments will involve writing code in R and interpreting R output. R software is free! It is also an open source programming language. What this means is you can download R, take this course, and start programming in R after just a few lectures.

经济与计量、金融计算自然也不能免俗...Introduction to Computational Finance and Financial Econometrics

Learn mathematical and statistical tools and techniques used in quantitative and computational finance. Use the open source R statistical programming language to analyze financial data, estimate statistical models, and construct optimized portfolios.

自然还有类似的金融课程:Financial Engineering and Risk Management

With regards to programming, we have designed the course so that all required "programming" questions can be completed within Excel. However some questions may be easier to complete using Matlab, R, Python etc.

然后居然还看到社会网络分析也是用R来辅助的:Social Network Analysis,这里有我最喜欢的Gephi和R,咿呀呀,不奇怪的嘛,好歹我也是研究了SNA那么久了呢。

We will be using Gephi for visualization and analysis. The interactive demonstrations will be primarily in NetLogo, which you will be able to access through your web browser. If you would like to complete the programming assignments, which will be done in NetLogo and R, NetLogo is freely available here and R is freely available here.

显然这远远没有结束...生统方面,Mathematical Biostatistics Boot Camp自然也是用R的:

临床上也是...Data Management for Clinical Research

  • What resources will I need for this class?
    For this course, you will need: 1) an Internet connection; 2) software package capable of generating table-based CSV files (e.g. Microsoft Excel, Google Docs, Numbers); and 3) and an install of the open-source R programming platform.

好吧,我已经不奇怪会看到更多的了。这两年业界对于R的需求井喷,绝对跟学校里面的教育脱不开干系。不过如果我只是一味列举而不是比较,那岂不是有违统计学风范?R跟C或者Java比显然没有意思...不是做一件事儿的嘛。SAS的结果只有一门,Passion Driven Statistics,我猜很大程度上是这东西不免费,不能让每个学生都有的用;Matlab就多很多了,主要是Matlab的计算确实强大,其他的跟它确实没法比啊(至今我写模拟还是喜欢在Matlab里面写矩阵运算...R里面涉及到无路可逃的循环真的是让人忍无可忍),目测有9门课程使用。注:Stata结果为0,哎,真的是打不过免费软件啊。

赘述完毕...

 

Categories
事儿关经济

中文文本聚类小尝试(Text Clustering in R)

众所周知的,我会经常百无聊赖的玩一些比较好玩的东西。比如画画旅行地图啦,恶搞一下COS的版猪啦,抓抓新浪围脖啦。这不R大会又要开始了么,有一点点小数据也要玩玩啦。比如,呃,君不见周六上午三场演讲都是文本挖掘的,那我不研究一下文本挖掘怎么去混演讲听啊~自己动手先。

A nearby galaxy cluster about 65 million light years from Earth.
文本挖掘自然也有有个情景嘛。这不正好会议要排日程表嘛,那得把我们16个讲座分成四个半天,每天大约4场。这个应该怎么分呢?从直觉上来说,听众肯定是希望相关的话题放在相邻的时间,这样他们就可以选择自己感兴趣的时间段去听啦,不用在那里一坐两天。同时也便于之后的集中讨论嘛。于是这个目的就是:根据演讲的题目、摘要和关键字,进行聚类。这显然是一个无监督的学习嘛,我又没有一个特定的结果变量。

那么首先,自然是要对中文文本进行分词啦。这个嘛就可以偷个懒,直接用现成的R包rmmseg4j。(中间鼓捣若干编码问题,不赘述...)

然后就是聚类。这里继续偷懒,调用现成的文本处理包tm,可以直接生成文本词对应的矩阵。比如,一个编号为1的句子是 “我 在 中国”,编号为2的句子是“我 爱 中国” 那么生成的矩阵就是:

句子 我 在 中国 爱

1 1 1 1 0

2 1 0 1 1

就是说,把每个词都作为一个变量,然后统计它在每个句子出现的次数作为变量值。这样一来,如果总共有10个句子,有不重复的100个词,那么就会给出一个10×100的矩阵了。

有了这个矩阵之后,我们就相当于知道了每个个体的观测特征,那么就可以聚类了。比较简单的,可以直接算余弦相似度(比如google识别相似新闻的做法);也可以调用kmeans聚类。这里我们的摘要直接不会有特别多的相似,所以余弦相似度的区分度可能会不好。那么就先试试kmeans吧。

到这里,代码如下:

#读数据
library(xlsx)

presentations <- read.xlsx("r-presentations.xlsx", sheetName="Sheet1") #读excel数据

summary(presentations)

presentations$Title <- as.character(presentations$Title) #转文本

Encoding(presentations$Title) <- "UTF-8" #转换编码

presentations$Title

presentations$Abstract <- as.character(presentations$Abstract)

Encoding(presentations$Abstract) <- "UTF-8"

presentations$Abstract

presentations$KeyWords <- as.character(presentations$KeyWords)

Encoding(presentations$KeyWords) <- "UTF-8"

#分词

library("rmmseg4j")

presentations$raw_word <-with(presentations,paste0(KeyWords,Abstract, sep=";")) #连接所有标题、摘要、关键字

presentations$raw_word <- with(presentations, str_replace_all(raw_word, "R","")) #去掉r

presentations$seg <- mmseg4j(presentations$raw_word) #分词

#kmeans聚类

library("tm")

presebtation_seg <- Corpus(DataframeSource(presentations[,c("Title","seg")])) #转换到tm专用格式

presebtation_term <- TermDocumentMatrix(presebtation_seg, control = list(stopwords = TRUE)) #生成词频矩阵

presebtation_term <- t(as.matrix(presebtation_term)) #转换为matrix并转置

summary(presebtation_term)

presebtation_kmeans <- kmeans(presebtation_term, 7) #kmeans聚为7类

为什么我会在kmeans里面聚成7类呢?理论上只是要聚4类嘛。可是直接聚四类的话,区分度没那么好,一半多的演讲都聚到一类去了,没法安排嘛~所以只能增加聚类的个数,看看到时候再把小类合并。

聚成7类的结果如下:

Title cluster_result
R语言在eBay搜索引擎反馈与测试中的应用 1
营销分析模型及其在广告界的应用 2
系统生物学和转换医学中的R语言 + R in Systems Biology and Translational Medicine 3
R/Bioconductor在生物多维组学数据整合中的应用 3
R Case Study from EBAY APD 4
网络用户浏览路径分析 4
啤酒与尿布的当代版--商品分析在电子商务中的应用 4
基于RHadoop的关联规则挖掘 5
模型预测的利器——随机森林 5
基于R的地理信息系统 (R-based GIS) 6
R语言和其他计算机语言的混合编程 6
ggplot和knitr包简介 6
R与面向对象统计分析 6
twitteR包入门和应用 6
短文本分类器与电商品类数据挖掘 7
R语言环境下的文本挖掘 7

比较理想的是,聚类之后识别出来了两个文本挖掘的演讲...还有一堆R包的演讲。但是还是没法安排演讲嘛。看到这里,大家有没有发现,这样做最大的问题就是,聚类的时候把一些没有实际意义的虚词也聚类进来了,比如“的”;还有一些几乎所有演讲都会涉及的词,比如“R”和“分析”。这些词在其中是没有意义的,也会影响我们算dissimilarity的结果——这到底是按内容聚类啊,还是按作者的行文风格聚类啊?此外,虽然我们规定演讲摘要大都在100-200字,但还是有长有短,到目前我还没有对文本的出现频率用语句长度来加权...这也是不科学的嘛。那些原来在Google搜搜里面排名作弊的,不就是同样的内容复制10几次,来提高关键词出现频数(而不是频率)嘛。

为了解决这些问题,首先就是要去掉没有意义的虚词。这个不算太麻烦,把一些常用的虚词和转折词连接词之类去掉就可以了。其次,要去掉每个演讲都有的词。这里虽然可以一个个去看,不过简单一点,我们先统计一下词频嘛:

#高频词统计

presentations$seg2 <- unique((strsplit(presentations$seg,split=" "))) #断词

all_key_words <- iconv(unlist(presentations$seg2), from="UTF-8", to="GBK") #转换到GBK编码

all_key_words_fre <- as.data.frame(table(all_key_words)) #统计词频

names(all_key_words_fre)

all_key_words_fre <- arrange(all_key_words_fre,desc(Freq)) #按词频排序

all_key_words_fre[1:20,]$all_key_words #100个高频词

然后看一下TOP 20高频词:

1 的 105

2 数据 27

3 分析 24

4 和 21

5 图 18

6 在 17

7 挖掘 15

8 用户 15

9 应用 14

10 分类 13

11 了 13

12 语言 13

13 介绍 11

14 是 11

15 文本 11

16 试验 10

17 平台 9

18 ebay 9

19 案例 8

20 模型 8

所以看来,“挖掘”,“用户”,“文本”,“试验”,“平台”,“ebay”,“案例”,“模型”等等还是比较有区分度的词。按照这个思路,选择有限的几十个词重新分类,效果可能会有所改善。

此外,鉴于样本量不大(16个),所以可以人工的去看每个简介,手动标注tag作为聚类的变量。事实上,最后我还是这么做了一下,来在上述原始聚类结果上进行了一下重新的分组处理,形成了4大类。但是这个东西也不完全是可以直接用的,总要考虑时间之类的其他因素。最终的结果更多是人工思考的排序,估计李舰哥在确定顺序的时候更多的是按照经验和以往R会议的风俗。算法虽然好玩,但毕竟捕捉的还是人的思维模式,暂时没办法完美的取代吧。不过其实也差的不远呢。

最终人工结果:

冯兴东:R语言和其他计算机语言的混合编程

刘思喆:R语言环境下的文本挖掘

张翔:短文本分类器与电商品类数据挖掘

沈羽、周春英:R语言在eBay搜索引擎反馈与测试中的应用

周扬:基于R的地理信息系统

肖凯:twitteR包入门和应用

陈钢:系统生物学和转换医学中的R语言

杭兴宜:R / Bioconductor在生物多维组学数据整合中的应用

陈逸波:基于RHadoop的关联规则挖掘

李忠:R Case Study from EBAY APD

洪健飞:啤酒与尿布的当代版——商品分析在电子商务中的应用

廖明:营销分析模型及其在广告界的应用

肖嘉敏:网络用户浏览路径分析

刘成昊:模型预测的利器——随机森林

王雨晨:R与面向对象统计分析

魏太云:R基础作图与可重复研究

纯属好玩而已~不过R会议也举行了整整五届了,每次15个演讲的话也有15*9=135个演讲了。在这个样本量下,如果我们要出个论文集什么的,倒是可以直接用聚类的办法划分chapter了...嘻嘻。

Categories
事儿关经济

R并行做大数据时间序列分析与bootstrap

好久没写关于经济学的文章了...今天看到两个搞matching的人拿到诺奖,瞬间想起当年和一个学心理学出身的童鞋一起搞的考虑到心理学因素的matching game...可惜最后没有时间完善,就放在那里当雏形了。我们的idea应该还是蛮有新意的呢~matching需要设计机制,然后用博弈来解,哎,很好玩的...

很多人都知道我是不搞时间序列分析的,尤其不喜欢基于时间序列的因果推断(格兰杰因果检验几乎是被我打入黑名单的一个词)。但这次为什么专门写这么一篇blog post呢?其实,我不反对时间序列作预测嘛~而google这篇paper,目的在于预测,又是latex排版的(搞不好还有sweave或者knitr的功劳),读起来赏心悦目的多。

这篇paper题目是:

Large-Scale Parallel Statistical Forecasting Computations in R

地址见:http://research.google.com/pubs/pub37483.html

亮点自然是:大数据计算、map reduce、时间序列分析、bootstrap (所以说google是一家让人尊敬的公司,http://t.cn/zlbAYJY这里总结了他将学术成果悄无声息的服务于大众的案例)。虽说大半夜的,但看到这个东西,再也睡不着了,索性写完一吐为快。(那次沙龙还说到来着,看吧,google早把这些好玩的东西都搞定了,只是不公开拿出来给大家用罢了~)

1. 并行算法

并行算法这里,主要是map reduce。他们的任务主要有:

  • Facilitate parallelism of computations on up to thousands of machines without access to shared NFS filesystems(上千个节点的并行计算,无共享的NFS文件系统).
  • Make distribution of code and required resources as seamless as possible for analysts to minimize code modifications required to enable parallelism.(无缝衔接各环节,减少分析师写并行算法工作量)
  • No setup or pre-installation of R or specific libraries should be required on the machines in the cluster. A virtual machine for the workers should be created dynamically based on the global environment and available libraries of the caller(节点上无需事先安装R,虚拟机会自动按需构建).
  • Return results of parallel computations in list form directly back to the calling interactive session, as with lapply in R.(lapply函数直接返回并行结果)
  • Allow the parallel functions to be used recursively, so that MapReduce workers can in turn spawn additional MapReduces.(并行函数可以循环调用)

系统构架如下:

2013-12-09 16_48_39-37483.pdf - Adobe Reader

然后若干技术细节还包括,搞定R函数复制和计算的环境并行同步,搞定data.frame和list存储格式与在map函数中直接调用。搞定后,基本底层就搭建好了,剩下的就是调用R了。

2.时间序列预测

最典型的就是google trend的预测了...

2013-12-09 16_49_18-37483.pdf - Adobe Reader

这里他们直接用R包googleparallelism,然后希望用一些时序模型都尝试,做预测,然后取他们的均值(剔除上下总计20%)作为估计值以减少误差。示意图如下:

2013-12-09 16_49_35-37483.pdf - Adobe Reader
我只能说,果然是做机器学习的人啊,和random forest思路一致,弱的分类器结合起来,可能有意想不到的结果。同样的,每个模型都是多少有效的话,平均一下就更稳健啦,尤其是在大数据支撑下...

这样的平均之后已经无法直接推导方差和置信区间,所以他们采取了更依赖机器计算的bootstrap方法,直接强行算出来置信区间...喵的,我只能说谁让当年高斯那群天才整出来大数定律和中心极限定理呢?推不出来估量量方差不要紧嘛,直接重抽样模拟就好了...汗。

3.训练集

果然是做机器学习的人,接下来的思路就是直接一期期训练模拟呗。这个没啥说的了,见下图。

2013-12-09 16_49_51-37483.pdf - Adobe Reader

模拟出来的结果示意图:

2013-12-09 16_50_05-37483.pdf - Adobe Reader

至此,整个问题解决完毕。细节还请大家直接去看原文paper。我的几点感触吧:

  • 1. 机器学习之所以在业界这么受欢迎,主要是其确实能够解决问题。迅速、有效,这个是其他方法比不上的。
  • 2. 大数据、大规模计算,使得一些很简单的idea借助模拟和重抽样方法,大放异彩。
  • 3. 预测,有时候不比因果推断次要。
  • 4. 传统模型,需要适应大数据。
  • 5. 说到底,理论体系还是有待完善的。希望这类方法是下一个微积分,可以先用,然后慢慢补充完相应理论体系。这样,我们才知道什么时候,需要勒贝格测度来取代原有牛顿积分。

总之,虽然无奈,但是有用之物必有有用的道理。期待对理论研究的冲击和激发。

 

Categories
日常应用

探索R包reshape2:揉数据的最佳伴侣

前几天放出来的那个R的展示中,有说到其实学R的过程更多的就是熟悉各种函数的过程(学习统计模型不在此列...我个人还是倾向于不要借助软件来学习理论知识,虽然可以直接看codes...笔和纸上的推导还是不可或缺的基本功),然后各种基础函数熟悉了之后很多被打包好的函数就是缩短代码长度的利器了。

excel里面有神奇的“数据透视表(pivot table)”,其实很多时候真的已经很神奇了....不过我还是喜欢R,喜欢R直接输出csv或者xlsx的简洁。揉数据呢(学名貌似叫数据整理),我也还是喜欢写出来代码的形式,而不是直接向excel那样面对结果。只是感觉更加不容易出错吧。

揉数据,顾名思义,就是在原有的数据格式基础上,变化出来其他的形式。比如,长长的时间序列变成宽一点的~当然这个可以简单的借助reshape()函数了。可惜我还是不死心,想找一个更好用的,于是就自然而然的看到了reshape2这个包。

这个包里面函数精华在melt()*cast()。说实话melt()耗了我一段时间来理解,尤其是为什么需要先melt再cast...后来发现这个步骤简直是无敌啊,什么样的形状都变得更加容易揉了,大赞。

warm-up完毕,还是回到正题吧,怎么用reshape2揉数据呢?虽然reshape2支持array, list和data.frame,但是我一般还是习惯于用data.frame,所以还是说说这东西怎么揉吧。揉数据的第一步就是调用melt()函数,不用担心你的input是什么格式,这个函数array, list和data.frame通吃。然后,要告诉他哪些变量是(唯一)识别一个个体的,这句话是什么意思呢?我们先看melt()的参数:

 melt(data, id.vars, measure.vars,
    variable.name = "variable", ..., na.rm = FALSE,
    value.name = "value")

其中id.vars可以指定一系列变量,然后measure.vars就可以留空了,这样生成的新数据会保留id.vars的所有列,然后增加两个新列:variable和value,一个存储变量的名称一个存储变量值。这样就相当于面板数据的长格式了。直接拷一个作者给出的例子:

原数据:

head(airquality)
  ozone solar.r wind temp month day
1    41     190  7.4   67     5   1
2    36     118  8.0   72     5   2
3    12     149 12.6   74     5   3
4    18     313 11.5   62     5   4
5    NA      NA 14.3   56     5   5
6    28      NA 14.9   66     5   6
dim(airquality)
[1] 153   6

然后我们将month和day作为识别个体记录的变量,调用melt(airquality, id=c("month", "day"))

head(melt(airquality, id=c("month", "day")))
  month day variable value
1     5   1    ozone    41
2     5   2    ozone    36
3     5   3    ozone    12
4     5   4    ozone    18
5     5   5    ozone    NA
6     5   6    ozone    28
dim(melt(airquality, id=c("month", "day")))
[1] 612   4

嗯,这样数据就变长了~然后,就可以随意的cast了...dcast()会给出宽格式的数据,比如我们想把day作为唯一的识别,那么:

names(airquality) <- tolower(names(airquality))
aqm <- melt(airquality, id=c("month", "day"), na.rm=TRUE)
head(dcast(aqm, day ~ variable+month))
  day ozone_5 ozone_6 ozone_7 ozone_8 ozone_9 solar.r_5 solar.r_6 solar.r_7 solar.r_8 solar.r_9 wind_5 wind_6 wind_7 wind_8 wind_9 temp_5 temp_6
1   1      41      NA     135      39      96       190       286       269        83       167    7.4    8.6    4.1    6.9    6.9     67     78
2   2      36      NA      49       9      78       118       287       248        24       197    8.0    9.7    9.2   13.8    5.1     72     74
3   3      12      NA      32      16      73       149       242       236        77       183   12.6   16.1    9.2    7.4    2.8     74     67
4   4      18      NA      NA      78      91       313       186       101        NA       189   11.5    9.2   10.9    6.9    4.6     62     84
5   5      NA      NA      64      35      47        NA       220       175        NA        95   14.3    8.6    4.6    7.4    7.4     56     85
6   6      28      NA      40      66      32        NA       264       314        NA        92   14.9   14.3   10.9    4.6   15.5     66     79
  temp_7 temp_8 temp_9
1     84     81     91
2     85     81     92
3     81     82     93
4     84     86     93
5     83     85     87
6     83     87     84

或者对于每个月,求平均数:

 head(dcast(aqm, month ~ variable, mean, margins = c("month", "variable")))
  month    ozone  solar.r      wind     temp    (all)
1     5 23.61538 181.2963 11.622581 65.54839 68.70696
2     6 29.44444 190.1667 10.266667 79.10000 87.38384
3     7 59.11538 216.4839  8.941935 83.90323 93.49748
4     8 59.96154 171.8571  8.793548 83.96774 79.71207
5     9 31.44828 167.4333 10.180000 76.90000 71.82689
6 (all) 42.12931 185.9315  9.957516 77.88235 80.05722

当然还有更强大的acast(),配合.函数:

library(plyr) # needed to access . function
acast(aqm, variable ~ month, mean, subset = .(variable == "ozone"))
             5        6        7        8        9
ozone 23.61538 29.44444 59.11538 59.96154 31.44828

嗯,基本上数据就可以这么揉来揉去了...哈哈。怎么感觉有点像数据透视表捏?只是更加灵活,还可以自定义函数。

此外还有recast()可以一步到位,只是返回的是list;colsplit()可以分割变量名...函数不多,却精华的很啊。

--------------------
题外废话:我的小册子哎~只能这样零零碎碎的写一些了,事后再统一整理进去好了。不要鄙视...