crosshyou

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

都道府県別の教育費のデータの分析2 - Rでデータをグラフにする。ggplot2パッケージのgeom_point(), geom_boxplot(), geom_line(), geom_histogram()

Photo by Antonio Sessa on Unsplash 

www.crosshyou.info

の続きです。

前回はデータをRに読み込みました。

今回はそのデータをいろいろグラフにしてみます。

まずは、都道府県別のed: 教育費(都道府県財政と市町村財政の合計)です。

mutate()関数の中でreorder()関数を使い、prefをedの値、中央値が大きい順に並び替えています。そして、ggplot()関数とgeom_pouint()関数でグラフを描き、coord_flip()関数でX軸とY軸を反転させ、theme_bw()関数でグラフの見た目を整えています。

東京都がダントツの金額規模です。鳥取県が一番少ないです。

edをlog()関数で対数変換してみます。

鳥取県はダントツで少ないです。

次は都道府県別のp_ratio: 都道府県財政の教育費の占める割合を見てみます。

東京都が一番p_ratioが低いのですね。高知県や和歌山県などがp_ratioが高いです。

次は、年別のedを見てみます。

ggplot(aes(year, ed, group = year))と yearでグループ化してから、geom_boxplot()関数で年別の箱ひげ図を作りました。全体としては横ばい傾向なのでしょうか?

折れ線グラフで年ごとの推移を見てみましょう。

ggpot(aes(year, ed, group = pref)) とgroup = pref で都道府県別にグループ化してから、

geom_line()関数で折れ線グラフにしました。やはり、横ばいという感じですね。

年別のp_ratioも見てみます。

p_ratioのほうが年ごとに違いがありそうですね。

折れ線グラフも描いてみましょう。

やはりp_ratioのほうが動きが激しい感じですね。

こんどはヒストグラムを描いてみます。

geom_histogram()関数でヒストグラムを描きます。上に青い丸の追加したいますが、3つの山があるようです。

対数に変換したヒストグラムも描いてみます。

対数変換した ed: 教育費でも山は3つあります。

p_ratioのヒストグラムを描きます。

p_ratioは左右対称の正規分布のようなヒストグラムですね。

今回は以上です。

今回はggplot2パッケージのgeom_point(), geom_line(), geom_boxplot(), geom_histogram()を使ってデータをグラフにしました。

初めから読むには、

 

www.crosshyou.info

です。

都道府県別の教育費のデータの分析1 - Rにデータを読み込ませる。

Photo by Vadym Chumak on Unsplash 

今回は、都道府県別の教育費のデータを分析しようと思います。

データは、政府統計の総合窓口(www.e-stat.go.jp)のウェブサイトから取得します。

まずは、47の都道府県を選択します。

教育費は都道府県財政の値と市町村財政の2つありました。

その他に人口、可住地面積、県内総生産額も取得します。

このようなCSVファイルがサイトからダウンロードできます。

このデータをRに読み込んで分析してみます。

まず、tidyverseというパッケージの読み込みをしておきます。

そうしたら、read_csv()関数を使ってデータを読み込みます。

head()関数を使って読み込んだデータのはじめの数行を表示させてみます。

nenとyearは、実質同じ情報なので、nenを削除します。

na.omit()関数をつかって NA のある行を削除します。

yearをはじめの4文字だけにします。str_sub()関数で1文字目から4文字目までを抽出して、as.numeric()関数で数値型にします。

ed_p: 教育費(都道府県財政)とed_c: 教育費(市町村財政)を合計した変数、edを作成します。

ed_p/edでp_ratioという変数を作ります。

最小が0.5869ということは、どの都道府県でも都道府県財政の教育費のほうが市町村財政の教育費よりも大きいのですね。

gdp: 県内総生産額は百万円単位、教育費は千円単位と単位が違いますので、どちらも1億円単位にして統一します。

pop: 人口は1人単位です。少しgdpやedなどと比べると値が大きいので、千人単位になおします。

人口の一番少ないところは、57万7千人、一番多いところ、東京都ですが1339万9千人です。

area: 可住面積はha単位です。これも値が大きいので、百ha単位に変換します。

一番小さいところは大阪府で、一番大きなところは北海道です。

prefのデータ型をファクター型に変換します。

summary()関数でデータフレームのサマリーを表示してみます。

summary()関数で表示すると、pref: 都道府県名が文字化けしていますが、head()関数では問題なく表示されていたので、よしとしましょう。

今回は以上です。

次回は

 

www.crosshyou.info

です。

 

都道府県別の建設総合統計のデータの分析2 - Rの pivot_longer()関数で 横長型のデータフレームを縦長型に変換

Photo by Drew Bae on Unsplash 

www.crosshyou.info

前回はCSVファイルのデータを取り込みました。

今回は取り込んだデータをいろいろと調べてみようと思います。

はじめに前回作成したデータフレームがどんなものだったか確認します。

横にhokkaido, aomori, iwate, miyagi ~~~ と並んでいるデータですね。

これを、都道府県名が縦に並ぶように、

現在の姿

これを、下の図のように変換したいです。

これは、pivot_longer()という関数でできます。

pivot_longer(データフレーム, 変換したい列名, names_to = 変換したい列名が値に入る新しい列名, values_to = 変換したい列名にあった値が入る新しい列名)

という構文です。ちょっとわかりにくいですが、実際にやってみると簡単です。

変換したい列名は、hokkaidoからokinawaで、これはprefという新しい列に入れて、それぞれの値をkingakuという新しい列名にいれます。

できました!

minkouがboth, shuruiがboth_totalのデータだけでとりあえず分析してみましょう。

今回は以上です。

初めから読むには、

 

www.crosshyou.info

です。

 

 

 

都道府県別の建設総合統計のデータの分析1 - CSVのデータファイルをRに取り込む

Photo by Mikita Yo on Unsplash 

いつものように、政府統計の総合窓口(www.e-stat.go.jp)を閲覧すると、建設総合統計のデータベースが更新されているようでした。

建設総合統計は、国内の建設活動を出来高ベースで把握することを目的とした加工統計とのことです。

クリックすると、下の画像になりました。

一番下にある、年度次[9件]をクリックしてみます。

年度表 <<表-21>> 総合表 都道府県別・種類別-年度別工事費

というのをクリックしてみます。

このようなデータでした。これをダウンロードします。

このようなCSVファイルでした。

これを眺めて、少し整理整頓したのが下の画像です。

このCSVファイルをRに読み込んで、分析しようと思います。

まず、tidyverseパッケージの読み込みをしておきます。

read_csv()関数でCSVファイルを読み込みます。

minkou: 民間か公共か

shurui: 建設の種類

の2つの変数のデータ型をファクター型に変換します。

summary()関数で読み込んだデータを確認してみます。

こんな感じです。yearを見るとデータは、2019年、2020年、2021年の3年間のようです。

各変数にNAが含まれているかどうか調べます。

na_countという名前のNAを数える関数を作って、apply()関数で各列について調べました。どうやらNAは無いようです。

今回は以上です。

次回は

 

www.crosshyou.info

です。

 

OECD Young self-employed data analysis 6 - panel data analysis using R - first differenced estimator using plm package

Photo by Daniel Olah on Unsplash 

www.crosshyou.info

This post is following of above post.

In this post, I will use First Differenced Estimator to estimate capi10K effect for men.

The background model is below

men = beta_0 + beta_1 * capi10K + beta_2 * y15 + a + u

a is LOCATION specific error it does not change in 2012 and 2015.

u is unobserved factor.

To remove a, we make 2015 - 2102

2015 model is men = beta_0 + beta_1 * capi10K + beta_2 + a + u

2012 model is men = beta_0 + beta_1 * capi10K +                a + u

So, 2015 - 2012 is

men2015 - men2012 = beta_0_new + beta_1 * (capi10K2015 - capi10K2012) + u_new.

beta_0_new and u_new can be expressed beta_0 and u because we only care about beta_1 and we can remove a; LOCATION specific factor.

So, First Differenced Estimator is

men2015 - men2012 = beta_0 + beta_1 *(capi10K2015 - capi10K2012) + u

Let's do it with R.

Firstly, I makde 2012 data frame and 2015 data frame.

Then, I merge the two data frame and make men2015 - men2012 and capi10K2015 - capi10K2012.

All right, let's make regression analysis.

coefficient of capi10K_diff is -0.7854, it has smaller impact than previous post's model_1 and model_2. And p-value is greater than 0.05.

So, I cannot say capi10K is significant for men.

Let's draw a sactter plot.

Next. I use plm package. It is much convenient for panel data.

Then, I make panel data frame using pdata.frame() function.

And I use plm() function with model = "fd" for First Differenced Estimator.

capi10K coefficient -0.78541 is the exactly same as model_diff coefficient.

That's it. Thank you!
To read the 1st post,

 

www.crosshyou.info

 








 

OECD Young self-employed data analysis 5 - panel data analysis using R - basic pooling cross section regression

Photo by Colin Watts on Unsplash 

www.crosshyou.info

This post is following of above post.

In the above post, I made panel data set. 
Let's analyze with the panel data.

Firstly, I male year dummy variable.

y15 is 1 when TIME is 2015 and 0 when TIME is 2012.

Then, I made pooling cross section regression.

My first model is men = beta_0 + beta_1 * capita + beta_2 * y15 + u

capita is in 1 USD, so, model_1 suggest that 1 USD increase of capita makes men lower by -0.0002222655%. It is a bit difficult to understand.
So, I made capi10K, that is capita in 10K, 10,000 USD.

The summary() function result is more easy to understand.

10,000 USD per capita increase means -2.1117 percent point decrese for men young self-empolyed.

Let's make a scatter plot.

I don't see large difference between 2012 and 2015.

Next, let's add interaction term.

men = beta_0 + beta_1 * capi10K + beta_2 * y15 + beta_3 * capi10K * y15 + u

So, when TIME is 2015, y15 = 1, so the model is

men = 16.6451 - 0.4777 + (-2.2107 + 0.1966) * capi10K + u

When TIME is 2012, t15 = 0, the model is

men = 16.6451 -2.2107 * capi10K + u

Let's draw a scatter plot.

Next, I use 1/capi10K instead of capi10K. 

men = beta_0 + beta_1 * (1/capi10K) + beta_2 * y15 + u

I(1/capi10K) coefficient is 18.7. so 1 inclease of 1/capi10K makes 18.7 increase of men.

Let's make a sacatter plot

Finally, let's compare three models with stargazer package's stargazer() function.

That's it. Thank you!

Next post is

 

www.crosshyou.info

 


To read the 1st post,

 

www.crosshyou.info

 

読書記録 - 「日本の品種はすごい - うまい植物をめぐる物語」 竹下大学 著 中公新書

 

ジャガイモ、ナシ、リンゴ、ダイズ、カブ、ダイコン、ワサビについての品種のお話。

美味しい品種を作ろうと頑張る人たちに感謝です。