Rで外れ値検定 | グラブス検定を行い箱髭図と蜂群図を描く

remove(list=ls())
imported.filename <- "testdata.csv" #読み込みたいファイル名に変える
test.all <- TRUE #全部のカラムに外れ値検定を行うならTRUE、そうでなければFALSE

#どのカラムに対して検定を行うか？の情報を配列text.collumに格納する
imported.table <- read.csv(imported.filename)
test.collum <- logical(ncol(imported.table))
if (test.all == TRUE) {
	for (i in 1:ncol(imported.table)) {
		test.collum[i] <- TRUE
		i <- i + 1
	}
} else {  
	test.collum <- c(FALSE, TRUE, TRUE) #行いたくないカラムをFALSEにする
}

#gurubbs.testを実行した結果を、配列result[]とp.value[]に格納する。この時点では文字列型
result <- character(ncol(imported.table))
p.value <- numeric(ncol(imported.table))
i <- 1
for (i in 1:ncol(imported.table)) {
	result.temp <- grubbs.test(imported.table[[i]])
	result[i] <- result.temp[[2]]
	p.value[i] <- result.temp[[3]]
	i <- i + 1
}

#今後の処理を行いやすくするため、外れ値が何行目かをまとめておく
#resultは文字列型なので、数値型の配列outlier[]を作って格納しておく
outlier.row <- numeric(ncol(imported.table))
outlier <- numeric(ncol(imported.table))
i <- 1
for (i in 1:ncol(imported.table)) {
	#grep()は文字列検索して存在したらTRUEを返す
	#外れ値が最大値の場合はgrep()で1が返るため、if関数の引数が1になる
	#外れ値が最小値の場合はgrep()で長さ0のベクトルが返るため、if関数の引数が0になる
	#max()中のna.rm=TRUEは欠損値を除外するオプション
	if (length(grep("highest", result[i]))) {
		outlier.row[i] <- which.max(imported.table[[i]])
		outlier[i] <- max(imported.table[[i]], na.rm = TRUE)
	} else {
		outlier.row[i] <- which.min(imported.table[[i]])
		outlier[i] <- min(imported.table[[i]], na.rm = TRUE)
	}
}

#p.valueが0.05未満かどうかの情報を、配列p.value.lowに格納する
p.value.low <- logical(ncol(imported.table))
i <- 1
for (i in 1:ncol(imported.table)) {
	p.value.low[i] <- p.value[i] < 0.05
	i <- i + 1
}

#除外したくない列として手動で指定したカラムがあれば反映する
exclude.ok <- logical(ncol(imported.table))
if (test.all == TRUE) {
	exclude.ok <- p.value.low
} else {
	exclude.ok <- p.value.low & test.collum
}

#外れ値を除外した表を作る
excluded.table <- imported.table
i <- 1
for (i in 1:ncol(imported.table)) {
	if (exclude.ok[i]) {
		excluded.table[outlier.row[i], i] <- NA
	}
	i <- i + 1
}
View(excluded.table)

#外れ値を格納するための表を作っておく
outlier.table <- imported.table
i <- 1
for (i in 1:ncol(imported.table)) {
	outlier.table[i] <- NA
	i <- i + 1
} 

#outlier.tableに外れ値を格納する
i <- 1
for (i in 1:ncol(imported.table)) {
	if (exclude.ok[i]) {
		outlier.table[outlier.row[i], i] <- outlier[i]
	}
	i <- i + 1
}

#result.tableに結果をまとめて表示する
result.table <- rbind(result, outlier, p.value, outlier.row, p.value.low, test.collum, exclude.ok)
colnames(result.table) <- colnames(imported.table)
View(result.table)

win.graph(10, 7) #プロットエリアのサイズ、単位はインチ
	#win.graph()が無いとRStudio右下にグラフが表示され、縦横サイズを保存前に聞かれる
	par(tcl = -0.5, #目盛の長さ：負の数だと外向き
		las = 1, #目盛の数字の向き
		lwd = 2, #線の太さ
		cex.axis=2 #文字の大きさ
	)

	yMax <- 7 #y軸の範囲
	box.color <- "white" #多色にする場合はc("","","")と表記
	beeswarm.color <- c("black", "blue", "red")
	
	#これより下は変更する必要無し
	boxplot(excluded.table,
			range = 0,#0は外れ値無しという指定
			ylim = c(-0.01,yMax), #最小値が0だとプロットが潰れてしまう
			medlwd = 3, #中央値の太さ
			col = box.color
	)

	beeswarm(excluded.table,
			method = "center",
			col = beeswarm.color,
			pch = 16, #プロットの種類：16は黒丸
			cex = 1.5, #プロットの大きさ
			ylim = c(-0.01,yMax),
			add = TRUE #グラフを上書きする
	)

	beeswarm(outlier.table,#外れ値を描画する
			pch = 1, #プロットの種類：1は白丸
			cex = 1.5, #プロットの大きさ
			ylim = c(-0.01,yMax),
			add = TRUE #グラフを上書きする
	)
	
	axis(side = 2,lwd = 2) #目盛線の太さを上書きする,side2は左