crosshyou

主にクロス表(分割表)分析をしようかなと思いはじめましたが、あまりクロス表の分析はできず。R言語の練習ブログになっています。

国民経済計算四半期GDP速報データの分析5 - R言語で株価、株価の騰落率とGDPデータの相関を調べる。

 

www.crosshyou.info

 の続きです。

今回は、株価、株価の騰落率とGDPデータの相関係数を調べてみようと思います。

相関係数はcor関数でわかります。

f:id:cross_hyou:20190420155849j:plain

soukan <- の行で、株価とGDPの各データの相関係数を計算しています。round関数で小数点以下3桁までにしています。

names(soukan) <- の行で算出した各データに名前を付けています。

rev(sort の行で相関係数の大きい順に表示しています。Kabukaと一番相関の高いのは、MinSetsu, つまり民間企業設備で、0.526ですね。いちおう確認しましょう。

f:id:cross_hyou:20190420160328j:plain

cor.test関数を使いました。一番下の行の0.5263771とあります。これが相関係数です。さきの計算結果と同じですね。p-value = 1.857e-08なので有意です。

1年後の株価騰落率、Kabuka.1Yとの相関も同じようにやってみましょう。

f:id:cross_hyou:20190420161615j:plain

Kabuka.1YにはNAの行があります。NAがあると相関係数がうまく計算できないので、na.omit関数でNAのないデータフレーム、data2をはじめにつくり、あとは株価との相関係数を算出するのと同じ手順で相関係数を算出しています。絶対値ベースでNetExpo(純輸出)が、-0.306と一番相関がありますね。確認します。

f:id:cross_hyou:20190420162049j:plain

cor.testの結果です。一番下の行の相関係数が -0.305555で一致しています。p-value = 0.003058ですから有意な相関ということです。マイナス相関ですから、純輸出がプラスのときは1年後は株価は安くなっている、純輸出がマイナスのときは1年後は株価は高くなっている、ということですね。

グラフで確認します。散布図を書いてみます。

f:id:cross_hyou:20190420162622j:plain

 

f:id:cross_hyou:20190420162634j:plain

lm関数で線形単回帰モデルを作り、そのモデルをabline関数で散布図の上に重ねてみました。右肩下がりですが、NetExpoだけで1年後の株価の騰落率を説明するのは難しいようですね。

今回は以上です。

今回のR言語のコードです。

# KabukaとGDP各データの相関係数
soukan <- round(cor(data$Kabuka, data[ , 1:12]), 3)
names(soukan) <- colnames(data)[1:12]
rev(sort(soukan))

# KabukaとMinsetsuの相関係数
cor.test(data$Kabuka, data$MinSetsu)

# Kabuka.1YとGDP各データの相関係数
data2 <- na.omit(data)
soukan2 <- round(cor(data2$Kabuka.1Y, data2[ , 1:12]), 3)
names(soukan2) <- colnames(data2)[1:12]
rev(sort(soukan2))

# Kabuka.1YとNetExpoの相関関係
cor.test(data2$Kabuka.1Y, data2$NetExpo)

# Kabuka.1YとNetExpoの散布図
plot(data2$NetExpo, data2$Kabuka.1Y)
abline(lm(Kabuka.1Y ~ NetExpo, data = data2), col = 2)

 

国民経済計算四半期GDP速報データの分析4 - R言語で各変数のヒストグラムとグラフを作成する。

 

www.crosshyou.info

 の続きです。

前回まではGDPと株価しかみていませんでしたが、データには、民間最終消費支出や、民間企業設備など様々な変数がありますから、それらを見ていきましょう。

まずは、names関数で変数名を確認します。

f:id:cross_hyou:20190420103931j:plain

MinShouhiは民間最終消費支出です。

hist関数とplot関数でグラフを描きましょう。

f:id:cross_hyou:20190420104727j:plain

f:id:cross_hyou:20190420104746j:plain

民間最終消費支出は山型の分布ですね。

次は、MinJuu(民間住宅)です。でもその前に、ヒストグラムと推移のグラフを作る関数を作ってしまいましょう。function関数です。

f:id:cross_hyou:20190420105545j:plain

こうして、ghp関数というヒストグラムと推移のグラフを描く関数を作成しました。ghp はgraph hist plot という意味です。この関数でMinJuu(民間住宅)を実行します。

f:id:cross_hyou:20190420105735j:plain

f:id:cross_hyou:20190420105748j:plain

うまくできました。民間住宅は右下がりですね。

次は、MinSetsu(民間企業設備)です。

f:id:cross_hyou:20190420110052j:plain

f:id:cross_hyou:20190420110103j:plain

民間企業設備は変動が大きいようです。

MinZaiko(民間在庫変動)はどうでしょうか?

f:id:cross_hyou:20190420110329j:plain

f:id:cross_hyou:20190420110342j:plain

双峰型のヒストグラムです。変動はさらに激しいのかトレンドが無い感じです。

SeiShouhi(政府最終消費支出)を見てみましょう。

f:id:cross_hyou:20190420110737j:plain

f:id:cross_hyou:20190420110754j:plain

政府最終消費支出は右肩上がりのトレンドです。

KoteiShihon(公的固定資本形成)はどうでしょうか?

f:id:cross_hyou:20190420111113j:plain

f:id:cross_hyou:20190420111124j:plain

公的固定資本形成は右肩下がりです。

KouZaiko(公的在庫変動)はどんな感じでしょうか?

f:id:cross_hyou:20190420111455j:plain

f:id:cross_hyou:20190420111523j:plain

公的在庫変動は、2005年ぐらいまでは大きく変動していましたが、以降は変動幅は小さくなっています。

NetExpo(財貨・サービス_純輸出)を見てみましょう。

f:id:cross_hyou:20190420111906j:plain

f:id:cross_hyou:20190420111918j:plain

2010年ぐらいから急激に落ち込んでいる時期があります。

次は、GrossExpo(財貨・サービス_輸出)です。

f:id:cross_hyou:20190420112256j:plain

f:id:cross_hyou:20190420112314j:plain

全体としては右肩上がりですが、2008年ぐらいに急激な落ち込みがありますね。リーマンショックのときでしょう。

GrossImpo(財貨・サービス_輸入)が最後です。

f:id:cross_hyou:20190420112712j:plain

f:id:cross_hyou:20190420112724j:plain

輸入も全体としては右肩上がりですが、輸入と同じように大きな落ち込みがあります。

今回は以上です。

 次回は

 

www.crosshyou.info

 

です。

今回のR言語のコードです。

# names関数で各変数の名前を確認
names(data)

# MinShouhi(民間最終消費支出)のグラフ
par(mfrow = c(1, 2))
hist(data$MinShouhi)
plot(data$Time, data$MinShouhi, type = "l")

# ヒストグラムと推移のチャートを描く関数の作成
ghp <- function(x) {
par(mfrow = c(1, 2))
hist(x)
plot(data$Time, x, type = "l")
}

# MinJuu(民間住宅)のグラフ
ghp(data$MinJuu)

# MinSetsu(民間企業設備)のグラフ
ghp(data$MinSetsu)

# MinZaiko(民間在庫変動)のグラフ
ghp(data$MinZaiko)

# SeiShouhi(政府最終消費支出)のグラフ
ghp(data$SeiShouhi)

# KoteiShihon(公的固定資本形成)のグラフ
ghp(data$KoteiShihon)

# KouZaiko(公的在庫変動)のグラフ
ghp(data$KouZaiko)

# NetExpo(財貨・サービス_純輸出)のグラフ
ghp(data$NetExpo)

# GrossExpo(財貨・サービス_輸出)のグラフ
ghp(data$GrossExpo)

# GrossImpo(財貨・サービス_輸入)のグラフ
ghp(data$GrossImpo)

 

 

国民経済計算四半期GDP速報データの分析3 - R言語でGDPと株価のチャートをつくる(その2)

 

www.crosshyou.info

 の続きです。

四半期ごとの株価のデータを景気動向指数の先行系列のデータから取得できたので、その株価データとGDPのデータをくっつけてチャートを作成したいと思います。

f:id:cross_hyou:20190418203554j:plain

このような株価のCSVファイルです。Kbuka-1Yは、1年後にどれだけ株価が上昇しているかを表しでいます。「来年の同じ四半期の株価 / 今年の四半期の株価」で計算しています。

まず、この株価のCSVファイルをread.csv関数で読込みます。skip = 4というオプションで5行目からデータを読込みます。

f:id:cross_hyou:20190418204246j:plain

str関数でデータ構造を、head関数で始めの6行を表示しました。

前回作成したデータフレーム、dataを同じようにstr関数とhead関数で確認しましょう。

f:id:cross_hyou:20190418204449j:plain

Timeの列で一緒にすればいいですね。merge関数を使います。

f:id:cross_hyou:20190418204745j:plain

dataにもkabukaにもKabukaという変数があったので統合後はKabuka.xとKabuka.yになっています。Kabuka.yのほうが正しいほうなので、Kabuka.xを削除します。subset関数を使います。

f:id:cross_hyou:20190418205157j:plain

Kabuka.xがなくなっていますね。

Kabuka.xをKabukaにしましょう。

f:id:cross_hyou:20190418210113j:plain

このように、Kabuka.yをKabukaという名前の新しい変数にデータをコピーして、subset関数でKabuka.yを削除しました。もっと簡単な方法がる気がしますが。

これでようやく四半期の株価がGDPデータに加わりました。plot関数でチャートを作ります。

f:id:cross_hyou:20190418210520j:plain

f:id:cross_hyou:20190418210532j:plain

前回のよりも株価が正確にチャートになりました。

株価をY軸、GDPをX軸に散布図を描き、回帰直線を重ねましょう。

f:id:cross_hyou:20190418210921j:plain

f:id:cross_hyou:20190418210937j:plain

今回は以上です。

 次回は

 

www.crosshyou.info

 

です。

 

今回のR言語のコードです。

# 株価(四半期)の株価を読み込む
kabuka <- read.csv("KABUKA_QUARTERY.csv", skip = 4)
str(kabuka)
head(kabuka)

# dataの確認
str(data)
head(data)

# merge関数でdataとkabukaを統合
data <- merge(data, kabuka, by = "Time")
str(data)
head(data)

# Kabuka.xの削除
data <- subset(data, select = -Kabuka.x)
str(data)

# Kabuka.yをKabukaに変更
data$Kabuka <- data$Kabuka.y
data <- subset(data, select = -Kabuka.y)
str(data)

# GDPと株価のチャート
par(mfrow = c(2,1))
plot(data$Time, data$GDP, col = "red", type = "l")
plot(data$Time, data$Kabuka, col = "blue", type = "l")

# 株価とGDPの散布図と回帰直線
plot(data$GDP, data$Kabuka, pch = 21, col = "red")
abline(lm(data$Kabuka ~ data$GDP), col = "blue")

 

国民経済計算四半期GDP速報データの分析2 - R言語で株価のデータをくっつけて、GDPと株価のチャートを作る。

 

www.crosshyou.info

 の続きです。

株価のデータを前回作成したデータフレームにくっつけようと思います。

1994年から2018年末までの株価のデータを探したら、日本取引所グループのウェブサイトに株価の長期時系列のファイルがありました。

f:id:cross_hyou:20190413154614j:plain

 

こういうファイルです。年初値、年末値、最低値、最高値があるのでこの4つの平均値をその年の平均値にしてしまいましょう。

f:id:cross_hyou:20190413155810j:plain

こんなふうに加工したCSVファイルをR言語のread.csv関数で読み込ませます。

f:id:cross_hyou:20190413160022j:plain

 

GDPのデータフレームにもこの株価のデータフレームにもTimeという名前の変数があり、これが共通です。このTime変数を基準にして、merge関数で二つのデータフレームを統合します。

f:id:cross_hyou:20190413160438j:plain

できました。

plot関数でGDPとKabukaの相関を見てみましょう。

f:id:cross_hyou:20190413160823j:plain

 

f:id:cross_hyou:20190413160839j:plain

均等に散らばっている感じですが、cor.test関数で相関関係を調べてみましょう。

f:id:cross_hyou:20190413161122j:plain

 

相関係数は0.503958でp-valueが9.033e-08ですから正の相関関係で有意ですね。

lm関数で線形単回帰分析をしてみます。

f:id:cross_hyou:20190413161447j:plain

p-valueは9.033e-08なので有意なモデルです。切片、GDPの係数ともにp値は0.05以下なので有意です。R-squareは0.254です。このモデルでは、株価をGDPで25%ぐらい説明できるということですかね。

回帰直線を散布図に重ねてみます。

f:id:cross_hyou:20190413161920j:plain

f:id:cross_hyou:20190413161952j:plain

残差プロットも描きます

f:id:cross_hyou:20190413162127j:plain

f:id:cross_hyou:20190413162138j:plain

横軸をTimeにして推移を描いてみましょう。

f:id:cross_hyou:20190413162516j:plain

f:id:cross_hyou:20190413162526j:plain

GDPはギザギザですね。Kabukaは4つ同じ値が続くので値動きがなめらかというか実際の株価の値動きよりもおだやかです。

今回は以上です。

株価の四半期のデータをどこかからとってきたいですね。

 次回は

 

www.crosshyou.info

 

です。

 

今回のR言語のコードです。

# 株価のCSVファイルの読込み
kabuka <- read.csv("kabuka.csv")
str(kabuka)

# dfとkabukaをTime列で統合
data <- merge(df, kabuka, by = "Time")
head(data)

# GDPとKabukaの散布図
with(data,
plot(GDP, Kabuka)
)

# GDPと株価の相関関係
cor.test(data$GDP, data$Kabuka)

# 株価 = a + b * GDP のモデル
model1 <- lm(Kabuka ~ GDP, data = data)
summary(model1)

# 回帰直線を散布図に重ねる
plot(data$GDP, data$Kabuka, pch = 21, bg = "red")
abline(model1, col = "blue")

# 残差プロット
plot(model1, which = 1)

# 横軸をTimeにしてGDPと株価の推移
par(mfrow = c(2, 1))
plot(data$Time, data$GDP, type = "l", col = "blue")
plot(data$Time, data$Kabuka, type = "l", col = "red")

 

国民経済計算四半期GDP速報データの分析1 - R言語のread.csv関数でデータを読み込み、分析しやすいようにデータを整える。

今回は、国民経済計算四半期GDP速報のデータをR言語のread.csv関数で読み込んでみたいと思います。

f:id:cross_hyou:20190413102954j:plain

https://www.e-stat.go.jp/dbview?sid=0003109741

このサイトからファイルをダウンロードしました。

f:id:cross_hyou:20190413103147j:plain

こういうファイルです。13行目だけ私が英語の変数名を付け足しています。

このファイルをread.csv関数で読込みます。

read.csv("ファイル名.csv")だと1行目から読み込まれてしまうので、skip = 12というオプションを付けて13行目から読込みます。

また、***や - や ーは数字が得られないものということなので、na.strings = c("***", "-", "ー")というオプションを加えてこれらの文字が入っているセルはNAになるように読込みます。

f:id:cross_hyou:20190413104449j:plain

このように始めの12行をスキップして、13行目から読み込むことができました。

Skip1とSkip2の変数はいらないので削除します。

f:id:cross_hyou:20190413104819j:plain

Timeをもう少し簡単にしたいですね。1994000103は1994年の1月から3月期、つまり1994年の第1四半期ですから、1994年のQ1ですね。数字にすると、1994.125が1994年Q1で、1994.375が1994年Q2で、1994.625が1994年Q3で1994.875が1994年Q4です。

とりあえず、このTimeからYearを作成します。

f:id:cross_hyou:20190413105745j:plain

お、うまくいきました。あとは、0.125, 0.375, 0.625, 0.875を足していけばいいですね。R言語はこういうのは簡単です。

f:id:cross_hyou:20190413110133j:plain

R言語は違う長さのベクトルの足し算は、短いベクトルを繰り返して長いベクトルに合わせて計算してくれますからね。

こうして作成したYearをdf$Timeと置き換えてしまいましょう。

f:id:cross_hyou:20190413110747j:plain

TimeJはもういらないので、削除してしまいましょう。

f:id:cross_hyou:20190413111013j:plain

はい、できました。

今回は以上です。

 次回は

 

www.crosshyou.info

 

です。

 

今回のR言語のコードです。


# データの読込み
df <- read.csv("gdp.csv", skip = 12, na.strings = c("***", "-", "-"))
head(df)

# Skip1とSkip2の削除
df <- df[ , -c(1, 2)]
head(df)

# Yearの作成
Year <- round(df$Time / 1000000, 0)
head(Year)

# Yearに小数点以下を足し算
Year <- Year + c(0.125, 0.375, 0.625, 0.875)
head(Year)
tail(Year)

# TimeをYearと置き換え
df$Time <- Year
head(df)

# TimeJを削除
df <- df[ , -2]
head(df)

 

 

商業統計調査データの分析5 - R言語のstep関数で回帰モデルの最適な変数を選ぶ。

 

www.crosshyou.info

 の続きです。

まずは、str関数でデータを確認しましょう。

f:id:cross_hyou:20190410190932j:plain

まず、Yearが何年からあるか確認します。table関数を使います。

f:id:cross_hyou:20190410191331j:plain

1988年から始まって、3年ごとに2012年までの7つの年のデータがあります。

一つの年で、77のデータがあるということは、Sectorが77種類あるということですね。

今回は、Corp_TotalをYear, Staff, Revenue, Inventory, Spaceで回帰分析してみて、どの変数がCorp_Totalに一番影響を与えているのかみてみましょう。

線形回帰分析のlm関数でやってみます。

f:id:cross_hyou:20190410192308j:plain

どの係数も有意でよくわからないですね。

step関数で効果の低い変数を段階的に削除できるらしいです。やってみます。

f:id:cross_hyou:20190410192736j:plain

結局、すべての変数が入っているモデルが一番いいということなのですね。

step関数では、何も変数のないモデルからはじめて、一つずつ変数を追加するという方法もあります。これをやってみます。

まずは変数のないゼロモデルを作成します。

f:id:cross_hyou:20190410193200j:plain

切片は、86112とありますが、これはCorp_Totalの平均値です。

f:id:cross_hyou:20190410193336j:plain

さて、このmodel_zeroからstep関数で変数を追加していきます。

f:id:cross_hyou:20190410193735j:plain

前進ステップワイズ回帰でもやっぱり全部の変数が含まれますね。

今回は以上です。

 

今回のR言語のコードです。


# データの確認
str(df)

# Yearの確認
table(df$Year)

# Corp_TotalをYear, Staff, Revenue, Inventory, Spaceで回帰
model1 <- lm(Corp_Total ~ Year + Staff + Revenue + Inventory + Space, data = df)
summary(model1)

# step関数
reduced_model <- step(model1, direction = "backward")

# ゼロモデルの作成
model_zero <- lm(Corp_Total ~ 1, data = df)
summary(model_zero)

# Corp_Totalの平均値
mean(df$Corp_Total)

# step関数(前進ステップワイズ回帰)
model_zenshin <- step(model_zero, direction = "forward",
scope = (~ Year + Space + Inventory + Revenue + Staff),
trace = 0)
summary(model_zenshin)

商業統計調査データの分析4 - R言語で法人と個人の事業所の割合を比較する。

 

www.crosshyou.info

 の続きです。

今回は、R言語で法人の事業所の割合と個人の事業所の割合を出して、経年変化をみてみたいと思います。

まずは、もとになるデータをhead関数とsummary関数で確認しましょう。

f:id:cross_hyou:20190406145917j:plain

f:id:cross_hyou:20190406145929j:plain

このようなデータです。Corp_Totalが法人と個人を合わせた事業所の数、Corp_Corpが法人の事業所の数、Corp_Indiが個人の事業所の数でした。

Corp_Corp / Corp_Total で法人の割合、Corp_Indi / Corp_Totalで個人の割合になります。早速作成しましょう。

f:id:cross_hyou:20190406150706j:plain

まずは法人の割合を作成しました。平均値と中央値が同じ値ですね。0.4082です。最大値は0.8766、最小値は0.0244です。

個人の割合はこうなりました。

f:id:cross_hyou:20190406150834j:plain

個人の割合も平均値と中央値が同じになりました。

こうしてみると、だいたい6割が個人の事業所、4割が法人の事業所だとわかりました。

hist関数でヒストグラムを作成してみます。

f:id:cross_hyou:20190406151213j:plain

f:id:cross_hyou:20190406151226j:plain

前回の一人当り売上高、売場面積当り売上高と比べると左右対称の度合いが高いですね。

年毎の比率を見てみましょう。1 - 法人の割合 = 個人の割合ですから、ここからは法人の割合だけを見ていきます。tapply関数を使います。

f:id:cross_hyou:20190406151614j:plain

1988年は0.32ぐらいでしたが、年が経つたびに上昇し、2012年には0.499と半分近くにまで高まっています。

平均値だけでなくて、年ごとの箱ひげ図を作成します。Yearをfactor関数でファクタにしてX軸にしてplot関数です。

f:id:cross_hyou:20190406152048j:plain

f:id:cross_hyou:20190406152100j:plain

法人の割合 = a + b * 年
という線形モデルを作成して年が法人の割合に有意に影響しているかを調べましょう。

lm関数を使います。

f:id:cross_hyou:20190406152421j:plain

p-value: < 2.2e-16とありますので、このモデルは有意です。Yearの行のPr(>|t|)を見ると、<2e-16とありますので、Yearは有意な影響を及ぼしています。

plot関数で散布図を描いて、abline関数でモデルの回帰直線を重ねてみます。

f:id:cross_hyou:20190406152933j:plain

f:id:cross_hyou:20190406152945j:plain

こんな感じです。モデルの残差をグラフにしてみましょう。

f:id:cross_hyou:20190406153206j:plain

f:id:cross_hyou:20190406153218j:plain

このようになりました。

ところで、このYearだけで法人の割合を回帰するmodel1は、Adjusted R-squaredが0.1301とかなり低いです。他の変数も加えてもう少しAdjusted R-squaredの数値をあげてみたいです。法人の割合が多いのはどんな場合でしょうか?一事業所当りの売上高が多い業種は法人の割合が多いような気がします。やってみましょう。まず、一事業所当りの売上高の変数を作成します。

f:id:cross_hyou:20190406153828j:plain

一事業所当りの売上高の平均値は役1億2100万円です。中央値は約5600万円です。

この変数を加えて回帰分析をしてみましょう。

f:id:cross_hyou:20190406154239j:plain

一番下の行のp-valueは2.2e-16よりも小さいの有意です。Yearの係数もRev_Corの係数もp値は2e-16以下なので有意です。Adjusted R-Squaredは0.2554ですので、model1よりも高くなりました。

回帰残差をプロットしてみます。

f:id:cross_hyou:20190406154917j:plain

こうなりました。

anova関数を使って、model1とmodel2が有意に異なるかチェックします。

f:id:cross_hyou:20190406155029j:plain

p値が2.2e-16以下なので有意に異なります。model1に一事業所当り売上高を加えたことは意味があったということです。

一事業所当りの従業者数も加えてみましょう。一事業者当りの従業者数が大きいということは大きな事業所ということを示唆しますし、大きな事業所ということは法人の運営する事業所だと思われます。

f:id:cross_hyou:20190406155527j:plain

平均すると一事業所当り約6人いますね。

それではこのStf_Corを加えたmodel3を作成してみましょう。

f:id:cross_hyou:20190406155926j:plain

モデル全体はp-value* < 2.2e-16とありますので有意ですが、Rev_Cor, Stf_Corの係数が有意ではないですね。Rev_Cor(一事業所当り売上高)を除いて、Stf_Cor(一事業所当り従業者数)だけを残したmodel4を試してみます。

f:id:cross_hyou:20190406160406j:plain

model4は、Adjusted R-squaredが0.2605と一番いいです。

回帰残差プロットを描いてみます。

f:id:cross_hyou:20190406160709j:plain

f:id:cross_hyou:20190406160723j:plain

こうなりました。

回帰分析はここまでにして、法人の割合の高い業種と個人の割合の高い業種を調べてみましょう。tapply関数とsort関数とrev関数とhead関数です。

f:id:cross_hyou:20190406161054j:plain

法人の割合の高い業種は上の業種でした。各種商品小売業が約8割です。以下、燃料小売業、その他の各種商品小売業(従業者が常時50人未満のもの)、楽器小売業、農業用機械器具小売業、その他の機械器具小売業と続きます。

個人の割合の高い業種は

f:id:cross_hyou:20190406161454j:plain

こうなりました。たばこ・喫煙具専門小売業が約96%です。以下、自転車小売業、履物小売業(靴を除く)、骨とう品小売業、卵・鳥肉小売業、鮮魚小売業と続きました。

今回は以上です。

 次回は

 

www.crosshyou.info

 

です。

 

今回のR言語のコードです。


# head関数
head(df)

# summary関数
summary(df)

# 法人の割合
df4 <- df
df4$R_Corp <- df$Corp_Corp / df$Corp_Total
summary(df4$R_Corp)

# 個人の割合
df4$R_Indi <- df$Corp_Indi / df$Corp_Total
summary(df4$R_Indi)

# ヒストグラム
par(mfrow = c(1, 2))
hist(df4$R_Corp)
hist(df4$R_Indi)

# 年毎の法人割合の平均値
tapply(df4$R_Corp, df4$Year, mean)

# 年ごとの箱ひげ図
plot(factor(df4$Year), df4$R_Corp)

# 法人の割合 = a + b * 年
model1 <- lm(R_Corp ~ Year, data = df4)
summary(model1)

# 散布図と回帰直線
plot(df4$Year, df4$R_Corp)
abline(model1, col = "red", lty = 2)

# model1の残差プロット
plot(model1, which = 1)

# 一事業所当りの売上高
df4$Rev_Cor <- df4$Revenue / df4$Corp_Total
summary(df4$Rev_Cor)

# 法人の割合 = a + b1 * 年 + b2 * 一事業所当りの売上高
model2 <- lm(R_Corp ~ Year + Rev_Cor, data = df4)
summary(model2)

# 回帰残差のプロット
plot(model2, which = 1)

# model1とmodel2が有意に異なるか
anova(model1, model2)

# 一事業所当りの従業者数
df4$Stf_Cor <- df4$Staff / df4$Corp_Total
summary(df4$Stf_Cor)

# model3(model2 + 一事業所当りの従業者数)
model3 <- lm(R_Corp ~ Year + Rev_Cor + Stf_Cor, data = df4)
summary(model3)

# model4(model1 + 一事業所当りの従業者数)
model4 <- lm(R_Corp ~ Year + Stf_Cor, data = df4)
summary(model4)

# model4の回帰残差プロット
plot(model4, which = 1)

# 法人の割合の高い業種
head(rev(sort(tapply(df4$R_Corp, df4$Sector, mean))))

# 個人の割合の高い業種
head(rev(sort(tapply(df4$R_Indi, df4$Sector, mean))))