smoothing是一种很有力的数据分析手段，通常应用于评估分布特征未知但假定光滑的数据。一般的做法是先将具有相似期望值的数据点进行分组，然后计算均值或者拟合一个简单的参数模型。本章节，作者还是以一个基因表达数据为例，对bin smoothing和loess这两种方法进行讲解。
讲解这两种方法之前，作者以线性回归作为切入点进行讲解:
library(Biobase)
library("SpikeIn")
library(hgu95acdf)
data(SpikeIn95)
library(rafalib)
##Example with two columns
i=10;j=9
##remove the spiked in genes and take random sample
siNames<-colnames(pData(SpikeIn95))
ind <- which(!probeNames(SpikeIn95)%in%siNames)
pms <- pm(SpikeIn95)[ ind ,c(i,j)]
##pick a representative sample for A and order A
Y=log2(pms[,1])-log2(pms[,2])
X=(log2(pms[,1])+log2(pms[,2]))/2
set.seed(4)
ind <- tapply(seq(along=X),round(X*5),function(i)
if(length(i)>20) return(sample(i,20)) else return(NULL))
ind <- unlist(ind)
X <- X[ind]
Y <- Y[ind]
o <-order(X)
X <- X[o]
Y <- Y[o]
mypar()
plot(X,Y)
202001071726.png
作者选择了两个重复样本，然后绘制MA图（横坐标是均数，纵坐标是差值）。因此理论上，无论取何值
都应该是0才对，但看图显然不是这个样子,换句话说由于
对
存在依赖性，所以使得分析产生了偏差。
我们可以通过预测给定值的条件下
的期望值这种手段去除这种偏差(条件期望)，也就是：
，而线性回归在这种情况下并不适用：
mypar()
plot(X,Y)
fit <- lm(Y~X)
points(X,Y,pch=21,bg=ifelse(Y>fit$fitted,1,3))
abline(fit,col=2,lwd=4,lty=2)
202001071739.png
很显然，位于回归直线上面的绿色数据点和下面的紫色数据点的分布并不均匀，因此需要使用更加灵活地方法去分析。
Bin Smoothing
所谓Bin Smoothing其实就是将数据分层然后计算均值，这种方法的大致逻辑就是如果曲线足够光滑，那么在一个很小的区域内的那一段曲线大致是恒定的（原文用的词是constant），也就是那一段区间内的值均相等，举例说明：
windowSize<-0.5
centers <- seq(min(X),max(X),0.1)
smooth<-rep(NA,length(centers))
mypar (4,3)
for(i in seq(along=centers)){
center<-centers[I]
ind=which(X>center-windowSize & X<center+windowSize)
smooth[i]<-mean(Y[ind])
if(i%%round(length(centers)/12)==1){ ##we show 12
plot(X,Y,col="grey",pch=16)
points(X[ind],Y[ind],bg=3,pch=21)
lines(c(min(X[ind]),max(X[ind])),c(smooth[i],smooth[i]),col=2,lwd=2)
lines(centers[1:i],smooth[1:i],col="black")
points(centers[i],smooth[i],col="green",pch=16,cex=1.5)
}
}
202001081429.png
Bin smoothing的具体做法是，首先设置bin的大小，本例是1，然后均匀选定中心点，最后以每个中心点左右0.5范围内的点为单位计算的拟合值（也就是均值），以新的Y值重新描点就能得到上图所示的黑色曲线。
这种通过拟合一个常数去画曲线的方法应该算是最简单的一种处理，出来的曲线也通常是弯弯曲曲的，如loess的方法则可对这一问题做出改善。
Loess（local weighted regression）
Loess与bin smoothing的最大区别就是它在局部用的是直线或者抛物线，所以这种方法选定的bin宽度肯定要比bin smoothing大，还是通过举例进行说明：
mypar (4,3)
windowSize<-1.25
smooth<-rep(NA,length(centers))
for(i in seq(along=centers)){
center<-centers[I]
ind=which(X>center-windowSize & X<center+windowSize)
fit<-lm(Y~X,subset=ind)
smooth[i]<-fit$coef[1]+fit$coef[2]*center
if(i%%round(length(centers)/12)==1){ ##we show 12
plot(X,Y,col="grey",pch=16)
points(X[ind],Y[ind],bg=3,pch=21)
a <- min(X[ind]);b <- max(X[ind])
lines(c(a,b),fit$coef[1]+fit$coef[2]*c(a,b),col=2,lwd=2)
lines(centers[1:i],smooth[1:i],col="black")
points(centers[i],smooth[i],col="green",pch=16,cex=1.5)
}
}
202001081459.png
仔细阅读代码后你就可以发现，其实这两趴最大的区别就是对于Y值的拟合方法上，一种是以bin中点的Y值的均值作为拟合值，一种是用线性回归的方法去获取Y的拟合值，显然后一种方法得到的曲线更加平滑。
这种曲线平滑的方法直接用R的一个函数loess()就可以实现：
fit <- loess(Y~X,degree = 1,span=1/3)
newx <- seq(min(X),max(X),len=100)
##predict函数可以给定x的值然后预测y值
smooth <- predict(fit,newdata=data.frame(X=newx))
mypar()
plot(X,Y,col="darkgrey",pch=16)
lines(newx,smooth,col="black",lwd=3)
202001081512.png
最后，作者讲解了几点loess与bin smoothing之间的几个主要区别：
- bin smoothing是保持窗口大小一定，而loess则是保持窗口内用于拟合的点的数目一定；
- 当使用参数模型拟合Y值时，loess会使用加权最小二乘法，也就是说离中心点越近的点所占权重越大；
- loess在局部拟合这方面可以实现鲁棒性，只要设置参数
family="symmetirc"即可，这一参数的作用就是在每次拟合时会对离群值进行检测，然后在下一次拟合时降低该点的权重。
阅读原文请戳
网友评论