保護犬・保護猫・保護動物のデータ解析その5ーどの要素が譲渡に重要か？ー
前回（保護犬・保護猫・保護動物のデータ解析その5）の続きです。 今回は、logistic regression、random forest、gradient boost machineを使って各要素の重要度を算出したいと思います。
logistic regressionは、理解しやすさが高いので採用します。 ramdom forestとgradient boost machineは、高精度かつ各要素の重要度を算出できるので採用です。 重要度を算出するためにvarImpという関数を使いますが、これの説明見ると他のモデルでも要素の重要度を計測できるようですね。
まずは、ライブラリとデータを読み込みます。
library(tidyverse)df <- read.csv("C:\\Users\\shelter\\Austin_Animal_Center_Outcomes.csv")
#譲渡を１，そのほかを0とする。
df$Outcome.Type <- ifelse(df$Outcome.Type=="Adoption","1","0")
#時間の処理
library(stringr)
df_year <- df %>% filter(str_detect(df$Age.upon.Outcome,"year"))
df_year$Age.upon.Outcome <- str_remove(df_year$Age.upon.Outcome, "years")
df_year$Age.upon.Outcome <- str_remove(df_year$Age.upon.Outcome, "year")
df_year$Age.upon.Outcome <- str_trim(df_year$Age.upon.Outcome, side = "both")
df_year <- df_year %>% filter(Age.upon.Outcome > 0)
df_week <- df %>% filter(str_detect(df$Age.upon.Outcome,"week"))
df_week$Age.upon.Outcome <- 0.5
df_day <- df %>% filter(str_detect(df$Age.upon.Outcome,"day"))
df_day$Age.upon.Outcome <- 0.5
df_age <- rbind(df_year,df_week,df_day)
#必要な変数だけ取り出す。
df2 <- df_age %>% select(Outcome.Type, #譲渡か否か
Animal.Type, #動物種（犬猫その他）
Sex.upon.Outcome, #性別・避妊去勢
Age.upon.Outcome, #顛末時の年齢
Breed) #ブリード
#犬と猫に分けます。
df_dog <- df2 %>% filter(Animal.Type=="Dog") %>% select(-Animal.Type)
df_cat <- df2 %>% filter(Animal.Type=="Cat") %>% select(-Animal.Type)多い10犬種に絞って解析します。
Breeddf_dog <- as.data.frame(table(df_dog$Breed))
Breeddf_dog <- Breeddf_dog %>% arrange(-Freq)
top10_dog <- Breeddf_dog[1:10,]
dog <- df_dog %>%
filter(Breed %in% top10_dog[,1])
#解析に時間がかかるため、データ数を30952から10000に減らします。
dog_10000 <- dog[1:10000,]
colnames(dog_10000) <- c("adopted","sex","years_old","breed")これからいよいよ解析です。まずCPUの並列計算を有効にして、計算を早くします。ふつうは、CPUを１つしか使いませんが、これを使うとCPUを複数使えて効率が良くなります。私のCPUはコアが４つあるので、makePSOCKclusterを4にしました。
library(doParallel)cl <- makePSOCKcluster(4)
registerDoParallel(cl)logistic regression
logistci regressionをする際には、outcome variable は、factorにしないといけないので、as.factorを使います。glmの機能でlogistic regressionをを実施するため、family=binomialにします。
dog_10000_glm <- dog_10000
dog_10000_glm$adopted <- as.factor(dog_10000_glm$adopted)
str(dog_10000_glm)## 'data.frame': 10000 obs. of 4 variables:
## $ adopted : Factor w/ 2 levels "0","1": 2 2 2 2 1 1 2 1 2 2 ...
## $ sex : chr "Neutered Male" "Spayed Female" "Neutered Male" "Neutered Male" ...
## $ years_old: chr "1" "7" "9" "10" ...
## $ breed : chr "Chihuahua Shorthair Mix" "Chihuahua Shorthair Mix" "Chihuahua Shorthair Mix" "Chihuahua Shorthair" ...factorがうまくできているか確認するために、元のデータを見てみます。
factorが、
元々1 → 2
元々2 → 1
となっていて気持ちがわるいですが、仕方ありません。
str(dog_10000)## 'data.frame': 10000 obs. of 4 variables:
## $ adopted : chr "1" "1" "1" "1" ...
## $ sex : chr "Neutered Male" "Spayed Female" "Neutered Male" "Neutered Male" ...
## $ years_old: chr "1" "7" "9" "10" ...
## $ breed : chr "Chihuahua Shorthair Mix" "Chihuahua Shorthair Mix" "Chihuahua Shorthair Mix" "Chihuahua Shorthair" ...lg <- glm(adopted~.,
data = dog_10000_glm,
family = binomial)
summary(lg)##
## Call:
## glm(formula = adopted ~ ., family = binomial, data = dog_10000_glm)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6754 -1.0384 -0.2869 1.0223 3.1983
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.988496 0.180524 -11.015 < 2e-16 ***
## sexIntact Male -0.228728 0.208955 -1.095 0.273679
## sexNeutered Male 2.893016 0.146005 19.815 < 2e-16 ***
## sexSpayed Female 3.109894 0.146896 21.171 < 2e-16 ***
## sexUnknown -10.597071 128.892699 -0.082 0.934475
## years_old10 -1.177969 0.157211 -7.493 6.74e-14 ***
## years_old11 -1.603682 0.267321 -5.999 1.98e-09 ***
## years_old12 -1.581002 0.237628 -6.653 2.87e-11 ***
## years_old13 -1.399969 0.290736 -4.815 1.47e-06 ***
## years_old14 -1.505191 0.422715 -3.561 0.000370 ***
## years_old15 -1.621673 0.425063 -3.815 0.000136 ***
## years_old16 -2.578575 0.760982 -3.388 0.000703 ***
## years_old17 -2.813136 1.060666 -2.652 0.007996 **
## years_old18 -1.431239 1.227776 -1.166 0.243729
## years_old19 -14.157790 535.411172 -0.026 0.978904
## years_old2 -0.430819 0.061358 -7.021 2.20e-12 ***
## years_old24 -14.301775 535.411176 -0.027 0.978690
## years_old3 -0.809599 0.077003 -10.514 < 2e-16 ***
## years_old4 -0.976409 0.094819 -10.298 < 2e-16 ***
## years_old5 -0.977289 0.096212 -10.158 < 2e-16 ***
## years_old6 -1.205447 0.110551 -10.904 < 2e-16 ***
## years_old7 -1.262485 0.127869 -9.873 < 2e-16 ***
## years_old8 -0.829825 0.128573 -6.454 1.09e-10 ***
## years_old9 -1.237510 0.168833 -7.330 2.30e-13 ***
## breedBoxer Mix -0.500577 0.165874 -3.018 0.002546 **
## breedChihuahua Shorthair -0.003021 0.170801 -0.018 0.985890
## breedChihuahua Shorthair Mix -0.312797 0.120086 -2.605 0.009193 **
## breedDachshund Mix -0.370085 0.161104 -2.297 0.021608 *
## breedGerman Shepherd Mix -0.168811 0.133930 -1.260 0.207510
## breedLabrador Retriever Mix -0.425458 0.121774 -3.494 0.000476 ***
## breedMiniature Poodle Mix -0.571329 0.161441 -3.539 0.000402 ***
## breedPit Bull -0.427695 0.159248 -2.686 0.007237 **
## breedPit Bull Mix -0.528144 0.117433 -4.497 6.88e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13560 on 9999 degrees of freedom
## Residual deviance: 11365 on 9967 degrees of freedom
## AIC: 11431
##
## Number of Fisher Scoring iterations: 12各係数と各変数の合計が、logoddsです。
log odds = ax1 + ax2 + ax3 …
factorで大きい方となる（今回は１-譲渡-になる）確率を計算するには、logoddsを変換します。
\[ 確率 = \frac{1}{1+e^-logodds}\ \] これはsigmoid functionで、０から１の間の確率を生み出します。以下を見てみましょう。
logodds <- seq(-10,10,length=1000)
probability <- 1/(1+exp(-logodds))
as.data.frame(cbind(logodds, probability)) %>%
ggplot(aes(x=logodds,y=probability))+
geom_line()
logoddsが大きくなれば１（譲渡）となる確率も上がるということです。
係数が大きいものを見てみます。
coefficient <- as.data.frame(coef(lg))
coefficient <- coefficient %>% arrange(-`coef(lg)`)
head(coefficient)## coef(lg)
## sexSpayed Female 3.109894093
## sexNeutered Male 2.893016073
## breedChihuahua Shorthair -0.003020574
## breedGerman Shepherd Mix -0.168811494
## sexIntact Male -0.228727903
## breedChihuahua Shorthair Mix -0.3127967501位Spayed female 3.1
2位Neutered male 2.9
3位チワワショートヘア -0.003
です。1位と2位の要素は、譲渡に成功するファクターとして重要だと分かりました。
同じ変数内では、dummy codingによって1つの基準カテゴリーが決められています。そのため、１つの基準カテゴリーのcoeffientが表示されていないことに注意です。
今度は、係数が低い方から見てみます。
coefficient <- coefficient %>% arrange(`coef(lg)`)
head(coefficient)## coef(lg)
## years_old24 -14.301775
## years_old19 -14.157790
## sexUnknown -10.597071
## years_old17 -2.813136
## years_old16 -2.578575
## (Intercept) -1.988496年齢が高齢というのは、逆にマイナスファクターだということが分かりました。
Random Forest
random forestはcaretというパッケージで実施できます。このパッケージは色々な寄せ集めで、gradient boost machineもできます。
caretというパッケージを使う際に、色々きかれます。random forestをする際には、e1071をインストールするように促されたので追加でインストールしました。gbmの際には、パッケージがパッケージが必要だけどインストールしますか、と聞かれたのでyesにしたら自動でインストールしてくれました。
今回は予測するわけではないので、training, validation, test setにデータを分けずに全部モデルに読み込ませます。10000データセットだと、時間がかかりすぎてしまうので、サイズを減らします。
library(caret)dog_1000 <- dog_10000[1:1000,]
train.rf <- train(adopted~., dog_1000, method="rf", importance=TRUE)重要度を計算します。
varImp.rf <- varImp(train.rf)
plot(varImp.rf, top=20)
避妊去勢が重要だと分かりました。年齢も、重要だと分かりましたが、logistic regressionの結果から、譲渡にはマイナスの方向で働くと思われます。
Gradient Boost Machine
random forestに比べて圧倒的に早い！
train.gbm <- train(adopted~., dog_10000, method="gbm")## Iter TrainDeviance ValidDeviance StepSize Improve
## 1 1.3287 nan 0.1000 0.0132
## 2 1.3065 nan 0.1000 0.0109
## 3 1.2886 nan 0.1000 0.0090
## 4 1.2733 nan 0.1000 0.0075
## 5 1.2601 nan 0.1000 0.0066
## 6 1.2489 nan 0.1000 0.0056
## 7 1.2429 nan 0.1000 0.0031
## 8 1.2339 nan 0.1000 0.0042
## 9 1.2292 nan 0.1000 0.0021
## 10 1.2251 nan 0.1000 0.0018
## 20 1.1990 nan 0.1000 0.0002
## 40 1.1704 nan 0.1000 0.0000
## 60 1.1586 nan 0.1000 0.0001
## 80 1.1508 nan 0.1000 0.0000
## 100 1.1460 nan 0.1000 -0.0000
## 120 1.1432 nan 0.1000 -0.0002
## 140 1.1403 nan 0.1000 -0.0001
## 150 1.1391 nan 0.1000 -0.0001重要度を計測します。
varImp.gbm <- varImp(train.gbm)varImp.gbm <- varImp(train.gbm)
のところで、以下のエラーがでてしまう。
Error in relative.influence(object, n.trees = numTrees) : could not find function “relative.influence”
library(gbm)をその前に実行することで解決。
library(gbm)varImp.gbm <- varImp(train.gbm)plot(varImp.rf, top=20)
randomforest とほぼ同じです。避妊メスが一番重要のようです。
ロジスティック回帰分析でもvarImpが使えるので、重要度を計算したいと思います。重要度は、t値の絶対値とされています。t値は、coefficient/Std. Errorで計算されます。ｔ値は、summary(glmモデル)で出てくる、ｚ valueと同じ値のようです。なぜこれが重要度かというと、そのcoefficientが大きく揺らぐかどうかの指標になるからです。ゆらぎ（standard error）が係数と同じくらいの数値だと、z valueは１に近づきます。そうすると、P valueが大きくなります。p valueは、係数が０ではないことの確率です。z値が大きければどんなケースでも、その係数が効いているという指標になるわけです。
summary(lg)##
## Call:
## glm(formula = adopted ~ ., family = binomial, data = dog_10000_glm)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6754 -1.0384 -0.2869 1.0223 3.1983
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.988496 0.180524 -11.015 < 2e-16 ***
## sexIntact Male -0.228728 0.208955 -1.095 0.273679
## sexNeutered Male 2.893016 0.146005 19.815 < 2e-16 ***
## sexSpayed Female 3.109894 0.146896 21.171 < 2e-16 ***
## sexUnknown -10.597071 128.892699 -0.082 0.934475
## years_old10 -1.177969 0.157211 -7.493 6.74e-14 ***
## years_old11 -1.603682 0.267321 -5.999 1.98e-09 ***
## years_old12 -1.581002 0.237628 -6.653 2.87e-11 ***
## years_old13 -1.399969 0.290736 -4.815 1.47e-06 ***
## years_old14 -1.505191 0.422715 -3.561 0.000370 ***
## years_old15 -1.621673 0.425063 -3.815 0.000136 ***
## years_old16 -2.578575 0.760982 -3.388 0.000703 ***
## years_old17 -2.813136 1.060666 -2.652 0.007996 **
## years_old18 -1.431239 1.227776 -1.166 0.243729
## years_old19 -14.157790 535.411172 -0.026 0.978904
## years_old2 -0.430819 0.061358 -7.021 2.20e-12 ***
## years_old24 -14.301775 535.411176 -0.027 0.978690
## years_old3 -0.809599 0.077003 -10.514 < 2e-16 ***
## years_old4 -0.976409 0.094819 -10.298 < 2e-16 ***
## years_old5 -0.977289 0.096212 -10.158 < 2e-16 ***
## years_old6 -1.205447 0.110551 -10.904 < 2e-16 ***
## years_old7 -1.262485 0.127869 -9.873 < 2e-16 ***
## years_old8 -0.829825 0.128573 -6.454 1.09e-10 ***
## years_old9 -1.237510 0.168833 -7.330 2.30e-13 ***
## breedBoxer Mix -0.500577 0.165874 -3.018 0.002546 **
## breedChihuahua Shorthair -0.003021 0.170801 -0.018 0.985890
## breedChihuahua Shorthair Mix -0.312797 0.120086 -2.605 0.009193 **
## breedDachshund Mix -0.370085 0.161104 -2.297 0.021608 *
## breedGerman Shepherd Mix -0.168811 0.133930 -1.260 0.207510
## breedLabrador Retriever Mix -0.425458 0.121774 -3.494 0.000476 ***
## breedMiniature Poodle Mix -0.571329 0.161441 -3.539 0.000402 ***
## breedPit Bull -0.427695 0.159248 -2.686 0.007237 **
## breedPit Bull Mix -0.528144 0.117433 -4.497 6.88e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13560 on 9999 degrees of freedom
## Residual deviance: 11365 on 9967 degrees of freedom
## AIC: 11431
##
## Number of Fisher Scoring iterations: 12同じようにプロットします。
varImp.lg <- varImp(lg)
plot(varImp.lg)
うまくプロットできないです。stack overflowでしらべたところggplot2を使えばプロットできるようです。
ggplot(data= varImp.lg, aes(x=rownames(varImp.lg),y=Overall)) +
geom_bar(position="dodge",stat="identity",width = 0, color = "black") +
coord_flip() + geom_point(color='skyblue') + xlab(" Importance Score")+
ggtitle("Variable Importance") +
theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.background = element_rect(fill = 'white', colour = 'black'))
避妊去勢が譲渡成功に非常に有効だと分かりました。
猫の解析
多い10猫種に絞って解析します。
Breeddf_cat <- as.data.frame(table(df_cat$Breed))
Breeddf_cat <- Breeddf_cat %>% arrange(-Freq)
top10_cat <- Breeddf_cat[1:10,]
cat <- df_cat %>%
filter(Breed %in% top10_cat[,1])
#解析に時間がかかるため、データ数を30952から10000に減らします。
cat_10000 <- cat[1:10000,]
colnames(cat_10000) <- c("adopted","sex","years_old","breed")logistic regression
cat_10000_glm <- cat_10000
cat_10000_glm$adopted <- as.factor(cat_10000_glm$adopted)
lg_cat <- glm(adopted~.,
data = cat_10000_glm,
family = binomial)係数が大きい方から見ます。
coefficient_cat <- as.data.frame(coef(lg_cat))
coefficient_cat <- coefficient_cat %>% arrange(-`coef(lg_cat)`)
head(coefficient_cat)## coef(lg_cat)
## sexSpayed Female 4.5255238
## sexNeutered Male 4.1514348
## breedSnowshoe Mix 0.5764071
## breedSiamese 0.4003164
## breedDomestic Longhair 0.3720200
## breedDomestic Medium Hair 0.26954111位Spayed female 4.5
2位Neutered male 4.1
3位Snow shoe mix 0.5
今度は、係数が低い方から見てみます。
coefficient_cat <- coefficient_cat %>% arrange(`coef(lg_cat)`)
head(coefficient_cat)## coef(lg_cat)
## years_old22 -17.560209
## years_old20 -16.661112
## sexUnknown -13.166468
## (Intercept) -4.179621
## years_old19 -1.802686
## years_old17 -1.321207年齢が高齢がマイナスファクター。
varImp <- varImp(lg_cat)
ggplot(data= varImp, aes(x=rownames(varImp),y=Overall)) +
geom_bar(position="dodge",stat="identity",width = 0, color = "black") +
coord_flip() + geom_point(color='skyblue') + xlab(" Importance Score")+
ggtitle("Variable Importance") +
theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.background = element_rect(fill = 'white', colour = 'black'))
Random Forest
cat_1000 <- cat_10000[1:1000,]
train.rf <- train(adopted~., cat_1000, method="rf", importance=TRUE)重要度を計算します。
varImp.rf <- varImp(train.rf)
plot(varImp.rf, top=20)
避妊去勢が重要だと分かりました。
Gradient Boost Machine
train.gbm <- train(adopted~., cat_10000, method="gbm")## Iter TrainDeviance ValidDeviance StepSize Improve
## 1 1.2358 nan 0.1000 0.0202
## 2 1.2031 nan 0.1000 0.0170
## 3 1.1737 nan 0.1000 0.0142
## 4 1.1493 nan 0.1000 0.0120
## 5 1.1288 nan 0.1000 0.0100
## 6 1.1109 nan 0.1000 0.0085
## 7 1.0954 nan 0.1000 0.0075
## 8 1.0825 nan 0.1000 0.0065
## 9 1.0709 nan 0.1000 0.0056
## 10 1.0609 nan 0.1000 0.0048
## 20 1.0103 nan 0.1000 0.0013
## 40 0.9956 nan 0.1000 0.0005
## 60 0.9912 nan 0.1000 -0.0002
## 80 0.9884 nan 0.1000 -0.0002
## 100 0.9868 nan 0.1000 -0.0003varImp.gbm <- varImp(train.gbm)plot(varImp.rf, top=20)
避妊去勢が大事ということが分かりました。
すでに日本の動物保護センターにおいても避妊去勢を実施する活動がなされていますね。