統計解析R演習4

注意事項

• 慶應義塾大学SFCで開講している「統計解析」の授業履修者向けの演習用ページです。
• 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.6.0を使っています。
• Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

パッケージのインストール

• 今回の演習では、以下のパッケージを使います。

install.packages("dplyr")
install.packages("magrittr")
install.packages("ggplot2")
install.packages("car")

library(dplyr)
library(magrittr)
library(ggplot2)
library(car)

データのダウンロード

• 今回の演習で用いるデータは、Titanicの生存者情報データです。
• kaggleのTitanic: Machine Learning from Disasterからダウンロードできます。
• Rのtitanicライブラリからも利用できます。

install.packages("titanic")
library(titanic)

• effectsライブラリのTitanicSurvivalデータも利用できます。
• 演習用データにも含まれています。

データの読み込み

• 使用するデータには、欠損値が含まれます。欠損値は、NAで定義されている場合と、空白となっている場合の両方があるので、na.strings=(c("NA", ""))により欠損値を定義します。

kakei <- read.csv("./フォルダを指定/train.csv", stringsAsFactors = F,na.strings=(c("NA", "")))

head(train)

summary(train)

  PassengerId       Survived          Pclass          Name               Sex           
 Min.   :  1.0   Min.   :0.0000   Min.   :1.000   Length:891         Length:891        
 1st Qu.:223.5   1st Qu.:0.0000   1st Qu.:2.000   Class :character   Class :character  
 Median :446.0   Median :0.0000   Median :3.000   Mode  :character   Mode  :character  
 Mean   :446.0   Mean   :0.3838   Mean   :2.309                                        
 3rd Qu.:668.5   3rd Qu.:1.0000   3rd Qu.:3.000                                        
 Max.   :891.0   Max.   :1.0000   Max.   :3.000                                        
                                                                                       
      Age            SibSp           Parch           Ticket               Fare       
 Min.   : 0.42   Min.   :0.000   Min.   :0.0000   Length:891         Min.   :  0.00  
 1st Qu.:20.12   1st Qu.:0.000   1st Qu.:0.0000   Class :character   1st Qu.:  7.91  
 Median :28.00   Median :0.000   Median :0.0000   Mode  :character   Median : 14.45  
 Mean   :29.70   Mean   :0.523   Mean   :0.3816                      Mean   : 32.20  
 3rd Qu.:38.00   3rd Qu.:1.000   3rd Qu.:0.0000                      3rd Qu.: 31.00  
 Max.   :80.00   Max.   :8.000   Max.   :6.0000                      Max.   :512.33  
 NA's   :177                                                                         
    Cabin             Embarked        
 Length:891         Length:891        
 Class :character   Class :character  
 Mode  :character   Mode  :character  
                                      
                                      
                                      
                                      

欠損値の処理

• summary()関数だけでは欠損値の状況がわかりにくいので、以下の方法で欠損値を把握します。

na_count_train <- sapply(train, function(y) sum(is.na(y)))
na_count_train

PassengerId    Survived      Pclass        Name         Sex         Age       SibSp 
          0           0           0           0           0         177           0 
      Parch      Ticket        Fare       Cabin    Embarked 
          0           0           0         687           2 

• 欠損値の処理には、主に欠損値を除去する方法と、何らかの方法で欠損値に値を代入する方法とがあります。ここでは、欠損値を除去する方法を採用することにします。
• na.omit()関数で欠損値を除去します。欠損値を除去したデータを、train2とします。

train2 <- na.omit(train)

• 欠損値除去後のデータは以下のようになります。

summary(train2)

  PassengerId    Survived Pclass      Name               Sex          Age       
 Min.   :  2.0   0: 60    1:158   Length:183         female:88   Min.   : 0.92  
 1st Qu.:263.5   1:123    2: 15   Class :character   male  :95   1st Qu.:24.00  
 Median :457.0            3: 10   Mode  :character               Median :36.00  
 Mean   :455.4                                                   Mean   :35.67  
 3rd Qu.:676.0                                                   3rd Qu.:47.50  
 Max.   :890.0                                                   Max.   :80.00  
                                                                                
     SibSp            Parch           Ticket               Fare       
 Min.   :0.0000   Min.   :0.0000   Length:183         Min.   :  0.00  
 1st Qu.:0.0000   1st Qu.:0.0000   Class :character   1st Qu.: 29.70  
 Median :0.0000   Median :0.0000   Mode  :character   Median : 57.00  
 Mean   :0.4645   Mean   :0.4754                      Mean   : 78.68  
 3rd Qu.:1.0000   3rd Qu.:1.0000                      3rd Qu.: 90.00  
 Max.   :3.0000   Max.   :4.0000                      Max.   :512.33  
                                                                      
    Cabin             Embarked          CabinLetter  CabinNumber      Title   
 Length:183         Length:183         C      :51   6      :  6   Mr     :81  
 Class :character   Class :character   B      :43   33     :  6   Miss   :44  
 Mode  :character   Mode  :character   D      :31   2      :  5   Mrs    :38  
                                       E      :30   20     :  5   Master : 7  
                                       A      :12   22     :  5   Dr     : 3  
                                       F      :11   (Other):149   Major  : 2  
                                       (Other): 5   NA's   :  7   (Other): 8  
     Surname   
 Carter  :  4  
 Fortune :  4  
 Allison :  3  
 Graham  :  3  
 Navratil:  3  
 Newell  :  3  
 (Other) :163  

head(train2)

変数の変換

• 以下の要領で、いくつかの変数をファクター変数に変換します。
• 変数の変換は、このサイトを参考にしました。

train2$Survived <- factor(train2$Survived)
train2$Pclass <- factor(train2$Pclass)
train2$Sex <- factor(train2$Sex)
train2$CabinLetter = factor(substr(train2$Cabin, 1, 1))
train2$CabinNumber = factor(as.numeric(substr(train2$Cabin, 2, 4)))
train2$Title = factor(gsub('(.*, )|(\\..*)', '', train2$Name))
train2$Surname <- factor(sapply(train2$Name, function(x){x=as.character(x); strsplit(x, split = '[,.]')[[1]][1];}))
train2$CabinNumber = factor(train2$CabinNumber)

ロジスティック回帰モデルの推定

• glm()関数を用いてロジスティック回帰モデルを推定します。リンク関数をfamily=binomial(link='logit')と指定します。

glm.logit1 <- glm(Survived~Pclass+Sex+Age+SibSp+Parch+Fare+Embarked+CabinLetter+Title,
            data=train2, family=binomial(link='logit'))
summary(glm.logit1)

Call:
glm(formula = Survived ~ Pclass + Sex + Age + SibSp + Parch + 
    Fare + Embarked + CabinLetter + Title, family = binomial(link = "logit"), 
    data = train2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3499  -0.6664   0.1865   0.3934   1.9563  

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)
(Intercept)        2.044e+00  9.224e+03   0.000    1.000
Pclass2           -6.710e-01  1.260e+00  -0.533    0.594
Pclass3           -1.754e-01  1.582e+00  -0.111    0.912
Sexmale           -1.781e+01  6.523e+03  -0.003    0.998
Age               -2.859e-02  1.826e-02  -1.566    0.117
SibSp              2.954e-01  4.173e-01   0.708    0.479
Parch             -5.683e-01  4.055e-01  -1.401    0.161
Fare               3.542e-03  3.559e-03   0.995    0.320
EmbarkedQ         -2.064e+00  2.109e+00  -0.979    0.328
EmbarkedS         -5.330e-01  5.433e-01  -0.981    0.327
CabinLetterB      -2.400e-01  9.574e-01  -0.251    0.802
CabinLetterC      -1.025e+00  9.397e-01  -1.091    0.275
CabinLetterD       3.513e-01  9.649e-01   0.364    0.716
CabinLetterE       6.872e-01  9.481e-01   0.725    0.469
CabinLetterF      -1.054e+00  1.775e+00  -0.594    0.552
CabinLetterG      -3.070e+00  2.246e+00  -1.367    0.172
CabinLetterT      -1.783e+01  6.523e+03  -0.003    0.998
TitleCol           3.581e+01  9.224e+03   0.004    0.997
TitleDr            1.801e+01  6.523e+03   0.003    0.998
TitleLady          1.746e+01  1.130e+04   0.002    0.999
TitleMajor         1.822e+01  6.523e+03   0.003    0.998
TitleMaster        3.627e+01  6.955e+03   0.005    0.996
TitleMiss          2.570e+00  9.224e+03   0.000    1.000
TitleMlle          1.769e+01  1.028e+04   0.002    0.999
TitleMme           1.720e+01  1.130e+04   0.002    0.999
TitleMr            1.673e+01  6.523e+03   0.003    0.998
TitleMrs           2.761e+00  9.224e+03   0.000    1.000
TitleSir           3.524e+01  9.224e+03   0.004    0.997
Titlethe Countess  1.793e+01  1.130e+04   0.002    0.999

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 231.55  on 182  degrees of freedom
Residual deviance: 132.99  on 154  degrees of freedom
AIC: 190.99

Number of Fisher Scoring iterations: 17

• carパッケージのvif()関数を使ってVIFを確認します。

car::vif(glm.logit1)

                    GVIF Df GVIF^(1/(2*Df))
Pclass      5.672965e+00  2        1.543308
Sex         1.533092e+08  1    12381.807908
Age         1.431840e+00  1        1.196595
SibSp       1.350778e+00  1        1.162230
Parch       1.523750e+00  1        1.234403
Fare        1.647715e+00  1        1.283634
Embarked    2.086932e+00  2        1.201924
CabinLetter 1.096882e+01  7        1.186580
Title       2.666478e+08 12        2.244299

• そこで次に、次式のモデルを推定します。

glm.logit2 <- glm(Survived~Pclass+Sex+Age+SibSp+Parch+Fare+Embarked,
            data=train2, family=binomial(link='logit'))
summary(glm.logit2)

Call:
glm(formula = Survived ~ Pclass + Sex + Age + SibSp + Parch + 
    Fare + Embarked, family = binomial(link = "logit"), data = train2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6730  -0.7700   0.3052   0.5401   2.0412  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.451473   0.941136   4.730 2.25e-06 ***
Pclass2      0.059891   0.845412   0.071  0.94352    
Pclass3     -1.512335   0.925014  -1.635  0.10206    
Sexmale     -2.904210   0.503842  -5.764 8.21e-09 ***
Age         -0.039658   0.014788  -2.682  0.00732 ** 
SibSp        0.176325   0.356422   0.495  0.62081    
Parch       -0.438347   0.323717  -1.354  0.17570    
Fare         0.001621   0.003024   0.536  0.59188    
EmbarkedQ   -1.882895   1.972139  -0.955  0.33971    
EmbarkedS   -0.373572   0.448718  -0.833  0.40511    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 231.55  on 182  degrees of freedom
Residual deviance: 159.76  on 173  degrees of freedom
AIC: 179.76

Number of Fisher Scoring iterations: 5

ロジスティック回帰モデルのオッズ比と95%信頼区間

• オッズ比は以下のようにして出力できます。

exp(coef(glm.logit2))

(Intercept)     Pclass2     Pclass3     Sexmale         Age       SibSp       Parch 
85.75313507  1.06172061  0.22039485  0.05479204  0.96111781  1.19282621  0.64510159 
       Fare   EmbarkedQ   EmbarkedS 
 1.00162224  0.15214902  0.68827120 

• 95%信頼区間は以下のようにして出力できます。

exp(confint(glm.logit2))

Waiting for profiling to be done...

                   2.5 %      97.5 %
(Intercept) 15.265845913 624.2027347
Pclass2      0.214503960   6.3519834
Pclass3      0.033190848   1.2950601
Sexmale      0.018429484   0.1366699
Age          0.932319986   0.9883203
SibSp        0.599045274   2.4548105
Parch        0.332704463   1.1927129
Fare         0.995873242   1.0086291
EmbarkedQ    0.002980473   7.9613402
EmbarkedS    0.281694081   1.6530164

• 以下の方法でも出力できます。

library(broom)
broom::tidy(glm.logit2, conf.int = TRUE, exponentiate = TRUE)

モデルを用いた予測と的中率

• predict()関数を用いてモデル上の（理論的な）生存確率を予測できます。

predict(glm.logit2, type="response")

         2          4          7         11         12         22         24         28 
0.96217871 0.95037004 0.29239251 0.89770494 0.86065670 0.47658486 0.53016350 0.62216968 
        53         55         63         67         76         89         93         97 
0.94315564 0.20287680 0.42572841 0.95278692 0.21120090 0.96246937 0.40730786 0.22929355 
        98        103        111        119        124        125        137        138 
0.57431375 0.50694849 0.35298097 0.63606581 0.94634176 0.21739232 0.92345909 0.49218310 
       140        149        152        171        175        178        184        194 
0.67345220 0.25951368 0.97039320 0.23305068 0.34890716 0.92519356 0.76340557 0.70986489 
       195        196        206        210        216        219        225        231 
0.93999469 0.91596858 0.88744061 0.50279481 0.97293304 0.96463648 0.49722097 0.95264151 
       246        249        252        253        258        263        264        269 
0.17050350 0.38449735 0.76321952 0.22407266 0.95383910 0.26474381 0.39827886 0.83035193 
       270        274        276        292        293        298        300        306 
0.94832444 0.42304140 0.86785423 0.98239695 0.54991420 0.97191044 0.91919235 0.66429339 
       308        310        311        312        319        320        326        328 
0.98417365 0.96623564 0.97428094 0.97439103 0.90370053 0.94380722 0.96244339 0.93884075 
       330        332        333        337        338        340        341        342 
0.96989722 0.35788266 0.37216944 0.57636739 0.95450096 0.36507820 0.71796431 0.96101026 
       346        357        367        370        371        378        391        394 
0.96109925 0.94563939 0.91453239 0.97371202 0.69466017 0.48565882 0.31867519 0.98014545 
       395        413        430        431        435        436        439        446 
0.68225218 0.82949650 0.16872907 0.52654836 0.36765518 0.95332851 0.07479936 0.56735021 
       450        453        454        457        461        463        474        485 
0.30172375 0.59928678 0.48119958 0.20406753 0.33472918 0.34799955 0.97395598 0.70677314 
       487        488        493        497        499        505        506        513 
0.95311640 0.33075394 0.27726322 0.93170088 0.93287565 0.97297452 0.76607104 0.44734633 
       516        517        521        524        537        540        541        545 
0.34635401 0.94302225 0.95433612 0.91389439 0.36172209 0.94172451 0.86859499 0.47830871 
       551        557        559        572        573        578        582        584 
0.54391409 0.94204514 0.91669848 0.91774579 0.44738641 0.94257643 0.94388022 0.54602025 
       586        588        592        600        610        619        622        626 
0.93191362 0.27570251 0.93657857 0.46819538 0.93936434 0.98123240 0.44264203 0.23270922 
       628        631        633        642        646        648        660        663 
0.96679997 0.12451568 0.58118388 0.97371202 0.48609413 0.35067827 0.19061746 0.34326580 
       672        680        682        690        691        699        700        701 
0.55104568 0.62519662 0.64585323 0.96730340 0.55304988 0.38265782 0.12005676 0.98638333 
       702        708        711        713        716        717        718        725 
0.45717018 0.38951575 0.97287808 0.38478327 0.25355294 0.96488088 0.95622946 0.59032589 
       731        738        742        743        746        749        752        760 
0.96339831 0.72902281 0.51252913 0.97124858 0.14813215 0.66431408 0.26998064 0.94830852 
       764        766        773        780        782        783        790        797 
0.89513902 0.91356941 0.86924370 0.90693807 0.97521520 0.51804898 0.46290784 0.89813513 
       803        807        810        821        824        836        854        858 
0.55764023 0.40782032 0.95398030 0.87047937 0.74586244 0.94145169 0.95559087 0.30877612 
       863        868        872        873        880        888        890 
0.90170655 0.50653286 0.88464021 0.46831027 0.87292303 0.96685064 0.63755655 

• 生存者の的中率は、以下のようにして計算できます。モデルglm.logit2を用いた場合には、的中率が77.596%となりました。

train3 <- train2 %>%
  dplyr::mutate(predict=predict(glm.logit2, type="response"),
         survive=dplyr::if_else(predict > 0.5, 1, 0))
sum(train3$Survived == train3$survive)/nrow(train3)*100

[1] 77.59563

プロビット回帰モデルの推定

• glm()関数を用いてロジスティック回帰モデルを推定します。リンク関数をfamily=binomial(link='probit')と指定します。

glm.probit2 <- glm(Survived~Pclass+Sex+Age+SibSp+Parch+Fare+Embarked,
                  data=train2, family=binomial(link='probit'))
summary(glm.probit2)

Call:
glm(formula = Survived ~ Pclass + Sex + Age + SibSp + Parch + 
    Fare + Embarked, family = binomial(link = "probit"), data = train2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.7068  -0.7768   0.2880   0.5752   1.9984  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.530143   0.512473   4.937 7.93e-07 ***
Pclass2      0.066280   0.481734   0.138  0.89057    
Pclass3     -0.846859   0.520452  -1.627  0.10370    
Sexmale     -1.668207   0.261009  -6.391 1.64e-10 ***
Age         -0.022104   0.008376  -2.639  0.00831 ** 
SibSp        0.108097   0.203716   0.531  0.59568    
Parch       -0.292705   0.183126  -1.598  0.10996    
Fare         0.001098   0.001785   0.615  0.53827    
EmbarkedQ   -1.106024   1.082945  -1.021  0.30711    
EmbarkedS   -0.226086   0.260391  -0.868  0.38526    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 231.55  on 182  degrees of freedom
Residual deviance: 159.74  on 173  degrees of freedom
AIC: 179.74

Number of Fisher Scoring iterations: 6

補対数対数モデルの推定

• glm()関数を用いてロジスティック回帰モデルを推定します。リンク関数をfamily = binomial('cloglog')と指定します。

glm.cloglog2 <- glm(Survived~Pclass+Sex+Age+SibSp+Parch+Fare+Embarked,
                   data=train2, family=binomial(link='cloglog'))
summary(glm.cloglog2)

Call:
glm(formula = Survived ~ Pclass + Sex + Age + SibSp + Parch + 
    Fare + Embarked, family = binomial(link = "cloglog"), data = train2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.7479  -0.8394   0.2412   0.6222   1.7647  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.988022   0.499628   3.979 6.92e-05 ***
Pclass2      0.103065   0.451284   0.228   0.8193    
Pclass3     -0.912063   0.562415  -1.622   0.1049    
Sexmale     -1.591601   0.237599  -6.699 2.10e-11 ***
Age         -0.020518   0.008776  -2.338   0.0194 *  
SibSp        0.106497   0.191274   0.557   0.5777    
Parch       -0.351570   0.186834  -1.882   0.0599 .  
Fare         0.001433   0.001677   0.854   0.3929    
EmbarkedQ   -0.916746   1.082665  -0.847   0.3971    
EmbarkedS   -0.238959   0.260193  -0.918   0.3584    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 231.55  on 182  degrees of freedom
Residual deviance: 159.98  on 173  degrees of freedom
AIC: 179.98

Number of Fisher Scoring iterations: 9