「統計学関連なんでもあり」の過去ログ--- 047

No.21406　総当たり法による多重ロジスティック回帰分析で、モデル内の説明変数を制限　　【鈴木康弘】　2014/10/24(Fri) 06:59

　青木先生の「総当たり法によるロジスティック回帰分析」をありがたく使わせていただいています。
　ところで，このプログラムは，説明（独立）変数が増えると指数関数的に時間がかかるから，10個以上組み込む時には熟慮して使うように，ということです。

　一方，多重ロジスティックモデル内に，説明変数は5つも6つもいらない，2つか3つで十分です，ただし検討に加えたい説明変数は20以上あるのですが，という事態はよく起こると思います。

　それで，説明変数が2つのモデル式だけを出力するように，プログラムの変更を考えました。

　....考える事，数週間，11行目あたりからの次の部分

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

　を，下の4行に変えてみました。

kumiawase <- t(combn(nv, 2)) #nv個の独立変数から2つ取り出す組み合わせ行列
n <- choose(nv, 2) #組み合わせはn通り
bincomb <- matrix(FALSE, nrow=n, ncol=nv) #n行nv列のFALSE行列を作成
for (j in 1:n){ bincomb[j, kumiawase[j,]] <- TRUE}
#FALSE行列のj行目の，"kumiawase行列のj行目の要素" 列目を，TRUEに変更

　これで説明変数が2つのモデルだけが，ずらり出力されました。
　データに説明変数が20以上あっても一瞬に出力。うれしいな。

　最後に青木先生におねだりです。
こんな機能を正式に取り入れてくださると，うれしいのですが..？

No.21407　Re: 総当たり法による多重ロジスティック回帰分析で，モデル内の説明変数を制限　　【青木繁伸】　2014/10/24(Fri) 12:19

モデルに取り込む独立変数の個数（の範囲）を指定できるように拡張しました。

No.21408　Re: 総当たり法による多重ロジスティック回帰分析で，モデル内の説明変数を制限　　【鈴木康弘】　2014/10/24(Fri) 17:36

　やったー！，ありがとうございます。へたなプログラムでお目汚しした甲斐がありました。

bc <- combn(nv, np, function(x) {y <- rep(FALSE, nv); y[x] <- TRUE; y})

　という1行の意味をこれからゆっくり考えてみます。

No.21409　Re: 総当たり法による多重ロジスティック回帰分析で，モデル内の説明変数を制限　　【青木繁伸】　2014/10/25(Sat) 04:01

bc <- combn(nv, np, function(x) 1:nv %in% x) の方がよかったですね。直しておきました。

No.21410　Re: 総当たり法による多重ロジスティック回帰分析で，モデル内の説明変数を制限　　【鈴木康弘】　2014/10/25(Sat) 12:53

　%in% って？
　ググってみると，集合演算子で，%in%の前にある集合要素が，%in%の後ろにある集合に含まれていればTRUEを，含まれていなければFALSEを返す．．
　な，なるほど，美しい！

No.21657　Re: 総当たり法による多重ロジスティック回帰分析で，モデル内の説明変数を制限　　【鈴木康弘】　2015/05/16(Sat) 10:43

　半年ぶりに追加発言です。（青木先生にではなく，読んでる方々で，Rのスクリプトが読めない方向き）

　総当たり法によるロジスティック回帰分析に追加いただいた，独立変数の個数を制限するスクリプト部品は，もちろん他の分析にも使えます。

　「総当たり法による線形判別分析」ならば，

1．最初の
all.disc <- function ( data, group)

を

all.disc <- function ( data, group, max.p=NULL, min.p=max.p )

に修正，

2．4行目から13行目までのBinConv関数の定義を

BinConv <- function(nv)
{
if(is.null(max.p)){
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))
}
return(bincomb)
}
}

に差し替えると，モデルに組み込む独立変数の個数を制限した総当たり線形判別分析ができます。

　こちらもありがたく使わせていただいてます。

（線形判別分析は，普通はロジスティック判別分析より正判別率が落ちるらしいけれど，正判別率が簡単に出るのが利点ですね。）

No.21695　Re: 総当たり法による多重ロジスティック回帰分析で，モデル内の説明変数を制限　　【鈴木康弘】　2015/06/06(Sat) 12:26

　一ヶ月ぶりの追加発言です。青木先生にではなく，Rのスクリプトが苦手な人向け。

　青木先生の「総当たり法によるロジスティック回帰分析」，devianceとAICだけでなく，回帰係数とそのp値も出ればいいなあ...と思い，

　モデルの独立変数が2個だけなら，こうすればいいのでした。

　スクリプト35行目あたり
result3 <- character(n) # 数値型ベクトル確保
result1 <- result2 <- numeric(n) # 数値型ベクトル確保
　の2行を
result3 <- resulta <- resultb <- character(n)
result1 <- result2 <- result4 <- result5 <- result6 <- result7 <- numeric(n)
　に変更。

　41行目あたり
result1[i] <- result$deviance # deviance
result2[i] <- result$aic # AIC
　の次に
resulta[i] <- rownames(result$coefficients)[2] #1変数め
result4[i] <- result$coefficients[2,1] #その係数
result5[i] <- result$coefficients[2,4] #そのp値
resultb[i] <- rownames(result$coefficients)[3] #2変数め
result6[i] <- result$coefficients[3,1] #その係数
result7[i] <- result$coefficients[3,4] #そのp値
　を追加。

　46行目あたり
return(structure(list(deviance=result1, AIC=result2, formula=result3),
　を
return(structure(list(deviance=result1, AIC=result2, rowname1=resulta, coefficients1=result4, pvalue1=result5, rowname2=resultb, coefficients2=result6, pvalue2=result7, formula=result3),
　に変更。

　print メソッドに行って，その5行目
result <- data.frame(obj$deviance, obj$AIC, obj$formula)
を
result <- data.frame(obj$deviance, obj$AIC, obj$rowname1, obj$coefficients1, obj$pvalue1, obj$rowname2, obj$coefficients2, obj$pvalue2, obj$formula)
　に変更。

　9行目
cat(sprintf("\n%12s %12s %s\n", "deviance", "AIC", "Formula"))
　を
cat(sprintf("\n%12s %12s %12s %12s %12s %12s %12s %12s %s\n", "deviance", "AIC", "1st var", "coefficients", "p-value", "2nd var", "coefficients", "p-value", "Formula"))
　に変更。

　12，13行目
cat(sprintf("%12.5f %12.5f %s\n",
as.double(x[1]), as.double(x[2]), x[3])))
　を
cat(sprintf("%12.5f %12.5f %.12s %12.5f %12.5f %.12s %12.5f %12.5f %s\n",
as.double(x[1]), as.double(x[2]), x[3], as.double(x[4]), as.double(x[5]), x[6], as.double(x[7]), as.double(x[8]), x[9])))

　に変更。これで各変数の回帰係数とp値もいっしょに出るようになりました。

　モデルの独立変数が3つ以上の一般的な場合は，パス。(^_^;)
　あらためて青木先生に感謝します。

No.21706　Re: 総当たり法による多重ロジスティック回帰分析で，モデル内の説明変数を制限　　【鈴木康弘】　2015/06/21(Sun) 09:20

　また追加発言です。このスレッドを自分だけで使ってるみたいですが．．．

　「総当たり法によるロジスティック回帰分析」を判別分析に使った場合，正判別率もいっしょに出ればいいなあ...と思い，
　データの従属変数が0か1の形式になってるなら，こうすればいいのでした。

　スクリプト35行目あたり
result1 <- result2 <- numeric(n) # 数値型ベクトル確保
　を
result1 <- result2 <- atari <- numeric(n) # 数値型ベクトル確保
　に変更。

　41行目あたり
result1[i] <- result$deviance # deviance
result2[i] <- result$aic # AIC
　の次に
kakeru <- (ans$fitted.values -0.5)*(df[,ncol(df)]-0.5) #回帰値と従属変数値から0.5ずつ引いてかけ算。正判別なら正，誤判別なら負になる。
atari[i] <- length(kakeru[kakeru >0]) #kakeruが正の例数が正判別個数
　を追加。

　46行目あたり
return(structure(list(deviance=result1, AIC=result2, formula=result3),
　を
return(structure(list(deviance=result1, AIC=result2,
atari=atari, cases=nrow(df), formula=result3),
　に変更。

　print メソッドに行って，その5行目
result <- data.frame(obj$deviance, obj$AIC, obj$formula)
　を
result <- data.frame(obj$deviance, obj$AIC, obj$atari, obj$cases, obj$formula)
　に変更，

　9行目
cat(sprintf("\n%12s %12s %s\n", "deviance", "AIC", "Formula")) # 表頭
　を
cat(sprintf("\n%12s %12s %12s %12s %s\n", "deviance", "AIC", "correct rate", "correct cases", "Formula")) # 表頭
　に変更，

　12，13行目
cat(sprintf("%12.5f %12.5f %s\n",
as.double(x[1]), as.double(x[2]), x[3])))
　を
cat(sprintf("%12.5f %12.5f %12.5f %12.1f %s\n",
as.double(x[1]), as.double(x[2]), as.double(x[3])/as.double(x[4]), as.double(x[3]), x[5])))
　に変更。

　これで正判別率もいっしょに出るようになりました。
　2つ前に，部分的総当たりを線形判別分析に応用する発言をしましたが，そちらは僕は不要になりました。(^_^;）

No.21909　Re: 総当たり法による多重ロジスティック回帰分析で，モデル内の説明変数を制限　　【鈴木康弘】　2016/02/14(Sun) 11:18

　半年ぶりに，2つ前の発言に追加です。

　「総当たり法によるロジスティック回帰分析」で，偏回帰係数だけでなく，標準化偏回帰係数も出るといいなあ，と思い．．，
　モデルの独立変数が2つだけならば，こうすればいいのでした。

　2つ前の発言の変数をそのまま使うなら，40数行目に
stdresult4[i] <- sd (df[ , resulta[i] ]) * result4[i]
stdresult6[i] <- sd (df[ , resultb[i] ]) * result6[i]

　を追加。
　その前後に，ベクトル確保とか表示のための処理とかを，前発言に準じてやれば，できあがり。

（できてみると簡単だけれど，R初心者としては，df[ , resulta[i] ] にたどり着く前に，df$resulta[i] とか，df[ , 'resulta[i]'] とかやって，計算してくれない．．．と半日悩みました。）