目的

総当たり法によるロジスティック回帰分析を行う。
モデルに組み込む独立変数の個数を制限することもできる（部分的総当たり法）。

使用法

all.logistic(df, limit=10, max.p=NULL, min.p=max.p)
print.all.logistic(obj, sort.by=c("deviance", "AIC"), models=20)

引数

df        データフレーム（従属変数を最後の列に置く）
limit     独立変数の個数の制限値。標準では 10 個
          注意：独立変数が多いと莫大な計算時間がかかり，出力行数も増える
          　　　独立変数が1個増えるだけで，計算時間はほぼ 2 倍になり，出力行数は 2 倍になる
max.p     完全な総当たり探索を行うとき（デフォルト）は指定しない
          モデルに組み込む独立変数の数を制限するときは，その値
min.p     モデルに組み込む独立変数の数の下限を指定するときは，その値
          min.p を指定しないときには，モデルに取り込む独立変数の数は max.p 個に限られる
obj       all.logistic が返すオブジェクト
sort.by   モデルの適否を表す2種の統計量のどれでソートして結果を表示するかの選択
          デフォルト（"deviance"）では 結果を deviance の小さい順にソートして出力する。
          "AIC" を指定すると，結果を AIC の小さい順にソートして出力する。
models    よいモデルから mmodels 番目までを表示する（デフォルトは 20）

ソース

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

# 総当たり法によるロジスティック回帰を行う
# 　データフレームには，分析に使用する独立変数と従属変数のみを含むこと。
# 　また，従属変数は最終列に置くこと。
#
all.logistic <- function(    df,                             # データフレーム（独立変数）
                                limit=10,                       # 独立変数の個数の上限（数が多いと計算時間が指数的に増える）
                                max.p=NULL,                     # モデルに組み入れる独立変数の数の上限
                                min.p=max.p)                    # モデルに組み入れる独立変数の数の下限
{
        df <- subset(df, complete.cases(df))                 # 欠損値を持つケースを除く
        nv <- ncol(df)-1                                     # 独立変数の個数
        if (is.null(max.p)) {                                   # 完全な総当たり
                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,]
        }
        else {                                                  # モデルに組み入れる独立変数の数を制限する
                if (max.p < 1 || max.p > nv || min.p < 1 || min.p > nv || min.p > max.p) {
                        stop("モデルに取り込む独立変数の数の指定が不正")
                }
                bincomb <- NULL
                for (np in min.p:max.p) {
                        bc <- combn(nv, np, function(x) 1:nv %in% x)
                        bincomb <- rbind(bincomb, t(bc))
                }
        }
        n <- nrow(bincomb)
        name <- names(df)                                    # 変数名を取り出す
        depname <- name[nv+1]
        name <- name[1:nv]
        result3 <- character(n)                                      # 数値型ベクトル確保
        result1 <- result2 <- numeric(n)                  # 数値型ベクトル確保
        for (i in 1:n) {                                        # 独立変数の全ての組み合わせについて，
                str <- name[bincomb[i,]]                     # どの独立変数が使われるかを割り出す
                form <- reformulate(str, depname)            # モデル式を作る（"formula" クラス）
                ans <- glm(form, df, family=binomial)                # ロジスティック回帰分析の結果
                result <- summary(ans)
                result1[i] <- result$deviance                        # deviance
                result2[i] <- result$aic                     # AIC
                temp <- as.character(form)                   # モデル式を文字列に変換
                result3[i] <- paste(temp[2], "~", temp[3])   # モデル式を記録
        }
        return(structure(list(deviance=result1, AIC=result2, formula=result3),
        class="all.logistic"))
}
# print メソッド
print.all.logistic <- function( obj,                         # "all.logistic" クラスのオブジェクトをプリント
                                sort.by=c("deviance", "AIC"),   # 結果を何で並べ替えるかを指示
                                models=20)                      # 良い方から何番目まで出力するか
{
        result <- data.frame(obj$deviance, obj$AIC, obj$formula)
        sort.by <- match.arg(sort.by)
        o <- order(switch (sort.by, "deviance"=result[,1], "AIC"=result[,2]))
        result <- result[o,]
        cat(sprintf("\n%12s %12s   %s\n", "deviance", "AIC", "Formula"))        # 表頭
        models <- min(models, nrow(result))
        apply(result[1:models,], 1, function(x)                 # 各行の出力
                cat(sprintf("%12.5f %12.5f   %s\n",
                as.double(x[1]), as.double(x[2]), x[3])))
        invisible(result)
}


使用例

> lr <- read.table("lr.data", header=TRUE)
> all.logistic(lr)	# lr.data を分析してみる

    deviance          AIC   Formula
    72.18397     78.18397   y ~ x1 + x2
    72.33492     76.33492   y ~ x2
    76.25592     80.25592   y ~ x1

> all.logistic(lr, max.p=2)	# 独立変数を 2 個だけ含むモデルを探索

    deviance          AIC   Formula
    72.18397     78.18397   y ~ x1 + x2

> all.logistic(lr, max.p=1)	# 独立変数を 1 個だけ含むモデルを探索

    deviance          AIC   Formula
    72.33492     76.33492   y ~ x2
    76.25592     80.25592   y ~ x1

> all.logistic(lr, max.p=2, min.p=1)	# 独立変数を 1〜2 個だけ含むモデルを探索

    deviance          AIC   Formula
    72.18397     78.18397   y ~ x1 + x2
    72.33492     76.33492   y ~ x2
    76.25592     80.25592   y ~ x1