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)
としますが、この２変数の分析をRで行う際に、xおよびyについて、どのような演算処理が行われるのでしょうか？わかりやすく、教えていただければ、幸いです。
みんなの回答
-
dilutionist
満足200pt
ロジスティック回帰は一般化線形モデル（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] 7glm()の出力と比較すると同じ係数と切片に収束していることがわかります。
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 -
dilutionistこの質問へのコメント
質問の情報
- 登録日時
- 2015-10-19 14:16:46
- 終了日時
- 2015-10-26 14:20:03
- 回答条件
- 1人2回まで
この質問のカテゴリ
この質問に含まれるキーワード
人気の質問
-
この英文どうなってますか？ The boy seemed to be waiting for his mother to come.（その少年は母親が来…
匿名質問者
3
1
-
不定詞の形容詞的用法らしいんですけど、副詞的用法と書いたら✖️ですか？ I want to buy a book to read …匿名質問者11
-
日向陽葵という○V女優が好きなのですが似ている○V女優はいますか匿名質問者
-
共通テストで使われるプログラミング言語は共通テスト特有の言語と聞いたのですが、受験勉強で身につけたプ…匿名質問者21
-
内定してほしいポケモンはいますか？匿名質問者