R言語で時間軸を小刻みに分割して回帰分析を実行し「勾配」の変化から傾向を探ってみる。「新型コロナウイルスの感染者数推移（全国）データ」を利用

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/cf62bb3c-609c-382e-9817-2356b2574d28.png)


##厚労省のオープンデータを利用

例えばn次曲線グラフのが急に変化する場合、その動きを確認（察知）するには、シンプルな（回帰）直線の勾配（傾き）を観察する方がわかりやすいのでは？、と考えました。
一定の期間ごとに分割し、その期間毎に回帰分析を実行することを試してみることにしました。
新型コロナウイルスの感染者数の推移（全国）のデータを使います。
データは[厚生労働省のサイト](https://www.mhlw.go.jp/stf/covid-19/open-data.html)から取得しました。
データ取得先はこちらです。

https://covid19.mhlw.go.jp/public/opendata/newly_confirmed_cases_daily.csv

Rで取得します。

```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
```

列名を変更します。

```R
colnames(covid)[2]<- "Pref"
colnames(covid)[3] <- "Cases"
colnames(covid)
[1] "Date"  "Pref"  "Cases"
```

##group_byとsummariseを使用

少し横道に逸れますが、「R」の復習がてら、県ごとにグルーピングし、県ごとのCasesの合計値を求めてみます。
まずグルーピングです。
その前に（重複計算するので）Pref=ALLを外しておきます。

```R
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の和を求めます。

```R
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
```

目的の日単位集計を実行します。

```R
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型だったので日付型に変換しています。
変換後、数字型に変えます。
数字の型でないと、回帰分析がかけられないからです。
（方法があるのかもしれませんが、今回は調べられませんでした）


```R
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）を使い、回帰分析をかけます。
散布図と近似直線を描きます。

```R
plot(DailyDataN, DailyDataY, cex=0.8, col="grey") #図
DailyDataR <- lm(DailyDataY~DailyDataN)
summary(DailyDataR)
abline(DailyDataR, col="grey", lty=2) #図
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/d6bbc292-797b-f1d3-4be4-4f836ee5515e.png)


今回の目的は、x軸を細分割して、回帰分析をかけ、分割した範囲ごとに近似直線を描くことです。
もう少しデータを整えます。
元のデータフレーム （DailyData）に数字にした日付の列（DailyDataN）を加えます。
曜日列、曜日番号（色の指定に使用）も加えます。

```R
#数字にした日付の列追加
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）を付るように列を追加します。

```R
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を使って日付順にソートします。

```R
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軸は指定しておきます。
（目盛りが上書きされると見にくくなったりするので）

```R
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) #図
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/a5bbffba-ff17-571f-91bd-e84bcd7a6fdf.png)


散布図の月曜のプロットの色（6：紫）を変えておきます。

```R
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) #図
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/bb3e3ba9-a878-a6ea-383c-b5c8d3d71d25.png)


##whileを使いx軸を分割して回帰分析

この散布図に期間を分割した近似曲線を描きます。
集計範囲（数字にした日付）を確認します。

```R
#集計データ開始日：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（日曜日）と週単位で分割できます。


```R
###（ここから）
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  #近似曲線の始まりを重ならせるための処理

}
###（ここまで）
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/22c434c5-bb76-9f01-22e5-c84bd32e7acc.png)


##ｘ軸を書き込む

x軸を書き込みます。


```R
#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") #図（月単位）
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/655002e3-78cb-dadf-e3b1-21956777e542.png)


##描画期間を変更する

描画範囲を2021年に絞ってみました。
21年の最初の月曜が「2021-01-04」（18631）でした。

```R
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以降のデータをデータフレーム として残します。

```R
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)
```

同様に可視化します。


```R
#初期設定
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  #近似曲線の始まりを重ならせるための処理

}
###（ここまで）
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/3a15f574-5202-f16c-256f-362706e007fb.png)


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/74b3a5ac-ae10-1c19-8032-e0b3eb9e4dc4.png)


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/f026f026-9ade-8349-90e7-4fec7a522cb1.png)


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/cf62bb3c-609c-382e-9817-2356b2574d28.png)


複数のグラフを重ね書きする際に、軸を一致させるのに思った以上に手こずりました。
少々スレが残っているかもしれませんが、確認と調整は後日に行うつもりです。

**了**