crosshyou

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

OECD Social spending data analysis 4 - Calculating Confidence Interval using R

UnsplashArda Demirkaynakが撮影した写真 

www.crosshyou.info

This post is following of above post. In the previous post, I made some visualizations with R ggplot2 package. In this post. In this post I will calculate confidence intervals.

Firstly, let's see average value of variables.

I see average priv_pc_gdp is 2.53, average pub_pc_gdp is 18.7 and average pub_usd_cap is 6848.

Let's see histograms with average value.

priv_pc_gdp

pub_pc_gdp

pub_usd_cap

I see priv_pc_gdp is the most skewed variable among those three.

To calclate confidence interval, we need standard error, which is standard deviation / sqrt(number of observations).

Let's caluclate them.

So, standard error of priv_pc_gdp is 0.0778, standard error of pub_pc_gdp is 0.172 and standard error of pub_usd_cap is 112.

All right, I can calculate confidence intarval. 95% confidence interval is calculated as followings. average value -/+ 1.96*standard error.

So, 95% confidence interval for priv_pc_gdp is [2.378, 2.682], 95% confidence interval for pub_pc_gdp is [18.36, 10.04] and 95% confidence interval for pub_usd _cap is [6628, 7067].

Let's add confidence intervals to previous histograms.

priv_pc_gdp:

pub_pc_gdp:

pub_usd_cap:

That's it. Thank you!

 

To read from the 1st post,

www.crosshyou.info

 

OECD Social spending data analysis 3 - Data Visualization with 5 Named Graphs (5NG) using R

UnsplashAlicia Steelsが撮影した写真 

 

www.crosshyou.info

This post is following of above post.
In the previous post, I made a dataframe for data analysis, named 'df2'.
Now, let's start data analysis with data visualization.
I will make 5 Named Graphs (5NG) using R.
My text book is

Chapter 2 Data Visualization | Statistical Inference via Data Science (moderndive.com)

Let's start with 5NG #1, Scatter Plot.

geom_point() makes scatter plot. I see relatively positive relationship between priv_pc_gdp and pub_pc_gdp.

I see relatively positive relationship between priv_pc_gdp and pub_usd_cap too.

I see there is strong correlation between pub_pc_gdp and pub_usd_cap.

Next, Line Graph.

geom_line() makes line graphs. prive_pc_gdp are increasing in general.

pub_pc_gdp has more variation than priv_pc_gdp.

pub_usd_cap are increasing for almost all countires.

 

#3 5NG is histogram.

geom_histogram() makes histograms. I use facet_wrap(~ continent) to make histograms for each continent.

I see Europe has higher percentage than other continents about pub_pc_gdp.

Sounth America has the lowest pub_usd_cap distribution.

The 4th 5NG is Boxplot.

geom_boxplot() makes a boxplot. I use mutate() and reorder() to reorder continent by average value. So, I can easily see Notrh America has the highest average priv_pc_gdp.

I see Europe has the higest average pub_pc_gdp and the highest medain too.

Europe has the highest average and median value for pub_usd_cap.

The last 5NG is Barplot.
To make a barplot, we can use geom_bar() and geom_col() function.

geom_bar() automatically count number of observations and make barplot.

Before using geom_col() function, we need to calculate summary statistics data such as number of observations.

That's it. Thank you!

Next post is

www.crosshyou.info

 

To read the 1st post,

www.crosshyou.info

 









 



















 

OECD Social spending data analysis 2 - Using filter(), select(), inner_join(), rename() function with R to make a dataframe to analyze.

UnsplashMilos Prelevicが撮影した写真 

www.crosshyou.info

This post is following of the above post. In the previous post, I load OECD Social spending data into R. I also load country ISO code and continent name data like below CSV file.

I get above data from List of Countries by Continent - StatisticsTimes.com

Let's use read_csv() function to load the CSV file data.

Using glimpse() function, I display continent dataframe.

Then, I use inner_join() function to merge 'df' and 'continent'.

Let' see if it is all right or not.

Good! I successfully merged 'df' and 'continent' dataframes.

To make a dataframe to analyze, I will check how many obervations each SUBJECT and MEASURE has by using table() function.

From the above table, I decide I will use PRIV & PC_GDP, PUB & PC_GDP and PUB & USD_CAP observations.

I make PRIVE & PC_GDP datagrame using filter() function, select() function and rename() function.

Let's see how 'priv_pc_gdp' looks like with glimpse() function.

I see AUS means Australia and it is in Oceania continent.

Next, I make PUB & PC_GDP datagrame.

Let's see how 'pub_pc_gdp' dataframe looks like with glimpse() functiion.

And I make PUB & USD_GDP dataframe.

Let's see 'pub_usd_cap' with glimpse() function.

Now, I have three dataframes, 'priv_pc_gdp', 'pub_pc_gdp' and 'pub_usd_cap'.
I will merge those datafames into one dataframe using inner_join() function.
LOCATION, TIME, name and Continent are common variable names.

Let' see how 'df2' looks like with glimpse() function.

Well done! 
I modifiy 'df2' dataframe with remane() function and select() function.

Let's see how 'df2' dataframe looks like with glimpse() funtion.

I change continent to factor.

Lastly, let's see summary statistics for 'df2' dataframe with summary() function.

I see there is no NA observations in 'df2' dataframe, there are 5 continent, Asia, Europe, North America, Oceania and South America. So there is not African countries in this dataframe.

That's it. Thank you!

Next post is

www.crosshyou.info

 

 

To read from the 1st post,

www.crosshyou.info

 

OECD Social spending data analysis 1 - Load CSV file data using R, read_csv() function.

UnsplashAlexander Schimmeckが撮影した写真 

In this post I will analyze OECD Social spending data using R.

OECD (2022), Social spending (indicator). doi: 10.1787/7497563b-en (Accessed on 26 November 2022)

This indicator is measured as a percentage of GDP or USD per capita.

From the web site, I downloaded CSV file as below.

Before loading the CSV file, I load some packages such as tidyverse and so on.

I loaded tidyvese, moderndive, infer and gridExtra packages. tidyverse is de-facto standard package in R for data science. I load moderndive package and infer package because I am learning Statistical Inference via Data Science (moderndive.com)

I loaded gridExtra package to display multiple plots at one screen.

Now, let's load the CSV file data with read_csv() function.

Then, using glimpse() function, let's check if data is correctly loaded.

Nice! The CSV file data is loaded succesfully into R.

There are 7 variables in the dataframe, let's see each variables one by one.

Location:

LOCATION is ISO country code, I see USA(United States of America) has the most observations.

 

INDICATOR:

INDICATOR has no variation, has just one value, SOCEXP. So, I will remove INDICATOR.

 

SUBJECT:

According to OECD web site below, PUB means Public, PRIV means Privae, PUBNET means Public net and TOTNET means Total net. We see PUB has the most observations.

 

MEASURE:

PC_GDP means a percentage of GDP, USD_CAP means ISD per capita.

 

FREQUENCY:

There is no variation in FREQUENCY, then I remove it.

 

TIME:

TIME means year, I see 2017 has the most observations, 190.

 

Value:

There is no NA in Value.

That's it. Thank you!

Next post is

www.crosshyou.info

 





 

都道府県別の被服及び履物費のデータの分析6 - R言語のplmパッケージでパネルデータ分析。First Difference, Fixed Effect, Random Effect Estimator.

UnsplashZoltan Tasiが撮影した写真 

www.crosshyou.info

の続きです。

前回はR言語のplmパッケージを使ってパネルデータフレームを作成し、model = "pooling" にして普通のクロスセクションで回帰分析をしました

今回は、First Differenced Estimator, Fixed Effect Estimator, Random Effect Estimatorの3つの方法でパネルデータを回帰分析してみます。

まずは、First Differenced Estimatorです。

First Differenced Estimator というのは、もともとの回帰分析モデルが、

yit = β0 + β1xit1 + β2xit2 + ... + βkxitk + ai + uit, 

t = 1, ... , T; i = 1, ..., n (tは時間、nは個体)

aiが個体i特有の要因

というモデルだという前提で、このaiを削除するために、

Δyit = yit - yit-1

       = Δβ1xit1 + Δβ2xit2 + ... + Δβkxitk + Δuit

t = 2, ... T ; i = 1, ... n 

と一つ前の時間との差分をとって(First Difference)、aiを削除して式で回帰分析をするものです。

plm()関数で回帰分析をするときに、model = "fd" とします。

結果をみてみます。

mitsudoの係数が-0.645と負の符号になっていますが、標準誤差が2.513ととても大きく、First Difference Estimatorでは、mitsudoの係数は有意な値は得られませんでした。

今度は、Fixed Effect Estimator でやります。

これは、前提となる回帰式モデルは、

yit = β0 + β1xit1 + β2xit2 + ... + βkxitk + ai + uit, 

t = 1, ... , T; i = 1, ..., n (tは時間、nは個体)、aiが個体i特有の要因

で同じなのですが、aiを取り除くために、それぞれの個体ごとに平均値をとります。

ybar = β0 + βbar1xbari1 + βbar2xbari2 + ,,, + βbarkxbarik + ai + ubari

そしてこの平均値と元の式の差分をとって、aiを除外するというものです。

plm()関数で実行するときは、model = "within" で実行できます。

結果をみてみましょう。

mitsudoの係数は1.551で有意な値です。

3番目は Random Effect Model です。これは勉強不足でまだ仕組みが良く分かっていないので説明は割愛します。plm()関数では、model = "random" で実行できます。

Pooling Cross Section, First Difference, Fixed Effect, Random Effect の4つの推計結果を並べて表示します。

Random Effect Estimator でのmitsudoの係数の推計結果は0.428で統計的に有意です。

どうやら、misudoが大きいほどwear_shoeが増えるのは間違いないようですね。まわりの人の目が大きいほど、被服及び履物には気を使うということでしょうか?笑

今回は以上です。

初めから読むには、

www.crosshyou.info

です。

都道府県別の被服及び履物費のデータの分析5 - R言語でパネルデータ分析。plmパッケージのplm()関数を使う

UnsplashAndrey Andreyevが撮影した写真 

www.crosshyou.info

の続きです。

今回は、R言語でパネルデータの分析をします。

こちらの本を参考にしてやってみます。

まず、パネルデータ分析用のデータのパッケージ、plmパッケージの読み込みをします。

パネルデータ分析用のデータパネル、pdata.frameを作ります。NAのあるwariai, percapita17, percapita23は除外します。

pdim()関数で作成したパネルデータの構造を確認します。

Balanced Panelとなっていますので、どの年も47都道府県のデータが揃っていることがわかります。n = 47となっています。47都道府県だとわかります。T = 33となっています。33年間のデータがあるとわかります。 N = 1551となっています。47*33 = 1551です。

sample_n()関数でランダムにデータを眺めてみます。

ibaaraki-1984は、1984年の茨城県、Oita-1990は1990年の大分県という意味ですね。

まずは、pooling cross section ということでパネルデータと意識せず、普通の回帰分析として分析してみます。被説明変数はwear_shoeで、説明変数は、mitsudoとyearをファクター型に変換したものとbig6にしてみます。plm()関数をつかってmodelをpoolingにすると、普通の回帰分析になります。

結果を表示するのに便利なstargazerパッケージの読み込みをしておきます。

stargazer()関数で結果を表示します。yearのファクター型に変換したものは32もありますので、mitsudoとbig6の係数だけを表示します。keeep = のところで表示したい変数だけにすることができます。big6, 東京都、大阪府、愛知県、神奈川県、千葉県、埼玉県だと他の都道府県よりも1172ぐらいwear_shoeが多いです。

mitsudoの係数は、0.139で統計的に有意に0ではないです。mistudoが1増えると、wear_shoeが0.139増えるということですね。

今回は以上です。

次回は

www.crosshyou.info

です。

 

初めから読むには、

www.crosshyou.info

です。

都道府県別の被服及び履物費のデータの分析4 - R言語で回帰分析。カテゴリーデータを説明変数に加える。重回帰分析。

UnsplashGwen Weustinkが撮影した写真 

www.crosshyou.info

の続きです。

前回はwear_shoeを被説明変数、wariaiを説明変数にして単回帰分析をしました。

今回はもう一つ説明変数を加えてみます。カテゴリーデータを加えてみましょう。

まず、カテゴリーデータを作成します。

nosea: 海が無い県は1、ある県は0というダミー変数からカテゴリー変数を作成します。

mutate()関数の中でifelse()関数とas.factor()関数を使って、no_sea, has_seaという2つの値を取るカテゴリカル変数を作成しました。

glimpse()関数でちゃんとできているか確認します。

select()関数、pull()関数とtable()関数を使って数を確認します。

has_seaが1287個、no_seaが264個です。

それでは、このseaを説明変数に加えて回帰分析をしてみます。

wariai*sea と * で結んでいるので、このモデルの式は、

wear_shoe = beta0 + beta1*wariai + beta2*sea + beta3*wariai*sea + u

という式です。

結果をget_regression_table()関数でみてみます。

これは、どういう式になるかというと、

seaがno_seaのときは、

wear_shoe = 1800 + 181*wariai -23626 + 380*wariai + u

                  = -26826 + 561*wariai + u

seaがhas_seaのときは、

wear_shoe = 1800 + 181*wariai + u

ということです。no_seaのときのほうが切片は小さく、傾きは大きいです。

散布図にしてみます。

散布図は確かにno_seaの緑色の直線のほうが切片が小さく、傾きは大きいです。

もうひとつ、交差項の無い回帰モデルを分析をしてみます。

こんどは、wariaiとseaを*ではなくて+で結んでいます。

これは、式で表すと、

wear_shoe = beta0 + beta1*wariai + beta2*sea + u

です。

結果みてみます。

式にしてみます。

seaがno_seaのときは、

wear_shoe = -1262 + 230*wariai + 967 + u

                  = -295 + 230*wariai + u

seaがhas_seaのときは、

wear_shoe = -1262 + 230*wariai + u

です。傾きは230で同じ、切片だけが違います。

moderndiveパッケージのgeom_parallel_slopes()関数をつかって散布図にこの回帰直線を重ねてみます。

前回の単回帰分析モデルとあわせて3つのモデルを作成しました。

この3つを比較してみます。moderndiveパッケージのget_regression_summaries()関数を使います。

adj_r_squaredの値を比較すると、lm_wariai_interactionモデルが一番、値が大きいです。

この値が大きいほうがよいので、交差項の入ったモデルが一番いいことがわかりました。

今回は以上です。

 

次回は

www.crosshyou.info

です。

 

初めから読むには、

www.crosshyou.info

です。