当我们想研究不同 sample 的某个变量 A 之间的差异时,往往会因为其它一些变量 B 对该变量的固有影响,而影响不同 sample 变量 A 的比较,这个时候需要对 sample 变量 A 进行标准化之后才能进行比较。标准化的方法是对 sample 的 A 变量和 B 变量进行 loess 回归,拟合变量 A 关于变量 B 的函数 f(b),f(b) 则表示在 B 的影响下 A 的理论取值,A-f(B)(A 对 f(b)残差)就可以去掉 B 变量对 A 变量的影响, 此时残差值就可以作为标准化的 A 值在不同 sample 之间进行比较。
LOWESS 最初由 Cleveland 提出, 后又被 Cleveland&Devlin 及其他许多人发展。在 R 中 loess 函数是以 lowess 函数为基础的更复杂功能更强大的函数。主要思想为: 在数据集合的每一点用低维多项式拟合数据点的一个子集, 并估计该点附近自变量数据点所对应的因变量值, 该多项式是用加权最小二乘法来拟合; 离该点越远, 权重越小, 该点的回归函数值就是这个局部多项式来得到, 而用于加权最小二乘回归的数据子集是由最近邻方法确定。 最大优点: 不需要事先设定一个函数来对所有数据拟合一个模型。并且可以对同一数据进行多次不同的拟合,先对某个变量进行拟合,再对另一变量进行拟合,以探索数据中可能存在的某种关系,这是普通的回归拟合无法做到的。
1. 以 x0 为中心确定一个区间,区间的宽度可以灵活掌握。具体来说,区间的宽度取决于 q=fn。其中 q 是参与局部回归观察值的个数,f 是参加局部回归观察值的个数占观察值个数的比例, n 是观察值的个数。在实际应用中,往往先选定 f 值,再根据 f 和 n 确定 q 的取值,一般情况下 f 的取值在 1/3 到 2/3 之间。q 与 f 的取值一般没有确定的准则。增大 q 值或 f 值,会导致平滑值平滑程度增加,对于数据中前在的细微变化模式则分辨率低,但噪声小,而对数据中大的变化模式的表现则比较好;小的 q 值或 f 值,曲线粗糙,分辨率高,但噪声大。没有一个标准的 f 值,比较明智的做法是不断的调试比较。 2. 定义区间内所有点的权数,权数由权数函数来确定,比如立方加权函数 weight = (1 - (dist/maxdist)^3)^3),dist 为距离 x 的距离,maxdist 为区间内距离 x 的最大距离。任一点 (x0,y0) 的权数是权数函数曲线的高度。权数函数应包括以下三个方面特性:(1) 加权函数上的点 (x0,y0) 具有最大权数。(2) 当 x 离开 x0(时,权数逐渐减少。(3) 加权函数以 x0 为中心对称。 3. 对区间内的散点拟合一条曲线 y=f(x)。拟合的直线反映直线关系,接近 x0 的点在直线的拟合中起到主要的作用,区间外的点它们的权数为零。 4. x0 的平滑点就是 x0 在拟合出来的直线上的拟合点 (y0,f( x0))。 5. 对所有的点求出平滑点,将平滑点连接就得到 Loess 回归曲线。
- loess(formula, data, weights, subset, na.action, model = FALSE,
- span = 0.75, enp.target, degree = 2,
- parametric = FALSE, drop.square = FALSE, normalize = TRUE,
- family = c("gaussian", "symmetric"),
- method = c("loess", "model.frame"),
- control = loess.control(...), ...)
formula 是公式,比如 y~x, 可以输入 1 到 4 个变量; data 是放着变量的数据框,如果 data 为空,则在环境中寻找; na.action 指定对 NA 数据的处理,默认是 getOption("na.action"); model 是否返回模型框; span 是 alpha 参数,可以控制平滑度, 相当于上面所述的 f,对于 alpha 小于 1 的时候,区间包含 alpha 的点,加权函数为立方加权,大于 1 时,使用所有的点,最大距离为 alpha^(1/p),p 为解释变量; anp.target,定义 span 的备选方法; normalize,对多变量 normalize 到同一 scale; family,如果是 gaussian 则使用最小二乘法,如果是 symmetric 则使用双权函数进行再下降的 M 估计; method,是适应模型或者仅仅提取模型框架; control 进一步更高级的控制,使用 loess.control 的参数; 其它参数请自己参见 manual 并且查找资料
- loess.control(surface = c("interpolate", "direct"),
- statistics = c("approximate", "exact"),
- trace.hat = c("exact", "approximate"),
- cell = 0.2, iterations = 4, ...)
surface,拟合表面是从 kd 数进行插值还是进行精确计算; statistics, 统计数据是精确计算还是近似,精确计算很慢 trace.hat, 要跟踪的平滑的矩阵精确计算或近似?建议使用超过 1000 个数据点逼近, cell, 如果通过 kd 树最大的点进行插值的近似。大于 cell floor(nspancell) 的点被细分。 robust fitting 使用的迭代次数。
- predict(object, newdata = NULL, se = FALSE,
- na.action = na.pass, ...)
object,使用 loess 拟合出来的对象; newdata, 可选数据框,在里面寻找变量并进行预测; se, 是否计算标准误差; 对 NA 值的处理
生物数据分析中,我们想查看 PCR 扩增出来的扩增子的测序深度之间的差异,但不同的扩增子的扩增效率受到 GC 含量的影响,因此我们首先应该排除掉 GC 含量对扩增子深度的影响。
amplicon 测序数据,处理后得到的每个 amplicon 的深度,每个 amplicon 的 GC 含量,每个 amplicon 的长度 先用 loess 进行曲线的拟合
- gcCount.loess <- loess(log(RC+0.01)~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
画出拟合出来的曲线
- predictions1<- predict (gcCount.loess,RC_DT$GC)
- #plot scatter and line
- plot(RC_DT$GC,log(RC_DT$RC+0.01),cex=0.1,xlab="GC Content",ylab=expression(paste("log(NRC"["lib"],"+0.01)",sep="")))
- lines(RC_DT$GC,predictions1,col = "red")
取残差,去除 GC 含量对深度的影响
- #sustract the influence of GC
- resi <- log(RC_DT$RC+0.01)-predictions1
- RC_DT$RC <- resi
- setkey(RC_DT,GC)
此时 RC_DT$RC 就是 normalize 之后的 RC 画图显示 nomalize 之后的 RC, 并将拟合的 loess 曲线和 normalize 之后的数据保存
- #plot scatter and line using Norm GC data
- plot(RC_DT$GC,RC_DT$RC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))
- gcCount.loess <- loess(RC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
- save(gcCount.loess,file="/home/ywliao/project/Gengyan/gcCount.loess.Robject")
- predictions2 <- predict(gcCount.loess,RC_DT$GC)
- lines(RC_DT$GC,predictions2,col="red")
- save(RC_DT,file="/home/ywliao/project/Gengyan/RC_DT.Rdata")
当然,也想看一下 amplicon 长度 len 对 RC 的影响,不过影响不大
全部代码如下:
- library(data.table)
- load("/home/ywliao/project/Gengyan/RC_DT.Rdata")
- ####loess GC vs RC####
- gcCount.loess <- loess(log(RC+0.01)~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
- predictions1<- predict (gcCount.loess,RC_DT$GC)
- #plot scatter and line
- plot(RC_DT$GC,log(RC_DT$RC+0.01),cex=0.1,xlab="GC Content",ylab=expression(paste("log(NRC"["lib"],"+0.01)",sep="")))
- lines(RC_DT$GC,predictions1,col = "red")
- #sustract the influence of GC
- resi <- log(RC_DT$RC+0.01)-predictions1
- RC_DT$RC <- resi
- setkey(RC_DT,GC)
- #plot scatter and line using Norm GC data
- plot(RC_DT$GC,RC_DT$RC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))
- gcCount.loess <- loess(RC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
- save(gcCount.loess,file="/home/ywliao/project/Gengyan/gcCount.loess.Robject")
- predictions2 <- predict(gcCount.loess,RC_DT$GC)
- lines(RC_DT$GC,predictions2,col="red")
- save(RC_DT,file="/home/ywliao/project/Gengyan/RC_DT.Rdata")
- ####loess len vs RC###
- setkey(RC_DT,Len)
- len.loess <- loess(RC_DT$RC~RC_DT$Len, control = loess.control(surface = "direct"),degree=2)
- predictions2<- predict (len.loess,RC_DT$Len)
- #plot scatter and line
- plot(RC_DT$Len,RC_DT$RC,cex=0.1,xlab="Length",ylab=expression(paste("log(RC"["GC"],"+0.01)",sep="")))
- lines(RC_DT$Len,predictions2,col = "red")
来源: http://www.cnblogs.com/ywliao/p/6863734.html