目的線形判別分析を行う関数 candisを下請けとして,総当たりによる線形判別分析を行い,正判別率の高い判別関数を探索する。
使用法 all.disc(data, group) print.all.disc(obj) 引数 data データ行列(行がケース,列が変数) group 各ケースがどの群であるかを表す変数 データフレーム中の変数の場合には,iris[, 5] ではなく iris[5] のように (1 列だけのデータフレームになるように)指定すること(変数名を関数に渡せるように) ソース 以下のほかに,disc が必要。 インストールは,以下の 1 行をコピーし,R コンソールにペーストする source("http://aoki2.si.gunma-u.ac.jp/R/src/all.disc.R", encoding="euc-jp") # 総当たり法による線形判別分析 all.disc <- function( data, # 説明変数のデータフレーム(データ行列) group) # 群を表す変数(データフレームから指定するときには iris[,5] ではなく,iris[5] のように) { BinConv <- function(nv) { n <- 2^nv # 独立変数を取り出す取り出し方 bincomb <- matrix(FALSE, nrow=n, ncol=nv) # e1071 パッケージの bincombinations より for (j in 1:nv) { bincomb[, j] <- rep(c(rep(FALSE, n/2^j), rep(TRUE, n/2^j)), length = n) } bincomb <- bincomb[-1,] return(bincomb) } nv <- ncol(data) vname <- colnames(data) # 変数名(なければ作る) if (is.null(vname)) { vname <- colnames(data) <- paste("x", 1:nv, sep="") } gname <- names(group) if (is.null(gname)) { # group を,データフレームから iris[, 5] のようにすると,名前がなくなる gname <- "" # iris[5] のように 1 列のデータフレームとして指定すること } group <- factor(as.matrix(group)) ok <- complete.cases(data, group) # 欠損値を含まないケースだけを対象にする data <- data[ok,] group <- group[ok] n <- nrow(data) bincomb <- BinConv(nv) nr <- nrow(bincomb) correct <- numeric(nr) # 正準判別関数の正判別率を基準とする for (i in 1:nr) { dat <- data[, bincomb[i,], drop=FALSE] a <- disc(dat, group) # 別途用意してある disc.R で定義 correct[i] <- sum(a$correct)/n # 正判別率を記録しておく } ans <- data.frame(correct, bincomb) colnames(ans) <- c("correct rate", vname) return(structure(list(ans=ans, name=vname, gname=gname), class="all.disc")) } # print メソッド print.all.disc <- function(obj) { ans <- obj$ans name <- obj$name gname <- obj$gname o <- order(ans[, 1], decreasing=TRUE) ans <- ans[o,] nc <- ncol(ans) cat("\ncorrect rate formula\n") for (i in 1:nrow(ans)) { cat(sprintf("%10.5f %s ~ %s\n",ans[i, 1], gname, paste(name[as.matrix(ans[i, 2:nc])], collapse=" + "))) } invisible(ans) # 結果をソートしただけのものを返す } 使用例 > all.disc(as.matrix(iris[1:4]), iris[5]) correct rate formula 0.98000 Species ~ Sepal.Width + Petal.Length + Petal.Width 0.98000 Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width 0.97333 Species ~ Sepal.Length + Petal.Length + Petal.Width 0.96667 Species ~ Sepal.Width + Petal.Width 0.96667 Species ~ Sepal.Length + Petal.Length 0.96667 Species ~ Sepal.Length + Sepal.Width + Petal.Length 0.96000 Species ~ Petal.Width 0.96000 Species ~ Petal.Length + Petal.Width 0.96000 Species ~ Sepal.Length + Petal.Width 0.96000 Species ~ Sepal.Length + Sepal.Width + Petal.Width 0.95333 Species ~ Sepal.Width + Petal.Length 0.94667 Species ~ Petal.Length 0.80000 Species ~ Sepal.Length + Sepal.Width 0.74667 Species ~ Sepal.Length 0.55333 Species ~ Sepal.Width