目的

総当たり法による重回帰分析を行う。

使用法

All.possible.subset.selection(df, limit)
print.All.possible.subset.selection(obj, sort.by=c("adj", "rsq", "AIC"), models=20)

引数

df        データフレーム（従属変数を最後の列に置く）
limit     独立変数の個数の制限値。標準では 10 個
          注意：独立変数が多いと莫大な計算時間がかかり，出力行数も増える
          　　　独立変数が1個増えるだけで，計算時間はほぼ 2 倍になり，出力行数は 2 倍になる
obj       
sort.by   モデルの適否を表す3種の統計量のどれでソートして結果を表示するかの選択
          デフォルト（"adj"）では 結果を自由度調整済み重相関係数の二乗（adj.r.square）の大きい順にソートして出力する。
          "rsq" を指定すると，結果を重相関係数の二乗（決定係数 r.square）の大きい順にソートして出力する。
          "AIC" を指定すると，結果を AIC の小さい順にソートして出力する。
models    よいモデルから mmodels 番目までを表示する（デフォルトは 20）

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/All_possible_subset_selection.R", encoding="euc-jp")

# 総当たり法による重回帰分析を行う
# 　データフレームには，分析に使用する独立変数と従属変数のみを含むこと。
# 　また，従属変数は最終列に置くこと。
#
All.possible.subset.selection <- function(df,                                # データフレーム（独立変数，従属変数）
                                          limit=10)                     # 独立変数の個数の上限（数が多いと計算時間が指数的に増える）
{
        df <- subset(df, complete.cases(df))                         # 欠損値を持つケースを除く
        nv <- ncol(df)-1                                             # 独立変数の個数
        if (nv > limit) {                                               # limit より多いと分析を中断する
                stop(paste("独立変数が", limit,
                "個以上である（多すぎる）。\n",
                "limit 引数で変更できる", paste=""))
        }
        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,]
        n <- n-1
        name <- names(df)                                            # 変数名を取り出す
        depname <- name[nv+1]
        name <- name[1:nv]
        result4 <- character(n)                                              # 数値型ベクトル確保
        result1 <- result2 <- result3 <- numeric(n)                    # 数値型ベクトル確保
        for (i in 1:n) {                                                # 独立変数の全ての組み合わせについて，
                str <- name[bincomb[i,]]                             # どの独立変数が使われるかを割り出す
                form <- reformulate(str, depname)                    # モデル式を作る（"formula" クラス）
                ans <- lm(form, df)                                  # 重回帰分析の結果
                result <- summary(ans)
                result1[i] <- result$r.square                                # 重相関係数の二乗（決定係数）
                result2[i] <- result$adj.r.square                    # 自由度調整済み重相関係数の二乗
                result3[i] <- AIC(ans)                                       # AIC
                temp <- as.character(form)                           # モデル式を文字列に変換
                result4[i] <- paste(temp[2], "~", temp[3])           # モデル式を記録
        }
        return(structure(list(rsq=result1, adj=result2, AIC=result3, form=result4),
        class="all.possible.subset.selection"))
}
print.all.possible.subset.selection <- function( obj,                        # "all.possible.subset.selection" クラスのオブジェクトをプリント
                                                 sort.by=c("adj", "rsq", "AIC"), # 結果を何で並べ替えるかを指示
                                                 models=20)             # 良い方から何番目まで出力するか
{
        result <- data.frame(obj$rsq, obj$adj, obj$AIC, obj$form)
        sort.by <- match.arg(sort.by)
        o <- order(switch (sort.by, "rsq"=result[,1], "adj"=result[,2], "AIC"=result[,3]), decreasing = sort.by != "AIC")
        result <- result[o,]
        cat("\nR square  Adjusted R square       AIC       Formula\n")  # 表頭
        models <- min(models, nrow(result))
        apply(result[1:models,], 1, function(x)                         # 各行の出力
                cat(sprintf("%8.5f %13.5f  %13.5f       %s\n",
                as.double(x[1]), as.double(x[2]), as.double(x[3]), x[4])))
        invisible(result)
}


使用例

All.possible.subset.selection(iris[1:4])	# 1〜3列目を独立変数，4列目（Petal.Width）を従属変数とする

> str(iris)	# iris の構成は以下の通り
`data.frame':	150 obs. of  5 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Species     : Factor w/ 3 levels "setosa","versic..",..: 1 1 1 1 1 1 1 1 1 1 ...

出力結果例

> All.possible.subset.selection(iris[1:4]) # デフォルトでは自由度調整済みの重相関係数の二乗で並べ替える

R square  Adjusted R square       AIC       Formula
 0.93785       0.93657      -63.50218       Petal.Width ~ Sepal.Length + Sepal.Width + Petal.Length
 0.92975       0.92879      -47.11940       Petal.Width ~ Sepal.Width + Petal.Length
 0.92902       0.92806      -45.58470       Petal.Width ~ Sepal.Length + Petal.Length
 0.92711       0.92662      -43.59109       Petal.Width ~ Petal.Length
 0.74293       0.73943      147.46929       Petal.Width ~ Sepal.Length + Sepal.Width
 0.66903       0.66679      183.37107       Petal.Width ~ Sepal.Length
 0.13405       0.12820      327.64025       Petal.Width ~ Sepal.Width

> a <- All.possible.subset.selection(iris[1:4]) # 他の基準で並べ替えるときは，付値してから，
> print(a, sort.by="AIC")                          # print メソッドの sort.by 引数で指定する

R square  Adjusted R square       AIC       Formula
 0.93785       0.93657      -63.50218       Petal.Width ~ Sepal.Length + Sepal.Width + Petal.Length
 0.92975       0.92879      -47.11940       Petal.Width ~ Sepal.Width + Petal.Length
 0.92902       0.92806      -45.58470       Petal.Width ~ Sepal.Length + Petal.Length
 0.92711       0.92662      -43.59109       Petal.Width ~ Petal.Length
 0.74293       0.73943      147.46929       Petal.Width ~ Sepal.Length + Sepal.Width
 0.66903       0.66679      183.37107       Petal.Width ~ Sepal.Length
 0.13405       0.12820      327.64025       Petal.Width ~ Sepal.Width