「ネットショッピングの支出額」を使い「R言語」で線形単回帰分析を再学習

# R言語を使用

* 今回のポイントとなった**「R」**アイテム
    * **read.csv()** ：csv読み込み
    * **names()** ：列名変更
    * **par(family = "HiraKakuProN-W3")**   ：　文字化け回避
    * **library(zoo)** ：時間軸への変更
    * **gsub** ：文字列変換
    * **mutate(index = row_number())**  ：行の数だけ連番発行
    * **lm(ye~x, data=xxx)**  ：　線形回帰分析
    * **dev.new()**  ：新しい描画ウインドウ
    * **predict()**  ：予測関数、信頼区間算出
    * **par(new=TRUE)**   ：　グラフ上書き
    * **par(xaxt="n"), par(xaxt ="s")**   ：　x軸描かない、描く

「R言語」を使って試したいことができたので、その準備のために「線形単回帰分析」の方法を再確認しました。
使ったデータは、[以前も使った](https://qiita.com/text_owleyes_0214/items/ee67ac14c3f25d953be9)「家計消費状況調査」（ネットショッピングの支出額：[5/11総務省統計局公表](https://www.stat.go.jp/data/joukyou/12.html)）です。
データは[「e-Stat」（政府統計の総合窓口）](https://www.e-stat.go.jp/)から取り出します。
2015年1月からデータが収録されています。
まず、「e-Stat」上で月ごとの推移グラフを確認します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/1f449aed-9fe5-2ba9-677e-864922b32e71.png)

全体的に増加傾向です。
回帰分析をかけます

```r
#"CP932”は、Macで読み込む際に文字化けを防ぐために使用
nets1 <- read.csv("FEH_netshop.csv", header=T,fileEncoding="CP932")
head(nets1)

     時間軸 インターネットを利用した支出総額_全国
1 2015年1月                                 9,531
2 2015年2月                                 7,835
3 2015年3月                                 8,982
4 2015年4月                                 8,449
5 2015年5月                                 8,727
6 2015年6月                                 8,831

str(nets1)
'data.frame':	75 obs. of  2 variables:
 $ 時間軸                               : Factor w/ 75 levels "2015年10月","2015年11月",..: 4 5 6 7 8 9 10 11 12 1 ...
 $ インターネットを利用した支出総額_全国: Factor w/ 75 levels "10,049","10,360",..: 71 50 65 59 62 63 68 67 51 48 ...
```
よくやってしまうのですが、支出額が”10,049”と「,」を含む値として認識されています。
「as.numeric(nets1[,2])」と数値型にしてもFactorのレベル値が数値化されます。
このため、元のcsvデータからカンマをなくす書式にし、再度Rで読み込みます。

```r
nets1 <- read.csv("FEH_netshop.csv", header=T,fileEncoding="CP932")
head(nets1)
 
     時間軸 インターネットを利用した支出総額_全国
1 2015年1月                                  9531
2 2015年2月                                  7835
3 2015年3月                                  8982
4 2015年4月                                  8449
5 2015年5月                                  8727
6 2015年6月                                  8831
str(nets1)

'data.frame':	75 obs. of  2 variables:
 $ 時間軸                               : Factor w/ 75 levels "2015年10月","2015年11月",..: 4 5 6 7 8 9 10 11 12 1 ...
 $ インターネットを利用した支出総額_全国: num  9531 7835 8982 8449 8727 ...
```
# 時系列型のx軸に

これもこなれた扱いができないのが、「YYYY年MM月」をx軸として扱う場合です。
取り込んだ時点では「時間軸」もファクターなので、こちらは時系列型にします。
その前に、扱いやすいように列名を変更します。

```r
names(nets1) <- c("yyyymm",  "expense")
head(nets1)

     yyyymm  expense
1 2015年1月     9531
2 2015年2月     7835
3 2015年3月     8982
4 2015年4月     8449
5 2015年5月     8727
6 2015年6月     8831
```
時系列型（という言い方は適切ではないかもしれませんが）に変換します。
「YYYY-MM」という形式に変換してから、「zoo」ライブラリを使って時系列のnum型にします。

```r
#「YYYY-MM」という形式に変換
nets2 <- nets1
nets2["yyyymm"] <- lapply(nets2["yyyymm"],  gsub, pattern="年",  replacement = "-")
head(nets2)
nets2["yyyymm"] <- lapply(nets2["yyyymm"],  gsub, pattern="月",  replacement = "")
head(nets2)

  yyyymm  expense
1 2015-1     9531
2 2015-2     7835
3 2015-3     8982
4 2015-4     8449
5 2015-5     8727
6 2015-6     8831

install.packages("zoo")
library(zoo)

nets2[,1] <- as.yearmon(nets2[,1])
str(nets2[,1])
 'yearmon' num [1:75]  1 2015  2 2015  3 2015  4 2015 ...

head(nets2)

   yyyymm  expense
1  1 2015     9531
2  2 2015     7835
3  3 2015     8982
4  4 2015     8449
5  5 2015     8727
6  6 2015     8831
```
散布図でプロットします。

```r
par(family = "HiraKakuProN-W3")
plot(nets2)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/a778efd0-db97-55ba-7342-f953b21602de.png)

# 線形単回帰分析を実行

まず、「cor」を使い相関係数を求めてみます。
先ほどx軸を時系列で扱えるようにしましたが、相関係数を求めるには、x軸は数値である必要がありました。

```r
#連番の列を追加します
library(dplyr)
nets3 <- nets2 %>% mutate(index = row_number()) 

#相関係数を求める
cor(nets3$index, nets3$expense)
[1] 0.9217676
```
「0.92」なので相関関係があることが読み取れます。
線形単回帰分析を実行します

>*y=ax+b*

の近似直線を求めます
「expense」を目的変数、「index」を説明変数とします。

```r
nets3.lm <- lm(nets3$expense~nets3$index, data=nets3)
summary(nets3.lm)

Call:
lm(formula = nets3$expense ~ nets3$index, data = nets3)

Residuals:
    Min      1Q  Median      3Q     Max 
-2587.9  -857.7  -217.5   682.7  4734.4 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6695.23     303.52   22.06   <2e-16 ***
nets3$index   140.96       6.94   20.31   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1301 on 73 degrees of freedom
Multiple R-squared:  0.8497,	Adjusted R-squared:  0.8476 
F-statistic: 412.6 on 1 and 73 DF,  p-value: < 2.2e-16

#回帰係数（aとb）は上記にも含まれているが、下記でも出力できる
round(coefficients(nets3.lm) ,2) 
(Intercept) nets3$index 
    6695.23      140.96 
```
p値（Pr(>|t|) ）が有意水準0,005（0.5%）より小さく、その横に「*」が3つ並んでいるので、係数が回帰モデルに役立っていると言えます。　※参考「Rによるデータサイエンス（金明哲著、森北出版）」
>*y=140.96x+6695.23*

という近似直線が得られました。
indexを軸にした散布図を描き、その上のにこの直線を描きます。

```r
plot(nets3$index, nets3$expense)
abline(nets3.lm, lwd=2, col="grey") 
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/57893be3-afc1-fd69-84e2-d5f44274c6dc.png)

# 残差と予測値

実際の値と近似直線の差を求め、散布図にします。

```r
#残差を求める
nets3_re <- residuals(nets3.lm)
dev.new() #新しいウインドウで描画
plot(nets3_re)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/a28684ef-f54d-8ac9-44e1-e9508d167d0e.png)

予測値は「pedict」を使い求ました。
「予測値」とは、回帰式で計算される値です。

```r
nets3_pr <- predict(nets3.lm)
dev.new() #新しいウインドウで描画
plot(nets3_pr)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/6a583866-d910-643a-7e5b-c048c6417f72.png)

# 信頼区間を求める

近似（回帰）直線の「95%信頼区間」を求め、描画します。

```r
#信頼区間の値を求める
nets3.pr95 <- predict(nets3.lm, interval="confidence", level=0.95)
#データフレームにする
nets3.pr95 <- as.data.frame(nets3.pr95) 
#元の実測値データと結合し1つのデータフレーム とする
nets4 <- data.frame(nets3, nets3.pr95)
head(nets4)

   yyyymm expense index      fit      lwr      upr
1  1 2015    9531     1 6836.194 6243.260 7429.128
2  2 2015    7835     2 6977.157 6396.119 7558.195
3  3 2015    8982     3 7118.120 6548.891 7687.349
4  4 2015    8449     4 7259.083 6701.570 7816.596
5  5 2015    8727     5 7400.046 6854.150 7945.943
6  6 2015    8831     6 7541.009 7006.624 8075.394
```
実測値の散布図と、近似直線を再度描き、95%信頼区間も加えます。
複数のグラフを重ねるので、y軸の範囲を指定します。

```r
par(family = "HiraKakuProN-W3")
plot(nets4$index, nets4$expense,, ylim=c(6000,22000))
abline(nets3.lm, lwd=2, col="grey") 

par(new=TRUE)  #グラフを上書き
plot(nets4$index, nets4$lwr, ylim=c(6000,22000), col="blue", type="l", lwd=1, ylab="", xlab="")
par(new=TRUE)  #グラフを上書き
plot(nets4$index, nets4$upr, ylim=c(6000,22000), col="blue", type="l", lwd=1, ylab="", xlab="")
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/e53f23fb-8b45-9a4e-4d4a-5cdbe32ed20a.png)

# x軸を時系列に戻す

回帰分析をかけるために、x軸を数値型にしています。
時系列（日付型）に戻そうと思いましたが、効率的な方法がわかりませんでした。
yyyymm列を数値型にして集計しておけば、（確か）「axis.Date(1, at = seq(min(x),max(x)," days”))」を使えば可能だと思いました。
今回は使わず、まず、x軸目盛りを表示しないで作成し、その後でyyyymmとexpeseの散布図を重ね書きすることにしました。

```r
#par(xaxt=“n”)はx軸を表示しない指定
par(family = "HiraKakuProN-W3")
par(xaxt="n")
plot(nets4$index, nets4$expense, ylim=c(6000,22000))
abline(nets3.lm, lwd=2, col="grey") 

par(new=TRUE)  #グラフを上書き
par(xaxt="n")
plot(nets4$index, nets4$lwr, ylim=c(6000,22000), col="blue", type="l", lwd=1, ylab="", xlab="")
par(new=TRUE)  #グラフを上書き
par(xaxt="n")
plot(nets4$index, nets4$upr, ylim=c(6000,22000), col="blue", type="l", lwd=1, ylab="", xlab="")

#par(xaxt=“s”)はx軸を表示する指定
par(new=TRUE)  #グラフを上書き
par(xaxt ="s")
plot(nets4$yyyymm, nets4$expense, ylim=c(6000,22000), ylab="", xlab="")
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1143963/ce514c1e-8fcf-19bc-e6dd-05cdf18d7b84.png)

>**参考させて頂いたサイト**

https://kitamix.net/archives/save-cdv-file-in-r/1096

https://www.karada-good.net/analyticsr/r-58

http://cse.naro.affrc.go.jp/takezawa/r-tips/r/48.html

https://www.it-swarm-ja.com/ja/r/%E5%B9%B4%E3%81%A8%E6%9C%88%EF%BC%88-yyyymm%E5%BD%A2%E5%BC%8F%EF%BC%89%E3%82%92%E6%97%A5%E4%BB%98%E3%81%AB%E5%A4%89%E6%8F%9B%E3%81%97%E3%81%BE%E3%81%99%E3%81%8B%EF%BC%9F/972925475/

https://fujinitaka.hatenablog.com/entry/2018/07/08/022044

http://by-oneself.com/r_contents_4/

https://qiita.com/weda_654/items/c480d6d45344fd947e6b

https://qiita.com/HinaKimu/items/cde3f828aa7e6c88db4b

https://blog.statsbeginner.net/entry/2014/10/18/130504

https://myopomme.hatenablog.com/entry/20120123/1327331897

http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html

http://blog.livedoor.jp/mondayish/archives/1581974.html

https://teramonagi.hatenablog.com/entry/20110520/1305891395

https://www.mitsumatado.com/zen/r-plot2/

http://takenaka-akio.org/doc/r_auto/chapter_13.html

https://qiita.com/tbpgr/items/989c6badefff69377da7

**了**