Data
d0 <- read.csv(file = 'https://stats.dip.jp/01_ds/data/titanic_data_jp.csv')
library(DT)
datatable(d0)summary(d0)## 乗船番号 生死 客室等級 名前
## Length:714 Length:714 Length:714 Length:714
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## 性別 年齢 兄弟配偶者数 親子数
## 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.33head(d0) ## 乗船番号 生死 客室等級 名前 性別 年齢 兄弟配偶者数 親子数 運賃 乗船地
## 1 No.1 死亡 3等 Braund 男性 22 1 0 7.2500 S
## 2 No.2 生存 1等 Cumings 女性 38 1 0 71.2833 C
## 3 No.3 生存 3等 Heikkinen 女性 26 0 0 7.9250 S
## 4 No.4 生存 1等 Futrelle 女性 35 1 0 53.1000 S
## 5 No.5 死亡 3等 Allen 男性 35 0 0 8.0500 S
## 6 No.7 死亡 1等 McCarthy 男性 54 0 0 51.8625 Stail(d0)## 乗船番号 生死 客室等級 名前 性別 年齢 兄弟配偶者数 親子数 運賃 乗船地
## 709 No.885 死亡 3等 Sutehall 男性 25 0 0 7.050 S
## 710 No.886 死亡 3等 Rice 女性 39 0 5 29.125 Q
## 711 No.887 死亡 2等 Montvila 男性 27 0 0 13.000 S
## 712 No.888 生存 1等 Graham 女性 19 0 0 30.000 S
## 713 No.890 生存 1等 Behr 男性 26 0 0 30.000 C
## 714 No.891 死亡 3等 Dooley 男性 32 0 0 7.750 Qd1 <- d0[,-1][,-3][,-8]
datatable(d1)summary(d1)## 生死 客室等級 性別 年齢
## Length:714 Length:714 Length:714 Min. : 0.42
## Class :character Class :character Class :character 1st Qu.:20.12
## Mode :character Mode :character Mode :character Median :28.00
## Mean :29.70
## 3rd Qu.:38.00
## Max. :80.00
## 兄弟配偶者数 親子数 運賃
## Min. :0.0000 Min. :0.0000 Min. : 0.00
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 8.05
## Median :0.0000 Median :0.0000 Median : 15.74
## Mean :0.5126 Mean :0.4314 Mean : 34.69
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.: 33.38
## Max. :5.0000 Max. :6.0000 Max. :512.33d1$生死 <- ifelse(d1$生死 == '生存', 1, 0)
datatable(d1)summary(d1)## 生死 客室等級 性別 年齢
## Min. :0.0000 Length:714 Length:714 Min. : 0.42
## 1st Qu.:0.0000 Class :character Class :character 1st Qu.:20.12
## Median :0.0000 Mode :character Mode :character Median :28.00
## Mean :0.4062 Mean :29.70
## 3rd Qu.:1.0000 3rd Qu.:38.00
## Max. :1.0000 Max. :80.00
## 兄弟配偶者数 親子数 運賃
## Min. :0.0000 Min. :0.0000 Min. : 0.00
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 8.05
## Median :0.0000 Median :0.0000 Median : 15.74
## Mean :0.5126 Mean :0.4314 Mean : 34.69
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.: 33.38
## Max. :5.0000 Max. :6.0000 Max. :512.33d1$性別 <- ifelse(d1$性別 == '男性', 1, 0)
d1$客室等級 <- ifelse(d1$客室等級 == '1等', 1,ifelse(d1$客室等級=='2等', 2, 3))
d11 <- na.omit(d1)Clean Data
datatable(d11)summary(d11)## 生死 客室等級 性別 年齢
## Min. :0.0000 Min. :1.000 Min. :0.0000 Min. : 0.42
## 1st Qu.:0.0000 1st Qu.:1.000 1st Qu.:0.0000 1st Qu.:20.12
## Median :0.0000 Median :2.000 Median :1.0000 Median :28.00
## Mean :0.4062 Mean :2.237 Mean :0.6345 Mean :29.70
## 3rd Qu.:1.0000 3rd Qu.:3.000 3rd Qu.:1.0000 3rd Qu.:38.00
## Max. :1.0000 Max. :3.000 Max. :1.0000 Max. :80.00
## 兄弟配偶者数 親子数 運賃
## Min. :0.0000 Min. :0.0000 Min. : 0.00
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 8.05
## Median :0.0000 Median :0.0000 Median : 15.74
## Mean :0.5126 Mean :0.4314 Mean : 34.69
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.: 33.38
## Max. :5.0000 Max. :6.0000 Max. :512.33(n <- nrow(d0))## [1] 714Data2
set.seed(5)
library(rsample)
d.trte <- initial_split(d11, prop = 4/5, strata = 生死)
d.trte## <Training/Testing/Total>
## <571/143/714>d.tr <- training(d.trte)
d.te <- testing (d.trte)
nrow(d.tr)## [1] 571nrow(d.te)## [1] 143Fitting
fit.all <- glm(生死 ~ ., data = d.tr, family = 'binomial')
summary(fit.all)##
## Call:
## glm(formula = 生死 ~ ., family = "binomial", data = d.tr)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.4276535 0.6619771 8.199 2.42e-16 ***
## 客室等級 -1.2569372 0.1801565 -6.977 3.02e-12 ***
## 性別 -2.7319807 0.2495250 -10.949 < 2e-16 ***
## 年齢 -0.0413017 0.0089734 -4.603 4.17e-06 ***
## 兄弟配偶者数 -0.3312736 0.1392555 -2.379 0.0174 *
## 親子数 -0.0696349 0.1360334 -0.512 0.6087
## 運賃 0.0008852 0.0024054 0.368 0.7129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 771.40 on 570 degrees of freedom
## Residual deviance: 500.79 on 564 degrees of freedom
## AIC: 514.79
##
## Number of Fisher Scoring iterations: 51.
運賃は生死を分ける要素ではない
Model upgrade
fit <- glm(生死 ~ 客室等級 + 性別 + 年齢, data = d.tr, family = 'binomial')
summary(fit)##
## Call:
## glm(formula = 生死 ~ 客室等級 + 性別 + 年齢, family = "binomial",
## data = d.tr)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.053870 0.564456 8.954 < 2e-16 ***
## 客室等級 -1.279224 0.156612 -8.168 3.13e-16 ***
## 性別 -2.622552 0.235270 -11.147 < 2e-16 ***
## 年齢 -0.035297 0.008464 -4.170 3.04e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 771.40 on 570 degrees of freedom
## Residual deviance: 508.68 on 567 degrees of freedom
## AIC: 516.68
##
## Number of Fisher Scoring iterations: 5fit1 <- glm(生死 ~ 客室等級 + 性別 + 年齢 + 兄弟配偶者数, data = d.tr, family = 'binomial')
summary(fit)##
## Call:
## glm(formula = 生死 ~ 客室等級 + 性別 + 年齢, family = "binomial",
## data = d.tr)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.053870 0.564456 8.954 < 2e-16 ***
## 客室等級 -1.279224 0.156612 -8.168 3.13e-16 ***
## 性別 -2.622552 0.235270 -11.147 < 2e-16 ***
## 年齢 -0.035297 0.008464 -4.170 3.04e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 771.40 on 570 degrees of freedom
## Residual deviance: 508.68 on 567 degrees of freedom
## AIC: 516.68
##
## Number of Fisher Scoring iterations: 5#大差ないのでどちらとも保留2.
2等客室にいる45才の女性
d.new <- data.frame(年齢 = 45, 客室等級 = 2,性別 = 0)
p.hat <- predict(fit, type = 'response', newdata = d.new) # 予測確率
sprintf('生存確率：%2.1f％', p.hat * 100)## [1] "生存確率：71.2％"confusion matrix
fit2 <- glm(生死 ~ 客室等級 + 年齢, data = d.tr, family = 'binomial')
p.hat <- predict(fit2, type = 'response', newdata = d.te)
threshold <- 0.5 # 閾値（しきいち），カットオフ値（cut-off）とも呼ばれる。
is.pred <- p.hat > threshold
is.ref <- d.te$生死 == 1
table(予測値 = is.pred, 真値 = is.ref)## 真値
## 予測値 FALSE TRUE
## FALSE 68 28
## TRUE 17 30p.hat <- predict(fit, type = 'response', newdata = d.te)
threshold <- 0.5 # 閾値（しきいち），カットオフ値（cut-off）とも呼ばれる。
is.pred <- p.hat > threshold
is.ref <- d.te$生死 == 1
table(予測値 = is.pred, 真値 = is.ref)## 真値
## 予測値 FALSE TRUE
## FALSE 70 16
## TRUE 15 42p.hat <- predict(fit1, type = 'response', newdata = d.te)
threshold <- 0.5 # 閾値（しきいち），カットオフ値（cut-off）とも呼ばれる。
is.pred <- p.hat > threshold
is.ref <- d.te$生死 == 1
table(予測値 = is.pred, 真値 = is.ref)## 真値
## 予測値 FALSE TRUE
## FALSE 71 15
## TRUE 14 43これを使う
is.ok <- is.pred == is.ref
n.ok <- sum(is.ok)
sprintf('生死予測精度：%2.1f％', n.ok / nrow(d.te) * 100)## [1] "生死予測精度：79.7％"library(caret)
confusionMatrix(data = as.factor(is.pred),
reference = as.factor(is.ref))## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 71 15
## TRUE 14 43
##
## Accuracy : 0.7972
## 95% CI : (0.7219, 0.8598)
## No Information Rate : 0.5944
## P-Value [Acc > NIR] : 2.035e-07
##
## Kappa : 0.5783
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.8353
## Specificity : 0.7414
## Pos Pred Value : 0.8256
## Neg Pred Value : 0.7544
## Prevalence : 0.5944
## Detection Rate : 0.4965
## Detection Prevalence : 0.6014
## Balanced Accuracy : 0.7883
##
## 'Positive' Class : FALSE
## ROC curve
library(pROC)
roc1 <- roc(response = d.te$生死, predict = p.hat,
of = 'thresholds', thresholds = 'best', print.thres = 'best',
percent = F, plot = T, print.auc = T, grid = T, ci = T, auc.polygon=T)
# 閾値
coords(roc1, 'best')## threshold specificity sensitivity
## 1 0.6164488 0.9294118 0.67241383.
thresholdbest <- coords(roc1, 'best', ret = 'threshold')
thresholdbest## threshold
## 1 0.6164488c <- coords(roc1)
library(DT)
datatable(round(c, 3))