d0 <- read.csv(file = 'https://stats.dip.jp/01_ds/data/titanic_data_jp.csv')
library(DT)
datatable(d0)
d0$生死 <- ifelse (d0$生死== '生存',1 ,0 )
summary(d0)
## 乗船番号 生死 客室等級 名前
## Length:714 Min. :0.0000 Length:714 Length:714
## Class :character 1st Qu.:0.0000 Class :character Class :character
## Mode :character Median :0.0000 Mode :character Mode :character
## Mean :0.4062
## 3rd Qu.:1.0000
## Max. :1.0000
## 性別 年齢 兄弟配偶者数 親子数
## Length:714 Min. : 0.42 Min. :0.0000 Min. :0.0000
## Class :character 1st Qu.:20.12 1st Qu.:0.0000 1st Qu.:0.0000
## Mode :character Median :28.00 Median :0.0000 Median :0.0000
## Mean :29.70 Mean :0.5126 Mean :0.4314
## 3rd Qu.:38.00 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :80.00 Max. :5.0000 Max. :6.0000
## 運賃 乗船地
## Min. : 0.00 Length:714
## 1st Qu.: 8.05 Class :character
## Median : 15.74 Mode :character
## Mean : 34.69
## 3rd Qu.: 33.38
## Max. :512.33
head(d0)
## 乗船番号 生死 客室等級 名前 性別 年齢 兄弟配偶者数 親子数 運賃 乗船地
## 1 No.1 0 3等 Braund 男性 22 1 0 7.2500 S
## 2 No.2 1 1等 Cumings 女性 38 1 0 71.2833 C
## 3 No.3 1 3等 Heikkinen 女性 26 0 0 7.9250 S
## 4 No.4 1 1等 Futrelle 女性 35 1 0 53.1000 S
## 5 No.5 0 3等 Allen 男性 35 0 0 8.0500 S
## 6 No.7 0 1等 McCarthy 男性 54 0 0 51.8625 S
tail(d0)
## 乗船番号 生死 客室等級 名前 性別 年齢 兄弟配偶者数 親子数 運賃 乗船地
## 709 No.885 0 3等 Sutehall 男性 25 0 0 7.050 S
## 710 No.886 0 3等 Rice 女性 39 0 5 29.125 Q
## 711 No.887 0 2等 Montvila 男性 27 0 0 13.000 S
## 712 No.888 1 1等 Graham 女性 19 0 0 30.000 S
## 713 No.890 1 1等 Behr 男性 26 0 0 30.000 C
## 714 No.891 0 3等 Dooley 男性 32 0 0 7.750 Q
(n <- nrow(d0))
## [1] 714
set.seed(5)
ii <- sample(x = 1:n, size = floor(0.8*n))
d.tr <- d0[ ii, ]
d.te <- d0[-ii, ]
(n.tr <- nrow(d.tr))
## [1] 571
(n.te <- nrow(d.te))
## [1] 143
fit.all <- glm(生死~客室等級+ 性別+ 年齢+ 兄弟配偶者数+ 親子数+ 運賃+ 乗船地, data = d.tr, family = 'binomial')
summary(fit.all)
##
## Call:
## glm(formula = 生死 ~ 客室等級 + 性別 + 年齢 + 兄弟配偶者数 +
## 親子数 + 運賃 + 乗船地, family = "binomial", data = d.tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6437 -0.6306 -0.3804 0.6662 2.4121
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 17.215715 605.555989 0.028 0.9773
## 客室等級2等 -1.540699 0.381403 -4.040 5.36e-05 ***
## 客室等級3等 -2.721618 0.398577 -6.828 8.59e-12 ***
## 性別男性 -2.512823 0.245347 -10.242 < 2e-16 ***
## 年齢 -0.046152 0.009290 -4.968 6.76e-07 ***
## 兄弟配偶者数 -0.353650 0.138036 -2.562 0.0104 *
## 親子数 -0.044062 0.139460 -0.316 0.7520
## 運賃 -0.003207 0.003342 -0.960 0.3372
## 乗船地C -12.272960 605.555797 -0.020 0.9838
## 乗船地Q -13.069235 605.556026 -0.022 0.9828
## 乗船地S -12.731745 605.555772 -0.021 0.9832
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 774.32 on 570 degrees of freedom
## Residual deviance: 519.84 on 560 degrees of freedom
## AIC: 541.84
##
## Number of Fisher Scoring iterations: 13
fit <- glm(生死 ~ 客室等級 + 性別 + 年齢 + 兄弟配偶者数, data = d.tr, family = 'binomial')
summary(fit)
##
## Call:
## glm(formula = 生死 ~ 客室等級 + 性別 + 年齢 + 兄弟配偶者数,
## family = "binomial", data = d.tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7973 -0.6277 -0.3904 0.6724 2.3959
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.389219 0.504021 8.708 < 2e-16 ***
## 客室等級2等 -1.508696 0.318108 -4.743 2.11e-06 ***
## 客室等級3等 -2.651199 0.313919 -8.445 < 2e-16 ***
## 性別男性 -2.480583 0.235670 -10.526 < 2e-16 ***
## 年齢 -0.045982 0.009104 -5.051 4.40e-07 ***
## 兄弟配偶者数 -0.405085 0.129990 -3.116 0.00183 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 774.32 on 570 degrees of freedom
## Residual deviance: 523.85 on 565 degrees of freedom
## AIC: 535.85
##
## Number of Fisher Scoring iterations: 5
fit <- glm(生死 ~ 客室等級 + 性別 + 年齢, data = d.tr, family = 'binomial')
summary(fit)
##
## Call:
## glm(formula = 生死 ~ 客室等級 + 性別 + 年齢, family = "binomial",
## data = d.tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7089 -0.6917 -0.4177 0.7013 2.4047
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.716434 0.437531 8.494 < 2e-16 ***
## 客室等級2等 -1.360400 0.307460 -4.425 9.66e-06 ***
## 客室等級3等 -2.557113 0.307238 -8.323 < 2e-16 ***
## 性別男性 -2.347847 0.225497 -10.412 < 2e-16 ***
## 年齢 -0.036569 0.008345 -4.382 1.18e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 774.32 on 570 degrees of freedom
## Residual deviance: 534.59 on 566 degrees of freedom
## AIC: 544.59
##
## Number of Fisher Scoring iterations: 4
d.new <- data.frame(客室等級 = '2等', 年齢 = 45, 性別 = '女性')
p.hat <- predict(fit, type = 'response', newdata = d.new)
sprintf('生存確率：%2.1f％', p.hat * 100)
## [1] "生存確率：67.0％"
p.hat <- predict(fit, type = 'response', newdata = d.te)
threshold <- 0.5
is.pred <- p.hat > threshold
is.ref <- d.te$生死 == 1
table(予測値 = is.pred, 真値 = is.ref)
## 真値
## 予測値 FALSE TRUE
## FALSE 77 9
## TRUE 12 45
is.ok <- is.pred == is.ref
n.ok <- sum(is.ok)
sprintf('生存予測精度：%2.1f％', n.ok / nrow(d.te) * 100)
## [1] "生存予測精度：85.3％"
library(caret)
## 要求されたパッケージ ggplot2 をロード中です
## 要求されたパッケージ lattice をロード中です
confusionMatrix(data = as.factor(is.pred),
reference = as.factor(is.ref))
## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 77 9
## TRUE 12 45
##
## Accuracy : 0.8531
## 95% CI : (0.7843, 0.9067)
## No Information Rate : 0.6224
## P-Value [Acc > NIR] : 1.037e-09
##
## Kappa : 0.691
##
## Mcnemar's Test P-Value : 0.6625
##
## Sensitivity : 0.8652
## Specificity : 0.8333
## Pos Pred Value : 0.8953
## Neg Pred Value : 0.7895
## Prevalence : 0.6224
## Detection Rate : 0.5385
## Detection Prevalence : 0.6014
## Balanced Accuracy : 0.8493
##
## 'Positive' Class : FALSE
##
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## 次のパッケージを付け加えます: 'pROC'
## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## cov, smooth, var
roc1 <- roc(response = d.te$生死, predict = p.hat,
plot = T, print.auc = T, grid = T, ci = T, auc.polygon=T)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
coords(roc1, 'best')
## threshold specificity sensitivity
## 1 0.4505387 0.8539326 0.8518519
c <- coords(roc1)
library(DT)
datatable(round(c, 3))