原文鏈接：http://tecdat.cn/?p=18266

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

?

本文通過(guò)R語(yǔ)言建立廣義線性模型(GLM)、多項(xiàng)式回歸和廣義可加模型（GAM）來(lái)預(yù)測(cè)誰(shuí)在1912年的泰坦尼克號(hào)沉沒(méi)中幸存下來(lái)。

1. 


2. str(titanic)

數(shù)據(jù)變量為：

• Survived：乘客存活指標(biāo)（如果存活則為1）

• Pclass：旅客艙位等級(jí)

• Sex：乘客性別

• Age：乘客年齡

• SibSp：兄弟姐妹/配偶人數(shù)

• Parch：父母/子女人數(shù)

• Embarked： 登船港口

• Name：旅客姓名

最后一個(gè)變量使用不多，因此我們將其刪除，

titanic = titanic[,1:7]

現(xiàn)在，我們回答問(wèn)題：

幸存的旅客比例是多少？

簡(jiǎn)單的答案是

?

1. mean(titanic$Survived)

2. [1] 0.3838384

可以在下面的列聯(lián)表中找到

1. table(titanic$Survived)/nrow(titanic)

2. 0 ? ? ? ? 1

3. 0.6161616 0.3838384

或此處幸存者的38.38 ％。也就是說(shuō)，也可以通過(guò)不對(duì)任何解釋變量進(jìn)行邏輯回歸來(lái)獲得（換句話說(shuō)，僅對(duì)常數(shù)進(jìn)行回歸）?；貧w給出了：

1. 


2. Coefficients:

3. (Intercept)

4. -0.4733

5. 


6. Degrees of Freedom: 890 Total (i.e. Null); ?890 Residual

7. Null Deviance: ? ?1187

8. Residual Deviance: 1187 AIC: 1189

給出β0的值，并且由于生存概率為

我們通過(guò)考慮

?

1. exp(-0.4733)/(1+exp(-0.4733))

2. [1] 0.3838355

我們也可以使用predict函數(shù)?

1. predict(glm(Survived~1, family=binomial,type="response")[1]

2. 1

3. 0.3838384

此外，在概率回歸中也適用，

1. reg=glm(Survived~1, family=binomial(link="probit"),data=titanic)

2. predict(reg,type="response")[1]

3. 1

4. 0.3838384

幸存的頭等艙乘客的比例是多少？

我們只看頭等艙的人，

?

[1] 0.6296296

約63％存活。我們可以進(jìn)行邏輯回歸

1. 


2. Coefficients:

3. (Intercept) ? ? ?Pclass2 ? ? ?Pclass3

4. 0.5306 ? ? ?-0.6394 ? ? ?-1.6704

5. 


6. Degrees of Freedom: 890 Total (i.e. Null); ?888 Residual

7. Null Deviance: ? ?1187

8. Residual Deviance: 1083 AIC: 1089

由于第1類(lèi)是參考類(lèi)，因此我們照舊考慮

1. exp(0.5306)/(1+exp(0.5306))

2. [1] 0.629623

predict預(yù)測(cè)：

1. 


2. predict(reg,newdata=data.frame(Pclass="1"),type="response")

3. 1

4. 0.6296296

我們可以嘗試概率回歸，我們得到的結(jié)果是一樣的，

1. 


2. predict(reg,newdata=data.frame(Pclass="1"),type="response")

3. 1

4. 0.6296296

卡方獨(dú)立性測(cè)試 ：在生存與否之間的檢驗(yàn)統(tǒng)計(jì)量是多少？

卡方檢驗(yàn)的命令如下

1. chisq.test(table( Survived, Pclass))

2. 


3. Pearson's Chi-squared test

4. 


5. data: ?table( Survived, ?Pclass)

6. X-squared = 102.89, df = 2, p-value < 2.2e-16

我們有一個(gè)列聯(lián)表，如果變量是獨(dú)立的，我們有

，然后是統(tǒng)計(jì)量

，我們可以看到，對(duì)測(cè)試的貢獻(xiàn)

1. 


2. 1 ? ? ? ? 2 ? ? ? ? 3

3. 0 -4.601993 -1.537771 ?3.993703

4. 1 ?5.830678 ?1.948340 -5.059981

這給了我們很多信息：我們觀察到兩個(gè)正值，分別對(duì)應(yīng)于“幸存”和“頭等艙”與“無(wú)法幸存”和“三等艙”之間的強(qiáng)（正）關(guān)聯(lián)，以及兩個(gè)很強(qiáng)的負(fù)值，對(duì)應(yīng)于“生存”和“第三等”之間的強(qiáng)烈負(fù)相關(guān)，以及“無(wú)法幸存”和“頭等艙”。我們可以在下圖上可視化這些值

?

?

1. 


2. ass(table( Survived, Pclass), shade = TRUE, las=3)

?

然后我們必須進(jìn)行邏輯回歸，并預(yù)測(cè)兩名模擬乘客的生存概率

假設(shè)我們有兩名乘客

1. newbase = data.frame(

2. Pclass = as.factor(c(1,3)),

3. Sex = as.factor(c("female","male")),

4. Age = c(17,20),

5. SibSp = c(1,0),

6. Parch = c(2,0),

讓我們對(duì)所有變量進(jìn)行簡(jiǎn)單回歸，

1. 


2. 


3. Coefficients:

4. Estimate Std. Error z value Pr(>|z|)

5. (Intercept) ?16.830381 607.655774 ? 0.028 ?0.97790

6. Pclass2 ? ? ?-1.268362 ? 0.298428 ?-4.250 2.14e-05 ***

7. Pclass3 ? ? ?-2.493756 ? 0.296219 ?-8.419 ?< 2e-16 ***

8. Sexmale ? ? ?-2.641145 ? 0.222801 -11.854 ?< 2e-16 ***

9. Age ? ? ? ? ?-0.043725 ? 0.008294 ?-5.272 1.35e-07 ***

10. SibSp ? ? ? ?-0.355755 ? 0.128529 ?-2.768 ?0.00564 **

11. Parch ? ? ? ?-0.044628 ? 0.120705 ?-0.370 ?0.71159

12. EmbarkedC ? -12.260112 607.655693 ?-0.020 ?0.98390

13. EmbarkedQ ? -13.104581 607.655894 ?-0.022 ?0.98279

14. EmbarkedS ? -12.687791 607.655674 ?-0.021 ?0.98334

15. ---

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

17. 


18. (Dispersion parameter for binomial family taken to be 1)

19. 


20. Null deviance: 964.52 ?on 713 ?degrees of freedom

21. Residual deviance: 632.67 ?on 704 ?degrees of freedom

22. (177 observations deleted due to missingness)

23. AIC: 652.67

24. 


25. Number of Fisher Scoring iterations: 13

兩個(gè)變量并不重要。我們刪除它們

1. 


2. 


3. Coefficients:

4. Estimate Std. Error z value Pr(>|z|)

5. (Intercept) ?4.334201 ? 0.450700 ? 9.617 ?< 2e-16 ***

6. Pclass2 ? ? -1.414360 ? 0.284727 ?-4.967 6.78e-07 ***

7. Pclass3 ? ? -2.652618 ? 0.285832 ?-9.280 ?< 2e-16 ***

8. Sexmale ? ? -2.627679 ? 0.214771 -12.235 ?< 2e-16 ***

9. Age ? ? ? ? -0.044760 ? 0.008225 ?-5.442 5.27e-08 ***

10. SibSp ? ? ? -0.380190 ? 0.121516 ?-3.129 ?0.00176 **

11. ---

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

13. 


14. (Dispersion parameter for binomial family taken to be 1)

15. 


16. Null deviance: 964.52 ?on 713 ?degrees of freedom

17. Residual deviance: 636.56 ?on 708 ?degrees of freedom

18. (177 observations deleted due to missingness)

19. AIC: 648.56

20. 


21. Number of Fisher Scoring iterations: 5

我們有年齡這樣的連續(xù)變量時(shí)，我們可以進(jìn)行多項(xiàng)式回歸

1. 


2. 


3. Coefficients:

4. Estimate Std. Error z value Pr(>|z|)

5. (Intercept) ? ? 3.0213 ? ? 0.2903 ?10.408 ?< 2e-16 ***

6. Pclass2 ? ? ? ?-1.3603 ? ? 0.2842 ?-4.786 1.70e-06 ***

7. Pclass3 ? ? ? ?-2.5569 ? ? 0.2853 ?-8.962 ?< 2e-16 ***

8. Sexmale ? ? ? ?-2.6582 ? ? 0.2176 -12.216 ?< 2e-16 ***

9. poly(Age, 3)1 -17.7668 ? ? 3.2583 ?-5.453 4.96e-08 ***

10. poly(Age, 3)2 ? 6.0044 ? ? 3.0021 ? 2.000 0.045491 *

11. poly(Age, 3)3 ?-5.9181 ? ? 3.0992 ?-1.910 0.056188 .

12. SibSp ? ? ? ? ?-0.5041 ? ? 0.1317 ?-3.828 0.000129 ***

13. ---

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

15. 


16. (Dispersion parameter for binomial family taken to be 1)

17. 


18. Null deviance: 964.52 ?on 713 ?degrees of freedom

19. Residual deviance: 627.55 ?on 706 ?degrees of freedom

20. AIC: 643.55

21. 


22. Number of Fisher Scoring iterations: 5

但是解釋參數(shù)變得很復(fù)雜。我們注意到三階項(xiàng)在這里很重要，因此我們將手動(dòng)進(jìn)行回歸

1. 


2. Coefficients:

3. Estimate Std. Error z value Pr(>|z|)

4. (Intercept) ?5.616e+00 ?6.565e-01 ? 8.554 ?< 2e-16 ***

5. Pclass2 ? ? -1.360e+00 ?2.842e-01 ?-4.786 ?1.7e-06 ***

6. Pclass3 ? ? -2.557e+00 ?2.853e-01 ?-8.962 ?< 2e-16 ***

7. Sexmale ? ? -2.658e+00 ?2.176e-01 -12.216 ?< 2e-16 ***

8. Age ? ? ? ? -1.905e-01 ?5.528e-02 ?-3.446 0.000569 ***

9. I(Age^2) ? ? 4.290e-03 ?1.854e-03 ? 2.314 0.020669 *

10. I(Age^3) ? ?-3.520e-05 ?1.843e-05 ?-1.910 0.056188 .

11. SibSp ? ? ? -5.041e-01 ?1.317e-01 ?-3.828 0.000129 ***

12. ---

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

14. 


15. (Dispersion parameter for binomial family taken to be 1)

16. 


17. Null deviance: 964.52 ?on 713 ?degrees of freedom

18. Residual deviance: 627.55 ?on 706 ?degrees of freedom

19. AIC: 643.55

20. 


21. Number of Fisher Scoring iterations: 5

可以看到，p值是相同的。簡(jiǎn)而言之，將年齡轉(zhuǎn)換為年齡的非線性函數(shù)是有意義的。可以可視化此函數(shù)

1. 


2. plot(xage,yage,xlab="Age",ylab="",type="l")

?

實(shí)際上，我們可以使用樣條曲線。廣義可加模型（?gam?）是完美的可視化工具

1. 


2. 


3. (Dispersion Parameter for binomial family taken to be 1)

4. 


5. Null Deviance: 964.516 on 713 degrees of freedom

6. Residual Deviance: 627.5525 on 706 degrees of freedom

7. AIC: 643.5525

8. 177 observations deleted due to missingness

9. 


10. Number of Local Scoring Iterations: 4

11. 


12. Anova for Parametric Effects

13. Df Sum Sq Mean Sq ?F value ? ?Pr(>F)

14. Pclass ? ? ?2 ?26.72 ?13.361 ?11.3500 1.407e-05 ***

15. Sex ? ? ? ? 1 131.57 131.573 111.7678 < 2.2e-16 ***

16. bs(Age) ? ? 3 ?22.76 ? 7.588 ? 6.4455 0.0002620 ***

17. SibSp ? ? ? 1 ?14.66 ?14.659 ?12.4525 0.0004445 ***

18. Residuals 706 831.10 ? 1.177

19. ---

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

我們可以看到年齡變量的變換，

?

并且我們發(fā)現(xiàn)變換接近于我們的3階多項(xiàng)式。我們可以添加置信帶，從而可以驗(yàn)證該函數(shù)不是真正的線性

?

我們現(xiàn)在有三個(gè)模型。最后給出了兩個(gè)模擬乘客的預(yù)測(cè)，

1. predict(reg,newdata=newbase,type="response")

2. 1 ? ? ? ? 2

3. 0.9605736 0.1368988

4. predict(reg3,newdata=newbase,type="response")

5. 1 ? ? ? ? 2

6. 0.9497834 0.1218426

7. predict(regam,newdata=newbase,type="response")

8. 1 ? ? ? ? 2

9. 0.9497834 0.1218426

可以看到萊昂納多·迪卡普里奧（??Leonardo DiCaprio）?有大約12％的幸存機(jī)會(huì)（考慮到他的年齡，他有三等票，而且船上沒(méi)有家人）。

最受歡迎的見(jiàn)解

1.用SPSS估計(jì)HLM層次線性模型模型

2.R語(yǔ)言線性判別分析（LDA），二次判別分析（QDA）和正則判別分析（RDA）

3.基于R語(yǔ)言的lmer混合線性回歸模型

4.R語(yǔ)言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真分析

5.在r語(yǔ)言中使用GAM（廣義相加模型）進(jìn)行電力負(fù)荷時(shí)間序列分析

6.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM

7.R語(yǔ)言中的嶺回歸、套索回歸、主成分回歸：線性模型選擇和正則化

8.R語(yǔ)言用線性回歸模型預(yù)測(cè)空氣質(zhì)量臭氧數(shù)據(jù)

9.R語(yǔ)言分層線性模型案例