前回の記事では熱中症搬送人員数に対する日最高気温の回帰曲線を描いた．今回はポアソン分布に基づく搬送人員数の95%信頼区間を求める．
データのインポート
この段階はユーザーインターフェースでの作業の記録である．
> library(readxl)
> HeatStroke <- read_excel("C:/Users/***/HeatStroke.xlsx",
+ col_types = c("date", "numeric", "text",
+ "numeric", "numeric", "numeric"))
ポアソン回帰モデルによる一般化線形回帰
日最高気温をx, 搬送人員数をyとしてポアソン回帰モデルにて一般化線形回帰を行う．回帰係数が求まる．
> y = HeatStroke$HeatStroke
> x = HeatStroke$Temp
> z = glm(y ~ x, family = poisson(link = "log"))
> 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というデータフレームを作成し，Tempには１から31までの整数を格納する．Numには前段で得られた回帰係数をもとに関数を作成し，回帰式から予想される搬送人員数を格納している．さらにLowには95%信頼下限を，Uppには95%信頼上限を格納している．ポアソン分布では分散は平均値に等しいこと，サンプルサイズが十分に大きい場合には2標準偏差が約95％をカバーすることを利用している．
> p = data.frame(Temp = seq(10, 40, length.out = 31)) > p$Num = exp(-7.475 + 0.306 * p$Temp) > p$Low = p$Num - 2 * sqrt(mean(y)) > p$Upp = p$Num + 2 * sqrt(mean(y))
pを展開する．p$Lowには負の数値が現れるが，ポアソン分布では非負数であることから，0で置換する必要がある．
> p Temp Num Low Upp 1 10 0.01209455 -5.7050657 5.729255 2 11 0.01642419 -5.7007361 5.733584 3 12 0.02230376 -5.6948565 5.739464 4 13 0.03028811 -5.6868721 5.747448 5 14 0.04113072 -5.6760295 5.758291 6 15 0.05585479 -5.6613055 5.773015 7 16 0.07584982 -5.6413104 5.793010 8 17 0.10300271 -5.6141575 5.820163 9 18 0.13987585 -5.5772844 5.857036 10 19 0.18994894 -5.5272113 5.907109 11 20 0.25794729 -5.4592130 5.975108 12 21 0.35028786 -5.3668724 6.067448 13 22 0.47568472 -5.2414755 6.192845 14 23 0.64597143 -5.0711888 6.363132 15 24 0.87721777 -4.8399425 6.594378 16 25 1.19124622 -4.5259140 6.908406 17 26 1.61769128 -4.0994690 7.334852 18 27 2.19679614 -3.5203641 7.913956 19 28 2.98321029 -2.7339500 8.700371 20 29 4.05114679 -1.6660135 9.768307 21 30 5.50138567 -0.2157746 11.218546 22 31 7.47078440 1.7536241 13.187945 23 32 10.14519303 4.4280328 15.862353 24 33 13.77699263 8.0598324 19.494153 25 34 18.70891223 12.9917520 24.426072 26 35 25.40637178 19.6892115 31.123532 27 36 34.50140335 28.7842431 40.218564 28 37 46.85229530 41.1351350 52.569456 29 38 63.62458804 57.9074278 69.341748 30 39 86.40106482 80.6839046 92.118225 31 40 117.33111729 111.6139570 123.048278
ループと条件分岐により置換
nというカウンタ変数をおいてデータフレームにループでアクセスし，p$Lowが負の場合は0に置換している．
> for (n in seq(1:31)) {
+ if (p$Low[n] < 0) {p$Low[n] = 0}
+ }
pを展開する．意図通り0に置換されている．
> p Temp Num Low Upp 1 10 0.01209455 0.000000 5.729255 2 11 0.01642419 0.000000 5.733584 3 12 0.02230376 0.000000 5.739464 4 13 0.03028811 0.000000 5.747448 5 14 0.04113072 0.000000 5.758291 6 15 0.05585479 0.000000 5.773015 7 16 0.07584982 0.000000 5.793010 8 17 0.10300271 0.000000 5.820163 9 18 0.13987585 0.000000 5.857036 10 19 0.18994894 0.000000 5.907109 11 20 0.25794729 0.000000 5.975108 12 21 0.35028786 0.000000 6.067448 13 22 0.47568472 0.000000 6.192845 14 23 0.64597143 0.000000 6.363132 15 24 0.87721777 0.000000 6.594378 16 25 1.19124622 0.000000 6.908406 17 26 1.61769128 0.000000 7.334852 18 27 2.19679614 0.000000 7.913956 19 28 2.98321029 0.000000 8.700371 20 29 4.05114679 0.000000 9.768307 21 30 5.50138567 0.000000 11.218546 22 31 7.47078440 1.753624 13.187945 23 32 10.14519303 4.428033 15.862353 24 33 13.77699263 8.059832 19.494153 25 34 18.70891223 12.991752 24.426072 26 35 25.40637178 19.689212 31.123532 27 36 34.50140335 28.784243 40.218564 28 37 46.85229530 41.135135 52.569456 29 38 63.62458804 57.907428 69.341748 30 39 86.40106482 80.683905 92.118225 31 40 117.33111729 111.613957 123.048278
ggplot2による描画
> w = ggplot(data = p) + geom_ribbon(aes(x = p$Temp, ymin = p$Low, ymax = p$Upp), alpha = 0.2, fill = "red") > w + geom_line(data = p, aes(x = p$Temp, y = p$Num), col = "red")
まとめ
回帰曲線の95%信頼区間をRで求め，ggplot2で描いた．もとの散布図とのレイヤー合成ができておらず，課題である．