Categorical Outcome Training
Hiroki Matsui
Thursday, January 15, 2015
はじめに
この資料は、カテゴリカルアウトカムを解析するうえで医療系で多用されるLogistic Model の使用方法についてRを用いて解説した資料です。
灰色の枠で囲われた部分がRのコードになっており、Rのコンソールにコピペすることで実行できます。
パッケージのインストール
今回使用するパッケージは以下の5つです。
install.packages(c("ggplot2","car","GGally","epicalc","pROC"))ライブラリの読み込み
まずは、ライブラリを読み込みます。
library(ggplot2)
library(car)
library(GGally)
library(dplyr)
library(epicalc)
library(pROC)演習に入る前に
- 皆さん、解析ノートをつけましょう。
- もしくは、自分の解析ログを適切に保管しましょう。
- データ
- 解析コード
- 解析結果
- 上記三つが１セットとなってフォルダ内に保存されるようにしましょう。
RStudio のご紹介
- Rstudio はR をより使いやすくするための統合開発環境です。
- http://www.rstudio.com/
- データ→解析→アウトプット過程をすべて保存する機能を有します。
- 今回のプレゼンなどは全てRstudio 上で作成しています。
データの用意
Web 上に仮想データを用意しました。 下記からダウンロードできます。参考資料[1]のサイトを参考に講演者が独自生成しました。
- gre テスト成績（800点満点）
- gpa 成績（1-4の連続数）
- rank 最終的な成績ランク（1-4グレード）
- gre_math 数学テストの成績（greの一部）
- household 家庭の年収
- sex 性別（1=female）
- y 大学合格（1=不合格）
# 元データの読み込み
train_data<- read.csv("https://raw.githubusercontent.com/Hiro-macchan/Pubs_20150121/master/train_data.csv")もしも、ネットワークがうまくつながらない場合は以下のコードをお試し下さい。
train_data <- read.csv("train_data.csv")因果推論
データの特性を確認
データの中身を確認します。
# カテゴリカルな変数を因子に変換
train_data$y<- as.factor(train_data$y)
train_data$rank<- as.factor(train_data$rank)
train_data$sex <- as.factor(train_data$sex)
# Training data の特性をチェック
summary(train_data)## gre gpa rank gre_math
## Min. :103.7 Min. :2.094 1:1552 Min. : 37.12
## 1st Qu.:510.6 1st Qu.:3.138 2:3529 1st Qu.:152.49
## Median :587.9 Median :3.392 3:3450 Median :175.93
## Mean :585.8 Mean :3.385 4:1469 Mean :175.71
## 3rd Qu.:665.6 3rd Qu.:3.649 3rd Qu.:200.42
## Max. :800.0 Max. :4.000 Max. :268.02
## household sex y
## Min. : 30.99 0:4873 0:6433
## 1st Qu.:112.91 1:5127 1:3567
## Median :149.08
## Mean :161.12
## 3rd Qu.:194.47
## Max. :604.54# データのはじめ10行を見る。
head(train_data)## gre gpa rank gre_math household sex y
## 1 521.5534 3.688633 2 153.8448 63.40774 0 0
## 2 637.6283 3.342274 3 217.0735 162.98985 1 1
## 3 482.1153 3.137813 1 140.7382 139.19278 1 1
## 4 606.0971 2.693034 2 186.3475 99.70843 1 0
## 5 625.8819 3.387371 2 198.6858 110.03578 0 0
## 6 492.6032 4.000000 4 149.7200 123.13360 1 0# 変数間の関係を見る。
GGally::ggpairs(train_data)
- 因子の中で世帯所得が対数正規分布している様子がうかがえます。
- 因子間の関連を見ると、gre とgre_mathが相関している様子が見えます。
因子ごとにイベントの出現頻度を確認
# イベントごとに因子を確認
filter(train_data, y==1) %>% summary## gre gpa rank gre_math
## Min. :172.2 Min. :2.135 1: 318 Min. : 45.61
## 1st Qu.:488.6 1st Qu.:3.062 2:1046 1st Qu.:145.81
## Median :564.2 Median :3.330 3:1485 Median :168.23
## Mean :563.6 Mean :3.318 4: 718 Mean :168.84
## 3rd Qu.:638.0 3rd Qu.:3.571 3rd Qu.:191.85
## Max. :800.0 Max. :4.000 Max. :266.73
## household sex y
## Min. : 30.99 0:1033 0: 0
## 1st Qu.:116.12 1:2534 1:3567
## Median :151.77
## Mean :165.08
## 3rd Qu.:201.44
## Max. :585.85filter(train_data, y==0) %>% summary## gre gpa rank gre_math household
## Min. :103.7 Min. :2.094 1:1234 Min. : 37.12 Min. : 32.2
## 1st Qu.:524.5 1st Qu.:3.189 2:2483 1st Qu.:156.93 1st Qu.:110.8
## Median :602.6 Median :3.430 3:1965 Median :180.20 Median :147.1
## Mean :598.1 Mean :3.422 4: 751 Mean :179.52 Mean :158.9
## 3rd Qu.:677.6 3rd Qu.:3.684 3rd Qu.:204.05 3rd Qu.:191.5
## Max. :800.0 Max. :4.000 Max. :268.02 Max. :604.5
## sex y
## 0:3840 0:6433
## 1:2593 1: 0
##
##
##
## イベントを0 に絞っても因子が大まかに出現していたので、イベントの明らかな偏りはない。
モデルの構築
まず、何も考えずモデルを組んで見たうえで、すべての\(\beta=0\)という仮説を検証
# model構築
# モデル全体の尤度比検定
model <- glm(y ~ 1, data = train_data, family = binomial(logit))
model_0 <- glm(y ~ ., data = train_data, family = binomial(logit))
anova(model,model_0,test= "Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ 1
## Model 2: y ~ gre + gpa + rank + gre_math + household + sex
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 9999 13030
## 2 9991 11196 8 1834.1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1少なくとも、すべての変数を強制投入したモデルはすべての変数を投入しなかった場合に比べ、データの関係性をよくあらわしている。
多重共線性のチェック
gre は有意だがgre_mathは有意ではない。 greとgre_math の間に関連が強く多重共線性を起こしている?
# 変数の有意差とオッズ比
logistic.display(model_0)##
## Logistic regression predicting y : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI)
## gre (cont. var.) 0.9972 (0.9968,0.9976) 0.9975 (0.9961,0.9989)
##
## gpa (cont. var.) 0.45 (0.4,0.5) 0.39 (0.35,0.45)
##
## rank: ref.=1
## 2 1.63 (1.42,1.89) 1.81 (1.56,2.11)
## 3 2.93 (2.55,3.37) 3.5 (3.01,4.07)
## 4 3.71 (3.16,4.35) 4.49 (3.78,5.34)
##
## gre_math (cont. var.) 0.9912 (0.99,0.9924) 0.9972 (0.9928,1.0017)
##
## household (cont. var.) 1.0013 (1.0007,1.002) 1.002 (1.0013,1.0026)
##
## sex: 1 vs 0 3.63 (3.33,3.97) 4.12 (3.76,4.53)
##
## P(Wald's test) P(LR-test)
## gre (cont. var.) < 0.001 < 0.001
##
## gpa (cont. var.) < 0.001 < 0.001
##
## rank: ref.=1 < 0.001
## 2 < 0.001
## 3 < 0.001
## 4 < 0.001
##
## gre_math (cont. var.) 0.219 0.219
##
## household (cont. var.) < 0.001 < 0.001
##
## sex: 1 vs 0 < 0.001 < 0.001
##
## Log-likelihood = -5597.888
## No. of observations = 10000
## AIC value = 11213.7759# マルチコチェック
car::vif(model_0)## GVIF Df GVIF^(1/(2*Df))
## gre 11.756115 1 3.428719
## gpa 1.012011 1 1.005987
## rank 1.022506 3 1.003716
## gre_math 11.735988 1 3.425783
## household 1.005007 1 1.002500
## sex 1.033329 1 1.016528VIF から多重共線性が疑われたため、gre_mathをモデルから取り除く。
model_1 <- glm(y ~ gre + gpa + rank + household + sex , data = train_data, family = binomial(logit))
logistic.display(model_1)##
## Logistic regression predicting y : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI)
## gre (cont. var.) 0.9972 (0.9968,0.9976) 0.9967 (0.9963,0.9971)
##
## gpa (cont. var.) 0.45 (0.4,0.5) 0.39 (0.35,0.45)
##
## rank: ref.=1
## 2 1.63 (1.42,1.89) 1.82 (1.56,2.11)
## 3 2.93 (2.55,3.37) 3.5 (3.01,4.07)
## 4 3.71 (3.16,4.35) 4.5 (3.79,5.35)
##
## household (cont. var.) 1.0013 (1.0007,1.002) 1.002 (1.0013,1.0026)
##
## sex: 1 vs 0 3.63 (3.33,3.97) 4.13 (3.76,4.53)
##
## P(Wald's test) P(LR-test)
## gre (cont. var.) < 0.001 < 0.001
##
## gpa (cont. var.) < 0.001 < 0.001
##
## rank: ref.=1 < 0.001
## 2 < 0.001
## 3 < 0.001
## 4 < 0.001
##
## household (cont. var.) < 0.001 < 0.001
##
## sex: 1 vs 0 < 0.001 < 0.001
##
## Log-likelihood = -5598.6422
## No. of observations = 10000
## AIC value = 11213.2844anova(model_0,model_1,test= "Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ gre + gpa + rank + gre_math + household + sex
## Model 2: y ~ gre + gpa + rank + household + sex
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 9991 11196
## 2 9992 11197 -1 -1.5085 0.2194gre_math を入れたモデルと比してもデータのあらわし具合は変わらない。
変数の線形性を確認
世帯所得は明らかな対数正規分布をしており、年収200万円の家庭と年収1000万円の仮定では合格率に対する効果は異なることが予想される。
# 線形性の確認
household_10 <- with(train_data,cut(household,breaks=quantile(x = household,probs = seq(0,1,0.05)),include.lowest = T))
y <- train_data$y
household_10 <- as.numeric(household_10)
p <- apply(table(household_10,y),1,function(x){x[2]/(x[1]+x[2])})
x_3 <- max(household_10)*seq(0,1,0.05)[-1]
plot(x_3,p, main = "x_1 vs Proportion", xlab = "x_1", ylab = "Propotion")
plot(function(x){log(x)},0,2)
わかりにくいが正の対数関係が見て取れる。
model_2 <- glm(y ~ gre + gpa + rank + log(household) + sex , data = train_data, family = binomial(logit))
logistic.display(model_2)##
## OR lower95ci upper95ci Pr(>|Z|)
## gre 0.9966930 0.9962844 0.9971018 1.850799e-56
## gpa 0.3929785 0.3463812 0.4458443 1.144910e-47
## rank2 1.8189414 1.5626380 2.1172837 1.160319e-14
## rank3 3.5063929 3.0165952 4.0757178 4.882035e-60
## rank4 4.5146753 3.7999851 5.3637824 7.117532e-66
## log(household) 1.3885948 1.2410563 1.5536728 1.015239e-08
## sex1 4.1208913 3.7529401 4.5249178 1.643016e-193変数の変換を実施。カテゴリー化してもよい。
交互作用項の検討
変数間の関連で効果がModifyされていないか確認する。 具体的には交互作用項をモデルに投入し、投入前と比較する。
# 交互作用項の検討
model_3 <- glm(y ~ gre + gpa + rank + log(household) + sex + sex:household, data = train_data, family = binomial(logit))
anova(model_2,model_3,test="Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ gre + gpa + rank + log(household) + sex
## Model 2: y ~ gre + gpa + rank + log(household) + sex + sex:household
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 9992 11198
## 2 9990 11193 2 4.855 0.08826 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1logistic.display(model_3)##
## OR lower95ci upper95ci Pr(>|Z|)
## gre 0.9966875 0.9962787 0.9970965 1.386319e-56
## gpa 0.3924575 0.3458980 0.4452842 9.525349e-48
## rank2 1.8152299 1.5595131 2.1128771 1.403521e-14
## rank3 3.5008592 3.0119651 4.0691094 6.333148e-60
## rank4 4.5123001 3.7978574 5.3611419 8.358799e-66
## log(household) 1.1995613 0.7955881 1.8086587 3.851270e-01
## sex1 5.1560935 4.0557353 6.5549888 6.741782e-41
## sex0:household 1.0016926 0.9992344 1.0041569 1.773383e-01
## sex1:household 1.0003316 0.9977720 1.0028978 7.997646e-01model_4 <- glm(y ~ gre + gpa + rank + log(household) + sex + sex:gre, data = train_data, family = binomial(logit))
anova(model_2,model_4,test="Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ gre + gpa + rank + log(household) + sex
## Model 2: y ~ gre + gpa + rank + log(household) + sex + sex:gre
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 9992 11198
## 2 9991 11161 1 36.604 1.447e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1logistic.display(model_4)##
## OR lower95ci upper95ci Pr(>|Z|)
## gre 0.9981501 0.9975290 0.9987716 5.544111e-09
## gpa 0.3902814 0.3438674 0.4429602 4.688223e-48
## rank2 1.8359946 1.5756004 2.1394231 6.923789e-15
## rank3 3.5440510 3.0459609 4.1235912 2.940311e-60
## rank4 4.5979175 3.8665335 5.4676484 9.489477e-67
## log(household) 1.3871463 1.2393262 1.5525976 1.254523e-08
## sex1 18.0498279 11.0472491 29.4911689 7.396789e-31
## gre:sex1 0.9974420 0.9966131 0.9982717 1.561698e-09sex:gre の交互作用項が有意であった。よって性別によってgreの不合格に対する効果は異なる。
予測モデルの構築
データから合格の可否を予測するモデルを作成する。
データの準備
説明のため、データをサンプリングする。
train_data$y<- as.factor(train_data$y)
train_data$rank<- as.factor(train_data$rank)
train_data$sex <- as.factor(train_data$sex)
train_data_num <- sample(1:nrow(train_data),size = nrow(train_data)*0.1)
train_data <- train_data[train_data_num,]
summary(train_data)## gre gpa rank gre_math household
## Min. :245.3 Min. :2.214 1:173 Min. : 77.87 Min. : 39.25
## 1st Qu.:514.9 1st Qu.:3.150 2:344 1st Qu.:153.60 1st Qu.:111.59
## Median :594.8 Median :3.389 3:334 Median :178.15 Median :147.39
## Mean :592.2 Mean :3.386 4:149 Mean :177.80 Mean :159.00
## 3rd Qu.:667.7 3rd Qu.:3.640 3rd Qu.:200.25 3rd Qu.:192.52
## Max. :800.0 Max. :4.000 Max. :266.73 Max. :604.54
## sex y
## 0:495 0:656
## 1:505 1:344
##
##
##
## head(train_data)## gre gpa rank gre_math household sex y
## 5029 613.6519 3.456330 2 187.1800 216.3354 1 1
## 1916 479.7427 3.365282 3 130.9771 218.3751 0 0
## 1851 800.0000 2.992912 3 228.4181 253.2760 1 0
## 6043 730.8965 3.510851 2 212.1587 129.8386 1 1
## 8916 596.9922 3.420729 1 190.3718 159.8472 1 0
## 4877 588.9583 3.197925 1 182.6535 174.6296 1 0model 構築と予測力の評価
モデルの構築は上記と同じ。 予測力の評価はROC曲線のAUCを用いる。
#モデル作成
model_4 <- glm(y ~ gre+ gre_math + gpa + rank + log(household) + sex + sex:gre + sex:log(household) + rank:gre +rank:log(household), data = train_data, family = binomial(logit))
#モデルをデータに当てはめ、予測値（合格割合）を得る。
pred <- predict(model_4,type = "response",newdata = train_data)
hist(pred)
#ROCを作成し格納する。
roc_1 <- pROC::roc(response = as.numeric(train_data$y)-1,predictor = pred)
#AUCとAUCの信頼区間、他の予測指標の算出
pROC::auc(roc_1)## Area under the curve: 0.7838pROC::ci(roc_1)## 95% CI: 0.7538-0.8138 (DeLong)pROC::coords(roc_1,"best", ret=c("threshold", "specificity", "sensitivity","accuracy","tn", "tp", "fn", "fp", "npv", "ppv", "1-specificity","1-sensitivity", "1-accuracy", "1-npv", "1-ppv"))## threshold specificity sensitivity accuracy tn
## 0.3503357 0.7347561 0.7296512 0.7330000 482.0000000
## tp fn fp npv ppv
## 251.0000000 93.0000000 174.0000000 0.8382609 0.5905882
## 1-specificity 1-sensitivity 1-accuracy 1-npv 1-ppv
## 0.2652439 0.2703488 0.2670000 0.1617391 0.4094118model に多項式を入れる。
ploy(x,degree=n)はx のn次多項式を生成する。つまり、poly(gre,degree = 15)で \[ x+x^2+x^3+\cdots + x^{15} \] を表現している。
model_5 <- glm(y ~ poly(gre,degree = 15)+ gre_math + gpa + rank + poly(log(household),15) + sex + sex:gre + sex:log(household) + rank:gre +rank:log(household) , data = train_data, family = binomial(logit))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredpred <- predict(model_5,type = "response",newdata = train_data)## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleadinghist(pred)
roc_2 <- pROC::roc(response = as.numeric(train_data$y)-1,predictor = pred)
pROC::auc(roc_2)## Area under the curve: 0.7985pROC::ci(roc_2)## 95% CI: 0.7697-0.8274 (DeLong)pROC::coords(roc_2,"best", ret=c("threshold", "specificity", "sensitivity", "accuracy","tn", "tp", "fn", "fp", "npv", "ppv", "1-specificity",
"1-sensitivity", "1-accuracy", "1-npv", "1-ppv"))## threshold specificity sensitivity accuracy tn
## 0.3784632 0.7728659 0.6831395 0.7420000 507.0000000
## tp fn fp npv ppv
## 235.0000000 109.0000000 149.0000000 0.8230519 0.6119792
## 1-specificity 1-sensitivity 1-accuracy 1-npv 1-ppv
## 0.2271341 0.3168605 0.2580000 0.1769481 0.3880208plot(roc_1,col="blue")##
## Call:
## roc.default(response = as.numeric(train_data$y) - 1, predictor = pred)
##
## Data: pred in 656 controls (as.numeric(train_data$y) - 1 0) < 344 cases (as.numeric(train_data$y) - 1 1).
## Area under the curve: 0.7838plot(roc_2,col="red",add=T)
##
## Call:
## roc.default(response = as.numeric(train_data$y) - 1, predictor = pred)
##
## Data: pred in 656 controls (as.numeric(train_data$y) - 1 0) < 344 cases (as.numeric(train_data$y) - 1 1).
## Area under the curve: 0.7985モデルの予測力を上げるために多項式を入れる手法はよく用いられる。 しかし、これには過剰適合のリスクがついて回る。
過剰適合の検証
モデル作成とモデル検証にデータを分ける。
model_data_num <- sample(1:nrow(train_data),size = nrow(train_data)*0.6)
model_data <- train_data[model_data_num,]
test_data <- train_data[-model_data_num,]多項式モデルを構築する。
model_6 <- glm(y ~ poly(gre,degree = 15)+ gre_math + gpa + rank + poly(log(household),15) + sex + sex:gre + sex:log(household) + rank:gre +rank:log(household), data = model_data, family = binomial(logit))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredpred <- predict(model_6,type = "response",newdata = test_data)## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleadinghist(pred)
roc_3 <- pROC::roc(response = as.numeric(test_data$y)-1,predictor = pred)
pROC::auc(roc_3)## Area under the curve: 0.7623pROC::ci(roc_3)## 95% CI: 0.7133-0.8114 (DeLong)pROC::coords(roc_3,"best", ret=c("threshold", "specificity", "sensitivity", "accuracy","tn", "tp", "fn", "fp", "npv", "ppv", "1-specificity",
"1-sensitivity", "1-accuracy", "1-npv", "1-ppv"))## threshold specificity sensitivity accuracy tn
## 0.3056309 0.6754717 0.7555556 0.7025000 179.0000000
## tp fn fp npv ppv
## 102.0000000 33.0000000 86.0000000 0.8443396 0.5425532
## 1-specificity 1-sensitivity 1-accuracy 1-npv 1-ppv
## 0.3245283 0.2444444 0.2975000 0.1556604 0.4574468ROCを比較する。
plot(roc_1,col="blue",lty=3)##
## Call:
## roc.default(response = as.numeric(train_data$y) - 1, predictor = pred)
##
## Data: pred in 656 controls (as.numeric(train_data$y) - 1 0) < 344 cases (as.numeric(train_data$y) - 1 1).
## Area under the curve: 0.7838plot(roc_2,col="blue", add=T)##
## Call:
## roc.default(response = as.numeric(train_data$y) - 1, predictor = pred)
##
## Data: pred in 656 controls (as.numeric(train_data$y) - 1 0) < 344 cases (as.numeric(train_data$y) - 1 1).
## Area under the curve: 0.7985plot(roc_3,col="red", add=T)
##
## Call:
## roc.default(response = as.numeric(test_data$y) - 1, predictor = pred)
##
## Data: pred in 265 controls (as.numeric(test_data$y) - 1 0) < 135 cases (as.numeric(test_data$y) - 1 1).
## Area under the curve: 0.7623多項式を入れたモデルは他のデータに対する汎用性がない。 これは、内的妥当性に問題があることを示しているし、外部データを用いて外的妥当性を検証しても同様の問題が生じる。