ポアソン回帰モデルの回帰曲線をRで描く

　前回の記事では熱中症搬送人員数に対する日最高気温と平均水蒸気圧の影響をポアソン回帰モデルまたは負の二項分布モデルを用いて回帰係数を求めた．今回はポアソン回帰モデルの回帰曲線をRで描く．

データを準備する

　下記コードのようにHeatStroke.xlsxをインポートする．

> library(readxl)
> HeatStroke <- read_excel("C:/Users/****/HeatStroke.xlsx", 
+     col_types = c("date", "numeric", "text", 
+         "numeric", "numeric", "numeric"))

　x, y にそれぞれ日最高気温(Temp), 平均水蒸気圧(Vapor)を格納する．

> y = HeatStroke$HeatStroke
> x = HeatStroke$Temp

　zにglm()関数を用いてポアソン分布を指定する．

> z = glm(y ~ x, family = poisson(link = "log"))

　zを展開して回帰係数を求める．

> z
Call:  glm(formula = y ~ x, family = poisson(link = "log"))
Coefficients:
(Intercept)            x  
     -7.475        0.306  
Degrees of Freedom: 82784 Total (i.e. Null);  82783 Residual
Null Deviance:	    1464000 
Residual Deviance: 698400 	AIC: 915500

データフレームにベクトルを格納する

　データフレームpにx軸のベクトルを格納する．

> p = data.frame(Temp = seq(0, 40, length.out = 100))

　式 y = exp(-7.475 + 0.306*x)の形になるようにy軸のベクトルをデータフレームに格納する．

> p$Num = exp(-7.475 + 0.306 * p$Temp)

プロットする

　下記コードでプロットする．

> plot(x, y) + lines(p$Temp, p$Num, type = "l", col = "red")

結果

　結果は以下の通りである．

ggplot2で描く

　下記コード参照．凡例の記述に問題がある．

> library(ggplot2)
> w = ggplot(data = HeatStroke, aes(x = x, y = y)) + geom_point(alpha = 0.1)
> w + geom_line(data = p, aes(x = p$Temp, y = p$Num, col = "red"))

まとめ

　散布図にポアソン回帰モデルの回帰曲線を追加した．