厚労省のオープンデータを利用
例えばn次曲線グラフのが急に変化する場合、その動きを確認（察知）するには、シンプルな（回帰）直線の勾配（傾き）を観察する方がわかりやすいのでは？、と考えました。
一定の期間ごとに分割し、その期間毎に回帰分析を実行することを試してみることにしました。
新型コロナウイルスの感染者数の推移（全国）のデータを使います。
データは厚生労働省のサイトから取得しました。
データ取得先はこちらです。
Rで取得します。
covid <- read.csv("https://covid19.mhlw.go.jp/public/opendata/newly_confirmed_cases_daily.csv")
head(covid)
Date Prefecture Newly.confirmed.cases
1 2020/1/26 ALL 1
2 2020/1/26 Hokkaido 0
3 2020/1/26 Aomori 0
4 2020/1/26 Iwate 0
5 2020/1/26 Miyagi 0
6 2020/1/26 Akita 0
列名を変更します。
colnames(covid)[2]<- "Pref"
colnames(covid)[3] <- "Cases"
colnames(covid)
[1] "Date" "Pref" "Cases"
group_byとsummariseを使用
少し横道に逸れますが、「R」の復習がてら、県ごとにグルーピングし、県ごとのCasesの合計値を求めてみます。
まずグルーピングです。
その前に（重複計算するので）Pref=ALLを外しておきます。
library(dplyr)
covid2 <- covid %>% filter(Pref != "ALL")
head(covid2)
Date Pref Cases
1 2020/1/26 Hokkaido 0
2 2020/1/26 Aomori 0
3 2020/1/26 Iwate 0
4 2020/1/26 Miyagi 0
5 2020/1/26 Akita 0
6 2020/1/26 Yamagata 0
groupby <- covid2 %>% group_by(Pref)
head(groupby)
# A tibble: 6 x 3
# Groups: Pref [6]
Date Pref Cases
<fct> <fct> <int>
1 2020/1/26 Hokkaido 0
2 2020/1/26 Aomori 0
3 2020/1/26 Iwate 0
4 2020/1/26 Miyagi 0
5 2020/1/26 Akita 0
6 2020/1/26 Yamagata 0
Casesの和を求めます。
summarise <- groupby %>%
summarise(DailyCases = sum(Cases))
summarise %>% arrange(DailyCases) %>%head(10)
# A tibble: 10 x 2
Pref DailyCases
<fct> <int>
1 Shimane 924
2 Tottori 1109
3 Akita 1185
4 Tokushima 1942
5 Fukui 2015
6 Kochi 2318
7 Iwate 2428
8 Yamagata 2553
9 Toyama 2969
10 Kagawa 2985
目的の日単位集計を実行します。
DailyData <- covid2 %>% group_by(Date) %>%
summarise(DailyCases = sum(Cases))
DailyData %>% arrange(DailyCases) %>%head(10)
DailyData %>% arrange(DailyCases) %>%tail(10)
# A tibble: 10 x 2
Date DailyCases
<fct> <int>
1 2021/8/4 14204
2 2021/8/8 14457
3 2021/8/5 15248
4 2021/8/6 15624
5 2021/8/7 15740
6 2021/8/11 15771
7 2021/8/15 17829
8 2021/8/12 18893
9 2021/8/14 20133
10 2021/8/13 2035
Factor型をDate型、Numeric型に変換
回帰分析をかけるために、x軸（Date）を数字に変更し、別の変数に収納します。
その前に、DateがFactor型だったので日付型に変換しています。
変換後、数字型に変えます。
数字の型でないと、回帰分析がかけられないからです。
（方法があるのかもしれませんが、今回は調べられませんでした）
str(DailyData)
Classes ‘tbl_df’, ‘tbl’ and 'data.frame': 568 obs. of 2 variables:
$ Date : Factor w/ 565 levels "2020/1/26","2020/1/27",..: 1 2 3 4 5 6 7 8 9 10 ...
$ DailyCases: int 1 0 3 1 3 0 619 666 431 272 ...
DailyData$Date <- as.Date(DailyData$Date)
str(DailyData)
Classes ‘tbl_df’, ‘tbl’ and 'data.frame': 568 obs. of 2 variables:
$ Date : Date, format: "2020-01-26" "2020-01-27" ...
$ DailyCases: int 1 0 3 1 3 0 619 666 431 272 ...
DailyDataN <- as.numeric(DailyData$Date)
str(DailyDataN)
num [1:568] 18287 18288 18289 18290 18291 ...
DailyDataY <- DailyData$DailyCases
str(DailyDataY)
int [1:568] 1 0 3 1 3 0 619 666 431 272 ...
回帰分析を実行
上記で収納したx軸要素（DailyDataN）とy軸要素（DailyDataY）を使い、回帰分析をかけます。
散布図と近似直線を描きます。
plot(DailyDataN, DailyDataY, cex=0.8, col="grey") #図
DailyDataR <- lm(DailyDataY~DailyDataN)
summary(DailyDataR)
abline(DailyDataR, col="grey", lty=2) #図
今回の目的は、x軸を細分割して、回帰分析をかけ、分割した範囲ごとに近似直線を描くことです。
もう少しデータを整えます。
元のデータフレーム （DailyData）に数字にした日付の列（DailyDataN）を加えます。
曜日列、曜日番号（色の指定に使用）も加えます。
# 数字にした日付の列追加
DailyData2 <- transform(DailyData, Datenum=DailyDataN)
head(DailyData2)
Date DailyCases Datenum
1 2020-01-26 1 18287
2 2020-01-27 0 18288
3 2020-01-28 3 18289
4 2020-01-29 1 18290
5 2020-01-30 3 18291
6 2020-01-31 0 18292
# 曜日列、曜日番号を追加
DailyData3 <- transform(DailyData2, Weekday=weekdays(DailyData2[,1]))
head(DailyData3)
library(lubridate)
DailyData4 <- transform(DailyData3, Weeknum=wday(DailyData2[,1]))
head(DailyData4)
Date DailyCases Datenum Weekday Weeknum
1 2020-01-26 1 18287 日曜日 1
2 2020-01-27 0 18288 月曜日 2
3 2020-01-28 3 18289 火曜日 3
4 2020-01-29 1 18290 水曜日 4
5 2020-01-30 3 18291 木曜日 5
6 2020-01-31 0 18292 金曜日 6
月曜のプロットに色を付ける
月曜のプロットに色（紫：6）を付るように列を追加します。
DailyData4$Colmon <- ifelse(DailyData4$Weeknum == 2, 6, 8) #6は紫、8は灰色
DailyData4 %>% head(20)
Date DailyCases Datenum Weekday Weeknum Colmon
1 2020-01-26 1 18287 日曜日 1 8
2 2020-01-27 0 18288 月曜日 2 6
3 2020-01-28 3 18289 火曜日 3 8
4 2020-01-29 1 18290 水曜日 4 8
5 2020-01-30 3 18291 木曜日 5 8
6 2020-01-31 0 18292 金曜日 6 8
7 2020-10-01 619 18536 木曜日 5 8
8 2020-10-10 666 18545 土曜日 7 8
9 2020-10-11 431 18546 日曜日 1 8
10 2020-10-12 272 18547 月曜日 2 6
念のため、Datenumを使って日付順にソートします。
DailyData5 <- DailyData4[order(DailyData4$Datenum),]
DailyData5 %>% head(20)
DailyData5 %>% tail(20)
Date DailyCases Datenum Weekday Weeknum Colmon
542 2021-07-27 7619 18835 火曜日 3 8
543 2021-07-28 9568 18836 水曜日 4 8
544 2021-07-29 10687 18837 木曜日 5 8
546 2021-07-30 10728 18838 金曜日 6 8
547 2021-07-31 12328 18839 土曜日 7 8
554 2021-08-01 10160 18840 日曜日 1 8
561 2021-08-02 8325 18841 月曜日 2 6
562 2021-08-03 12061 18842 火曜日 3 8
563 2021-08-04 14204 18843 水曜日 4 8
564 2021-08-05 15248 18844 木曜日 5 8
565 2021-08-06 15624 18845 金曜日 6 8
566 2021-08-07 15740 18846 土曜日 7 8
567 2021-08-08 14457 18847 日曜日 1 8
568 2021-08-09 12054 18848 月曜日 2 6
555 2021-08-10 10566 18849 火曜日 3 8
556 2021-08-11 15771 18850 水曜日 4 8
557 2021-08-12 18893 18851 木曜日 5 8
558 2021-08-13 20357 18852 金曜日 6 8
559 2021-08-14 20133 18853 土曜日 7 8
560 2021-08-15 17829 18854 日曜日 1 8
x、y軸を描かずプロット
散布図と全体の近似直線も描き直しておきます。
この後グラフの上書きを行うので、ｘ軸は表示させないように、y軸は指定しておきます。
（目盛りが上書きされると見にくくなったりするので）
dev.new() #は描画別ウインドウ開く
par(xaxt="n") #x軸を書かない
par(yaxt="n") #y軸を書かない
par(ann=F) #軸タイトルを書かない
DailyDataN <- DailyData5$Datenum
DailyDataY <- DailyData5$DailyCases
min <- min(DailyDataN)
max <- max(DailyDataN))
maxy <- max(DailyDataN)
# plot(DailyDataN, DailyDataY, xlim=c(min(DailyDataN),max(DailyDataN)), ylim=c(0,max(DailyDataN)), cex=0.8, col="grey", pch=19) #図
plot(DailyDataN, DailyDataY, xlim=c(min,max), ylim=c(0,maxy), cex=0.8, col="grey", pch=19) #図
DailyDataR <- lm(DailyDataY~DailyDataN)
abline(DailyDataR, col="grey", lty=2) #図
散布図の月曜のプロットの色（6：紫）を変えておきます。
DailyData5Mon <- DailyData5[DailyData5$Colmon== 6,]
head(DailyData5Mon)
tail(DailyData5Mon)
# par(new=T)でplot重ね書き
par(new=T)
par(yaxt="s")
plot(DailyData5Mon[,1], DailyData5Mon[,2], xlim=c(min,max), ylim=c(0,maxy), col=DailyData5Mon$Colmon, pch=19) #図
whileを使いx軸を分割して回帰分析
この散布図に期間を分割した近似曲線を描きます。
集計範囲（数字にした日付）を確認します。
# 集計データ開始日：2020-01-26（日曜日）
s <- min(DailyDataN)
s
[1] 18287
# 最新のデータ集計日：2020-08-15（日曜日）
max <- max(DailyDataN)
max
1] 18854
「s」から「max」の期間を7日間ずつ分割します。
各期間の始点を月曜日にしたいので、上記sの初期値を以下に設定します。
「while」を使って期間分割を進めつつ、回帰分析も実行し、実行毎に近似直線を描きます。
今回の軸の範囲は2020-01-27（月曜日）〜2020-08-15（日曜日）と週単位で分割できます。
### （ここから）
s <- 18287+1 #2020-01-27
col <- 2 #赤：色の初期値として設定
while (s < max) {
# 7日間ずつ刻む
DailyData5pks1 <- DailyData5[DailyData5$Datenum>=s,] #
s <- s + 6
DailyData5pks2<- DailyData5pks1[DailyData5pks1$Datenum<=s,] #
s <- s + 1
# 上記でピックアップした7日間のx軸値とy軸値を抽出
DailyData5pks2x <- DailyData5pks2$Datenum #
print(DailyData5pks2x) #確認
DailyData5pks2y <- DailyData5pks2$DailyCases #
print(DailyData5pks2y) #確認
while (col > 4) { #色番号2,3,4：近似直線に赤，緑，青を順番に使うための処理
col = 2
}
# 回帰分析を実行し、近似直線を描く
ResultR <- lm(DailyData5pks2y~DailyData5pks2x) #
# lines(DailyData5pks2x , predict(ResultR), lty = "dashed", lwd = 2, col="grey") #
lines(DailyData5pks2x , predict(ResultR), lty = "dotted", lwd = 2.2, col=col, xlim=c(min,max), ylim=c(0,maxy)) #図
col = col + 1 #近似曲線の始まりを重ならせるための処理
}
### （ここまで）
ｘ軸を書き込む
x軸を書き込みます。
# x軸に時系列目盛りを追加にする。「las=2」でx軸の目盛りを縦に。最初の「1」が下に配置。
par(family = "HiraKakuProN-W3") #Macでの日本語の文字化け回避
par(xaxt="s") #軸を書く
Date2 <- DailyData5$Date
axis.Date(1,at=seq(min(Date2),max(Date2), "2 week"),format="%m月%d日", las=2, col="grey", col.axis="grey", cex.axis=0.8) #図（２週間単位）
# axis.Date(1,at=seq(min(Date2),max(Date2), "month"),format="%y年%m月", las=2, col="black", col.axis="black") #図（月単位）
描画期間を変更する
描画範囲を2021年に絞ってみました。
21年の最初の月曜が「2021-01-04」（18631）でした。
DailyData5Mon %>% tail(35)
Date DailyCases Datenum Weekday Weeknum Colmon
73 2020-12-14 1681 18610 月曜日 2 6
81 2020-12-21 1791 18617 月曜日 2 6
88 2020-12-28 2397 18624 月曜日 2 6
367 2021-01-04 3330 18631 月曜日 2 6
344 2021-01-11 4901 18638 月曜日 2 6
351 2021-01-18 4947 18645 月曜日 2 6
359 2021-01-25 2760 18652 月曜日 2 6
373 2021-02-01 1782 18659 月曜日 2 6
2021-01-04以降のデータをデータフレーム として残します。
DailyData6 <-DailyData5[DailyData5$Datenum>=18631,]
head(DailyData6)
Date DailyCases Datenum Weekday Weeknum Colmon
367 2021-01-04 3330 18631 月曜日 2 6
368 2021-01-05 4958 18632 火曜日 3 8
369 2021-01-06 6066 18633 水曜日 4 8
370 2021-01-07 7793 18634 木曜日 5 8
371 2021-01-08 8045 18635 金曜日 6 8
372 2021-01-09 7528 18636 土曜日 7 8
tail(DailyData6)
同様に可視化します。
# 初期設定
DailyDataN2 <- DailyData6$Datenum
DailyDataY2 <- DailyData6$DailyCases
min2 <- min(DailyDataN2)
max2 <- max(DailyDataN2)
maxy2 <- max(DailyDataN2)
# 散布図プロット
par(xaxt="n") #x軸を書かない
par(yaxt="n") #y軸を書かない
par(ann=F) #軸タイトルを書かない
plot(DailyDataN2, DailyDataY2, xlim=c(min2,max2), ylim=c(0,maxy2), cex=0.8, col="grey", pch=19) #図
# 月曜に色つけ。y軸追加
DailyData6Mon <- DailyData6[DailyData6$Colmon== 6,]
head(DailyData6Mon)
tail(DailyData6Mon)
# par(new=T)でplot重ね書き
par(new=T)
par(yaxt="s")
plot(DailyData6Mon[,1], DailyData6Mon[,2], xlim=c(min2,max2), ylim=c(0,maxy2), col=DailyData6Mon$Colmon, pch=19) #図
# ｘ軸分割して回帰分析
### （ここから）
s2 <- 18631
col <- 2 #赤：色の初期値として設定
while (s2 < max2) {
# 7日間ずつ刻む
DailyData6pks1 <- DailyData6[DailyData6$Datenum>=s2,] #
s2 <- s2 + 6
DailyData6pks2<- DailyData6pks1[DailyData6pks1$Datenum<=s2,] #
s2 <- s2 + 1
# 上記でピックアップした7日間のx軸値とy軸値を抽出
DailyData6pks2x <- DailyData6pks2$Datenum #
print(DailyData6pks2x) #確認
DailyData6pks2y <- DailyData6pks2$DailyCases #
print(DailyData6pks2y) #確認
while (col > 4) { #色番号2,3,4：近似直線に赤，緑，青を順番に使うための処理
col = 2
}
# 回帰分析を実行し、近似直線を描く
ResultR2 <- lm(DailyData6pks2y~DailyData6pks2x) #
# lines(DailyData6pks2x , predict(ResultR2), lty = "dashed", lwd = 2, col="grey") #
lines(DailyData6pks2x , predict(ResultR2), lty = "dotted", lwd = 3.0, col=col, xlim=c(s2,max), ylim=c(0,maxy)) #図
col = col + 1 #近似曲線の始まりを重ならせるための処理
}
### （ここまで）
複数のグラフを重ね書きする際に、軸を一致させるのに思った以上に手こずりました。
少々スレが残っているかもしれませんが、確認と調整は後日に行うつもりです。
了