R語言通過loess去除某個變數對數據的影響

當我們想研究不同sample的某個變數A之間的差異時，往往會因為其它一些變數B對該變數的固有影響，而影響不同sample變數A的比較，這個時候需要對sample變數A進行標準化之後才能進行比較。標準化的方法是對sample 的 A變數和B變數進行loess回歸，擬合變數A關於 ...

當我們想研究不同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之間進行比較。

Loess局部加權多項式回歸

LOWESS最初由Cleveland 提出,後又被Cleveland&Devlin及其他許多人發展。在R中loess 函數是以lowess函數為基礎的更複雜功能更強大的函數。主要思想為:在數據集合的每一點用低維多項式擬合數據點的一個子集,並估計該點附近自變數數據點所對應的因變數值,該多項式是用加權最小二乘法來擬合;離該點越遠,權重越小,該點的回歸函數值就是這個局部多項式來得到,而用於加權最小二乘回歸的數據子集是由最近鄰方法確定。
最大優點:不需要事先設定一個函數來對所有數據擬合一個模型。並且可以對同一數據進行多次不同的擬合，先對某個變數進行擬合，再對另一變數進行擬合，以探索數據中可能存在的某種關係，這是普通的回歸擬合無法做到的。

LOESS平滑方法

1. 以x₀為中心確定一個區間，區間的寬度可以靈活掌握。具體來說，區間的寬度取決於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的最大距離。任一點(x₀,y₀)的權數是權數函數曲線的高度。權數函數應包括以下三個方面特性：(1)加權函數上的點(x₀,y₀)具有最大權數。(2)當x離開x₀₍時，權數逐漸減少。(3)加權函數以x₀為中心對稱。
3. 對區間內的散點擬合一條曲線y=f(x)。擬合的直線反映直線關係，接近x₀的點在直線的擬合中起到主要的作用，區間外的點它們的權數為零。
4. x₀的平滑點就是x₀在擬合出來的直線上的擬合點(y₀,f( x₀))。
5. 對所有的點求出平滑點，將平滑點連接就得到Loess回歸曲線。

R語言代碼

 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")