多元線性回歸分析 (Multiple Linear Regression)

套路43: 多元線性回歸分析 (Multiple Linear Regression)

1. 使用時機: 以多個獨立自變項預測一個應變項。
2. 分析類型: 母數分析(parametric analysis)。
3. 範例資料: 咪路測量菌菌的生長條件與菌數，資料如下表。咪路希望能得到一回歸方程式可用來預測菌菌生長。
item
Temp
pH
hour
ml
CFU(106)
1
6
5.7
1.6
2.12
9.9
2
1
6.4
3
3.39
9.3
3
-2
5.7
3.4
3.61
9.4
4
11
6.1
3.4
1.72
9.1
5
-1
6
3
1.8
6.9
6
2
5.7
4.4
3.21
9.3
7
5
5.9
2.2
2.59
7.9
8
1
6.2
2.2
3.25
7.4
9
1
5.5
1.9
2.86
7.3
10
3
5.2
0.2
2.32
8.8
11
11
5.7
4.2
1.57
9.8
12
9
6.1
2.4
1.5
10.5
13
5
6.4
3.4
2.69
9.1
14
-3
5.5
3
4.06
10.1
15
1
5.5
0.2
1.98
7.2
16
8
6
3.9
2.29
11.7
17
-2
5.5
2.2
3.55
8.7
18
3
6.2
4.4
3.31
7.6
19
6
5.9
0.2
1.83
8.6
20
10
5.6
2.4
1.69
10.9
21
4
5.8
2.4
2.42
7.6
22
5
5.8
4.4
2.98
7.3
23
5
5.2
1.6
1.84
9.2
24
3
6
1.9
2.48
7
25
8
5.5
1.6
2.83
7.2
26
8
6.4
4.1
2.41
7
27
6
6.2
1.9
1.78
8.8
28
6
5.4
2.2
2.22
10.1
29
3
5.4
4.1
2.72
12.1
30
5
6.2
1.6
2.36
7.7
31
1
6.8
2.4
2.81
7.8
32
8
6.2
1.9
1.64
11.5
33
10
6.4
2.2
1.82
10.4
4. 執行回歸。
第一步: 輸入建立資料。
  m <- read.table(header = T, text = "
Temp pH hour ml CFU
6 5.7 1.6 2.12 9.9
1 6.4 3 3.39 9.3
-2 5.7 3.4 3.61 9.4
11 6.1 3.4 1.72 9.1
-1 6 3 1.8 6.9
2 5.7 4.4 3.21 9.3
5 5.9 2.2 2.59 7.9
1 6.2 2.2 3.25 7.4
1 5.5 1.9 2.86 7.3
3 5.2 0.2 2.32 8.8
11 5.7 4.2 1.57 9.8
9 6.1 2.4 1.5 10.5
5 6.4 3.4 2.69 9.1
-3 5.5 3 4.06 10.1
1 5.5 0.2 1.98 7.2
8 6 3.9 2.29 11.7
-2 5.5 2.2 3.55 8.7
3 6.2 4.4 3.31 7.6
6 5.9 0.2 1.83 8.6
10 5.6 2.4 1.69 10.9
4 5.8 2.4 2.42 7.6
5 5.8 4.4 2.98 7.3
5 5.2 1.6 1.84 9.2
3 6 1.9 2.48 7
8 5.5 1.6 2.83 7.2
8 6.4 4.1 2.41 7
6 6.2 1.9 1.78 8.8
6 5.4 2.2 2.22 10.1
3 5.4 4.1 2.72 12.1
5 6.2 1.6 2.36 7.7
1 6.8 2.4 2.81 7.8
8 6.2 1.9 1.64 11.5
10 6.4 2.2 1.82 10.4")  # 資料間以空白分隔
第二步: 使用基本模組(base)的attach及names函數賦予m中的資料名稱(標題)。
  attach(m)  # 讓R能透過變數名稱搜尋資料。
  names(m)  # 賦予資料名稱(標題)。
第三步: 使用基本模組(base)的lm函數代入m資料求回歸方程式，結果儲存至變數x。
  x <- lm(CFU ~ Temp + pH + hour + ml, data = m)
  # CFU因變數，Temp、pH、hour及ml為自變數。
  # CFU ~ Temp + pH + hour + ml即回歸方程式的組成:
     CFU = 係數*Temp + 係數*pH + 係數*hour + 係數*ml
第四步: 顯示判讀結果。
  coefficients(x)
 # 以coefficients 函數顯示求得之回歸方程式係數(model coefficients)。
(Intercept)        Temp        pH        hour         ml
 13.9229254   0.1116533  -0.9924694   0.3241469  -0.2109881

  summary(x)
  # 以summary 函數顯示求得之回歸方程式係數及相關資料p值。
Call:
lm(formula = CFU ~ Temp + pH + hour + ml, data = m)
Residuals:
   Min    1Q  Median    3Q    Max
-2.2849 -1.0231  0.1202  0.9970  2.5673
Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)  13.9229     4.1799   3.331  0.00244 **
Temp          0.1117     0.1065   1.048  0.30364  
pH           -0.9925     0.6695  -1.482  0.14939  
hour          0.3241     0.2534   1.279  0.21129  
ml           -0.2110     0.6321  -0.334  0.74102  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.42 on 28 degrees of freedom
Multiple R-squared:  0.2001,    Adjusted R-squared:  0.08578
F-statistic: 1.751 on 4 and 28 DF,  p-value: 0.167

  confint(x, level = 0.95)
  # 以confint函數求信賴區間(CIs for model parameters)。
                   2.5 %     97.5 %
(Intercept)    5.3606934 22.4851573
Temp        -0.1066003  0.3299069
pH          -2.3638398  0.3789010
hour         -0.1948710  0.8431647
ml          -1.5057751  1.0837989
   
  anova(x)
  # 顯示anova table
Analysis of Variance Table
Response: CFU
              Df Sum Sq Mean Sq F value  Pr(>F) 
Temp       1  7.631  7.6311  3.7822 0.06191 .
pH         1  2.924  2.9245  1.4495 0.23870 
hour       1  3.348  3.3480  1.6594 0.20823 
ml         1  0.225  0.2248  0.1114 0.74102 
Residuals 28 56.494  2.0176                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

  vcov(x)
  # covariance matrix for model parameters
            (Intercept)         Temp           pH        hour          ml
(Intercept)  17.4719516 -0.120085075 -2.544985069  0.30277089 -1.08175295
Temp         -0.1200851  0.011352466 -0.005645347 -0.01089930  0.05271535
pH           -2.5449851 -0.005645347  0.448205285 -0.03677360  0.01105980
hour          0.3027709 -0.010899303 -0.036773602  0.06419957 -0.08129613
ml           -1.0817529  0.052715348  0.011059800 -0.08129613  0.39954353

  # 其他相關資料可以下列函數取得。
  fitted(x) # predicted values
  residuals(x) # residuals
  influence(x) # regression diagnostics

5. 挑選變數，重新求回歸方程式。
第一步: 建立初始模型CFU~1及全模型CFU~. (全模型就是用上所有的變數)。
  null <- lm(CFU~1, data = m)
  null
Call:
lm(formula = CFU ~ 1, data = m)
Coefficients:
(Intercept) 
      8.885
  full <- lm(CFU~., data = m)
  full
Call:
lm(formula = CFU ~ ., data = m)
Coefficients:
(Intercept)         Temp          pH         hour          ml 
    13.9229       0.1117      -0.9925       0.3241      -0.2110 
第二步: 使用基本模組(base)的step函數代入初始模型null及全模型full以forward addition方法挑選變數。
  step(null, scope = list(lower = null, upper = full), direction = "forward")
Start:  AIC=27.11
CFU ~ 1
       Df Sum of Sq    RSS    AIC
+ Temp  1    7.6311 62.991 25.334
<none>              70.622 27.108
+ ml    1    3.1752 67.447 27.590
+ hour  1    2.2965 68.326 28.017
+ pH    1    1.4950 69.127 28.402
Step:  AIC=25.33
CFU ~ Temp
       Df Sum of Sq    RSS    AIC
<none>              62.991 25.334
+ pH    1   2.92447 60.067 25.765
+ hour  1   1.88749 61.104 26.330
+ ml    1   0.12047 62.871 27.271
Call:
lm(formula = CFU ~ Temp, data = m)
Coefficients:
(Intercept)         Temp 
     8.3186       0.1271
第三步: 使用基本模組的step函數以backward elimination方法挑選變數。
  step(full, data = m, direction = "backward")
Start:  AIC=27.74
CFU ~ Temp + pH + hour + ml
       Df Sum of Sq    RSS    AIC
- ml    1    0.2248 56.719 25.873
- Temp  1    2.2156 58.710 27.011
- hour  1    3.3021 59.796 27.616
<none>              56.494 27.742
- pH    1    4.4341 60.928 28.235
Step:  AIC=25.87
CFU ~ Temp + pH + hour
       Df Sum of Sq    RSS    AIC
- hour  1     3.348 60.067 25.765
<none>              56.719 25.873
- pH    1     4.385 61.104 26.330
- Temp  1     8.928 65.647 28.697
Step:  AIC=25.77
CFU ~ Temp + pH
       Df Sum of Sq    RSS    AIC
- pH    1    2.9245 62.991 25.334
<none>              60.067 25.765
- Temp  1    9.0606 69.127 28.402
Step:  AIC=25.33
CFU ~ Temp
       Df Sum of Sq    RSS    AIC
<none>              62.991 25.334
- Temp  1    7.6311 70.622 27.108
Call:
lm(formula = CFU ~ Temp, data = m)
Coefficients:
(Intercept)         Temp 
     8.3186       0.1271 
第四步: 使用基本模組的step函數以stepwise regression方法挑選變數。
  step(null, scope = list(upper = full), data = m, direction = "both")
Start:  AIC=27.11
CFU ~ 1
       Df Sum of Sq    RSS    AIC
+ Temp  1    7.6311 62.991 25.334
<none>              70.622 27.108
+ ml    1    3.1752 67.447 27.590
+ hour  1    2.2965 68.326 28.017
+ pH    1    1.4950 69.127 28.402
Step:  AIC=25.33
CFU ~ Temp
       Df Sum of Sq    RSS    AIC
<none>              62.991 25.334
+ pH    1    2.9245 60.067 25.765
+ hour  1    1.8875 61.104 26.330
- Temp  1    7.6311 70.622 27.108
+ ml    1    0.1205 62.871 27.271
Call:
lm(formula = CFU ~ Temp, data = m)
Coefficients:
(Intercept)         Temp 
     8.3186       0.1271 
# 步驟二到四三種挑選模型的方法都挑出CFU ~ Temp為最佳模型。
第五步: 比較兩組回歸方程式統計上有無差別(就是省略某些變數不影響)。
  y1 <- lm(CFU ~ Temp + pH + hour + ml, data = m)
  y2 <- lm(CFU ~ Temp, data = m)
  anova(y1, y2)
第六步: 判讀結果。
Analysis of Variance Table
Model 1: CFU ~ Temp + pH + hour + ml
Model 2: CFU ~ Temp
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     28 56.494                          
2     31 62.991 -3   -6.4973 1.0734  0.3763
  # Pr < 0.05，H0: 兩個方程式效果相當，不成立。
  # Pr > 0.05，H0: 兩個方程式效果相當，成立。

來勁了嗎? 想知道更多?? 補充資料(連結):
1. 關於Regression analysis (https://en.wikipedia.org/wiki/Regression_analysis)
2. 關於Linear regression (https://en.wikipedia.org/wiki/Linear_regression)
3. STHDA Linear Regression Essentials in R (http://www.sthda.com/english/articles/40-regression-analysis/165-linear-regression-essentials-in-r/)
4. 關於R基礎，R繪圖及統計快速入門:
   a. R Tutorial: https://www.tutorialspoint.com/r/index.htm
   b. Cookbook for R: http://www.cookbook-r.com/
   c. Quick-R: https://www.statmethods.net/
   d. Statistical tools for high-throughput data analysis (STHDA): http://www.sthda.com/english/
e. The Handbook of Biological Statistics: http://www.biostathandbook.com/
f. An R Companion for the Handbook of Biological Statistics: http://rcompanion.org/rcompanion/index.html
5. Zar, JH. 2010. Biostatistical Analysis, Fifth Edition, Pearson.

留言

張貼留言

這個網誌中的熱門文章

統計不球人 目錄 (Table of Contents)

用 R 做生醫資料科學   (Biomedical Data Science using R) 如何選擇統計方法 1 如何選擇統計方法  2 ***************************   第一段 : 起手式   ****************************** 套路 0: 安裝 R 不球人 套路 1: R 算統計不球人 ， 使用 R 函數 套路 2: 使用 ggplot2 繪圖 (Graphics using ggplot2) 套路 3: 敘述性統計 (Descriptive Statistics) 套路 4: 如何檢查資料是否接近常態分布 (Check Normality using R) 套路 5: 機率分布 (Probability Distributions) ***************************   第二段 : 比較平均值、中數   *************************** 套路 6: 單一樣本 Z 檢定 (One-Sample Z test) 套路 7: 單一樣本 t 檢定 (One-Sample t test) 套路 8: 單一樣本 Wilcoxon signed-rank 檢定 ( One-Sample Wilcoxon signed-rank test ， non-parametric ) 套路 9: 單一樣本符號檢定 (One Sample Sign test) 套路 10: 兩組獨立樣本 Z 檢定 (Two-Sample Z Test) 套路 11: 兩組變異數相同獨立樣本 t 檢定 ( t Test for Two-Samples with Equal Variances) 套路 12: 兩組變異數不同樣本 t 檢定 ( t Test for Two-Samples with Unequal Variances ， Welch Two Sample t test) 套路 13: 成對樣本 t 檢定 (Paired Sample t test) 套路 14: 兩組獨立樣本曼恩 - 惠尼 U 檢定 (...
閱讀完整內容

比較二或多組變異數Levene’s 檢定 (Levene’s Test for Comparing Two or More Variances)

圖片
套路 27: 比較二或多組變異數 Levene’s 檢定   (Levene’s Test for Comparing Two or More Variances) 什麼是比較多組獨立樣本變異數檢定 ? 說白了就是多組分別獨立取樣的資料做比較的假設檢定。 ** 注意 ** ” 比較變異數 ” (comparing variances) 和 ” 變異數分析 ” (analysis of variance) 不同。變異數分析是多組資料比較平均值。統計假設檢定檢定什麼 ? 看 H 0 。例如多組獨立樣本假設檢定 H 0 : s 1 2 = s 2 2 = … = s k 2 ， H A : 至少有一組變異數不同 ， 是檢定多組資料的變異數是否相同。假設相等時為雙尾 (two-tailed test) 檢定。 1. 使用時機 : 用於比較觀測到的兩組或多組變異數 (variances) 。 2. 分析類型 : 母數分析 (parametric analysis) 。 直接使用資料數值算統計叫 parametric 方法，把資料排序之後用排序的名次算統計叫 non-parametric 方法。 3. 前提假設 : 資料偏離常態分布 (normal distribution) 仍可用。 4. 資料範例 1 ， 兩組變異數 : 咪路調查高一和大一學生體重 (kg) ， 資料如下 : 高一 41 35 33 36 40 46 31 37 34 30 38 大一 52 57 62 55 64 57 56 55 60 59 請問高一和大一學生體重變異數是否相同 ? H 0 : σ 2 1 = σ 2 2 ， H A : σ 2 1 ≠ σ 2 2 。 第一步 : 輸入建立資料 ， 用小 c 將資料放入名稱為 h1 及...
閱讀完整內容

三因子變異數分析 (Three Way ANOVA)

圖片
套路 24: 三因子變異數分析 (Three Way ANOVA) 應用三因子變異數分析的資料有三影響因子 ，就是有三個自變項 。此檢定有七組 H 0 和 H A 。 H 0 : m 因子 1-1 = m 因子 1-2 = … = m 因子 1-m ， H A : 至少有一組平均值不同。 H 0 : m 因子 2-1 = m 因子 2-2 = … = m 因子 2-n ， H A : 至少有一組平均值不同。 H 0 : m 因子 3-1 = m 因子 3-2 = … = m 因子 3-o ， H A : 至少有一組平均值不同。 H 0 : 因子 1 與因子 2 互不影響 ， H A : 因子 1 與因子 2 互相影響。 H 0 : 因子 1 與因子 3 互不影響 ， H A : 因子 1 與因子 3 互相影響。 H 0 : 因子 2 與因子 3 互不影響 ， H A : 因子 2 與因子 3 互相影響。 H 0 : 因子 1 、 因子 2 與因子 3 互不影響 ， H A : 因子 1 、 因子 2 與因子 3 互相影響。 1. 使用時機 : 變異數分析適用連續資料之差異分析，若自變項有 3 個，就是三因子變異數分析。 2. 分析類型 : 母數分析 (parametric analysis) 。 直接使用資料數值算統計叫 parametric 方法，把資料排序之後用排序的名次算統計叫 non-parametric 方法。 3. 假設前提 : 每組資料 (observations within each cell) 都是常態分布 (normally distributed) 且具相同變異數 (equal variances) 。 4. 範例資料 : 咪路研究溫度對孔雀魚 ( Poecilia reticulata ) 、寶蓮燈魚 ( Paracheirodon axelrodi ) 及斑馬魚 ( Danio rerio ) 生長的影響。在 20 ° C 、 24 ° C 及 28 ° C 培養下，三種魚身長變化 (cm) 資料如下： 孔雀魚 ( Poecilia reticulata ) 20 ° C 24 ° C ...
閱讀完整內容