目的

総当たり法によるロジスティック回帰分析を行う。

使用法

all.logistic(df, limit)
print.all.logistic(obj, sort.by=c("adj", "rsq", "AIC"))

引数

df        データフレーム（従属変数を最後の列に置く）
limit     独立変数の個数の制限値。標準では 10 個
          注意：独立変数が多いと莫大な計算時間がかかり，出力行数も増える
          　　　独立変数が1個増えるだけで，計算時間はほぼ 2 倍になり，出力行数は 2 倍になる
obj       
sort.by   モデルの適否を表す3種の統計量のどれでソートして結果を表示するかの選択
          デフォルト（"adj"）では 結果を 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)                       # 独立変数の個数の上限（数が多いと計算時間が指数的に増える）
{
        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]
        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