欢迎您光临bv1946伟德体育官方网站!

R语言通过loess去除某个变量对数据的影响

时间:2019-11-28 17:01

奥迪Q5语言代码

 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"卡塔尔(英语:State of Qatar);
  model是不是重临模型框;
  span是阿尔法参数,可以垄断平滑度,也就是地点所述的f,对于阿尔法小于1的时候,区间包含阿尔法的点,加权函数为立方加权,大于1时,使用具备的点,最大间隔为阿尔法^(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, ...)

  平板电脑,拟合表面是从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值的管理

Loess局地加权多项式回归

  LOWESS最早由Cleveland 提议,后又被Cleveland&Devlin及任何不菲人升高。在途睿欧中loess 函数是以lowess函数为根基的更头晕目眩功用更加强硬的函数。首要理念为:在数额群集的每一点用低维多项式拟合数分部的三个子集,并推断该点周围自变量数总部所对应的因变量值,该多项式是用加权最小二乘法来拟合;离该点越远,权重越小,该点的回归函数值正是那么些片段多项式来赢得,而用于加权最小二乘回归的多少子集是由多年来邻方法明确。
  最大优点:不需求事前设定二个函数来对具备数据拟合贰个模型。何况能够对同风流倜傥数据开展一再不等的拟合,先对有个别变量举行拟合,再对另生龙活虎变量实行拟合,以商讨数据中大概存在的某种关联,那是平淡无奇的回归拟合相当小概做到的。

数据

amplicon 测序数据,处理后获得的各样amplicon的纵深,各种amplicon的GC含量,各种amplicon的尺寸
图片 1
先用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")

图片 2

取残差,去除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")

图片 3

本来,也想看一下amplicon 长度len 对RC的熏陶,可是影响相当的小
图片 4

漫天代码如下(经过改变,可能与地方完全相称卡塔尔(قطر‎:

library(data.table)

load("/home/ywliao/project/Gengyan/RC_DT.Rdata")
RRC_DT <- RC_DT[Type=="WBC" & !is.na(RC),]

lst <- list()
for (Samp in unique(RC_DT$Sample)){
RC_DT <- RRC_DT[Sample==Samp]
####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$NRC <- resi
setkey(RC_DT,GC)
#plot scatter and line using Norm GC data
#plot(RC_DT$GC,RC_DT$NRC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))
gcCount.loess <- loess(NRC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
predictions2 <- predict(gcCount.loess,RC_DT$GC)
#lines(RC_DT$GC,predictions2,col="red")
lst[[Samp]] <- RC_DT
}
NRC_DT <- rbindlist(lst)
save(RC_DT,file="/home/ywliao/project/Gengyan/NRC_DT.Rdata")

####loess len vs RC###
setkey(RC_DT,Len)
len.loess <- loess(RC_DT$NRC~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$NRC,cex=0.1,xlab="Length",ylab=expression(paste("log(RC"["GC"],"+0.01)",sep="")))
lines(RC_DT$Len,predictions2,col = "red")

实例

  生物数据解析中,大家想查看PC奥迪Q3扩大与扩大出来的扩大与扩大子的测序深度以内的差别,但不一致的扩大与扩充子的扩大与扩张效用受到GC含量的熏陶,因而我们首先应该裁撤掉GC含量对扩大与扩大子深度的影响。

  当大家想切磋差异sample的某些变量A之间的差距时,往往会因为其余一些变量B对该变量的固有影响,而影响不一致sample变量A的相比较,这时需求对sample变量A实行规范之后技巧张开相比。标准化的章程是对sample 的 A变量和B变量举办loess回归,拟合变量A关于变量B的函数 f(b卡塔尔国,f(b卡塔尔(قطر‎则意味在B的熏陶下A的申辩取值,A-f(B卡塔尔(英语:State of Qatar)(A对f(b)残差)就能够去掉B变量对A变量的影响,这时候残差值就足以视作典型的A值在不一样sample之间开展比较。

LOESS平滑方法

  1. 以x0为着力分明多个间距,区间的大幅能够灵活领会。具体来讲,区间的宽度决议于q=fn。当中q是参预部分回归观望值的个数,f是在座一些回归观察值的个数占观望值个数的比例,n是阅览值的个数。在事实上行使中,往往先选定f值,再依照f和n明确q的取值,日常情状下f的取值在45%到2/3里面。q与f的取值日常从不规定的守则。增大q值或f值,会促成平滑值平滑程度增添,对于数据中前在的细微变化格局则分辨率低,但噪声小,而对数据中山高校的浮动情势的变现则相比好;小的q值或f值,曲线粗糙,分辨率高,但噪声大。未有叁个职业的f值,相比较明智的做法是再三的调理比较。
  2. 概念区间内全体一点点的权数,权数由权数函数来规定,比方立方加权函数weight = (1 - (dist/maxdist卡塔尔(英语:State of Qatar)^3卡塔尔(قطر‎^3卡塔尔国,dist为间隔x的偏离,maxdist为间隔内间隔x的最大间距。任一点(x0,y0卡塔尔国的权数是权数函数曲线的莫斯中国科学技术大学学。权数函数应满含以下三个地点特征:(1卡塔尔加权函数上的点(x0,y0卡塔尔具备最大权数。(2卡塔尔(英语:State of Qatar)当x离开x0(时,权数渐渐回退。(3卡塔尔(英语:State of Qatar)加权函数以x0为核心对称。
  3. 对区间内的散点拟合一条曲线y=f(x卡塔尔。拟合的直线反映直线关系,临近x0的点在直线的拟合中起到首要的效能,区间外的点它们的权数为零。
  4. x0的平滑点正是x0在拟合出来的直线上的拟合点(y0,f( x0))。
  5. 对全部的点求出平滑点,将平滑点连接就得到Loess回归曲线。

上一篇:九一八后张学良还有机会翻盘 不听蒋介石命令自毁长城
下一篇:没有了