今回は、『入門機械学習による異常検知-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))
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)
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