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について、どのような演算処理が行われるのでしょうか?わかりやすく、教えていただければ、幸いです。

回答の条件
  • 1人2回まで
  • 13歳以上
  • 登録:2015/10/19 14:16:46
  • 終了:2015/10/26 14:20:03

回答(1件)

id:dilutionist No.1

dilutionist回答回数154ベストアンサー獲得回数512015/10/20 18:59:05

ポイント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] 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
他8件のコメントを見る
id:kojiro_i619

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

2015/10/28 16:58:10
id:kojiro_i619

ちなみに、上の点について数学的には、http://takeno.iee.niit.ac.jp/~shige/math/lecture/graduate/corel1/node5.html

2015/10/29 15:11:03

コメントはまだありません

この質問への反応(ブックマークコメント)

「あの人に答えてほしい」「この質問はあの人が答えられそう」というときに、回答リクエストを送ってみてましょう。

これ以上回答リクエストを送信することはできません。制限について

絞り込み :
はてなココの「ともだち」を表示します。
回答リクエストを送信したユーザーはいません