こちらのサイトポワソン回帰を拝見させていただき、Rで実装しながら勉強しています。
等価な方法「その3」として、以下が紹介されています。これは、link関数がidentityの場合ですが、link関数がlogの場合も可能なのでしょうか。
x = c(1,2,3,4)
y = c(2,3,5,4)
w = c(1,1,1,1)
for (i in 1:10) {
r = lm(y ~ x, weights=w)
lambda = predict(r)
print(c(as.numeric(r$coef), -sum(y*log(lambda)-lambda)))
w = 1 / lambda
}
#省略
#[1] 1.2783467 0.8886613 -4.0609501
glm(y~x,family=poisson(link="identity"))
#(Intercept) x
# 1.2784 0.8887
以下、link関数をlogを想定して、単純にyをlog(y)に、wをlambdaにしてみましたが、違うのでしょうか。どなたかアドバイスをお願いしてもよろしいでしょうか。
x = c(1,2,3,4)
y = c(2,3,5,4)
w = c(1,1,1,1)
for (i in 1:10) {
r = lm(log(y) ~ x, weights=w)
lambda = predict(r)
print(c(as.numeric(r$coef), -sum(y*log(lambda)-lambda)))
w = lambda
}
#省略
#[1] 0.6195090 0.2334289 1.7332280
glm(y~x,family=poisson(link="log"))
#(Intercept) x
# 0.6393 0.2320
###################### 追加投稿 20200619 ######################
user12399423さん(回答日時: 6月13日13:19)
ご回答いただき、ありがとうございました。分かりやすい説明で、ずっと疑問に思っていたことが分かり、すっきりしました。勉強不足でお恥ずかしいですが、もう1点、同じようなことで分からないことがあり、もしまだこの質問が生きていたら、アドバイスを頂けると幸いです。
こちらのサイトIRLS で解く最尤推定量も拝見させていただき、勉強していました。
以下の加重最小二乗法では、link関数がlogの場合です。例えば、こちらのプログラムでidentityの場合(単純にlambdaの項からexpを削除した場合)は、可能なのでしょうか。実行してみると、うまくいきません(回帰分析と同じになる?)。なぜ、この場合は、link関数をlogにするとポワソン回帰になるのでしょうか。
混乱してしまいました。どなたかアドバイスをお願いしてもよろしいでしょうか。
#log
x<-c(1,2,3,4)
X<-cbind(rep(1, length(x)),x)
Y<-c(2,3,5,4)
beta <- matrix(0, ncol = 2, nrow = 100)
beta[1, ] <- c(1, 10)
for (m in 2:100) {
lambda <- exp(X %*% beta[m - 1, ])
W <- diag(lambda[, 1])
XtWX <- t(X) %*% W %*% X
U <- t(Y - lambda) %*% X
beta[m, ] <- beta[m - 1, ] + solve(XtWX) %*% t(U)
}
tail(beta)
# [100,] 0.6392647 0.2320399#OK
これは、どうなのでしょうか。。
#identity?
x<-c(1,2,3,4)
X<-cbind(rep(1, length(x)),x)
Y<-c(2,3,5,4)
beta <- matrix(0, ncol = 2, nrow = 100)
beta[1, ] <- c(1, 10)
for (m in 2:100) {
lambda <- (X %*% beta[m - 1, ])
W <- diag(lambda[, 1])
XtWX <- t(X) %*% W %*% X
U <- t(Y - lambda) %*% X
beta[m, ] <- beta[m - 1, ] + solve(XtWX) %*% t(U)
}
tail(beta)
# [100,] 1.499999 0.8000007
glmを実行する。
#glm
summary(glm(y~x))
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 1.5000 1.1619 1.291 0.326
#x 0.8000 0.4243 1.886 0.200