ポアソン回帰

余所のサイトで勉強させて頂きました。

# pois regression

library(dplyr)

Y <- c( 6,  9, 11, 11,
       21, 14, 13, 20,
       38, 36, 23, 29,
       30, 33, 46, 40,
       61, 48, 54, 44)

X <- cbind(rep(1, length(Y)),
       c(1, 1, 1, 1,
         2, 2, 2, 2,
         3, 3, 3, 3,
         4, 4, 4, 4, 
         5, 5, 5, 5)
       )

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)
}

params=beta[nrow(beta),]

predict=exp(X%*%params)

plot(c(1:length(Y)),(Y),type="p",col=2,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))

par(new=T)

plot(c(1:length(Y)),predict,type="p",col=3,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))

data=read.csv("~/診療サンプル.csv")

library(dplyr)

Y <- as.numeric(data$発生件数)

X <- cbind(rep(1, length(Y)),
as.numeric(data$診療科),as.numeric(data$処方薬剤数)
       )

beta <- matrix(0, ncol = ncol(X), nrow = 100)

beta[1, ] <- rep(10,ncol(X))

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)
}

params=beta[nrow(beta),]

predict=exp(X%*%params)

plot(c(1:length(Y)),(Y),type="p",col=2,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))

par(new=T)

plot(c(1:length(Y)),predict,type="p",col=3,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))

データ：