用局部加权回归散点平滑法观察二维变量之间的关系
局部加权回归散点平滑法
二维变量之间的关系研究是很多统计方法的基础，例如回归分析通常会从一元回归讲起，然后再扩展到多元情况。局部加权回归散点平滑法（locally weighted scatterplot smoothing，LOWESS或LOESS）是查看二维变量之间关系的一种有力工具。
LOWESS主要思想是取一定比例的局部数据，在这部分子集中拟合多项式回归曲线，这样我们便可以观察到数据在局部展现出来的规律和趋势；而通常的回归分析往往是根据全体数据建模，这样可以描述整体趋势，但现实生活中规律不总是（或者很少是）教科书上告诉我们的一条直线。我们将局部范围从左往右依次推进，最终一条连续的曲线就被计算出来了。显然，曲线的光滑程度与我们选取数据比例有关：比例越少，拟合越不光滑（因为过于看重局部性质），反之越光滑。
本文的数据文件：物种数目与海拔高度（感谢中科院植物所赖江山博士提供数据并授权使用）
R程序代码：
# 从本站counts.txt文件直接将数据读入R
x = read.csv("http://cos.name/wp-content/uploads/2008/11/counts.txt")
par(las = 1, mar = c(4, 4, 0.1, 0.1))
plot(x, pch = 20, col = rgb(0, 0, 0, 0.5))
# 取不同的f参数值
for (i in seq(0.01, 1, length = 100)) {
lines(lowess(x$altitude, x$counts, f = i), col = gray(i),
lwd = 1.5)
Sys.sleep(0.15)
}
以上Sys.sleep()语句只是为了让读者看清楚添加LOWESS曲线的过程，实际画图过程中可以去掉。以上代码生成的图形如下：
局部加权回归散点平滑法
上图中，曲线颜色越浅表示所取数据比例越大。不难看出白色的曲线几乎已呈直线状，而黑色的线则波动较大。总体看来，图中大致有四处海拔上的物种数目偏离回归直线较严重：450米、550米、650米和700米附近。若研究者的问题是，多高海拔处的物种数最多？那么答案应该是在650米附近。如果仅仅从回归直线来看，似乎是海拔越高，则物种数目越多。如此推断下去，恐怕月球或火星上该物种最多。以下是回归直线的图示：
par(las = 1, mar = c(4, 4, 0.1, 0.1), mgp = c(2.5,
1, 0))
plot(x, pch = 20, col = rgb(0, 0, 0, 0.5))
abline(lm(counts ~ altitude, x), lwd = 2, col = "red")
物种数目与海拔高度的关系：回归直线
为了确保我们用LOWESS方法得到的趋势是稳定的，我们可以进一步用Bootstrap的方法验证。因为Bootstrap方法是对原样本进行重抽样，根据抽得的不同样本可以得到不同的LOWESS曲线，最后我们把所有的曲线添加到图中，看所取样本不同是否会使得LOWESS有显著变化；以下是R代码：
set.seed(711) # 设定随机数种子，保证本图形可以重制
par(las = 1, mar = c(4, 4, 0.1, 0.1), mgp = c(2.5,
1, 0))
plot(x, pch = 20, col = rgb(0, 0, 0, 0.5))
for (i in 1:400) {
idx = sample(nrow(x), 300, TRUE) # 有放回抽取300个样本序号
lines(lowess(x$altitude[idx], x$counts[idx]), col = rgb(0,
0, 0, 0.05), lwd = 1.5) # 用半透明颜色，避免线条重叠使得图形看不清
Sys.sleep(0.05)
}
dev.off()
生成图形如下：
物种数目与海拔高度的关系：Bootstrap结合LOWESS查看
可以看出，经过400次重抽样并计算LOWESS曲线，刚才在第一幅图中观察到的趋势大致都还存在（因为默认取数据比例为2/3，因此拟合曲线都比较光滑），只是700米海拔附近物种数目减小的趋势并不明显了，这是因为这个海拔附近的观测样本量较少，在重抽样的时候不容易被抽到，因此在图中代表性不足，最后得到的拟合曲线分布稀疏。
作者注：只是一副散点图而已，能做的文章并不少。本文是基于赖博士的另外一个问题而引发出来的思考，供生物与生态专业的同仁参考。值此新建站点之际，谨以此文抛砖，望能引来更多高人对COS网站贡献的“美玉”。作者联系方式：xie[at]yihui.name
[...] 用局部加权回归散点平滑法观察二维变量之间的关系 [...]
局部加权回归散点平滑图示…
在Frank Harrell强调LOWESS的重要性之前，我没有意识到它的作用。从德国回来之后，在John Maindonald的书上再一次看到了用LOWESS观察散点图的例子，这下子才渐渐有了感觉。统计就是以糟蹋数据为…
就technique background 和 application来说讲得基本到位哈。但是我觉得对于统计方法的描述上，如果能结合背景知识来讲的话效果会更好。因为数据只是数据，而统计模型就是因数据而生，看上去任何一种方法都很好。但是实际上，统计模型具有很浓厚的背景性。比如lowess方法，它有一个hidden assumption就是，你必须知道你的这个smooth regression model是实际而且可以intepret的。（这个东西本身就是从生物，医学或者社会学领域发展来的，人家已经知道两个变量之间是有关系的，所以才发展了这些非线性方法去建模）。可想而知，这样的模型应该是比简单的线性要好。但是，如果在实际应用当中这个前提不能满足的话，这个模型的uncertainty就大大增加了。很显然，参数越多，模型看上去对数据的拟合越好，那么你怎么去计算你得到的这个relationship的可靠度呢？（type-i and type-ii errors）？ 这一直是统计学界最担心的问题，统计模型的解释和应用。比如，chaos理论里有大量的简单的dynamic system可以模拟出非常随机的数据，而很多统计模型可以把这些数据拟合得非常好，但是事实上这些模型都是错的。：（
嗯，确实应该请给我数据的人来评价一下，统计学方法只有和相关专业知识结合起来才有用，而不能闭门造车还觉得自己造的挺好看
谢老师，能否将LOESS拟合的脚本发一份给我？谢谢！
见正文R代码。数据、程序都有。
yihui，你好。
我看到有些文献中通过cox regression计算martingale residuals，随后通过绘制mr的scatter plot，计算和0交界的值来确定cut off。
能介绍一下原理吗？
抱歉我已经忘了martingale residuals的意思，但我感觉应该是观察残差的趋势（普通的散点图可能不容易观察）
谢老师，谢谢你的指点，还有一个问题想请教下，如果loess拟合只要得到第二个图那样的只有点位和曲线的图，该怎样写代码？请原谅我的无知。多谢……
看正文R代码，依葫芦画瓢吧。
谢老师，你好：
又来打搅你了，拟合曲线我画出来了，现在我需要在曲线上找到breakpoint，但不知道怎么用rpart函数，能不能就本例指点下？谢谢.
什么叫breakpoint？如何定义？回归树中自变量的拆分点？看
?rpart及其返回值及其相应的帮助文档吧，关键是了解返回结果的数据对象类型以及怎么从中取值出来。益辉，呵呵，我来转转，这些东西我都不懂，不过能看出来这个网站真的很棒，你很用心，辛苦了！yujie
哈哈，你是咱班第一个来这里吼一嗓子的，谢谢支持！