多重ロジスティックモデル(Walker-Duncan 法) Last modified: Nov 15, 2004
目的
ロジスティック回帰を 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)
***** 多重ロジスティック回帰
サンプルサイズ 98
生存（打ち切り） 85
死亡（故障） 13
平均値 標準偏差
x1 132.6735 14.50473
x2 223.2347 49.22509
対数尤度 = -36.09198547355
パラメータ推定値
偏回帰係数 標準偏差 t 値 P 値 標準化偏回帰係数
x1 0.008297108 0.02120826 0.3912206 0.69651 0.1203473
x2 0.01138648 0.005739863 1.983755 0.05017 0.5605007
定数項 -5.645581 3.048239 1.852079 0.06712
t の自由度 = 95
以上 以下 期待値 リスク 観察値 故障率 サンプル数
1 -3.090951 -2.6634718 0.5402969 0.05402969 0 0.0000000 10
2 -2.661530 -2.4903799 0.7291859 0.07291859 0 0.0000000 10
3 -2.472638 -2.4368895 0.7120887 0.07912096 0 0.0000000 9
4 -2.398670 -2.2024521 0.9167646 0.09167646 3 0.3000000 10
5 -2.160525 -2.0547805 1.0976992 0.10976992 1 0.1000000 10
6 -2.029713 -1.9455944 1.2061576 0.12061576 1 0.1000000 10
7 -1.944181 -1.7427561 1.3501939 0.13501939 1 0.1000000 10
8 -1.724572 -1.5495385 1.4273797 0.15859774 3 0.3333333 9
9 -1.522529 -1.2521662 1.9860571 0.19860571 0 0.0000000 10
10 -1.238838 0.5878548 3.0341763 0.30341763 4 0.4000000 10
左の2列は，各区間のλの値（最小値と最大値）
解説ページ
直前のページへ戻る E-mail to Shigenobu AOKI
