内容


使用 R 编写统计程序

第 1 部分. 初涉大量统计工具

整理原始的数据

Comments

系列内容:

此内容是该系列 # 部分中的第 # 部分: 使用 R 编写统计程序

敬请期待该系列的后续内容。

此内容是该系列的一部分:使用 R 编写统计程序

敬请期待该系列的后续内容。

R 环境 本身并不是要成为一门编程语言,而是要成为一个用于研究数据集(包括范围广泛 的数据属性图形化描述的生成)的交互式工具。您可以将在一个会话中生成的图形以及所采取的步骤都 保存下来以备后用,这对于重新获得上次退出时的工作环境(每一个项目的)来说尤其实用。默认情况 下,R 命令保存在一个会话历史记录中,不过,您也可以将特别有用的指令序列保存到 .R 文件中,这样在会话中您就可以通过 source() 来使用它。

R 的创建者在“An Introduction to R”(参阅下面 参考资料 中的完整链接) 中描述了他们的目标:

在使用 R 进行分析时,建议您应该使用单独的工作目录。在分析的过程中,创建名称为 x 和 y 的 对象是非常常见的。这种名称在单个分析环境中通常是有意义的,但是,当 [...] 在同一 个目录中进行多个分析时,就难以判定它们是什么了。

每个目录下都会生成包含有会话中工作对象二进制序列的隐藏文件。这样就使得您可以使用先前的活动 变量来重新开始一个会话。

(GPL'd)R 编程语言继承于两门语言,一门是私有的 S/S-PLUS 编程语言,R 的大部分语法 都是继承自它,另一门是 Scheme 编程语言,R 的语义方面很多(更微妙的)继承自它。S 的出现 可以追溯到 1984 年,在后来的版本(包括 S-PLUS)中有了很多改进。当然,Scheme(和 Lisp 一样)的 历史最为古老。1997 年,R 作为 S/S-frPLUS 的一个自由软件扩展集问世,从那时起它就拥有了一个 日益兴旺的用户和开发者社区。您可以受益于 R,但不必须担心它的继承权问题。Cameron Laird 在 他的 developerWorks 文章 R handy for crunching data 中介绍了与 R 相关的一些恰到好处的背景和参考资料。

在很多计算机平台上都可以使用编译过的 R:Linux、Windows、Mac OS X 以及 Mac OS Classic。 当然,也可以获得源代码,为其他平台进行编译(例如,合著者 Brad 为 FreeBSD 编译了 R)。 在多种平台上 R 会遇到一些小问题:例如,在 David 的 Mac OS X 机器上,使用 Quartz 绘制 输出时生成了一个没有反应的显示窗口;更糟糕的是,在 Brad 的 FreeBSD/AMD Athlon 机器上, 退出 R 竟然会导致重新引导(这可能与不正确的 SSE 内核选项有关,但是这种行为还是令人厌烦)。 虽然如此,但是,R 总的来讲是稳定的、快速的,而且附带了范围惊人的统计和数学函数。可选的 程序包为近乎无限的标准程序包和函数集合提供了更进一步的补充。

R 的数据模型

在 R 中,基本的数据对象是向量。向量的多个变种提供了另外的能力,比如(多维)数组(arrays)、 数据帧(data frames)、(不同种类的)列表(lists),以及矩阵(matrices)。与 NumPy/NumArray 或 Matlab 中非常类似,对向量等的操作会作用于成员数据的每一个元素。可以实际地通过 R 的一些生动示例来 感受它的语法(在这些初始的例子中包括了 shell 提示符和响应):

清单 1. 对向量及所有元素的操作
> a <- c(3.1, 4.2, 2.7, 4.1)  # Assign with "arrow" rather than "="
> c(3.3, 3.4, 3.8) -> b       # Can also assign pointing right
> assign("c", c(a, 4.0, b))      # Or explicitly to a variable name
> c                              # Concatenation "flattens" arguments
[1] 3.1 4.2 2.7 4.1 4.0 3.3 3.4 3.8
> 1/c                            # Operate on each element of vector
[1] 0.3225806 0.2380952 0.3703704 0.2439024 0.2500000 0.3030303 0.2941176
[8] 0.2631579
> a * b                          # Cycle shorter vector "b" (but warn)
[1] 10.23 14.28 10.26 13.53
Warning message:
longer object length
        is not a multiple of shorter object length in: a * b
> a+1                            # "1" is treated as vector length 1
[1] 4.1 5.2 3.7 5.1

下面有关于 R 语法的注标(indexing)、分割、可选参数以及其他内容的更多示例。R shell 提示符 —— 尤其是如果您安装了 GNU readlines —— 是一个极好的用于研究的界面。工作时别忘记了 通过 help(function) 命令来学习更多内容(您也可以使用 ?function)。Python shell 用户将发现 R shell 实在是很熟悉 —— 而且将会高度评价两者的实用程序。

温度数据集

近一年以来,Brad 一直在通过他房子里面以及附近的四个温度计来收集温度数据, 并使用 GnuPlot(参阅 参考资料 以获得关于 Gnuplot 更多资料的链接) 将不断变化的读数窗口自动编译为可以通过 Web 访问的图形。尽管此类 hackerish 的数据收集 不能真正地用于任何主要的科学用途,但是,它拥有很多 类似于 科学数据的 极好特性。

每三分钟收集一次数据,这样,一年中就会有非常多数据点(四个测量位置大约共有 750,000 个)。 由于温度计、传输通路或者进行记录的计算机的各种失误,有一些数据会丢失。在较为少见的一些情况下, 据了解,由于计时错误,会导致单线传输通路会发生同时读取的问题。换句话说,Brad 的温度数据看 起来与 相当好 的实际科学数据很类似,但是还是有一些小问题,也并不是理想的。

读取数据

温度数据收集到四个以收集位置命名的独立数据文件中,每个的格式如下:

清单 2. 初始温度数据文件的格式
2003 07 25 16 04 27.500000
2003 07 25 16 07 27.300000
2003 07 25 16 10 27.300000

读取数据的第一个步骤可能是这样:

清单 3. 读取温度数据的第一个步骤
> lab <- read.fwf('therm/lab', width=c(17,9))   # Fixed width format
> basement <- read.fwf('therm/basement', width=c(17,9))
> livingroom <- read.fwf('therm/livingroom', width=c(17,9))
> outside <- read.fwf('therm/outside', width=c(17,9))
> l_range <- range(lab[,2]) 	  # Vector of data frame: entire second column
> b_range <- range(basement[,2])    # range() gives min/max
> v_range <- range(livingroom[,2])
> o_range <- range(outside[,2])
> global <- range(b_range, l_range, v_range, o_range)
> global          			  # Temperature range across all sites
[1] -19.8  32.2

整理数据

未经处理的初始数据格式有一些问题。一方面,没有显式地标明丢失的数据,而只是通过 不存在的行和时间戳进行了简单的标记。另外,日期的存储格式也是不标准的(而不是 ISO8601/W3C),有内部的空格。有一个小问题,那就是四个文件中重复的时间戳浪费了空间。 当然,我们可以在 R 本身中来整理数据,不过我们采取的是“R Data Import/Export”一文中 R 作者的建议(参阅 参考资料 以获得链接)。文本处理最好使用 专门进行文本处理的语言来完成:在我们的例子中,我们编写了一个 Python 脚本来生成 R 直接 可读的统一标准的数据文件。例如,新的数据文件 glarp.temps 的前几行是:

清单 4. 统一标准的温度数据格式
timestamp      basement     lab     livingroom  outside
2003-07-25T16:04 24.000000      NA     29.800000  27.500000
2003-07-25T16:07 24.000000      NA     29.800000  27.300000
2003-07-25T16:10 24.000000      NA     29.800000  27.300000

让我们来使用改进的数据集:

清单 5. 使用统一标准的温度数据
> glarp <- read.table('glarp.temps', header=TRUE, as.is=TRUE)
> timestamps <- strptime(glarp[,1], format="%Y-%m-%dT%H:%M")
> names(glarp)              # What column names were detected?
[1] "timestamp"  "basement"   "lab"        "livingroom" "outside"
> class(glarp[,'basement'])  # Kind of data is in basement column?
[1] "numeric"
> basement <- glarp[,2]         # index by position
> lab <- glarp[,'lab']          # index by name
> outside <- glarp$outside      # equiv to prior indexing
> livingroom <- glarp$living    # name with unique initial name
> summary(glarp)  		# Handy built-in to describe most R objects
  timestamp            basement            lab            livingroom
 Length:171349      Min.   :   6.40   Min.   :  -6.40   Min.   :   7.20
 Class :character   1st Qu.:  17.00   1st Qu.:  16.60   1st Qu.:  18.10
 Mode  :character   Median :  19.10   Median :  17.90   Median :  20.30
                    Mean   :  18.88   Mean   :  18.12   Mean   :  20.17
                    3rd Qu.:  20.50   3rd Qu.:  19.50   3rd Qu.:  22.00
                    Max.   :  27.50   Max.   :  25.50   Max.   :  31.30
                    NA's   :1854.00   NA's   :2406.00   NA's   :1855.00
    outside
 Min.   : -19.800
 1st Qu.:   2.100
 Median :   9.800
 Mean   :   9.585
 3rd Qu.:  17.000
 Max.   :  32.200
 NA's   :1858.000

基本的统计分析

我们已经看到 range() 函数: min()max() 找到了一个数据范围的单独的极值。 summary() 显然 也会显示这一信息,不过不直接用于其他计算中。首先让我们来寻找此数据的另外一些非常基本的统计属性:

清单 6. 温度数据的基本统计计算
> mean(basement)        	# Mean fails if we include unavailable data
[1] NA
> mean(basement, na.rm=TRUE)
[1] 18.87542
> sd(basement, na.rm=TRUE)  	# Standard deviation must also exclude NA
[1] 2.472855
> cor(basement, livingroom, use="all.obs")   # All observations: no go
Error in cor(basement, livingroom, use = "all.obs") :
        missing observations in cov/cor
> cor(basement, livingroom, use="complete.obs")
[1] 0.9513366
> cor(outside, livingroom, use="complete.obs")
[1] 0.6446673

正如您直观上所期待的,两个室内温度的相关性比任意一个与室外温度的相关性更强。 这也很容易验证。

温度的分布

您已经看到均值和标准偏差,直观上您可能会觉得温度是正态分布的。让我们来看看:

清单 7. 生成一个短线形式的直方图
> hist(outside)

很多 R 命令都会弹出第二个窗口来显示数据集的 plot、chart 或者 diagram。由于平台 和个人配置的不同,如何这样做的细节也不相同。您也可以将这些图形重定向到外部文件 中以备后用。上面的 hist() 命令生成:

图 1. 室外温度的默认直方图
室外温度的默认直方图
室外温度的默认直方图

第一次尝试很不错。有一些参数可以限制舍入界限:

清单 8. 改变直方图舍入密度
hist(outside, breaks=1000, border="blue")
图 2. 室外温度的稠密直方图
室外温度的稠密直方图
室外温度的稠密直方图

注意,在稠密直方图中 7-12 度附近出现了偶然性的不平滑,有一些温度频率非常高, 也有一些频率出人意料地低。我们相信,这些非常明显的不连续表明存在取样偏差,也许要归因于 工具特性。另一方面,在 24 度附近有一个非常高而且窄的突起 —— 恰好是在室内温度自动 调节区域附近 —— 可能是我们前面提到度量变换的结果,与工具传输通路有关。无论如何,这 个图中包含有一些有趣的内容,值得去研究和分析。

两个更为生动的变化展示出室内温度的分布:

清单 9. 起居室温度直方图
> hist(livingroom, breaks=40, col="blue", border="red")
> hist(livingroom, breaks=400, border="red")
图 3. 起居室温度的 40-级 直方图
起居室温度的 40-级 直方图
起居室温度的 40-级 直方图
图 4. 起居室温度的 400-级 直方图
起居室温度的 400-级 直方图
起居室温度的 400-级 直方图

起居室温度分布看起来更合理。在高解析度的图中出现了一些不连续,似乎是小规模工具 取样的原因。不过,总体上的模式遵循了我们所期望三峰分布,这种分布基于 Brad 的 计时器控制的自动调温器(高峰值在 21 附近,稍低的峰值在 16 和 24 度附近)。

数据可视化进阶

每一个测量位置是温度值的一个线性向量。不过,直观地讲,我们会期望数据中有两个 主要的循环:每天和每年的(晚上和冬天冷)。

我们的第一个问题是将一维数据向量转化为数据点的二维矩阵。然后,我们将可视化这个二维 数据集:

清单 10. 重构向量并绘制温度图
> oarray <- outside[1:170880] # Need to truncate a few last-day readings
> dim(oarray) <- c(480,356)   # Re-dimension the vector to a 2-D array
> plot(oarray[1,], col="blue", type="l", main="4 p.m. outside temp",
+      xlab="Day of year (starting July 25, 2003)",
+      ylab="Temperature (Celsius)")
图 5. 下午 4 点的室外温度
4 p.m. outside temperature
4 p.m. outside temperature

一旦我们将向量转换为一个“Time X Day”矩阵(二维数组),就可以自然地抽取出每天的温度, 并绘出每年度的模式图。您可以通过其他方式完成 —— 从每个向量中提取出第 480 个点;R 的方式 比这好得多。

三维数据

如何描述全年的温度度量?一种方法是使用色彩编码的温度图:

清单 11. 创建一个温度图
> x <- 1:480                  # Create X axis indices
> y <- 1:356                  # Create Y axis indices
> z <- oarray[x,y]            # Define z-axis (really same as oarray)
> mycolors <- c(heat.colors(33),topo.colors(21))[54:1]
> image(x,y,z, col=mycolors,
+       main="Outside temperatures near Glarp", # Brad's name for his house
+       xlab="Minutes past 4 p.m",
+       ylab="Days past July 25, 2003" )
> dev2bitmap(file="outside-topo.pdf", type="pdfwrite")
图 6. 每日 x 每年 温度的温度图
每日 x 每年 温度的温度图
每日 x 每年 温度的温度图

稍微解释一下我们刚才的工作。定义数轴和端点,如果您识得创建一列数字的类 Python 片段符号, 这是相当显而易见的。在创建 z 的过程中通过使用 xy 进行注标创建了刻度宽和高的一个数组。在这个例子中, z 是微不足道的 —— oarray 也是如此;但是它可以足够简单地系统地修改值或者用于 访问它们的偏移量。颜色表 mycolor 是通过反复实验得到的:我们觉得红色和 黄色适于用来表示“热”的温度(也就是高于 0 摄氏度),而不适于表示冷的温度。所以,我们向颜色向量添加 一些蓝色/绿色。结果,我们还希望对标准颜色表函数生成的颜色进行倒序排列 —— 通过注标很容易实现。

您可能注意到了,温度图比先前的图形更为尖锐。使用 image() 命令在屏幕上绘制 出了一个适当的但给人印象并不深的图像。将“当前图像”导出到一个外部文件通常可以生成更好的结果,这里也是如此。

虽然我们个人喜欢先前的平面温度图,但是很多观察图形的人可能会通过伪透视三维数据而获得 更好的信息。在 R 中,要生成三维数据的极好的透视地形图(topographic maps),需要稍微多做一些工作。例如:

清单 12. 创建一个地形表面图
> persp(x,y,z, theta=10, phi=60, ltheta=40, lphi=30, shade=.1, border=NA,
+       col=mycolors[round(z+20)], d=.5,
+       main="Outside temperatures near Glarp",
+       xlab="Minutes past 4 p.m",
+       ylab="Days past July 25, 2003",
+       zlab="Temperature", )
图 7. 每日 x 每年 温度的地形图
每日 x 每年 温度的地形图
每日 x 每年 温度的地形图

结束语

我们在本文中进行的一些基本统计分析以及生成的基本图形,仅仅描述了 R 的统计能力的一个 小的子集。对于几乎科学的每一个领域,以及几乎所有众所周知的(甚至鲜为人知的)统计技术, 在 R 中都有支持相应的数学方法的标准函数或者扩展程序包。

我们希望通过本文能使您对 R 有一个初步了解,并喜欢使用它。不过,R 具有很多我们在这里没有 展示的功能。在这个由三部分构成的系列文章的第二部分,我们将深入研究 R 中的中级及高级技术 —— 首先是线性回归和非线性回归。


相关主题

  • 您可以参阅本文在 developerWorks 全球站点上的 英文原文
  • 查看 用于统计计算的 R 项目 主页。 这个 R Web 站点上有大量的包括从教程到完全的 API 描述在内的所有文档。对于 刚接触 R 的读者来说尤其会感兴趣的两个文档是:
  • 作为 Python 用户,David 的 Charming Python 专栏的一些读者对 R 的 Python 绑定表示出特别的喜好。实际上,有两种这样的绑定: RPy, 以及更老的 RSPython,这个也不错。(不过,David 的感觉是 RPy 绑定更好一些)。这两种绑定都让您可以通过 Python 代码透明地调用所有的 R 函数,使用 Python 对象作为函数的参数。
  • 不敢肯定的是,在大部分系统上,您都可以在 R 命令行中输入 help.start(), 以打开一个包含有为 R 文档生成的 HTML 页的浏览器。
  • 可以在线获得 S 和 S-Plus 历史摘要
  • Cameron Laird 的 服务器诊所:数据处理的利器 R (developerWorks,2003 年)给出了一个关于 R(以及 S)的很好概述, 并且有很多参考资料。
  • Gnuplot 也可以帮助您绘制和分析数据。Nishanth Sastry 在他的 自由控制高级图表和数据绘图 (developerWorks,2004 年)一文中介绍了 Gnuplot,并讨论了数据点丢失的问题。
  • 您可以练习转换和分析本文中 David 和 Brad 使用的温度数据(这里还有他们用来将初始的日志格式 转化为更好的制表格式的 Python 脚本):
  • developerWorks Linux 专区 可以找到 更多为 Linux 开发者准备的参考资料。
  • 在 Developer Bookstore Linux 专栏中定购 打折出售的 Linux 书籍

评论

添加或订阅评论,请先登录注册

static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=10
Zone=Linux
ArticleID=23200
ArticleTitle=使用 R 编写统计程序: 第 1 部分. 初涉大量统计工具
publish-date=10112004