

こちらのサイト[ポワソン回帰](https://oku.edu.mie-u.ac.jp/%7Eokumura/stat/poisson_regression.html)を拝見させていただき、Rで実装しながら勉強しています。

等価な方法「その3」として、以下が紹介されています。これは、link関数がidentityの場合ですが、link関数がlogの場合も可能なのでしょうか。

[code]

    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  
    
[/code]

以下、link関数をlogを想定して、単純にyをlog(y)に、wをlambdaにしてみましたが、違うのでしょうか。どなたかアドバイスをお願いしてもよろしいでしょうか。

[code]

    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 
    
[/code]

###################### 追加投稿 20200619 ######################

user12399423さん（回答日時: 6月13日13:19）  
ご回答いただき、ありがとうございました。分かりやすい説明で、ずっと疑問に思っていたことが分かり、すっきりしました。勉強不足でお恥ずかしいですが、もう1点、同じようなことで分からないことがあり、もしまだこの質問が生きていたら、アドバイスを頂けると幸いです。  
こちらのサイト[IRLS で解く最尤推定量](https://stats.biopapyrus.jp/glm/mle-
irls.html)も拝見させていただき、勉強していました。

以下の加重最小二乗法では、link関数がlogの場合です。例えば、こちらのプログラムでidentityの場合（単純にlambdaの項からexpを削除した場合）は、可能なのでしょうか。実行してみると、うまくいきません（回帰分析と同じになる？）。なぜ、この場合は、link関数をlogにするとポワソン回帰になるのでしょうか。  
混乱してしまいました。どなたかアドバイスをお願いしてもよろしいでしょうか。

[code]

    #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
    
[/code]

これは、どうなのでしょうか。。

[code]

    #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
    
[/code]

glmを実行する。

[code]

    #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
    
[/code]


CC BY-SA 4.0
2020-05-23T18:51:04.393
66928
2020-06-18T19:22:16.580
2020-06-18T19:22:16.580
37329
37329
question
r
最小二乗法の重付けでポワソン回帰を確認できるのでしょうか（Rプログラム）


違います. 元の方法はリンク関数が恒等写像のときに y = lambda = a * x + b になることを利用しているので,
リンク関数を対数に設定すると成り立ちません(その1, その2も良い近似にならないと思います). そしてお考えの方法では極値で成り立つのが
log(lambda) * (log(y) - 1) * lambda' * exp(lambda) = 0 となり,
ポアソン回帰の尤度の1階条件とはかけ離れてしまっています. ちなみに `lambda = exp(predict(r))` とすると比較的近い値が出ますが,
多分 lambda の値があまり変化しないことによる偶然だと思います


CC BY-SA 4.0
2020-06-13T13:19:49.647
67630
2020-06-13T13:19:49.647
40575
66928
answer