Rによるデータ分析

Rによる異常の定量化の分析

1つの変数だけで決まっている異常は、その変数をグラフにすれば簡単に分析できます。 「5以上は異常」といったデータの見方ができます。

ところが、2つ以上の変数で決まっている以上は、データを見ているだけでは、 「5以上は異常」といったシンプルなデータの見方ができません。 このページの方法は、多変量でも異常度をシンプルに表現する方法になります。

このページのコードの前提

このページのコードは、入力データに「Y」という名前の変数が入っていることを想定しています。

「Y」の変数は、「0」と「1」の2つの数値が入っていることを想定しています。 「0」のサンプルは、単位空間であることを想定しています。 「0」のサンプルだけでモデルが作られます。

「1」のサンプルは信号空間です。 単位空間のデータに対して、異常なのかを見たいデータになります。

Y以外の変数は、名前に決まりはありません。

計算結果はグラフになるようにしましたが、定量的な結果として、誤判別になるサンプル数がわかるようにもしてみました。 最後に出てくる数値がそれです。

MT法をベースにした方法

MT法

基本的な MT法 の場合です。 この場合は、量的変数しか使えません。 質的変数が含まれる場合は、下に差し替えの方法を書きました。

library(ggplot2) # ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- Data
# 質的変数がある場合にこの行を差し替え
Data2 <- subset(Data1, Data1$Y == 0)
# Yが0のサンプルのデータを別に作る
Data2$Y <- NULL
# Yの列を消す
n <- nrow(Data2)
# Yが0のサンプルのサンプル数を計算
Ave1 <- colMeans(Data2)
# Yが0のサンプルの各変数の平均値を計算
Var1 <- var(Data2)*(n-1)/n
# Yが0のサンプルの共分散行列を計算
k <- ncol(Data2)
# 変数の数を計算
Data1$Y <- NULL
# Yの列を消す
MD <- mahalanobis(Data1, Ave1, Var1)/k
# MD(の2乗)を計算
Data3 <- cbind(Data, MD)
# 元のデータとマハラノビス距離の計算結果を合わせる
Data3$Y <- factor(Data3$Y)
# Yの列を文字型にする
ggplot(Data3, aes(x=Y, y=MD)) + geom_jitter(size=1, position=position_jitter(0.1))
# 一次元ジター散布図を描く
MaxMDinUnit <- max(subset(Data3, Data3$Y == 0)$MD)
# Yが0のサンプルのMDのMaxを求める
Data6 <- subset(Data3, Data3$Y == 1)
# Yが1のサンプルを抜き出す
nrow(subset(Data6, Data6$MD < MaxMDinUnit))
# Yが1で、Yが0のサンプルのMDのMaxよりも小さなサンプル数を求める
MT法

質的変数が含まれる場合

MT法のコードは、量的変数のみにしていますが、質的変数が含まれる場合は、
Data1 <- Data# 質的変数がある場合にこの行を差し替え」
となっている1行を、下記の2行に書き変えます。

library(fastDummies) # ライブラリを読み込み
Data1 <- dummy_cols(Data,remove_first_dummy = TRUE,remove_selected_columns = TRUE)
# 質的変数はダミー変換

主成分MT法

主成分MT法 の場合です。 質的変数が混ざっていてもできるようにしてあります。

基本のMT法では、多重共線性で計算ができない場合に役立ちます。 「変数の数 > サンプルの数」でも役立つ場合もあります。

library(dummies) # ライブラリを読み込み
library(ggplot2)
# ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- Data
# 加工用のデータを作る
Y <- Data1$Y
# Yのデータを別に作る
Data1 <- dummy.data.frame(Data)
# 質的変数はダミー変換
Data1$Y <- NULL
# Yの列を消す
Data2 <- subset(Data1, Data$Y == 0)
# Yが0のサンプルのデータを別に作る
Data2$Y <- NULL
# Yの列を消す
model <- prcomp(Data2, scale=TRUE, tol=0.01)
# Yが0のサンプルのデータで主成分分析のモデルを作る。 寄与率の小さな主成分は除く
Data3 <- predict(model, Data1)
# 主成分得点を取る
Data4 <- subset(Data3, Data$Y == 0)
# 単位空間のデータを別に作る
n <- nrow(Data4)
# Yが0のサンプルのサンプル数を計算
Ave1 <- colMeans(Data4)
# Yが0のサンプルの各変数の平均値を計算
Var1 <- var(Data4)*(n-1)/n
# Yが0のサンプルの共分散行列を計算
k <- ncol(Data4)
# 変数の数を計算
MD <- mahalanobis(Data3, Ave1, Var1)/k
# MDの2乗を計算
Data5 <- cbind(Data, MD,Data3)
# 元のデータ、マハラノビス距離、主成分得点を合わせる
Data5$Y <- factor(Data5$Y)
# Yの列を文字型にする
ggplot(Data5, aes(x=Y, y=MD)) + geom_jitter(size=1, position=position_jitter(0.1))
# 一次元ジター散布図を描くく
MaxMDinUnit <- max(subset(Data5, Data5$Y == 0)$MD)
# Yが0のサンプルのMDのMaxを求める
Data6 <- subset(Data5, Data5$Y == 1)
# Yが1のサンプルを抜き出す
nrow(subset(Data6, Data6$MD < MaxMDinUnit))
# Yが1で、Yが0のサンプルのMDのMaxよりも小さなサンプル数を求める
MT法

カーネル主成分MT法

カーネル主成分MT法 の場合です。 質的変数が混ざっていてもできるようにしてあります。

library(dummies) # ライブラリを読み込み
library(ggplot2)
# ライブラリを読み込み
library(kernlab)
# ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- Data
# 加工用のデータを作る
Y <- Data1$Y
# Yのデータを別に作る
Data1 <- dummy.data.frame(Data)
# 質的変数はダミー変換
Data1$Y <- NULL
# Yの列を消す
Data2 <- subset(Data1, Data$Y == 0)
# Yが0のサンプルのデータを別に作る
Data2$Y <- NULL
# Yの列を消す
model <- kpca(as.matrix(Data2),kernel="rbfdot",kpar=list(sigma=0.1))
# Yが0のサンプルのデータでカーネル主成分分析のモデルを作る
Data3 <- predict(model, as.matrix(Data1))
# 主成分得点を取る
Data4 <- subset(Data3, Data$Y == 0)
# Yが0のサンプルのデータを別に作る
n <- nrow(Data4)
# Yが0のサンプルのサンプル数を計算
Ave1 <- colMeans(Data4)
# Yが0のサンプルの各変数の平均値を計算
Var1 <- var(Data4)*(n-1)/n
# Yが0のサンプルの共分散行列を計算
k <- ncol(Data4)
# 変数の数を計算
MD <- mahalanobis(Data3, Ave1, Var1)/k
# MDの2乗を計算
Data5 <- cbind(Data, MD,Data3)
# 元のデータ、マハラノビス距離、主成分得点を合わせる
Data5$Y <- factor(Data5$Y)
# Yの列を文字型にする
ggplot(Data5, aes(x=Y, y=MD)) + geom_jitter(size=1, position=position_jitter(0.1))
# 一次元ジター散布図を描く
MaxMDinUnit <- max(subset(Data5, Data5$Y == 0)$MD)
# Yが0のサンプルのMDのMaxを求める
Data6 <- subset(Data5, Data5$Y == 1)
# Yが1のサンプルを抜き出す
nrow(subset(Data6, Data6$MD < MaxMDinUnit))
# Yが1で、Yが0のサンプルのMDのMaxよりも小さなサンプル数を求める
MT法

混合分布MT法

混合分布MT法 は、多次元正規分布が複数あるデータが、単位空間になっている時の方法です。 クラスター分析 でグループに分けてから、それぞれに対してMDを計算して、各サンプルはMDの最小値(一番近いグループとのMD)を出力します。
MT

library(ggplot2) # ライブラリを読み込み
library(mclust)
# ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- Data
# 質的変数がある場合にこの行を差し替え
Data2 <- subset(Data1, Data1$Y == 0)
# Yが0のサンプルのデータを別に作る
Data2$Y <- NULL
# Yの列を消す
YData <- Data1$Y
# Yの列を別にしておく
Data4<- as.data.frame(YData)

Data1$Y <- NULL
# Yの列を消す
k <- ncol(Data2)
# 変数の数を計算
kNumber <- 2
#クラスターの数を指定
mc <- Mclust(Data2,kNumber)
# 混合分布で分類。これは3個のグループ分けの場合
output <- mc$classification
# 分類結果の抽出
Data3 <- cbind(Data2, output) 
# 最初のデータセットにグループ分けの結果を付ける
for (i2 in 1:kNumber) {
#
Data101 <- subset(Data3, Data3$output == i2)
# outputが1のサンプルのデータを別に作る
Data101$output <- NULL
# outputの列を消す
n1 <- nrow(Data101)
# outputが1のサンプルのサンプル数を計算
Ave1 <- colMeans(Data101)
# outputが1のサンプルの各変数の平均値を計算
Var1 <- var(Data101)*(n1 - 1)/n1
# outputが1のサンプルの共分散行列を計算
MD <- mahalanobis(Data1, Ave1, Var1)/k
# MD(の2乗)を計算
Data4 <- cbind(Data4,MD)

colnames(Data4)[i2+1]<-i2

}
#
Data5 <- Data4

Data4$YData <- NULL
# Yの列を消す
MinMD<- apply(Data4, 1, min)

MinMD<- as.data.frame(MinMD)

Data5 <- cbind(Data5,MinMD)

Data5$YData <- factor(Data5$YData)
# Yの列を文字型にする
ggplot(Data5, aes(x=YData, y=MinMD)) + geom_jitter(size=1, position=position_jitter(0.1))
# 一次元ジター散布図を描く
MT
2つのグループのどちらからも離れたサンプルは、MDが高くなっています。

近傍法をベースにした方法

近傍法 をベースにした方法は、特定の分布を仮定しないことが利点になります。 その代わり、単位空間のデータの個々の値に敏感に反応するので、 MT法をベースにした方法と比べると、単位空間のデータの内容のチェックがとても大事になります。

LOFをベースにした方法

LOF を使う方法です。

下記のコードを使う場合、「Y = 1」のサンプルは1個が基本です。 2個以上でも計算はできますが、正常値から見たら「異常」でも、異常値同士が近いと「正常領域」と計算上は判定される可能性があります。

「Y = 1」が複数の時は、1個ずつLOFの計算を回す方法ができなくもないのですが、 ただでさえ計算時間が長いのに、それが「Y = 1」のサンプル数の分、倍増することになります。 そこまでする利点がないように思いましたので、 コード例を掲載するのはやめました。

library(Rlof)# ライブラリを読み込み
library(dummies)
# ライブラリを読み込み
library(ggplot2)
# ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- Data
# 加工用のデータの作成
Y <- Data1$Y
# Yのデータを別に作る
Data1 <- dummy.data.frame(Data)
# 質的変数はダミー変換
Data1$Y <- NULL
# Yの列を消す
Data2 <- subset(Data1, Data$Y == 0)
# Yが0のサンプルのデータを別に作る
model <- prcomp(Data2, scale=TRUE, tol=0.01)
# Yが0のサンプルのデータで主成分分析のモデルを作る。
Data3 <- predict(model, as.matrix(Data1))
# 主成分得点を取る
std <- apply(model$x, 2, sd)
# 列ごとに標準偏差を計算
Data4 <- Data3
# 加工用のデータの作成
for (i in 1:ncol(Data3)) {
# ループの始まり。データの列数を数えて同じ回数繰り返す
Data4[,i] <- Data3[,i]/std[i]
# 標準化
}
# ループの終わり
Data4 <- cbind(Data4 ,Y) 
# データを合体
LOF <- lof(Data4,5)
# LOFを計算(5個の近傍点を使う)
Data8 <- as.data.frame(cbind(LOF ,Data))
# データを合体
Data8$Y <- factor(Data8$Y)
# Yの列を文字型にする
ggplot(Data8, aes(x=Y, y=LOF)) + geom_jitter(size=1, position=position_jitter(0.1))
# 一次元ジター散布図を描く

LOF

最小距離法をベースにした方法

最小距離法 を使う方法です。

library(dummies) # ライブラリを読み込み
library(ggplot2)
# ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- Data
# 加工用のデータの作成
Y <- Data1$Y
# Yのデータを別に作る
Data1 <- dummy.data.frame(Data)
# 質的変数はダミー変換
Data1$Y <- NULL
# Yの列を消す
Data2 <- subset(Data1, Data$Y == 0)
# Yが0のサンプルのデータを別に作る
model <- prcomp(Data2, scale=TRUE, tol=0.01)
# Yが0のサンプルのデータで主成分分析のモデルを作る。
Data3 <- predict(model, as.matrix(Data1))
# 主成分得点を取る
std <- apply(model$x, 2, sd)
# 列ごとに標準偏差を計算
Data4 <- Data3
# 加工用のデータの作成
for (i in 1:ncol(Data3)) {
# ループの始まり。データの列数を数えて同じ回数繰り返す
Data4[,i] <- Data3[,i]/std[i]
# 標準化
}
# ループの終わり
Data5 <- dist(Data4, diag = TRUE, upper = TRUE)
# サンプル間の距離を計算
Data5 <-as.matrix(Data5)
# 行列形式に変換
diag(Data5) <- 1000
# 対角成分をとても大きな値にする
Data5 <-as.data.frame(Data5)
# 形式を変換
Data5 <- cbind(Data5 ,Y) 
# データを合体
Data6 <- subset(Data5, Data5$Y == 0)# Y
# Yが0のサンプルのデータを別に作る
Data6$Y <- NULL
# データからYの列を消す
Data7 <- as.data.frame(t(Data6))
# データを転置
min0 <- apply(Data7, 1,min)
# 各行の最小値を計算
Data8 <- as.data.frame(cbind(min0 ,Y))
# データを合体
Data8$Y <- factor(Data8$Y)
# Yの列を文字型にする
ggplot(Data8, aes(x=Y, y=min0)) + geom_jitter(size=1, position=position_jitter(0.1))
# 一次元ジター散布図を描く
MaxMDinUnit <- max(subset(Data8, Data8$Y == 0)$min0)
# Yが0のサンプルのmin0のMaxを求める
Data9 <- subset(Data8, Data8$Y == 1)
# Yが1のサンプルを抜き出す
nrow(subset(Data9, Data9$min0 < MaxMDinUnit))
# Yが1で、Yが0のサンプルのmin0のMaxよりも小さなサンプル数を求める

sammon



参考文献

同志社大学 金明哲先生のページ
Rと判別分析 https://www1.doshisha.ac.jp/~mjin/R/Chap_18/18.html
マハラノビス距離を使う判別分析として、混合分布MT法とほぼ同じ方法が紹介されています。 このページのコードを作る時に参考にさせていただきました。





Tweet データサイエンス教室