微信公眾號：生物信息學習

筆者邀請您，先思考：
1 線性回歸是什么？
2 線性回歸怎么應(yīng)用？

本文解釋了如何在R中運行線性回歸。本教程將介紹線性回歸的假設(shè)以及如果假設(shè)不滿足如何處理。 它還包括擬合模型和計算模型性能指標以檢查線性回歸模型的性能。 線性回歸是最流行的統(tǒng)計技術(shù)之一。 它已被使用了三十多年。 它幾乎在每個領(lǐng)域都被廣泛接受，因為它很容易理解線性回歸的輸出。

線性回歸

它是一種發(fā)現(xiàn)稱為一個連續(xù)的因變量或者目標變量與一個或者多個（連續(xù)或者不連續(xù)）的自變量之間的關(guān)系的方法論。

image

這是一條直線曲線。 在上圖中，對角紅線是一條回歸線，也稱為最佳擬合直線。 點與回歸線之間的距離是誤差（errors）。 線性回歸旨在通過最小化點與回歸線之間的垂直距離的平方和來找到最佳擬合直線。

變量類型

線性回歸需要因變量是連續(xù)的，即數(shù)值（無類別或組）。

簡單與元線性回歸

當只有一個自變量時，線性回歸是簡單的線性回歸（一元線性回歸）。 而多元線性回歸會有多個自變量。

回歸方程

image

解釋：

當所有自變量變量（Xs）等于0且時，b0是因變量（Y）的期望平均值的截距。b1是斜率，它代表因變量（Y）的變化量，如果我們改變X1一個單位保持其他變量不變。

重要術(shù)語：殘差

觀察到的（實際）因變量值與從回歸線預(yù)測的因變量值之間的差異。

算法

線性回歸基于最小二乘估計，其表示回歸系數(shù)（估計值）應(yīng)該以使每個觀察到的響應(yīng)與其擬合值的平方距離的總和最小化的方式來選擇。

最小樣本量

線性回歸分析中需要每個自變量5個案例。

線性回歸分析的假設(shè)

1.線性關(guān)系：線性回歸需要依賴因變量和自變量之間的線性關(guān)系。

2.殘差的正態(tài)性：線性回歸需要殘差應(yīng)該正態(tài)分布。

3.同方差性：線性回歸假定所有預(yù)測的因變量值的殘差大致相等。換句話說，它意味著誤差方差的不變性。

4.沒有離群值問題

5.多重共線性：這意味著自變量之間有很高的相關(guān)性。線性回歸模型不得面臨多重共線性問題。

6.誤差的獨立性 - 無自相關(guān)

它指出，與一個觀測相關(guān)的誤差與任何其他觀測的誤差都不相關(guān)。當您使用時間序列數(shù)據(jù)時，這是一個問題。假設(shè)你從八個不同的地區(qū)收集了勞動者的數(shù)據(jù)。每個地區(qū)內(nèi)的勞動力可能更傾向于來自不同地區(qū)的勞動力，也就是說，他們的錯誤不是獨立的。

線性回歸的分布

線性回歸假設(shè)目標或因變量是正態(tài)分布的。 正態(tài)分布與高斯分布相同。 它使用高斯家族的連接函數(shù)。

標準化系數(shù)

標準化或標準化系數(shù)（估計）的概念出現(xiàn)在預(yù)測變量（又稱自變量）以不同單位表示時。 假設(shè)你有3個獨立變量 - 年齡，身高和體重。 變量“年齡”以年表示，身高以厘米為單位，體重以千克為單位。如果我們需要根據(jù)非標準化系數(shù)對這些預(yù)測因子進行排序，那么它們將不會是一個公平的比較，因為這些變量的單位并不相同。

標準化系數(shù)（或估計）主要用于排列預(yù)測變量（或自變量或解釋變量），因為它消除了自變量和因變量的測量單位）。 我們可以用標準化系數(shù)的絕對值對自變量進行排序。 最重要的變量將具有最大的標準化系數(shù)絕對值。

標準化系數(shù)的解釋

1.25的標準化系數(shù)值表示自變量中一個標準偏差的變化導(dǎo)致因變量增加1.25個標準偏差。

模型性能測量
1. R平方

它用模型中的所有自變量來解釋因變量的方差比例。 它假定模型中的每個獨立變量都有助于解釋因變量的變化。 事實上，有些變量不會影響因變量，它們也無助于建立一個好的模型。

image

上述方程式的分子，yi-hat是預(yù)測值。 Y的平均值出現(xiàn)在分母中。

規(guī)則：

R平方越高，模型越適合您的數(shù)據(jù)。 在心理調(diào)查或研究中，我們通常發(fā)現(xiàn)低R平方值低于0.5。 這是因為我們試圖預(yù)測人類行為，預(yù)測人類并不容易。 在這些情況下，如果您的R平方值很低，但您具有統(tǒng)計學上顯著的獨立變量（又稱預(yù)測變量），您仍然可以生成關(guān)于預(yù)測變量值中的變化如何與響應(yīng)值變化相關(guān)聯(lián)的見解。

R平方可否為負？

是的，當水平線比您的模型更好地解釋數(shù)據(jù)時。 它主要發(fā)生在你不包括截距的情況下。 沒有截距，在預(yù)測目標變量方面，回歸可能會比樣本均值差。 這不僅是因為沒有截距。 即使包含截距，它也可能是負的。

在數(shù)學上，當模型的誤差平方大于水平線上的總平方和時，這是可能的。

R-squared = 1 - [(Sum of Square Error)/(Total Sum of Square)]

2.調(diào)整R平方

它衡量的只是那些真正影響因變量的自變量所解釋的方差比例。 它懲罰你添加不影響因變量的自變量。

image

調(diào)整后的R-Squared比R-squared更重要

每次向模型添加自變量時，即使自變量不顯著，R平方也會增加。 它永不衰落。 而調(diào)整R平方僅在自變量顯著并且影響因變量時增加。

3. RMSE（均方根誤差）

它解釋了實際數(shù)據(jù)點與模型預(yù)測值的接近程度。 它測量殘差的標準差。

image

在上面的公式中，yi是因變量的實際值，yi-hat是預(yù)測值，n是樣本大小。

很重要的一點

RMSE與因變量具有相同的單位。 例如，因變量是以美元計量的銷售額。 假設(shè)這種銷售模式的RMSE出來21.我們可以說它是21美元。

什么是好的RMSE分數(shù)？

沒有關(guān)于好或壞RMSE評分的拇指規(guī)則。 這是因為它依賴于你的因變量。如果您的目標變量介于0到100之間。RMS的0.5可以認為是好的，但同樣的0.5 RMSE可以被認為是差分，如果因變量的范圍從0到10.因此，通過簡單地查看 在價值。

較低的RMSE值表示更合適。 RMSE是衡量模型如何準確預(yù)測響應(yīng)的一個很好的衡量標準，如果模型的主要目的是預(yù)測，那么它就是擬合最重要的標準。

如果您已經(jīng)建立了良好的模型，那么您的訓練和測試集的RMSE應(yīng)該非常相似。 如果測試集的RMSE遠高于訓練集的RMSE，那么很可能是您的數(shù)據(jù)嚴重過擬合，即您創(chuàng)建了一個在樣本中運行良好的模型，但在樣本外的測試集上幾乎沒有預(yù)測價值。

RMSE與MAE

與平均絕對誤差（MAE）相比，RMSE嚴格懲罰大的誤差。

R平方與RMSE

R平方是個比例并且沒有與目標變量相關(guān)聯(lián)的單位，而RMSE具有與目標變量相關(guān)的單位。 因此，R平方是擬合的相對度量，RMSE是擬合的絕對度量。

下面的代碼涵蓋了對模型性能的假設(shè)測試和評估：

1. 數(shù)據(jù)準備
2. 多重共線性檢驗
3. 多重共線性的處理
4. 檢查自相關(guān)
5. 檢查異常值
6. 檢查異方差
7. 殘差的正態(tài)性檢驗
8. 前進，后退和逐步選擇
9. 計算RMSE
10. 因變量的Box Cox變換
11. 手動計算R平方和調(diào)整R平方
12. 計算殘差和預(yù)測值
13. 計算標準化系數(shù)

R腳本：線性回歸

理論部分結(jié)束了。 讓我們用R -

加載所需的包

 library(ggplot2)
 library(car)
 library(caret)
 library(corrplot)

確保上面列出的軟件包已經(jīng)安裝并加載到R中。如果它們尚未安裝，則需要使用命令install.packages（“package-name”）進行安裝。

讀取數(shù)據(jù)

我們將使用cars包中的mtcars數(shù)據(jù)集。 這些數(shù)據(jù)摘自Motor Trend US雜志，包括32種汽車的燃料消耗和10個方面的汽車設(shè)計和性能。

 data(mtcars)
 str(mtcars)

image

以上變量的變量描述將在下面列出它們各自的變量名稱。

image

數(shù)據(jù)摘要

在這個數(shù)據(jù)集中，mpg是一個目標變量。 使用head（）函數(shù)查看前6行數(shù)據(jù)。

 head(mtcars)

image

要查看變量的分布情況，請使用summary（）函數(shù)。

  summary(mtcars)

image

數(shù)據(jù)準備

確保分類變量存儲為因子。 在下面的程序中，我們將變量轉(zhuǎn)換為因子。

  mtcars$am <- as.factor(mtcars$am)
  mtcars$cyl <- as.factor(mtcars$cyl)
  mtcars$vs <- as.factor(mtcars$vs)
  mtcars$gear <- as.factor(mtcars$gear)

識別和修正共線性

在這一步中，我們正在識別彼此高度相關(guān)的自變量。 由于mpg是一個因變量，我們將在下面的代碼中刪除它。

  mtcars_a <- subset(mtcars, select = -c(mpg))
  numericData <- mtcars_a[sapply(mtcars_a, is.numeric)]
  descrCor <- cor(numericData)
  print(descrCor)

image

 corrplot(
 descrCor,
 order = "FPC",
 method = "color",
 type = "lower",
 tl.cex = 0.7,
 tl.col = rgb(0,0,0)
 )

image

  highlyCorrelated <- findCorrelation(descrCor, cutoff=0.7)
  highlyCorCol <- colnames(numericData)[highlyCorrelated]
  dat3 <- mtcars[, -which(colnames(mtcars) %in% highlyCorCol)]
  dim(dat3)

有三個變量“hp”“disp”“wt”被發(fā)現(xiàn)是高度相關(guān)的。 我們已經(jīng)刪除它們以避免共線。 現(xiàn)在，我們有7個獨立變量和1個因變量。

開發(fā)回歸模型

在這一步，我們正在建立多元線性回歸模型。

  fit <- lm(mpg ~ ., data=dat3)
  summary(fit)
  summary(fit)$coeff
  anova(fit)
  par(mfrow=c(2,2))
  plot(fit)

查看線性回歸模型的系數(shù)和ANOVA表

性回歸模型測試估計等于零的零假設(shè)。 具有小于0.05的p值的獨立變量意味著我們拒絕5％顯著性水平的零假設(shè)。這意味著該變量的系數(shù)不等于0.大p值意味著變量對預(yù)測目標變量沒有意義。

  summary(fit)
  anova(fit)

線性回歸模型的關(guān)鍵圖如下所示 :-

• 殘留物與擬合值
• 正常Q-Q
• 縮放位置
• 殘差與杠桿

image

計算模型性能度量

  summary(fit)$r.squared
  summary(fit)$adj.r.squared
  AIC(fit)
  BIC(fit)

更高的R平方和調(diào)整的R平方值，更好的模型。 然而，更低的AIC和BIC得分，更好的模型。

理解AIC和BIC

AIC和BIC是擬合度的衡量標準。 他們懲罰復(fù)雜的模型。 換句話說，它會懲罰更多的估計參數(shù)。 它相信一個概念，即一個具有較少參數(shù)的模型將被一個具有更多參數(shù)的模型要好。 一般來說，BIC比AIC更為免費參數(shù)懲罰模型。 兩個標準都取決于估計模型的似然函數(shù)L的最大值。

AIC值大致等于參數(shù)的數(shù)量減去整個模型的似然性。 假設(shè)你有兩個模型，AIC和BIC分數(shù)較低的模型更好。

變量選擇方法

有三種變量選擇方法 - 向前，向后，逐步。

1.以單個變量開始，然后基于AIC（“前進”）一次添加一個變量

2.從所有變量開始，基于AIC（’后退’）迭代地去除那些重要性低的變量

3.雙向運行（’逐步’）

  library(MASS)
  step <- stepAIC(fit, direction="both")
  summary(step)
  step <- stepAIC(fit, direction="forward")
  summary(step)
  n <- dim(dat3)[1]
  stepBIC <- stepAIC(fit,k=log(n))
  summary(stepBIC)

在基于BIC執(zhí)行逐步選擇之后，查看以上估計值。 變量已經(jīng)減少，但調(diào)整后的R-Squared保持不變（稍微改善）。 AIC和BIC分數(shù)也下降，這表明一個更好的模型。

  AIC(stepBIC)
  BIC(stepBIC)

計算標準化系數(shù)

標準化系數(shù)有助于根據(jù)標準化估計值的絕對值排列預(yù)測值。 值越高，變量越重要。

  #使用QuantPsyc包的lm.beta函數(shù)計算
  library(QuantPsyc)
  lm.beta(stepBIC)

  #自定義函數(shù)計算
  stdz.coff <- function (regmodel)
  { b <- summary(regmodel)$coef[-1,1]
  sx <- sapply(regmodel$model[-1], sd)
  sy <- sapply(regmodel$model[1], sd)
  beta <- b * sx / sy
  return(beta)
  }

  std.Coeff <- data.frame(Standardized.Coeff = stdz.coff(stepBIC))
  std.Coeff <- cbind(Variable = row.names(std.Coeff), std.Coeff)
  row.names(std.Coeff) <- NULL

計算方差膨脹因子（VIF）

與獨立變量高度相關(guān)的情況相比，差異膨脹因子衡量的是系數(shù)的變化幅度。 它應(yīng)該小于5。

  vif(stepBIC)

測試其它假設(shè)

 #Autocorrelation Test
  durbinWatsonTest(stepBIC)

  #Normality Of Residuals (Should be > 0.05)
  res=residuals(stepBIC,type="pearson")
  shapiro.test(res)

  #Testing for heteroscedasticity (Should be > 0.05)
  ncvTest(stepBIC)

  #Outliers – Bonferonni test
  outlierTest(stepBIC)
  #See Residuals
  resid = residuals(stepBIC)

  #Relative Importance
  install.packages("relaimpo")
  library(relaimpo)
  calc.relimp(stepBIC)

查看實際值和預(yù)測值

 #See Predicted Value
  pred <- predict(stepBIC,dat3)
  #See Actual vs. Predicted Value
  finaldata <- cbind(mtcars,pred)
  print(head(subset(finaldata, select = c(mpg,pred))))

其它有用的函數(shù)

  #Calculating RMSE
  rmse = sqrt(mean((dat3$mpg - pred)^2))
  print(rmse)

  #Calculating Rsquared manually
  y = dat3[,c("mpg")]
  R.squared = 1 - sum((y-pred)^2)/sum((y-mean(y))^2)
  print(R.squared)

  #Calculating Adj. Rsquared manually
  n = dim(dat3)[1]
  p = dim(summary(stepBIC)$coeff)[1] - 1
  adj.r.squared = 1 - (1 - R.squared) * ((n - 1)/(n-p-1))
  print(adj.r.squared)

  #Box Cox Transformation
  library(lmSupport)
  modelBoxCox(stepBIC)

K-fold交叉驗證

在下面的程序中，我們正在進行5倍交叉驗證。 在5倍交叉驗證中，數(shù)據(jù)被隨機分成5個相同大小的樣本。 在5個樣本中，隨機20％數(shù)據(jù)的單個樣本保留為驗證數(shù)據(jù)，其余80％用作訓練數(shù)據(jù)。 然后這個過程重復(fù)5次，5個樣本中的每一個都只用作驗證數(shù)據(jù)一次。 稍后我們將結(jié)果平均。

  library(DAAG)
  kfold = cv.lm(data=dat3, stepBIC, m=5)

完整代碼

library(ggplot2)
library(car)
library(caret)
library(corrplot)
data(mtcars)
str(mtcars)
head(mtcars)
summary(mtcars)
mtcars$am <- as.factor(mtcars$am)
mtcars$cyl <- as.factor(mtcars$cyl)
mtcars$vs <- as.factor(mtcars$vs)
mtcars$gear <- as.factor(mtcars$gear)
mtcars_a <- subset(mtcars, select = -c(mpg))
numericData <- mtcars_a[sapply(mtcars_a, is.numeric)]
descrCor <- cor(numericData)
print(descrCor)
corrplot(
descrCor,
order = "FPC",
method = "color",
type = "lower",
tl.cex = 0.7,
tl.col = rgb(0,0,0)
)
highlyCorrelated <- findCorrelation(descrCor, cutoff=0.7)
highlyCorCol <- colnames(numericData)[highlyCorrelated]
dat3 <- mtcars[, -which(colnames(mtcars) %in% highlyCorCol)]
dim(dat3)
fit <- lm(mpg ~ ., data=dat3)
summary(fit)
summary(fit)$coeff
anova(fit)
par(mfrow=c(2,2))
plot(fit)
summary(fit)$r.squared
summary(fit)$adj.r.squared
AIC(fit)
BIC(fit)
library(MASS)
step <- stepAIC(fit, direction="both")
summary(step)
step <- stepAIC(fit, direction="forward")
summary(step)
n <- dim(dat3)[1]
stepBIC <- stepAIC(fit,k=log(n))
summary(stepBIC)
AIC(stepBIC)
BIC(stepBIC)
library(QuantPsyc)
lm.beta(stepBIC)
stdz.coff <- function (regmodel)
{ b <- summary(regmodel)$coef[-1,1]
sx <- sapply(regmodel$model[-1], sd)
sy <- sapply(regmodel$model[1], sd)
beta <- b * sx / sy
return(beta)
}
std.Coeff <- data.frame(Standardized.Coeff = stdz.coff(stepBIC))
std.Coeff <- cbind(Variable = row.names(std.Coeff), std.Coeff)
row.names(std.Coeff) <- NULL
vif(stepBIC)
#Autocorrelation Test
durbinWatsonTest(stepBIC)
#Normality Of Residuals (Should be > 0.05)
res=residuals(stepBIC,type="pearson")
shapiro.test(res)
#Testing for heteroscedasticity (Should be > 0.05)
ncvTest(stepBIC)
#Outliers – Bonferonni test
outlierTest(stepBIC)
#See Residuals
resid = residuals(stepBIC)
#Relative Importance
install.packages("relaimpo")
library(relaimpo)
calc.relimp(stepBIC)
#See Predicted Value
pred <- predict(stepBIC,dat3)
#See Actual vs. Predicted Value
finaldata <- cbind(mtcars,pred)
print(head(subset(finaldata, select = c(mpg,pred))))

來自搜狐網(wǎng)：https://www.sohu.com/a/230584172_274950

生物信息學習.jpg