Zenn
👻

ポアソン回帰のスクラッチ実装 - IWLSとFisher Scoring/Netwon Raphson

2024/05/06に公開
R
統計
glm
ポアソン回帰
tech

ゴール

  • ポアソン回帰のスクラッチ実装を理解する

やること

  • IWLSを実装する
  • Fisher scoring (Canonicalリンク関数を用いるのでNewton Raphsonと同じ)を実装する

関連記事

コード

データ生成

こちらはModigied Poissonの記事のn=50のデータと同じものになる。

# データを生成する関数 
datagen <-
  function(n) {
    N <- n # サンプルサイズ
    
    v1 <- rnorm(N, mean = 60, 10)
    v2 <- rbinom(N, size = 1, prob = 0.5)
    v3 <- rbinom(N, size = 1, prob = 0.5)
    v4 <- rbinom(N, size = 1, prob = 0.5)
    X <- do.call(cbind, list(1,v1,v2,v3,v4))
    colnames(X) <- c('const','v1', 'v2', 'v3', 'v4')
    P <- ncol(X) # 説明変数の数（定数項含む）
    
    true_beta <- c(-2.2, 0, 0.8, 0, 0) # betaの真値
    Y <- unlist(lapply(exp(X%*%true_beta), rbinom, n = 1, size = 1))
    
    X[,2] <- scale(X[,2])
    mydata <- data.frame(X[,-1],Y)
    
    return(list(X=X,mydata=mydata))
  }

# generate data
n <- 50
d <- datagen(n)
mydata <- d$mydata
X <- d$X

普通のポアソン回帰の実行

# Poisson regression model
myposfit <- glm(Y~v1+v2+v3+v4, data = mydata, family = poisson())

> summary(myposfit)$coefficients
               Estimate  Std. Error      z value   Pr(>|z|)
(Intercept) -21.9963357 4586.287488 -0.004796109 0.99617327
v1           -1.4633147    0.773754 -1.891188526 0.05859918
v2           18.5646941 4586.287385  0.004047870 0.99677028
v3           -0.5604871    1.444061 -0.388132442 0.69791802
v4            1.1496874    1.455447  0.789920475 0.42957421

QR分解で共分散行列を計算する関数の作成

vcov_qr <- function(mymat) {
  
  myqr <- qr(mymat, tol = 1e-10)
  R <- qr.R(myqr)
  Q <- qr.Q(myqr)
  
  inv_R <- backsolve(R, diag(ncol(R)))
  return(inv_R %*% t(inv_R))
}

IWLS

irwls_fun <-
  function(X,b_int,y,Iter) {
    
    N <- nrow(X)
    l <- matrix(NA, ncol=1, nrow=Iter+1)
    b <- matrix(NA, ncol=ncol(X), nrow=Iter+1)
    b[1,] <- b_int
    
    # log likelihood
    l[1] <- t(y)%*%X%*%b[1,] - sum(exp(X%*%b[1,]))
    
    for (i in 1:Iter){
      u <- exp(X%*%b[i,]) %>% as.vector()
      W <- diag(u, ncol = N)
      
      # Information matrix
      J_n_inv <- vcov_qr(sqrt(W)%*%X)
      
      # z
      z <- X%*%b[i,] + (y-u)*u^-1
      
      # b
      b[i+1,] <- J_n_inv%*%t(X)%*%W%*%z %>% t()
      
      # log likelihood
      l[i+1] <- t(y)%*%X%*%b[i+1,] - sum(exp(X%*%b[i+1,]))
      
      # same as glm.control with default epsilon
      l_diff <- abs(-2*l[i+1] - (-2*l[i]))/(-2*l[i+1] + 0.1) < 1e-8
      
      if (l_diff) {
        break
      }
    }
    
    res <- cbind(seq(0,Iter,1),b,l)
    colnames(res) <- c("Iter","b0","b1","b2","b3","b4","loglik")
    return(list(
      coef = res %>% 
        data.frame() %>%
        filter(!is.na(loglik)),
      vcov = J_n_inv)
    )
    
  }

out <-  irwls_fun(X, rep(0,ncol(X)), y = mydata$Y, Iter = 100) 

res <- 
out$coef %>%
  slice(n()) %>%
  select(-Iter, -loglik)

> data.frame(est = unlist(res),
+            se = sqrt(diag(out$vcov))
+            )
           est          se
b0 -21.8623896 4289.189577
b1  -1.4633147    0.773754
b2  18.4307480 4289.189467
b3  -0.5604871    1.444061
b4   1.1496874    1.455447

Fisher scoring

fisherscoring_fun <-
  function(X,b_int,y,Iter) {
    
    N <- nrow(X)
    l <- matrix(NA, ncol=1, nrow=Iter+1)
    b <- matrix(NA, ncol=ncol(X), nrow=Iter+1)
    b[1,] <- b_int
    
    # log likelihood
    l[1] <- t(y)%*%X%*%b[1,] - sum(exp(X%*%b[1,]))
    
    for (i in 1:Iter){
      u <- exp(X%*%b[i,]) %>% as.vector()
      W <- diag(u, ncol = N)
      
      # Observed Information matrix
      J_n_inv <- vcov_qr(sqrt(W)%*%X)
      
      # Observed score
      observed_score <- X*(y-u)
      observed_score_sum <- apply(observed_score, 2, sum)
      
      # b; the negative sign for Hessian is necessary
      b[i+1,] <- b[i,] - (-J_n_inv)%*%observed_score_sum %>% t()
      
      # log likelihood
      l[i+1] <- t(y)%*%X%*%b[i+1,] - sum(exp(X%*%b[i+1,]))
      
      # same as glm.control with default epsilon
      l_diff <- abs(-2*l[i+1] - (-2*l[i]))/(-2*l[i+1] + 0.1) < 1e-8
      
      if (l_diff) {
        break
      }
    }
    
    res <- cbind(seq(0,Iter,1),b,l)
    colnames(res) <- c("Iter","b0","b1","b2","b3","b4","loglik")
    return(list(
      coef = res %>% 
             data.frame() %>%
             filter(!is.na(loglik)),
      vcov = J_n_inv)
      
    )
    
  }


out2 <-  fisherscoring_fun(X, rep(0,ncol(X)), y = mydata$Y, Iter = 100)

res2 <-
out2$coef %>%
  slice(n()) %>%
  select(-Iter, -loglik)

> data.frame(est = unlist(res2),
+            se = sqrt(diag(out2$vcov))
+ )
           est          se
b0 -21.8623895 4289.189601
b1  -1.4633147    0.773754
b2  18.4307480 4289.189491
b3  -0.5604871    1.444061
b4   1.1496874    1.455447

感想

  • IWLSとFisher scoringの結果はほぼ一致した。
  • 自前のアルゴリズムでもglm()と推定値は近い値を得られた。b0とb2の標準誤差は2倍近くの大きさになってしまった。
    • 収束判定をglm.controlと同じにしたところ標準誤差もかなり近しくなった。
  • n=50で二項分布からYを生成したデータという特殊例なので一旦はこれでOKとしたいが、SEについても自前の2つはかなり近しいので、glm()との差分がどこから生まれるかは気になるところではある。

Discussion