読者です 読者をやめる 読者になる 読者になる

疫学と医療統計学と遺伝学と時々、大学院生活

疫学を専門とする大学院生の研究に関する備忘録的ページ。

R リッジ回帰モデルと異常検知

今回は、『入門機械学習による異常検知-Rによる実践ガイド-』(コロナ社、井出剛著、2015)の「6.3 リッジ回帰と異常検知(P165-172)」の内容をもとに政府・官公庁データをもとにして作成した都道府県別の10万人あたりの自殺者数を予測するモデルを作成し、実測値を比較した。さらに異常度を計算し、それを可視化した。

まずは、データの概要

#pref: 都道府県名
#suicide: 人口10万人あたりの自殺者数(平成27年人口動態調査)
#poverty: 子供の貧困率(平成24年就業構造基本調査をもとに戸室氏が作成)
#irregular: 若年者の非正規(平成24年就業構造基本調査)
#income: 県民所得(平成25年度県民経済)
#GDP: 県内総生産(平成25年度県民経済)
#E: 都道府県別1000児童あたりの不登校児童数・小学校(平成27年度「児童生徒の問題行動等生徒指導上の諸問題に関する調査」)
#J: 都道府県別1000児童あたりの不登校児童数・中学校(平成27年度「児童生徒の問題行動等生徒指導上の諸問題に関する調査」)
#H: 都道府県別1000児童あたりの不登校児童数・高校(平成27年度「児童生徒の問題行動等生徒指導上の諸問題に関する調査」)
#U: 完全失業率(平成27年度労働力調査)
#S_2024: 20〜24歳被生活保護者実数(平成27年度被保護者調査)
#S_2529: 25〜29歳被生活保護者実数(平成27年度被保護者調査)
#counselor: 都道府県別教育委員会が設置する「教育支援センター」の常勤職員数(平成27年度「児童生徒の問題行動等生徒指導上の諸問題に関する調査」)
#bullying: 都道府県別1000人あたりのいじめ認知件数(平成27年度「児童生徒の問題行動等生徒指導上の諸問題に関する調査」)

1. リッジ回帰による予測値の算出と予測値および実測値の図示

suicide <- read.csv("suicide.csv", header=T, sep=",")
head(suicide)
##       pref suicide poverty irregular   income      GDP   E    J    H   U
## 1 hokkaido    19.5    19.7      42.8 18757883 18268793 3.8 27.1  7.1 3.4
## 2   aomori    20.5    17.6      37.9  4493046  4411514 3.5 27.1  7.3 4.5
## 3    iwate    23.3    13.9      37.6  4637299  4516178 3.0 23.6 14.3 2.9
## 4   miyagi    17.4    15.3      39.3  8987335  8816646 4.7 35.3 21.9 3.7
## 5    akita    25.7     9.9      35.3  3564708  3477343 2.5 20.4 12.2 3.6
## 6 yamagata    21.7    12.0      35.8  4039485  3830374 3.0 22.2 15.5 2.7
##   S_2024 S_2529 counselor bullying
## 1    811   1094        52     11.3
## 2    172    202        10      8.8
## 3     92    104         7     24.5
## 4    108    150         6     70.8
## 5     66     76        13     17.8
## 6     70     81        15     48.4

X <- suicide[,-c(1,2)]; M <- ncol(X) #不要な変数の除去
y <- suicide[,2]; N <- length(y) #対象とするアウトカムの指定
lambdas <- seq(0, 5, length=50) #ラムダの候補生成
model <- lm.ridge(y ~ ., cbind(X,y), lambda=lambdas) #リッジ回帰の実行
bestIdx <- which.min(model$GCV) #一般化交差検証の評価値が最小となるものを指定
coefs <- coef(model)[bestIdx,] #回帰係数も同様
lam <- model$lambda[bestIdx] #一般化交差検証の評価値が最小になるラムダを選択 
ypred <- as.matrix(X) %*% as.matrix(coefs[2:13]) + coefs[1] #予測値の算出
#図示
par(pin=c(3,3))
plot(y, ypred, xlab="actual", ylab="predicted", xlim=c(15,25), ylim=c(15,25))

f:id:ryosukefujii0320:20170428173926p:plain

2. 推定値をもとに異常度の算出とそれを図示する(詳細は参考図書P169〜172のあたり参照)

sig2 <- (lam*sum(coefs[2:13]^2) + sum(as.numeric(ypred)-y)^2) / N
X_ <- t(scale(X, scale=F))
H <- t(X_) %*% solve(X_ %*% t(X_) + lam*diag(M), X_)
TrHN <- sum(diag(H)) / N
a <- (as.numeric(ypred) - y)^2 / ((1-TrHN)^2*sig2)
plot(a, xlab="index", ylab="anomaly score")
th <- sort(a)[N*(1-0.05)]
lines(0:50, rep(th, length(0:50)), col="red", lty=2)

f:id:ryosukefujii0320:20170428173959p:plain

3. 異常度が高いと検知した都道府県について、データを確認する

select <- which(a > th) #今回異常度の高いと検出された都道府県を抽出する
data.frame(PREF=suicide$pref[select], ACT=y[select], PRED=ypred[select]) #実測値と予測値を比較する
##      PREF  ACT     PRED
## 1   akita 25.7 21.61400
## 2   fukui 15.4 19.38527
## 3 shimane 22.9 18.60211

4. 今回の結果

抽出されたのは、秋田県福井県島根県秋田県島根県に関しては、ACTがPREDを上回っており、予測されるよりも自殺者数が多いことがこの異常度の高さに寄与していると考えられる。一方で福井県は、PREDがACTを上回っており、なんらかの活動の成果とも考えることができる。
ちょっとモデルの結果をみると、実測値が高い都道府県では、予測値が実測値よりは低い傾向(右肩上がりになっていない)にあるので、予測モデルの改善も必要かと考えられる。

参考図書

入門機械学習による異常検知-Rによる実践ガイド-』(コロナ社、井出剛著、2015)

RF
20170428