使用 R 编写统计程序: 第 2 部分. 函数式程序设计和数据探索

文档选项 未显示需要 JavaScript 的文档选项

级别： 初级

David Mertz (mertz@gnosis.cx), 开发人员, Gnosis Software
Brad Huntting (huntting@glarp.com), 数学家, University of Colorado

2004 年 10 月 18 日

本系列共分三个部分，这是第 2 部分，David 和 Brad 将基于他们的第一期文章继续对 R（R 是一个作为自由软件发布的功能丰富的统计 环境）进行论述。既然数据是井然有序的，所以我们将深入研究语言的功能性。

R 既是一门强大的函数式程序设计语言，也是一个用来对数据集进行统计探索的通用环境。 作为一个环境，R 允许您通过一个定制的命令行来创建数据的很多图形化表示。参阅本系列文章的 第 1 部分以及 Cameron Laird 的关于 R 的文章（都在 参考资料 中 列出了），以深入了解 R 及其背景，并获得关于可以运行 R 的平台和很多它可用的程序包的资料。

与 Python、Ruby、用于 Tcl/Tk 的 wish 以及很多 Lisp 环境的 交互式 shell 一样，R shell 是用于探索操作极好途径 —— 命令的重新调用和编辑让您可以尝试命令 的多种变体。与很多其他程序设计语言交互式 shell 相反，除了与 R 的以数据为中心的本质保持一致 以外，R shell 可以任意地保存一个完整的环境（每个工作目录一个）。命令历史是其中的一个实用 部分：可以在 .Rhistory 文件中对它进行检查和修改。不过，保存环境的主要方面是将工作数据以二进 制格式保存到 .RData 中。顺便说一句， rm() 命令的灵活运用可以停止 无限制地保存数据环境（ remove() 等同于 rm() ）。

本文是关于 R 的文章的第二期。 第 1 部分 介绍了 R 的数据类型的一些基本概念 —— 从向量开始，还包括多维数组（2-D 数组也称为“矩阵”）、用于“灵巧”表 的数据帧、用于异构集合的列表，等等。前一期还对合著者 Brad 手边的一个大规模数据集进行了一些基本的统计 分析，并绘制了图形，那个数据集是 Brad 住宅周围一年的温度历史，每三分钟采样一次。与现实世界中的很多数据 集类似，我们知道或者猜想 Brad 的温度数据中存在瑕疵和异常 —— 实际上，第二期的部分内容将尝试 去理解数据的可疑特性。

总体上讲，本文将做两件事情。首先，我们将继续使用第一期文章中引入的数据集来比上一次更深入地 研究 R 语言本身。然后，我们将研究数据中的一般模式，并向您展示如何使用 R 找到它们。

R 作为程序设计语言

（GPL'd）R 编程语言继承自两门语言：S/S-PLUS 和 Scheme（Lisp）。在这里值得强调它与 Scheme 的联系， 因为 R 的函数式程序设计方面继承自家族的 Scheme 一方。细心的读者可能已经注意到，在第一期的例子中 奇怪地缺少一些内容：流程控制！本文中同样仍然没有这些内容。实际上，R 有非常好的 if 、 else 、 while 和 for 命令，所有这些命令都 与 Python 中同样名称的命令非常类似。（其他语言对同一命令的拼写稍有不同。） 为了适度的冗余，另外还提供了 repeat 、 break next 和 switch 。命令式语言程序员将 对不使用流程控制就可以完成的工作而感到惊讶；如果不使用它，很多工作会更加简单，这将使他们更为惊讶!

声明您要做什么

按其函数式程序设计传统，使用普通声明语句，您就可以在 R 中完成几乎所有的事情。 大部分情况下，R 的两个特性使得没有必要再使用命令流程控制。首先，您已经看到， 大部分对集合对象的操作会按元素（elementwise）起作用。不需要手工地在一个向量的数据中进行循环 来对其元素进行某项处理，因为您可以简单地对向量中的全部元素进行某项处理：

> a = 1:10            # Range of numbers 1 through 10
> b = a+5             # Add 5 to each element
> b                   # Show b vector
 [1]  6  7  8  9 10 11 12 13 14 15

或者，通过使用一个“下标数组”，您可以对具有特定下标的元素有选择地进行操作：

> c = b               # Make a (lazy) copy of b
> pos = c(1,2,3,5,7)  # Change the prime indices
> c[pos] = c[pos]*10  # Reassign to indexed positions
> c                   # Show c vector
 [1]  60  70  80   9 100  11 120  13  14  15

或者，或许是最好的方式，您可以使用 Haskell 或 Python 中非常类似于 list comprehensions 的语法，只对具有期望属性的元素进行操作：

> d = c
> d[d %% 2 == 0] = -1 # Reassign -1 to all even elements
> d
 [1] -1 -1 -1  9 -1 11 -1 13 -1 15

非常敏锐的读者可能已经注意到，这些例子使用了 = ，而在 上一期文章中那些地方使用的是 <- 。等号可以完成与赋值箭头 （assignment arrow）相同的功能，但是只能用于最高层次范围，不能用于嵌套的表达式。 如果拿不准，那么使用箭头更安全。

回页首

探索温度数据

在上一期文章中，我们已经使用 read.table() 函数将四个位置每个的 大约前 170,000 个读数读入到一个数据帧中。不过，为了更方便地访问各个温度序列，我们将数据 列拷贝到以测试位置命名的向量中： outside 、 livingroom 、 basement 和 lab 。另外要记住，每个序列中都丢失了一些数据点：有时是因为负责 记录的计算机关闭了；有时是因为指令失败；四个温度计没有全部准确地在同一天联机；等等。换句话 说，我们的温度数据与现实生活中您可能要处理的数据非常相似。

在我们对温度数据集进行的最初的探索中，在绘制直方图时我们注意到有一个异常。室外温度 刚好在 24 摄氏度附近出现了一个大的尖峰。我们的第一个猜测是，这个突起代表温度读数的 偶然互换，即一个室内温度记录到了室外日志中（相应地，室外温度记录到了室内日志中）。 作为研究 R 的一种方法，来看我们是否能证明或者驳斥这一解释。

找出异常

探索异常数据的第一步只是 找到 它。具体讲，一个单点异常将由数据点每一侧温度 的突然波动标识出来。更具体地说，我们认为一个异常点的两个相邻点比参考点或者特别高， 或者特别低。例如，温度的抽样序列（以三分钟为间隔）“20, 16, 13”描述了温度的反常快速 下降，但是并不代表中间的读数是单点错误。当然，我们不会 先验地 排除 其他类型错误或异常（不是只与单个读数相关的）的存在。

我们最初识别异常读数的想法变得复杂了。或许我们可以寻找序列的 Fourier 变换中 的高频率成分。或者可能我们可以对温度向量进行解析求导（analytic derivatives， 要得到形变可能需要进行二次求导）。或者，我们还可以尝试在温度向量上做卷积 （convolutions）。所有这些选项在 R 中都已经具备。但是，反复考虑了 这些高级想法以后，我们决定回到实际中来。我们正在寻找的简单模式是与相邻读数 有大尺度差异的温度读数。换句话说，我们可以借助减法的魔力。

我们用来寻找相邻点之间的差异的技巧是，创建一个数据帧，它的列是给定连续数据点， 分别是先前数据点（prior）、当前数据点（current）、后继数据点（following）。 这就让我们可以将数据帧作为一个整体进行过滤，挑选出可能引起注意的连续点。

> len <- length(outside)    # Shorthand name for scalar length of vector
> i <- 1:(len-2)            # Shorthand vector of data frame rows nums
> # Create a data frame for window of measurements per row
> odf <- data.frame(lst=outside[1:(len-2)],
+                   now=outside[2:(len-1)],
+                   nxt=outside[3:len] )
> # Create vector of local temperature fluctuations, add to data frame
> odf$flux <- (odf[i,"now"]*2) - (odf[i,"lst"] + odf[i,"nxt"])
> odf2 <- odf[!is.na(odf$flux),]  # Exclude "Not Available" readings
> oddities <- odf2[abs(odf2$flux) > 6,] # It's odd if flux is over 6

那么，过滤后我们会得到什么呢？让我们来看：

> oddities
        lst  now  nxt flux
2866   21.3 15.0 14.9 -6.2
79501  -1.5 -6.2 -4.1 -6.8
117050 21.2 24.6 21.6  6.4
117059 20.6 23.4 20.1  6.1
127669 24.1 21.2 24.7 -6.4
127670 21.2 24.7 21.5  6.7

结果证明，我们有一些读数相对于其相邻点有大的波动。但是它们大部分看起来不太像是被互换 的读数。例如，在时间步 79501，温度 6.2 被两个明显温暖的温度环绕着。但是，这些 任何 一个都不可能是室外温度 —— 貌似更合理的解释是，有一阵特别寒冷的微风。

非常有可能在时间步 117059 或者 127670 附近发生了互换。中间的温度恰好与室内的 24 度 接近，而它相邻的读数虽然明显不至于结冰，但是显然比较低。可能是在春季发生了一些互换。

将实用的操作包装到函数中

现在我们想要知道，我们少量候选的互换是否能在其他读数中出现。如果在室内的某个位置 找到一个“丢失的”室外温度，那么这将有力地支持我们的假设。但是我们 实在 不愿意 重新输入全套步骤，而只是将各处 outside 替换为 lab 。我们真正应该做的是，将全部步骤序列参数化为一个 可重用的函数：

oddities <- function(temps, flux) {
  len <- length(temps)
  i <- 1:(len-2)
  df <- data.frame(lst=temps[1:(len-2)],
                   now=temps[2:(len-1)],
                   nxt=temps[3:len])
  df$flux <- (df[i,"now"]*2) - (df[i,"lst"] + df[i,"nxt"])
  df2 <- df[!is.na(df$flux),]
  oddities <- df2[abs(df2$flux) > flux,]
  return(oddities)
}

一旦我们拥有了一个函数，对不同的数据集和波动临界点进行过滤就简单多了。结果表明， 运行 oddities(lab,6) 不能产生任何与室外温度候选互换 相匹配的时间步（ livingroom 或 basement 同 样也不能）。但是对实验室温度的考察得到的结果本身有些令人惊讶：

> oddities(lab, 30)
       lst  now  nxt  flux
47063 19.9 -2.6 19.5 -44.6
47268 17.7 -2.6 18.2 -41.1
84847 17.1 -0.1 17.0 -34.3
86282 14.9 -1.0 14.8 -31.7
93307 14.2 -6.4 14.1 -41.1

这些类变化不是我们所期望的。或许最好的解释是 Brad 在极度寒冷的冬天有时 会打开实验室的窗户。如果是这样，合著者 David 会对 Brad 的炉子三分钟 排气就能恢复全部热量这种效率感到惊奇。

真实的解释

可能到根本的完全数据中去查找奇怪的时间步将会使事情明了：

> glarp[47062:47066,]
             timestamp basement  lab livingroom outside
47062 2003-10-31T17:07     21.5 20.3       21.8    -2.8
47063 2003-10-31T17:10     21.3 19.9       21.2    -2.7
47064 2003-10-31T17:13     20.9 -2.6       20.9    -2.6
47065 2003-10-31T17:16     20.8 19.5       20.8    -2.6
47066 2003-10-31T17:19     20.5 19.4       20.7    -2.8

首先，我们注意到 oddities() 返回的结果有差 1 （off-by-one）错误。噢，是的，我们使用一个偏移量来获得每个数据帧序列的一个窗口。 但是完全数据集本身好像是支持“打开门”的想法 —— 在科罗拉多州，万圣节前夜的晚上 5 点是非常冷的，实际上可能那时 trick-or-treaters（玩“不给糖，就捣蛋”游戏的人）来到 了 Brad 门口（用了三分钟时间收取糖果）。这样，这个谜团好像揭开了。

尽管如此，那么我们开始探索时的那个谜团呢？为什么室外温度中 24 度附近有那样一个尖峰？ 或许我们应该更仔细地看一下直方图。特别的，我们可以使用明确的标准来索引一个向量，将 我们的直方图收缩到只是我们感兴趣的范围：

hist(outside[outside < 26 & outside > 23],
     breaks=90, col="green" border="blue")

观察对温度尖峰的特写，我们可以看到恰好在尖峰附近有一个明显的波谷。看起来好像如果您 能平滑 24.7 度附近的十分之几度的区域，那么您将几乎得到预期的结果。很可能这让我们想到 温度计仪器内部的不正常的舍入行为。好了，可能 —— 我们还是不敢保证。而且，即使是经过了 一些平滑，那里还是会有微小的尖峰。

回页首

中间统计分析

计算线性与非线性回归模型的能力是 R 的强大功能之一。为了查看一个小例子， 让我们来创建两个向量。 x 是以天为单位的抽样时间， 自数据集的开始之处得到。 y 是相应的室外温度。

> y <- glarp$outside
> x <- 1:length(y)/(24*60/3)

我们可以用下面的方法给此数据匹配一条线：

> l1 = lm(y ~ x)
> summary(l1)
Call:
lm(formula = y ~ x)
Residuals:
     Min       1Q   Median       3Q      Max
-29.4402  -7.4330   0.2871   7.4971  23.1355
Coefficients:
        Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.2511014  0.0447748  228.95   <2e-16 ***
x           -0.0037324  0.0002172  -17.19   <2e-16 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 9.236 on 169489 degrees of freedom
Multiple R-Squared: 0.00174,    Adjusted R-squared: 0.001734
F-statistic: 295.4 on 1 and 169489 DF,  p-value: < 2.2e-16

“~”语法指的是一个公式对象。这样的结果是，要求 R 找到系数 A 和 B 使得 sum((y[i]-(A*x[i]+B))^2) 最小。当 A 为 -0.0037324 （非常平坦），B 为 10.2511014 时得到最佳匹配。注意，标准化残差（residual standard error） 为 9.236，这几乎与 y 初始时的标准偏差大小相同。这就 告诉我们，对室外温度来说，时间的简单线性函数是一个非常不好的模型。

更好的模型可能是使用以 1 天和 1 年为周期的正弦和余弦函数。要实现这个模型， 我们可以将公式改变为：

> l2 = lm(y ~ +I(sin(2*pi*x/365)) +I(cos(2*pi*x/365))
+             +I(sin(2*pi*x)) +I(cos(2*pi*x)) )

这个公式语法有些复杂：在 I() 调用内部，算术符号与它们通常的含义相同， 所以，例如第一个 I() 表示：2 乘 pi 乘 x，再除以 365。 I() 前的“+”号不是表示加，而是表示 + 后的 因式要包含在模型之中。与往常一样，可以使用 summary() 命令显示结果。

> summary(l2)
Call:
lm(formula = y ~ +I(sin(2 * pi * x/365)) +I(cos(2 * pi * x/365))
                 +I(sin(2 * pi * x)) +I(cos(2 * pi * x)))
Residuals:
     Min       1Q   Median       3Q      Max
-21.7522  -3.4440   0.1651   3.7004  17.0517
Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)             9.76817    0.01306  747.66   <2e-16 ***
I(sin(2 * pi * x/365)) -1.17171    0.01827  -64.13   <2e-16 ***
I(cos(2 * pi * x/365)) 10.04347    0.01869  537.46   <2e-16 ***
I(sin(2 * pi * x))     -0.58321    0.01846  -31.59   <2e-16 ***
I(cos(2 * pi * x))      3.64653    0.01848  197.30   <2e-16 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 5.377 on 169486 degrees of freedom
Multiple R-Squared: 0.6617,     Adjusted R-squared: 0.6617
F-statistic: 8.286e+04 on 4 and 169486 DF,  p-value: < 2.2e-16

残差还是比较大（5.377），但是我们可以聊以自慰的是，科罗拉多州以气候难以预料而闻名。

R 还具有更多用于时间序列分析的工具。例如，我们可以绘制出起居室温度的自相关函数：

> acf(ts(glarp$livingroom, frequency=(24*60/3)),
+     na.action=na.pass, lag.max=9*(24*60/3))

嵌入的 ts() 函数通过向量 glarp$livingroom 创建 一个时间序列对象。抽样频率用每天抽样次数指定。毫无意外，当偏移是天的整倍数时，温度之间相关性很强。另外， 注意 7 天时的稍大的尖峰。这是因为实际上 Brad 的自动调温器在周末时（当 Brad 通常白天不在房子里时）设置 了不同的温度，结果是七天窗口的相关性稍微更强。

回页首

结束语

到目前为止，我们已经使用 R 在一个数据集中对模式和异常进行了分析，这个数据集正好有 影响几乎所有现实世界科学数据的那种不连续和潜在问题。本文还研究了 R 的函数式程序设计 风格可以如何有助于加速模式、数据和分析情形的的探索。本系列的第 3 部分将使用更复杂的统计方法（不过还会总体上略微谈到 R 中的大量工具），继续在大数据集中 寻找模式。

参考资料

• 您可以参阅本文在 developerWorks 全球站点上的 英文原文。
  
  
• 本系列的 第一期 （developerWorks，2004 年 9 月）介绍了 R 以及它的一些特性和语言特征。
  
  
• Cameron Laird 的 服务器诊所：数据处理的利器 R （developerWorks，2003 年 7 月）给出了一个关于 R（以及 S）的很好概述， 并且有很多参考资料。
  
  
• 查看 用于统计计算的 R 项目 的主页。 这个 R Web 站点上有大量的包括从教程到完全的 API 描述在内的所有文档。对于 刚接触 R 的读者来说尤其会感兴趣的两个文档是：
  • An Introduction to R（PDF 格式）
  • R Data Import/Export（PDF 格式）
  
  
• 作为 Python 用户，David 的 Charming Python 专栏的一些读者对 R 的 Python 绑定表示出特别的喜好。实际上，有两种这样的绑定： RPy， 以及更老的 RSPython，这个也不错。（不过，David 的感觉是 RPy 绑定更好一些）。这两种绑定都让您可以通过 Python 代码透明地调用所有的 R 函数，使用 Python 对象作为函数的参数。
  
  
• 不敢肯定的是，在大部分系统上，您都可以在 R 命令行中输入 help.start() ， 以打开一个包含有为 R 文档生成的 HTML 页的浏览器。
  
  
• 可以在线获得 S 和 S-Plus 历史摘要。
  
  
• David 撰写了关于 可爱的 Python: Numerical Python 的一篇 Charming Python 连载文章，感觉上它与 R 类似并具有很多相同的能力（不过，R 的范围要广泛得多）。
  
  
• 可以下载温度数据以及用来将初始的日志格式转换为更好的制表格式的 Python 脚本：
  • basement.gz
  • lab.gz
  • livingroom.gz
  • outside.gz
  • measures.py
  
  
• 在 developerWorks Linux 专区 可以找到 更多为 Linux 开发者准备的参考资料。
  
  
• 自 developerWorks 的 Speed-start your Linux app 专区下载可以运行于 Linux 之上的经过挑选的 developerWorks Subscription 产品免费测试版本，包括 WebSphere Studio Site Developer、WebSphere SDK for Web services、WebSphere Application Server、DB2 Universal Database Personal Developers Edition、Tivoli Access Manager 和 Lotus Domino Server。要更快速地开始上手，请参阅针对各个产品的 how-to 文章和技术支持。
  
  
• 通过参与 developerWorks blogs 加入到 developerWorks 社区。
  
  
• 在 Developer Bookstore Linux 专栏中订购 打折出售的 Linux 书籍。
  

作者简介

		对 David Mertz 来说，世界就是一个舞台，他的职责就是专门提供临场指导。欢迎对过去或将来的文章提出您的意见和建议。 您可以通过 mertz@gnosis.cx 与 David 联系，通过他的 个人主页 来了解他生活的方方面面。阅读他所著的 Text Processing in Python 一书。

		14 年来，Brad 一直在为三个不同的公司从事 UNIX® 系统管理和网络工程方面的工作。 他他当前在 University of Colorado（位于 Boulder）攻读应用数学博士学位，通过为计算机科学系提供 UNIX 支持来获得学费。

对本文的评价

太差！ (1) 需提高 (2) 一般；尚可 (3) 好文章 (4) 真棒！(5) 建议？

回页首