全産業活動指数・建設業活動指数・鉱工業生産指数・第３次産業活動指数の分析４ - R言語のlm関数で回帰分析

www.crosshyou.info

 の続きです。

今回は、全産業活動指数を被説明変数、その他の三つの指数を説明変数にして回帰分析をしてみます。lm関数を使います。

まずは、建設業活動指数で回帰分析します。

p-valueは5.995e-15なので有意なモデルです。R-squaredは0.3641です。

次は、鉱工業生産指数で回帰分析します。

p-valueは2.2e-16よりも小さいので有意です。R-squaredは、0.7077です。

第３次産業活動指数はどうでしょうか？

p-valueは、2.2e-16よりも小さいので有意です。R-squaredは0.8826です。

第３次産業活動指数との回帰分析モデルがR-squaredの値が大きいので一番あてはまりが良いようです。

plot関数とabline関数でグラフにしてみます。

グラフにすると、建設業活動指数はあんまり全産業活動指数と連動していないことがわかりますね。

今度は、当月の全産業活動指数を、前月のそれぞれの指数で回帰分析してみます。

まずはこのように、全産業活動指数は、始めの値を削除し、その他の三つは最後の値を削除します。

まずは、建設業活動指数からです。

p-valueは7.898e-13と0.05よりも低いので有意です。R-squaredは0.3191です。

次は鉱工業生産指数です。

p-valueは2.2e-16よりも小さいので有意です。R-squaredは0.6053です。

第３次産業活動指数はどうでしょうか？

p-valueは2.2e-16よりも小さいので有意です。R-squaredは0.7802です。

いままで6つの回帰モデルを作成しました。この6つの回帰モデルのR-squaredを並べてみます。

建設業活動指数：0.3641

前月の建設業活動指数：0.3191

鉱工業生産指数：0.7077

前月の鉱工業生産指数：0.6054

第３次産業活動指数：0.8826

前月の第３次産業活動指数：0.7802

です。前月の指数よりも同じ月の指数のほうがR-squaredは大きいですね。

そして、第３次産業活動指数が一番、全産業活動指数をよく説明する指数であることがわかりました。

今回は以上です。

全産業活動指数・建設業活動指数・鉱工業生産指数・第３次産業活動指数の分析３ - plot関数・hist関数・boxplot関数で各種グラフを作成

www.crosshyou.info

 の続きです。

まずは、plot関数とlines関数で各指数の推移をグラフにしてみます。

黄色い線は鉱工業生産指数で大きく落ち込むときがあったことがわかります。

赤が全産業活動指数で、緑が第３次産業活動指数です。同じような動きだとわかります。

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

boxplot関数で箱ひげ図を作成します。

今回は以上です。

全産業活動指数・建設業活動指数・鉱工業生産指数・第３次産業活動指数の分析２ - Varianceの練習

www.crosshyou.info

 の続きです。

今回は

Statistics: An Introduction Using R by Michael J. Crawley(2014-11-24)

 の第4章、Varianceのところを練習してみようと思います。

varianceの定義ですが、まず、平均値を計算して、個々のデータと平均値の差を2乗し、それを合計します。その合計した値をデータの個数より一つ少ない数で割ります。

文章よりも数式のほうが簡単でしょう。

function関数で、my_varという名前でvarianceを計算する関数をつくりました。Rにはもともとvar関数というvarianceを計算する関数があります。全産業活動指数のvarianceは8.7ぐらいですね。他の三つも計算します。

建設業活動指数のvarianceは37.7ぐらい、鉱工業生産指数のvarianceは42.4ぐらい、第３次産業活動指数のvarianceは5.3ぐらいです。建設業と鉱工業生産のvarianceは大きく、全産業と第３次産業は小さいですね。

さて、varianceは何に使うかというと、信頼性の尺度で使うそうです。

standard errorsという尺度で、varianceをデータの個数で割って平方根を計算した値がstandard errorだそうです。早速計算してみます。

これはどういうことかと言うと、

全産業活動指数を例にすると、平均値とデータの個数も使って

全産業活動指数は、102.1496±0.252(1 s.e., n = 137)と書くそうです。

さて、もう一つvarianceから派生する値があります。これはconfidence intervalというもで、〇〇%の確率で平均値がconfidence intervalの中に入る、ということのようです。

例えば、95%の確率だと、qt(0.975, 自由度) x standard errors という計算式のようです。qt関数は自由度ど確率を与えてt値を計算する関数のようです。

早速計算してみましょう。

全産業活動指数のconfidence intervalは0.498ぐらいです。

これは、全産業活動指数は、102.1496±0.498(95% CI, n = 137)と書くようです。

今度は、平均値、データの個数、standard error、confidence intervalを一度に計算する関数を作成して、4つの変数で試して終わりにしましょう。

まずは、このように関数を作成しました。平均値、観測数、standard errors, confidence intervalを返す関数です。

さっそく４つの変数で試します。

こうなりました。文章で表現すると、

全産業活動指数は、102.15±0.25(1 s.e., n = 137), 102.15±0.50(95% CI, n = 137)

建設業活動指数は、106.43±0.52(1 s.e., n = 137), 106.43±1.04(95% CI, n = 137)

鉱工業生産指数は、99.16±0.56(1 s.e., n = 137), 99.16±1.10(95% CI, n = 137)

第３次産業活動指数は、102.66±0.20(1 s.e., n = 137), 102.66±0.39(95% CI, n = 137)

となります。

今回は以上です。

全産業活動指数・建設業活動指数・鉱工業生産指数・第３次産業活動指数の分析１ - 基本統計量

今回は全産業活動指数の分析をしようと思います。

政府統計の総合窓口(e-Stat)のウェブサイトからデータを取得しました。

今回は、

政府統計名：全産業活動指数

提供統計名：全産業活動指数(IAA)

提供分類1：長期時系列データ

提供分類2：平成22年(2010年)基準

提供分類3：時系列データ(平成20年01月~最新)

のデータの

表番号：全産業活動指数(IAA)

統計表：季節調整済指数の月次のファイルをダウンロードしました。

こういうファイルです。横長のファイルですね。2行目にCol1, Col2と付け足しました。

まず、このファイルをread.csv関数でR言語に読込ませます。

読込んだファイルは、観測数18 x 変数140 のファイルでした。

なので、本当は変数18、観測数140のデータですね。

読込んだデータフレームの左上の部分を見てみましょう。

2列目、Col2の値が本当の変数名ですね。

全産業活動指数から<参考系列>全産業活動指数(直接調整法)までの17種類ですね。

今回はこの中から、

全産業活動指数・建設業活動指数・鉱工業生産指数・第３次産業活動指数の４つのデータに注目してみます。

とりあえず、月、全産業活動指数、建設業活動指数、鉱工業生産指数、第３次産業活動指数の４つのデータをそれぞれ、独立したベクトルとして抽出しました。

１行目のコマンドで、それぞれのデータを抽出し、それをas.matrix関数でデータフレーム型からマトリックス型に変換し、as.vector関数でマトリックス型からベクトル型に変換しました。str関数でデータ構造を確認しています。すべて137個の要素のある数値型のベクトルになりました。

summary関数で月の基本統計量を算出しました、最小値は200801、最大値は201905なので、2008年01月から2019年05月までのデータですね。

sd関数で標準偏差を算出し、標準偏差をmean関数で算出する平均値で割って変動係数を算出しています。

他の指数も同じようにやってみます。

最小値に着目すると、鉱工業生産指数の76.60が一番低く、

最大値に着目すると、建設業活動指数の118.7が一番高く、

変動係数に着目すると、第３次産業活動指数の0.022が一番低く、鉱工業生産指数の0.066が一番高く、

中央値に着目すると、鉱工業生産指数の99.20が一番低く、建設業活動指数の107.8が一番高く、

平均値に着目すると、鉱工業生産指数の99.16が一番低く、建設業活動指数の106.4が一番高いことがわかります。

今回は以上です。

 次回は

www.crosshyou.info

 です。

都道府県別の財政力指数の分析２ - Central Tendancyの練習

今日は、

www.crosshyou.info

 の続き、というかこのデータを使って、

Statistics: An Introduction Using R

 の第3章、Central Tendencyの練習をしようと思います。

まず、ZaiAllというベクトルを作成し、as.character関数とpaste関数を使って、Yearとprefを連結した文字列ベクトルを作成し、names関数でZaiAllに名前の属性をつけて、head関数ではじめの6データを表示しました。

hist関数でヒストグラムを描きました。label = TRUEとして度数を表示しています。

mean関数で平均値、median関数で中央値を算出し、ablineでヒストグラムに重ねました。赤が平均値、青が中央値です。

ヒストグラムは、右の裾野が広い分布形状です。

上に掲載した本では、Central Tendancyとして、arithmetic mean, median, geometric mean, harmonic meanの4種類が記述されています。

arithmetic meanはR言語ではmean関数がはじめからあります。

sum(x) / length(x) 、つまりデータの合計値をデータの個数で割ったものです。

function関数でそのような関数をつくり計算しましたが、mean関数と一致します。

中央値はmedian関数という出来合いの関数がありますが、function関数で作成しました。

ベクトルの個数が奇数ならその真ん中の値になるので、sort(x)で小さい順に並べて、ceiling(ベクトルの個数 / 2)番目の値が中央値になります。

ベクトルの個数が偶数なら、ちょうど半分のところの上下の値の平均値になります。

%% 2 で2で割った余りを計算して、これが0なら偶数の個数、そうでないなら奇数と判断しています。

geometric mean, 幾何平均は全てのデータを掛けてそれをデータの個数分の平方根？だそうですが、exp(mean(log(x)))で計算できるそうです。

harmonic mean 調和平均はこのようにして計算します。

東京から大阪まで行きは時速100キロ、帰りは時速200キロで移動したときの平均速度は？という質問の答えがこのharmonic meanだそうです。

boxplot関数で箱ひげ図を描いてみます。

黄色の中央値と青の幾何平均がほとんど同じでわからないですね。

1.0以上の値をはずして箱ひげ図を描いてみます。

ZaiAll[ZaiAll < 1]とすると、1.0以上の値は除外できます。

今回は以上です。

都道府県別の財政力指数の分析 - 人口が多い都道府県、面積が狭い都道府県ほど指数は高い。東京都が一番高い

今回は都道府県別の財政力指数を分析しようと思います。

財政力指数とは、

政府統計の総合窓口(e-Stat)にありました。地方公共団体の財政力の強さを表す指標だそうです。

このデータをCSVファイルにダウンロードして、R言語で分析します。

read.csv関数でCSVにあるデータをR言語に読込みます。

na.omit関数でNAのある行を削除します。

as.character関数でファクターを文字列に直します。

as.character関数で文字列をファクターに直します。

summary関数で各変数の統計値を出します。Prefの各都道府県の度数が42とありますから、42年間のデータがあることがわかります。

Zaiが財政力指数です。47都道府県 X 42年間 = 1974個のデータの最小値が0.1970で最大値が1.6395です。平均値は0.4749, 中央値は0.4305です。

まずは、hist関数で財政力指数のヒストグラムを描いてみましょう。

右側の裾野が広い分布ですね。

年度ごとの箱ひげ図をplot関数で描きます。

箱ひげ図の上にある外れ値を見ると、1990年近辺が高いですね。

tapply関数で年度別の平均値を算出してみます。

tapply関数で年毎の財政力指数の平均値を出し、それをnengotoという名前のオブジェクトとして保存します。

as.data.frame関数でnengotoをデータフレーム型に変換します。

rownames関数を使って、Yearという名前で年度を表す変数を追加します。

as.character関数でYearを文字列型に変換します。

substr関数でYearのはじめの4文字、つまり、1975とか1976だけを抜き取ります。

as.numeric関数で文字列を数値に変換します。

head関数ではじめの6行だけを表示しました。

こうして、nengotoという名前のデータフレームができましたから、plot関数で年毎に平均値がどのように推移しているか見てみましょう。

上に行ったり、下に行ったりしていますね。

table関数でYearの度数を確認しました。一番最新の年度が2016年度ですね。2016年度のデータで分析しましょう。

df1$Year == "2016年度"という方法で2016年度のみのデータフレームにしました。

財政力指数の高いところ、低いところを見てみましょう。

まず、Zai2016という名前のベクトルをつくり、そこの財政力指数のデータを入れます。そして、names関数で都道府県名の属性を加えます。

そして、sort関数で小さい順に並び替えました。島根県、高知県、鳥取県などが低く、東京都、愛知県、神奈川県などが財政力指数が高いです。

この財政力指数を人口(Pop)と面積(Area)で回帰分析してみましょう。

左の散布図が人口と財政力指数、右が面積と財政力指数です。人口とは正の相関のようです。面積はあんまり関係なさそうです。

lm関数で線形回帰モデルを作って分析します。

モデルのp-valueは1.829e-14なので0.05よりも小さく有意です。

Pop2016:Area2016は不要のようです。

新しいモデルを調べてみます。

こちらのモデルのほうが、p-valueは低く、Adjusted R-squaredは高いのでこちらのモデルのほうがいいですね。

Pop2016はp値あ3.77e-16で有意、Area2016も0.0172で有意です。

人口が多いほど、面積が狭いほど、財政力指数は高まるということですね。

人口密度が高いほど財政力指数は高いのでしょうね。検証してみましょう。

このモデル(model2)も有意ですが、model2のほうがいいですね。

model2の残差グラフを描いてみます。

栃木県、群馬県、静岡県は、モデルの推定する財政力指数よりも実際の指数のほうが高いのですね。

今回は以上です。

都道府県別の商店主の数の分析 - 商店主がこ1980年度から2015年度で136万6540人の商店主が減少した。

今日は都道府県別の商店主の数を分析しようと思います。

政府統計の総合窓口、e-Statからデータを取得しました。

地域は47都道府県です。

項目は総人口(人)、総面積(ha)、商店主(人)です。

これをCSVファイルにします。

こんな感じです。

このCSVファイルをR言語のread.csv関数で読込み、分析してみます。

read.csv関数でCSVファイルを読込み、na.omit関数でNAの行を削除しました。str関数でデータ構造を確認しています。376の観測、5の変数です。Year(年度)とPref(都道府県)はファクターで、Pop(総人口)、Area(面積)、Tenshu(商店主)は整数型です。

as.character関数でYearを文字列型に戻し、それをas.factor関数でファクター型に変換しました。これで、度数が0の年度がファクターの水準からなくなります。

table関数でYearの度数を数えました。1980年から2015年まで5年ごとにあります。すべて47度数あるので、都道府県の欠損は無いです。

tapply関数で、年毎の商店主の合計値を計算しました。1980年度は182万1140人だったのが、2015年度には45万4600人と4分の1になっています。

人口1万人あたりの商店主数、面積1万haあたりの商店主数を算出しました。そのあとでsummary関数で基本統計値をだしています。

平均で人口1万人当り100人の商店主が、面積1万ha当り588人の商店主がいる計算です。ただ、この値は1980年度のように商店主が多くいた時期も併せての平均ですので、2015年度だけだともっと少なくなるでしょうね。

tapply関数とmean関数で年毎の人口1万人当り、面積1万ha当りの商店主数の平均を算出しています。1980年度は人口1万人当り160人、面積1万ha当り940人の商店主がいましたが、2015年度は人口1万人当り42人、面積1万ha当り221人になっています。35年間で4分の1ですね。

2015年度だけのデータフレームを作成しました。df1$Year == "2015年度"として2015年度だけを選択しています。

order関数、rev関数を使って、並びかえました。

人口1万人当りの商店主の少ない県は、神奈川県、千葉県、埼玉県、東京都、愛知県、北海道です。

その反対に多いのは、高知県、和歌山県、島根県、宮崎県、長崎県、徳島県です。

面積1万ha当りの商店主が少ない都道府県は、北海道、岩手県、秋田県、島根県、福島県、山形県です。

多いのは、大阪府、東京都、神奈川県、埼玉県、愛知県、福岡県です。

df2015$Popなどと、いちいちデータフレームを指定するのが面倒なので、人口、面積、商店主のデータをそれぞれ独立したベクトルにしました。

plot関数で散布図を描きました。

商店主の数を被説明変数、人口と面積を説明変数にして回帰分析をしてみます。

モデルのp値は2.2e-16より小さいので有意です。P2015だけが係数としては有意ですね。I(log(A2015))をはずしたモデルをチェックしましょう。

Pr(>F)が0.8824なので、modelとmodel1は有意な差は無いです。なので、説明変数が少ないmodel1のほうがいいです。

さらに、P2015:A2015の交差項をはずしたmodel2をチェックしましょう。

A2015もいらなそうです。

I(log(P2015))もはずしてみましょう。

商店主の数は人口でほぼ説明ができるようです。

散布図を回帰直線を描きます。

年毎の商店主の数もグラフにしましょう。

tapply関数で計算した結果をbarplot関数で棒グラフにしました。

商店主の減少度合いがよくわかりますね。

最後にtapply関数の計算結果をデータフレーム型になおして、商店主の減少数を計算しました。136万6540人減少しました。

今回は以上です。