人力検索はてな
モバイル版を表示しています。PC版はこちら
i-mobile

20個の2変数x,yデータで以下です。
39.2964, 0
43.2808, 1
39.1153, 0
38.9163, 0
38.2939, 0
45.0221, 1
38.5792, 0
35.7308, 0
41.1779, 1
39.0925, 0
39.0529, 0
35.9808, 0
38.5518, 0
39.507 , 0
40.0184, 1
39.4895, 0
42.835 , 1
38.7532, 0
37.0281, 0
41.0135, 1
で、Rでロジスティック分析
d <- data.frame(x = c(43.2808, 39.1153, 38.9163, 38.2939, 45.0221, 38.5792,
35.7308, 41.1779, 39.0925, 39.0529, 38.5518, 35.9808, 39.507, 40.0184, 39.4895,
42.835, 38.7532, 37.0281, 41.0135), y = c(1, 0, 0, 0, 1, 0, 0, 1, 0, 0,
0, 0, 0, 1, 0, 1, 0, 0, 0))
ans <- glm(y ~ x, d, family = binomial)
としますが、この2変数の分析をRで行う際に、xおよびyについて、どのような演算処理が行われるのでしょうか?わかりやすく、教えていただければ、幸いです。

●質問者: kojiro_i619
●カテゴリ:ウェブ制作
○ 状態 :終了
└ 回答数 : 1/1件

▽最新の回答へ

1 ● dilutionist
●200ポイント

ロジスティック回帰は一般化線形モデル(GLM)で扱うモデルのうち、従属変数が二項分布をとる場合についてロジット変換を使って回帰モデルを立てるもので、Rではglm()でfamilyをbinomial、リンク関数にlogitを使用します。
質問にあるモデル式ではリンク関数の記述が省略されているので、binomialの省略時デフォルトであるlogitが使用されます。

ans <- glm(y ~ x, d, family = binomial(logit))

GLMにおけるパラメータ推定は最尤法が用いられ、反復再重み付け最小二乗法(IRLS)が使用されます。最尤法の詳しい説明は参考書やGLMの解説をしたサイトなどいろいろあると思いますので、そういったものを参考にして下さい。
http://www7.atwiki.jp/hayatoiijima/sp/pages/18.html
乱暴に言いますと従属変数をロジット変換して線形回帰に持ち込むことで、データがロジスティック曲線にフィットするようなパラメータ推定を行えるようにしている、という感じでしょうか。

さてこのモデルにおいてどのような演算処理が行われているのかというご質問ですが、ロジスティック回帰の手順をRのコードで書いたものを付けておきますので参考にしてみて下さい。ここに書いたコードはこちらから拝借しました。
https://socserv.socsci.mcmaster.ca/jfox/Books/Companion-1E/Ch8-script.txt

まずデータフレームの読み込み

d <- data.frame(x = c(43.2808, 39.1153, 38.9163, 38.2939, 45.0221, 38.5792, 35.7308, 41.1779, 39.0925, 39.0529, 38.5518, 35.9808, 39.507, 40.0184, 39.4895, 42.835, 38.7532, 37.0281, 41.0135), y = c(1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0))

ロジスティック回帰(John Fox版)

lreg <- function(x, y, max.iter=10, tol=1E-6){
# x はモデル行列(説明変数のベクトル)
# y は反応ベクトル(0 または 1)
# max.iter は最大反復回数
# tol は収束基準
x <- cbind(1, x)
b <- b.last <- rep(0, ncol(x))
it <- 0
while (it <= max.iter){
p <- as.vector(1/(1 + exp(-x %*% b)))
V <- diag(p * (1 - p))
var.b <- solve(t(x) %*% V %*% x)
b <- b + var.b %*% t(x) %*% (y - p)
if (max(abs(b - b.last)/(abs(b.last) + 0.01*tol)) < tol) break
b.last <- b
it <- it + 1
}
if (it > max.iter) warning('maximum iterations exceeded')
list(coefficients=as.vector(b), var=var.b, iterations=it)
}
lreg(d$x, d$y)

上のコードの出力

$coefficients
[1] -82.219446 2.022379

$var
 x
 1812.2569 -45.043403
x -45.0434 1.120192

$iterations
[1] 7

glm()の出力と比較すると同じ係数と切片に収束していることがわかります。

ans <- glm(y ~ x, d, family = binomial(logit))
summary(ans)$coefficients
 Estimate Std. Error z value Pr(>|z|)
(Intercept) -82.219446 42.570413 -1.931375 0.05343665
x 2.022379 1.058386 1.910814 0.05602845

kojiro_i619さんのコメント
ご回答、ありがとうございます。Rのコードまで、提示して頂き、感謝します。検定にあたっての、回帰式を求める方法が、log(y/(1-y)=ax+bの式になると思いますが、yに比率があれば、式は求まると思いますが、1,0のデータだと、左辺が求まらないので、どうするのか、理解ができないわけです。ネットの中では、なにか適当にあたりをつけて、やっていくような感じに思えますが、そうなのでしょうか?whileのコードもあったので、最適値を、少しづつずらして求めるのでしょうか? 最小二乗法での解の式は無いのでしょうか?

dilutionistさんのコメント
式は確率を表したモデルであり、観測値を直接放り込んで使うことはしません。 パラメータ選択は少しずつ変化させて尤度が最大になるものを選ぶという形になります。 無理矢理0.5を足したり引いたりする「経験ロジット」というものを使う方法もあるようですがよく知りません。GLMで最尤法を使う方が一般的だと思います。 http://www.snap-tck.com/room04/c01/stat/stat10/stat1003.html#note01

kojiro_i619さんのコメント
色々、ありがとうございます。「GLMで最尤法を使う」の最尤法は、パラメータを少しずつ変えていくことでしょうか?a,bの2変数で、まずbを固定して、aの最小値をもとめ、aを固定して、bの最小値を求め,を繰り返すのでしょうか?

dilutionistさんのコメント
IRLS法ではa, bそれぞれを直接変化させるのではなく適当な初期値から尤度関数をもとに間接的にaとbのセットを求め、その値をもとにまた再計算して次のaとbのセットを...というような経過を辿ります。具体的にイメージするのは難しいですが。 [http://stat.biopapyrus.net/glm/mle-irls.html:title]

kojiro_i619さんのコメント
丁寧なご説明、ありがとうございます。数学の世界で、解はあるとおもったのですが、無いこともあるということでしょうか。

dilutionistさんのコメント
そうですね。私も数学に詳しいわけではないですが。 [http://masa-cbl.hatenadiary.jp/entry/20130429/1367247054:title]

kojiro_i619さんのコメント
だとすると、世紀の難問にならないのかなあ

kojiro_i619さんのコメント
くやしいですね。

kojiro_i619さんのコメント
平面での点を通る距離が最短な直線の式は、私、最小2乗法と偏微分で、半年がかり以上かで、求めました。その際の、PCで動く数式ソフト「数式展開」も作りました。どんな方法を使っても、解は求まらないのですね。不思議だ。

kojiro_i619さんのコメント
ちなみに、上の点について数学的には、http://takeno.iee.niit.ac.jp/~shige/math/lecture/graduate/corel1/node5.html
関連質問

●質問をもっと探す●



0.人力検索はてなトップ
8.このページを友達に紹介
9.このページの先頭へ
対応機種一覧
お問い合わせ
ヘルプ/お知らせ
ログイン
無料ユーザー登録
はてなトップ