目的

ロジスティック回帰を Walker-Duncan 法により行う。
R では，glm 関数を用いて分析を行うこともできる。
（多項ロジスティックモデルは nnet パッケージにある multinom 関数，
　順序ロジスティックモデルは MASS パッケージにある polr 関数。）

使用法

logistic.regression(data)

引数

data	データマトリックス。最終列が従属変数（0/1データであること）。

ソース

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

# 多重ロジスティックモデル(Walker-Duncan 法)
logistic.regression <- function(data)                                                   # データ行列（最右端の列が従属変数）
{

        printf <- function(fmt, ...)                                                    # 書式付き print 関数
        {
                cat(sprintf(fmt, ...))
        }

        diff <- function(x, m, n, coeff)                                                # 数値微分
        {
                mp1 <- m+1
                temp <- coeff[mp1]+coeff[1:m]%*%t(x[, 1:m, drop=FALSE])
                p0 <- 1/(1+exp(-temp))
                p1 <- 1-p0
                pp <- ifelse(x[, mp1] == 1, p1, -p0)
                diff1 <- numeric(mp1)
                diff1[mp1] <- sum(pp)
                diff1[1:m] <- diff1[1:m]+colSums(x[,1:m, drop=FALSE]*pp)
                temp <- -x[, 1:m, drop=FALSE]*p0*p1
                diff2 <- matrix(0, nrow=mp1, ncol=mp1)
                diff2[1:m, 1:m] <- t(x[,1:m, drop=FALSE])%*%as.matrix(temp)
                diff2[lower.tri(diff2, diag=FALSE)] <- 0
                diff2[mp1, mp1] <- -sum(p0*p1)
                diff2[1:m, mp1] <- diff2[1:m, mp1]+colSums(temp[,1:m, drop=FALSE])
                diff2 <- diff2+t(diff2)
                diag(diff2) <- diag(diff2)/2
                return(list(diff1=diff1, diff2=diff2))
        }

        llh <- function(x, m, coeff)                                                    # 対数尤度
        {
                temp <- colSums(t(x[,1:m])*coeff[1:m])+coeff[m+1]
                -sum(log(1+exp(ifelse(x[,m+1] == 1, -temp, temp))))
        }

        newton.logist <- function(x, m, n, sds)                                         # ニュートン法によるあてはめ
        {
                mp1 <- m+1
                coeff0 <- numeric(mp1)
                coeff <- rep(1e-14, mp1)
                for (itr in 1:500) {
                        temp <- diff(x, m, n, coeff)
                        coeff0 <- solve(temp$diff2, temp$diff1)
                        converge <- all(abs(coeff0/coeff) < 1e-10)
                        coeff <- coeff-coeff0
                        if (converge) {
                                break
                        }
                }
                stopifnot(converge == TRUE)
                se <- sqrt(-diag(solve(temp$diff2)))
                printf("\n対数尤度 = %.14g\n", llh(x, m, coeff))
                printf("\nパラメータ推定値\n\n")
                printf("         %14s %14s %12s %8s %14s\n", "偏回帰係数", "標準偏差", "t 値", "P 値", "標準化偏回帰係数")
                t <- abs(coeff[1:m]/se[1:m])
                p <- pt(t, n-mp1, lower.tail=FALSE)*2
                for (i in 1:m) {
                        printf("%8s %14.7g %14.7g %12.7g %8.5f %14.7g\n", vnames[i], coeff[i], se[i], t[i], p[i], coeff[i]*sds[i])
                }
                t <- abs(coeff[mp1]/se[mp1])
                p <- pt(t, n-mp1, lower.tail=FALSE)*2
                printf("%8s %14.7g %14.7g %12.7g %8.5f\n%60s%i\n\n", "　定数項", coeff[mp1], se[mp1], t, p, "t の自由度 = ", as.integer(n-mp1))
                return(coeff)
        }

        fitness <- function(data, m, n, coeff)                                          # 当てはまり具合の確認
        {
                lambda <- coeff[1:m]%*%t(data[,1:m])+coeff[m+1]                         # λ
                pred <- 1/(1+exp(-lambda))                                              # リスクの予測値
                y <- data[,m+1]                                                         # エンドポイント
                div <- round(seq(0, n-1, by=n/10),0)[-1]                                # 対象をほぼ同数に10区分するときのサンプル数
                xs <- sort(lambda)                                                      # λをソートして
                div2 <- (xs[div]+xs[div+1])/2                                           # 分割点を求める
                g <- findInterval(lambda, div2)+1                                       # λがどの区分にはいるか群分け
                from <- c(min(lambda), xs[div+1])                                       # λの区間 [from, to]
                to <- c(xs[div], max(lambda))
                mid <- (from+to)/2                                                      # 区間の級中心
                pred <- tapply(pred, g, sum)                                            # エンドポイントを持つものの期待値
                obs <- tapply(y, g, sum)                                                # エンドポイントを持つものの観察値
                cnt <- tapply(y, g, length)                                             # サンプル数
                valid <- as.integer(names(pred))                                        # 区間幅が正のもの
                from <- from[valid]                                                     # 調整
                to <- to[valid]                                                         # 調整
                table <- data.frame("以上"=from, "以下"=to,
                                "期待値"=pred, "リスク"=pred/cnt,
                                "観察値"=obs, "故障率"=obs/cnt, "サンプル数"=cnt)
                print(table)
                printf("%60s%s\n", "", "左の2列は，各区間のλの値（最小値と最大値）")
                plot(c(min(from), max(to)), c(0, max(ifelse(pred>obs, pred, obs)/cnt)), # リスク，故障率，理論曲線の描画
                        type="n", xlab="Lambda", ylab="Risk")
                for (i in 1:10) {
                        lines(c(from[i], to[i]), rep(pred[i]/cnt[i], 2), lty=3)
                        lines(c(from[i], to[i]), rep(obs[i]/cnt[i], 2))
                }
                x <- seq(min(lambda), max(lambda), length.out=100)
                lines(x, 1/(1+exp(-x)))
        }

# logistic 関数本体

        data <- subset(data, complete.cases(data))                                      # 欠損値を除く
        n <- nrow(data)                                                                 # サンプルサイズ
        mp1 <- ncol(data)                                                               # 列数
        m <- mp1-1                                                                      # 独立変数の個数
        vnames <- colnames(data)                                                        # 変数名を取り出す
        if (is.null(vnames)) vnames <- paste("変数", 1:mp1, sep="")                     # 変数名がなかったら定義する
        num <- table(data[, mp1])                                                       # サンプルの内訳

        printf("***** 多重ロジスティック回帰\n\n")
        printf("サンプルサイズ　　　　%5i\n", n)
        printf("　　生存（打ち切り）　%5i\n", as.integer(num[1]))
        printf("　　死亡（故障）　　　%5i\n", as.integer(num[2]))
        if (num[1] == 0 || num[2] == 0 || num[1]+num[2] == 2 || n <= mp1) {
                stop("有効ケース数が 1 以下です\n")
        }

        means <- colMeans(data)                                                         # 各変数の平均値
        sds <- apply(data, 2, sd)                                                       # 各変数の標準偏差
        printf("\n          %15s %15s\n", "平均値", "標準偏差") 
        for (i in 1:m) {
                printf("%8s  %15.7g %15.7g\n", vnames[i], means[i], sds[i])
        }

        coeff <- newton.logist(data, m, n, sds)                                         # パラメータ推定
        fitness(data, m, n, coeff)                                                      # 当てはまり具合の確認
        invisible(coeff)                                                                # 係数のみを返す（自動的に表示はされない）
}


使用例

分析に使うデータは lr.data という名前で用意されているとする（lr.data を見てみる）。

data <- read.table("lr.data", header=TRUE)
logistic.regression(data)

出力結果例

> data <- read.table("lr.data", header=TRUE)
> logistic.regression(data)
***** Multiple logistic regression

Number of cases    98
  Negarive cases                  85
  Positive cases                  13

                     mean            s.d.
    Var1         132.6735        14.50473
    Var2         223.2347        49.22509

log likelihood = -36.09198547355

estimates of coefficients

                      B        S.E.(B)      t value  P value           beta
    Var1    0.008297108     0.02120826    0.3912206  0.69651      0.1203473
    Var2     0.01138648    0.005739863     1.983755  0.05017      0.5605007
constant      -5.645581       3.048239     1.852079  0.06712

                                                d.f. of t = 95


 解説ページ