文章來源于： sci

這一系列里前面的三個部分都是用于比較組間差異的各種方法,，詳情點擊：R語言系列第四期：①R語言單樣本雙樣本差異性檢驗,、R語言系列第四期：②R語言多組樣本方差分析與KW檢驗、R語言系列第四期：③R語言表格數(shù)據(jù)率的比較

在這個部分里,，我們會為大家介紹如何使用R進行基礎(chǔ)回歸和相關(guān)分析,，以及模型作圖,、置信區(qū)間的預估和展示,。

我們使用數(shù)據(jù)集thuesen作為這一部分的例子，如下導入：

> library(ISwR)

> attach(thuesen)

我們使用函數(shù)lm( linear model,線性模型 )進行線性分析：

> lm(short.velocity~blood.glucose)

Call:

lm(formula = short.velocity ~ blood.glucose)

Coefficients:

(Intercept)  blood.glucose 

1.09781        0.02196 

#Tips：函數(shù)lm()里面的參數(shù)是模型方程,，波浪號(~)可以認為是“通過…來描述”,。它在之前出現(xiàn)過幾次，比如圖形展示部分箱式圖boxplot(),t檢驗,，anova檢驗里等等,。

#Tips：lm()函數(shù)的原始輸出格式非常簡單。你能看見的只有估計出來的截距α與斜率β,?？梢钥闯鲎顑?yōu)擬合直線為：short.velocity=1.098+0.0220*blood.glucose，但是沒有給出其他任何像顯著性檢驗之類的信息,。

#Tips：其實,，函數(shù)lm()可以處理比簡單線性回歸復雜很多的模型,。除了一個解釋變量與一個因變量之外，模型方程還能描述很多其他的情況,。比如，要在y上通過x1,，x2,，x3進行多元線性回歸分析（后文會介紹），可以通過y~x1+x2+x3來完成,。

如果要獲悉模型的其他信息,，可以使用summary()：

> summary(lm(short.velocity~blood.glucose))

Call:

lm(formula = short.velocity ~ blood.glucose)

Residuals:

Min       1Q   Median       3Q      Max

-0.40141 -0.14760 -0.02202  0.03001  0.43490

Coefficients:

Estimate Std. Error t value Pr(>|t|)   

(Intercept)    1.09781    0.11748   9.345 6.26e-09 ***

blood.glucose  0.02196    0.01045   2.101   0.0479 * 

—

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2167 on 21 degrees of freedom

(1 observation deleted due to missingness)

Multiple R-squared:  0.1737,    Adjusted R-squared:  0.1343 

F-statistic: 4.414 on 1 and 21 DF,  p-value: 0.0479

##Tips：首先，是

Call:

lm(formula = short.velocity ~ blood.glucose)

輸出的開頭本質(zhì)上在重復一個函數(shù)的調(diào)用,。當然如果你一開始就把模型賦值給了一個變量,，那么調(diào)用summary()之后這個部分就是那個變量了。

Residuals:

Min  1Q   Median   3Q    Max

-0.40141 -0.14760 -0.02202  0.03001  0.43490

這部分簡單地描述了殘差的分布,，可以幫助用戶對分布的假設(shè)做快速檢查,。

Coefficients:

Estimate  Std. Error  t value  Pr(>|t|)   

(Intercept)    1.09781    0.11748   9.345 6.26e-09 ***

blood.glucose  0.02196    0.01045   2.101   0.0479 * 

—

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

這里我們再次見到了回歸系數(shù)和截距，不過這次還伴有標準誤,，t檢驗和p值,，以及最右側(cè)的統(tǒng)計檢驗表示顯著性水平的圖像化標志。當然,，如果你不喜歡這些星號,，可以使用options(show.signif.stars=F)關(guān)閉它們。

Residual standard error: 0.2167 on 21 degrees of freedom

(1 observation deleted due to missingness)

殘差波動情況,。

Multiple R-squared:  0.1737,    Adjusted R-squared:  0.1343 

F-statistic: 4.414 on 1 and 21 DF,  p-value: 0.0479

上式第一項是R²,，在簡單線性回歸里可以被理解為Pearson相關(guān)系數(shù)的平方，另一個是修正后的R²,；第二行是對假設(shè)回歸系數(shù)是0進行的F檢驗,，對整體模型的檢驗。實際上,，F檢驗是t檢驗的平方：4.414=（2.101）²,。

我們先把數(shù)據(jù)點和回歸線畫出來：

> plot(blood.glucose,short.velocity)

> abline(lm(short.velocity~blood.glucose))

#Tips：abline()函數(shù)根據(jù)截距和斜率畫一條直線。它能夠接受數(shù)值參數(shù),，比如abline(1.1,0.022),；不過更方便的是，它也能夠從一個用lm擬合的線性回歸中直接提取相關(guān)信息,。

無論是否計算了置信帶和預測帶,，我們都能夠用函數(shù)predict析取出預測值，不加其他參數(shù),，它就只會輸出回歸值,。我們?yōu)榱朔奖阌米兞?/span>lm.velo代替模型：

> lm.velo<-lm(short.velocity~blood.glucose)

> predict(lm.velo)

1    2    3    4    5    6    7

1.433841 1.335010 1.275711 1.526084 1.255945 1.214216 1.302066

8    9    10   11  12  13  14

1.341599 1.262534 1.365758 1.244964 1.212020 1.515103 1.429449

15      17      18      19      20      21      22

1.244964 1.190057 1.324029 1.372346 1.451411 1.389916 1.205431

23      24

1.291085 1.306459

如果你加上了interval=”confidence”或者interval=”prediction”,，那么你就能在預測值向量的基礎(chǔ)上得到各類估計邊界的值。同時也可以縮寫這兩個單詞：

> predict(lm.velo,int=”c”)

fit      lwr      upr

1  1.433841 1.291371 1.576312

2  1.335010 1.240589 1.429431

3  1.275711 1.169536 1.381887

…..

22 1.205431 1.053805 1.357057

23 1.291085 1.191084 1.391086

24 1.306459 1.210592 1.402326

> predict(lm.velo,int=”p”)

fit       lwr      upr

1  1.433841 0.9612137 1.906469

2  1.335010 0.8745815 1.795439

3  1.275711 0.8127292 1.738693

…

22 1.205431 0.7299634 1.680899

23 1.291085 0.8294798 1.752690

24 1.306459 0.8457315 1.767186

Warning message:

In predict.lm(lm.velo, int = “p”) :

用當前數(shù)據(jù)得到的預測結(jié)果對_未來_響應有用

#Tips：前一個是置信帶,，后一個是預測帶,。lwr和upr分別是下界和上界。Warning信息里提醒我們：這個預測邊界不能用來考察我們做回歸線所使用的已觀測數(shù)據(jù),。

我們可以利用外部數(shù)據(jù)根據(jù)已有的模型做出它的預測情況圖,，圖形的展示可以通過下面的組合函數(shù)來實現(xiàn)：

> pred.frame<-data.frame(blood.glucose=4:20)

> pp<-predict(lm.velo,int=”p”,newdata=pred.frame)

> pc<-predict(lm.velo,int=”c”,newdata=pred.frame)

> plot(blood.glucose,short.velocity,ylim=range(short.velocity,pp,na.rm=T))

> pred.gluc<-pred.frame$blood.glucose

> matlines(pred.gluc,pc,lty=c(1,2,2),col=”black”)

> matlines(pred.gluc,pp,lty=c(1,3,3),col=”black”)

#Tips：我們采用由4到20這幾個數(shù)據(jù)作為新的x來做出預測情況的圖。Predict()函數(shù)里newdata=的參數(shù)就是調(diào)用新數(shù)據(jù)的參數(shù),；plot()函數(shù)里的ylim參數(shù)使用range()函數(shù)來保證圖形全部在范圍內(nèi),；matlines()函數(shù)里的lty是設(shè)置線型。

相關(guān)系數(shù)的計算可以使用cor()函數(shù),，但是如果對thuesen中的兩個向量也進行這樣簡單的操作,，就會發(fā)生下面狀況：

> cor(blood.glucose,short.velocity)

[1] NA

R中所有的基本統(tǒng)計函數(shù)都要求輸入的參數(shù)沒有缺失值，或者你明確指定如何處理缺失值,。對于函數(shù)mean（）,，var（），sd（）以及類似的單向量函數(shù),，你可以傳遞na.rm=T這個參數(shù)告訴它們在計算之前應該移除缺失值,。對于cor()函數(shù)，use=“complete.obs”代表提取所有變量全部非空的觀測,，你也可以這樣寫：

> cor(blood.glucose,short.velocity,use=”complete.obs”)

[1] 0.4167546

我們還可以通過如下的代碼得到一個數(shù)據(jù)框中多種變量的相關(guān)系數(shù)矩陣：

> cor(thuesen,use=”complete.obs”)

blood.glucose short.velocity

blood.glucose   1.0000000   0.4167546

short.velocity    0.4167546    1.0000000

#Tips：當然,，數(shù)據(jù)框中變量超過兩個結(jié)果會更有意思。

但是,，上面的結(jié)果只是跟我們展示了關(guān)聯(lián)系數(shù),，但是沒有做出檢驗，因此我們需要cor.test()函數(shù),，可以通過指定兩個變量來使用：

> cor.test(blood.glucose,short.velocity)

Pearson’s product-moment correlation

data:  blood.glucose and short.velocity

t = 2.101, df = 21, p-value = 0.0479

alternative hypothesis: true correlation is not equal to 0

95 percent confidence interval:

 0.005496682 0.707429479

sample estimates:

cor

0.4167546

#Tips：我們也能得到一個真實相關(guān)系數(shù)的置信區(qū)間,。注意，這里的p值和之前回歸分析的p值是一樣的,。同樣與之前回歸模型的anova表里的p值是一樣的,。

與前面的部分所講的單樣本和雙樣本問題一樣，相關(guān)問題也有非參數(shù)的方法,，這些方法的優(yōu)點在于不需要假設(shè)數(shù)據(jù)的正態(tài)分布性,，而且結(jié)果也不會受到單調(diào)變換的影響。

相關(guān)性檢驗的幾個方法都打包進了cor.test中,，沒有額外提供專門的spearman.test()函數(shù),。所以可以在cor.test()中指明：

> cor.test(blood.glucose,short.velocity,method=”spearman”)

Spearman’s rank correlation rho

data:  blood.glucose and short.velocity

S = 1380.4, p-value = 0.1392

alternative hypothesis: true rho is not equal to 0

sample estimates:

rho

0.318002

Warning message:

In cor.test.default(blood.glucose, short.velocity, method = “spearman”) :

無法給連結(jié)計算精確p值

而等級相關(guān)同理，只需要在method參數(shù)中做出修改就可以實現(xiàn)：

> cor.test(blood.glucose,short.velocity,method=”kendall”)

Kendall’s rank correlation tau

data:  blood.glucose and short.velocity

z = 1.5604, p-value = 0.1187

alternative hypothesis: true tau is not equal to 0

sample estimates:

tau

0.2350616

Warning message:

In cor.test.default(blood.glucose, short.velocity, method = “kendall”) :

無法給連結(jié)計算精確p值

#Tips：注意,，這兩個非參數(shù)方法的相關(guān)系數(shù)在5%置信水平下不顯著,，而pearson相關(guān)系數(shù)也僅僅是剛過線的顯著,。

到目前為止，簡單的統(tǒng)計分析方法就告一段落了,，而進一步的統(tǒng)計分析探討我們會在下面一個系列為大家介紹,，敬請期待。