R Possion Regression Analysis (ポアソン回帰)
・ポアソン分布の前提3つ
-時間を細かく分けていくと、各時間帯で発生するイベントは1回だけである
-細かく分けた時間帯でのイベントの発生確率は同じである
-他の時間帯のイベントの発生状況の影響は受けない
-----------------------------------------------------------------------------------------------------------------
P(Y=k) = λ^k / k! × e^-λ (k=0,1,2…)
λは単位時間あたりのイベントの平均発生回数
単位時間あたりのイベントの発生関数の分散もλ
-----------------------------------------------------------------------------------------------------------------
ポアソン分布をパラメータを変えて描画
X <- NULL
K <- c(0.5, 1, 2, 5, 10)
col <- c("Red", "blue", "green", "black", "yellow")
for (j in K) {
for (i in 1:8) {
X[i] <- j^i / gamma(i+1) * exp(-j)
}
par(new=T)
if (j == 10) {
plot(1:8, X[1:8], type="l", col=c(j), ylim=c(0,0.6),xlim=c(1,8))
} else {
plot(1:8, X[1:8], type="l", col=c(j), ylim=c(0,0.6),xlim=c(1,8),ann=FALSE)
}
}
#λが増加すると、正規分布に近づく
t <- c(0:5)
x <- c(8,6,8,3,4,1)
m <- sum(x)
s <- t %*% x
r <- s / m
print(c(s, r))
ci <- c(r-2 * sqrt(r/m), r+2*sqrt(r/m))
x <- c(8,6,8,3,4,1)
names(x) <- c(0:5)
y <- dpois(c(0:5), lambda = r) * sum(x)
lines(barplot(x, ylim=c(0,10), col=rainbow(6)), y, type="l")
#Poisson分布への適合度をgoodfitで測定
#準備として、データをmaxrixに修正
x <- matrix(c(8,6,8,3,4,1,c(0:5)),nrow=6)
print(x)
library(vcd)
result <- goodfit(x, type = "poisson", method = "MinChisq")
#Chisqを最小化するようにpoisson式を当てはめ
print(result)
summary(result)
result$par
plot(result)
#日付毎の件数を入力する場合には、ベクトル形式として入力して
#goodfitを利用可
x <- c(0,0,0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,4,4,4,4,5)
library(vcd)
result <- goodfit(x, type="poisson", method="MinChisq")
plot(result)
############Poisson Regression Analysis############
boxplot(count ~ spray, data = InsectSprays)
#グラフを見ると、平均値が大きいほど、分散が大きいという
result <- glm(count ~ spray, data=InsectSprays, family ="poisson")
#logλ = a + b1x1 + b2x2 + … + bnxnとして推定値(a,b1,…,bn)
#を計算するモデル。名義変数の場合は、1変数を基準として、
#他の変数との差をパラメータb1, b2, …, bnとして推定する。
summary(result)
exp(result$coefficients)
cbind(InsectSprays, SurviveRate = predict(TN4, Titanic, type ="response"))
############オフセットによる調整法###############
accident <- c(45703,8906,7938,12091,7327,9820,11526,6525)
pref <- as.factor(c('A','B','C','D','E','F','G','H'))
result1 <- glm(accident~pref,family ="poisson")
result1
exp(result1$coefficients)
#A=福岡県,B=佐賀県,C=長崎県,D=熊本県,E=大分県,F=宮崎県,
#G=鹿児島県,H=沖縄県
#Aは人口が多いために交通事故が多いことを考慮すべき
pop <- c(5056, 859, 1453, 1828, 1203, 1143, 1730, 1373)
#各県の人口の対数
result2 <- glm(accident ~ pref, offset = log(pop), family = "poisson")
result2
exp(result2$coefficients)
#佐賀の方が交通事故の頻度が多いことが分かる
#offsetを加えてもAICは不変でるあることに注意
############過分散である場合の解析方法###############
#ポアソン回帰分析においては、過分散が問題となることがある。
#n個の回数データがあり、各回数に対するλは一定ではなく、
#ガンマ分布に沿って変動するものと考える。
#f(λ) = λ^a-1 * exp(-λ/b) /Γ(a)*b^a aをshapeパラメータ,bをscaleパラメータという
#λがガンマ分布に従い、yが平均λのポアソン分布に従うとすると、
#λの変動も含めたyの平均はab,分散はab(1+b)
set.seed(1)
n <- 1000
lambda <- rgamma(n, shape=5, scale=3) #
y <- rpois(n, lambda)
mean(y)
var(y)
#Γ分布のプロット(on progess)
X <- NULL
K <- rbind(c(1,1), c(2,3), c(3,2), c(4,5), c(5,4))
col <- c("Red", "blue", "green", "black", "yellow")
for (l in 1:5) {
j = K[l,]
for (i in 1:8) {
X[i] <- i^(j[1]-1) * exp(-i/j[2]) / gamma(j[1])*j[2]^j[1]
}
par(new=T)
if (l == 5) {
plot(1:8, X[1:8], type="l", col=c(j), ylim=c(0,0.6),xlim=c(1,8))
} else {
plot(1:8, X[1:8], type="l", col=c(j), ylim=c(0,0.6),xlim=c(1,8),ann=FALSE)
}
}
#容量反応関係を表すデータの作成
set.seed(2)
n = 20
v = 3
dose <- c(0:4) * 5
x <- rep(dose, n)
a <- exp(3 + 0.05 * dose)/v
lambda <- rgamma(5 * n, shape = a, scale =v)
y <- rpois(5 *n, lambda)
dataNB <- data.frame(x,y)
boxplot(y ~ x, data = dataNB)
by(y, x ,mean)
by(y, x, var)
#上記のような、分散が平均に比べて大きなデータは、
#MASSパッケージのglm.nbを使う
glm.nb(y ~ x, data=dataNB)
##################練習問題######################
Y <- c(314,39,78,105,77,93,140,69)
X <- as.factor(c('福岡','佐賀','長崎','熊本','大分','宮崎','鹿児島','沖縄'))
result1 <- glm(formula = Y ~ X, family = 'poisson')
exp(result1$coefficients)
result2 <- glm(formula = Y ~ X, offset = log(pop), family = 'poisson')
exp(result2$coefficients)
