www.crosshyou.info

政府統計の総合窓口のデータや、OECDやUCIやのデータを使って、Rの練習をしています。ときどき、読書記録も載せています。

OECD Discriminatory family code data analysis 2 - Making a histogram, a boxplot and an ECDF plot with R

Photo by Kentaro Toma on Unsplash 

www.crosshyou.info

This post is following of the above post.
In this post, I will make histograms, boxplots and ECDF plots with R.

Before making those plots, I made some changes to the dataframe.

I changed variable names with colnames() function.

I changed subject values to "em" and "atwm" with ifelse() function.

"em" stands for "Early Marriage".

"atwm" stands for "Attitudes Towards Working Mothers".

So, when "em" is high, "Early Marriage" ratio is high, when "atwm" is high, a country is kind to working mothers.

Let's make histograms.

First, I use ggplot() function and geom_histogram() function.

I also use hist() function and dencity() function and lines() function.

Then, I made boxplots.

Above is using ggplot() and geom_boxplot() function.
Below iss using boxplot() function

I will make ECDF plot with ggplot() and stat_ecdf() function.

I can make ECDF plot with plot() and ecdf() function.

Let's see which country is higer/lower Early Marriage, Attitudes Towards Working Mothers.

JOR has the highest Attitudes Towards Working Mothers.

Let's do it with basic functions.

ISL has the lowest Attitudes Towards Working Mothers.

Let's do it with basic functions.

Which country has the higest Early Marriage?

NER has the higest Early Marriage.

Let's do it with basic functions.

Which country has the loweat Early Marriage?

LTU has 0 value. I assume LTU has ban for Early Marriage.

Let's do it with basic functions, order(), subset() function.

That's it. Thank you!

Next post is

 

www.crosshyou.info

 



To see the first post,

 

www.crosshyou.info

 






















 

読書記録 - 「物理学はいかに創られたか - 初期の観念から相対性理論及び量子論への思想の発展」 (上巻・下巻 ) アインシュタイン、インフェルト 著 岩波新書

 

抽象的な話が多く、よくわからなかったが、下巻の最後のほうに

「私たちの理論的の構成によってこの実在を把握することが可能であるという信念がなくては、また私たちの世界が内的調和をもっているという信念を欠いては、科学はまるであり得ないでしょう。この信念こそは、あらゆる科学的創造に対する根本的な動機なのでありますし、またいつになってもそうなのでありましょう。」

という記述があって、これは何となくわかりました。

物理学はこの世界を理論的の構成によって把握することの営みなのだと感じました。

 

 

OECD Discriminatory family code data analysis 1 - load CSV file with read_csv() fundtion and display dataframe summary with summary() function in R.

Photo by Redd on Unsplash 

I will analyize ODEC Discriminatory family code.

Inequality - Discriminatory family code - OECD Data

I downloaded CSV file likde below from aboce web site.

Let's analyze with R.

Before load the CSV file, I load tidyverse package.

Then, I use read_csv() function to load the CSV file data into R.

There are 7 variables. Let's check it one by one.

Let's begin with LOCATION. LOCATION is country code.

I see each location has 1 or 2 counts.

Next, INDICATOR

INDICATOR has only one value; DISCRIMFAMCODE, so we can remove this varibale from the data frame.

 

Next, SUBJECT

SUBJECT has two value, ATTWORKMUM and EARMARRIAGE.
SInce there are only 2 kinds of value, I change SUBJECT to factor type.

 

Next, MEASURE

MEASURE has only one value, PC. PC means percentage. I remove MEASURE.

 

Next, FREQUENCY

FREQUENCY has only one value, A. A means annual.

I remove FREQUENCY.

 

Next, TIME

TIME is year. I see min, max, median, max, and 1st quantile, 3rd quintile are all 2019.

I remove TIME.

Last variable is Value.

There is not NA in Value.

All right, let's use summary() function for the df.

I see the df have 243 observations, two kinds of subject, one is ATTWORKMUM and the other is EARMARRIAGE, Value measn those SUBJECT's value.

That's it.
Thank you!

Next post is

 

www.crosshyou.info

 

 

小売物価統計調査のデータ分析5 - R言語のlm()関数で重回帰分析をする。そして、scatterplot3d()関数で3D散布図を描く

f:id:cross_hyou:20220416072738j:plain

Photo by Ian Parker on Unsplash 

www.crosshyou.info

の続きです。

今回は、R言語のlm()関数で重回帰分析をしてみます。

前回までは、sougou: 総合をhouse: 住居、utility: 水道・光熱費とそれぞれ一つの説明変数で回帰分析していました。今回は、house: 住居と utility: 水道・光熱費を同時に説明変数にして回帰分析します。

f:id:cross_hyou:20220416073536p:plain

utility: 水道・光熱費は単回帰分析のときは有意な説明変数ではありませんでしたが、house: 住居と一緒に重回帰分析をすると有意な説明変数になるのですね。

前回と同じように、SST(Total Sum of Squares), SSE(Explained Sum of Squares), SSR(Residual Sum of Squares), R2(R-Squared)を計算しました。

f:id:cross_hyou:20220416074109p:plain

R2が0.688と算出されましたが、summary()関数での出力結果のR2と同じ値になりました。

f:id:cross_hyou:20220416074548p:plain

続いて、SER(Standard Error of the Regression)を計算しました。

f:id:cross_hyou:20220416074604p:plain

SERは0.884と計算されました。これもsummary()関数の出力結果と一致しました。

f:id:cross_hyou:20220416074746p:plain

Adjusted R-squaredを計算してみます。Adjusted R-squaredは説明変数の数を調整したR2です。

f:id:cross_hyou:20220416075439p:plain

0.6867と算出されました。これもsummary()関数の出力結果と一致しました。

f:id:cross_hyou:20220416075607p:plain

次は、それぞれの説明変数のStandard Error を計算してみます。

house: 住居のStandard Error から計算してみます。

f:id:cross_hyou:20220416075843p:plain

単回帰分析のStandard Error と違うのは、house: 住居をその他の説明変数(今回はutility: 水道・光熱費だけですが)で回帰分析をしてそのR2を利用することですね。分子のほうにこのR2が入っていて、1 - R2_house となっていますから、このR2が大きくなればなるほど、つまり対象の説明変数がその他の説明変数と相関が大きくなればなるほどStandard Error の値も大きくなる、ということですね。0.004817 と算出されました。

これも、summary()関数の出力結果と一致しました。

f:id:cross_hyou:20220416080507p:plain

utility: 水道・光熱費のStandard Error も計算してみます。

f:id:cross_hyou:20220416080640p:plain

0.00839 と算出されました。これもsummry()関数の出力結果と一致します。

f:id:cross_hyou:20220416080800p:plain

回帰分析の特徴である、残差の平均値が 0 であることを確認します。

f:id:cross_hyou:20220416081215p:plain

次は、残差が Homoskedasticity であるかどうかを確認します。残差の2乗が house: 住居、utility: 水道・光熱費と関係なければ Homoskedasticity です。日本語だと、等分散、分散の均一性といいますね。

f:id:cross_hyou:20220416081710p:plain

house: 住居のp値が0.0122と0.05 よりも小さく、回帰式全体のp値も0.04282 と0.05よりも小さいので、Homoskedasticity ではなくて、Heteroskedasticity ですね。

なので、通常のStandard Error よりも Heteroskedasticity-Robust Standard Error のほうがより適切です。carパッケージとlmtestパッケージを読み込んで、Heteroskedasticity Robust Standard Error を求めてみます。

f:id:cross_hyou:20220416082056p:plain

coeftest()関数で Standard Error を出力します。

f:id:cross_hyou:20220416082330p:plain

上のものが通常の Standard Error で、下のものが Heteroskedasticity-Robust Standard Error です。house: 住居のほうはHeteroskedasticity-Robust のほうが小さくなって、utility: 水道・光熱費のほうは大きくなっていますが、p値の結論はかわらず、どちらの説明変数も有意な説明変数です。

最後に3Dの散布図を描いてみましょう。scatterplot3dパッケージのscatterplot3d()関数を使うのが簡単です。

f:id:cross_hyou:20220416085819p:plain

f:id:cross_hyou:20220416085758p:plain

今回は以上です。

初めから読むには、

 

www.crosshyou.info

です。

HistDataパッケージのGuerry

f:id:cross_hyou:20220409220544j:plain

Photo by Boris Ho on Unsplash 

HistDataパッケージのGuerryのデータは、Guerryという人が1833年に世界で初めて犯罪や自殺など社会学的なデータを組織的に集めて分析したデータです。

データを呼び出してみます。

f:id:cross_hyou:20220409220927p:plain

86の観測と23の変数があります。

ヘルプのコードは、

data(Guerry)
## maybe str(Guerry) ; plot(Guerry) ... 
と書いてあるだけで特にないです。

なので、数値型の変数のヒストグラムを描いてみました。

f:id:cross_hyou:20220409221220p:plain

f:id:cross_hyou:20220409221231p:plain

ある決まった数の分布変数もあれば山型の分布の変数もあります。

下記が今回のコードです。

# HistData_Guerry
#
library(HistData)
data(Guerry)
str(Guerry)
#
par(mfrow = c(4, 5))
hist(Guerry$dept)
hist(Guerry$Crime_pers)
hist(Guerry$Crime_prop)
hist(Guerry$Literacy)
hist(Guerry$Donations)
hist(Guerry$Infants)
hist(Guerry$Suicides)
hist(Guerry$Wealth)
hist(Guerry$Commerce)
hist(Guerry$Clergy)
hist(Guerry$Crime_parents)
hist(Guerry$Infanticide)
hist(Guerry$Donation_clergy)
hist(Guerry$Lottery)
hist(Guerry$Instruction)
hist(Guerry$Prostitutes)
hist(Guerry$Distance)
hist(Guerry$Area)
hist(Guerry$Pop1831)
par(mfrow = c(1, 1))
#

小売物価統計調査のデータ分析4- OLSの3つの特性とSST(total sum of squares), SSE(explained sum of squares), SSR(residual sum of squares), SER(standard error of the regression), 説明変数の標準誤差

f:id:cross_hyou:20220409155549j:plain

Photo by Subtle Cinematics on Unsplash 

www.crosshyou.info

の続きです。

今回は、

Introductory Econometorics: A Modef\rn Approace 7e, by Jeffrey M. Wooldridgeを参考にしてOLSの特性を確認します。

2-3b Algebraic Properties of OLS Statisticsの部分で説明されている3つの特性を確認します。

一つ目は、

残差の合計が0、別の言い方で言えば、残差の平均が0ということです。

残差はresid()関数で求まりますので簡単に確認できます。

f:id:cross_hyou:20220409160439p:plain

0になっています。

2番目は、説明変数と残差の共分散は0、別の言い方で言えば、説明変数と残差を掛け算したものの合計は0、ということです。

f:id:cross_hyou:20220409161234p:plain

cov()関数で共分散を、sum()関数で合計を計算していますが、どちらの方法でも0になっています。

3番目は、被説明変数の平均値と説明変数の平均値は、OLSの回帰直線の上にある、ということです。

これはグラフにして確認しましょう。

f:id:cross_hyou:20220409162455p:plain

f:id:cross_hyou:20220409162537p:plain

赤丸、緑丸が平均値の位置です。回帰線の上に載っていることがわかります。

続いて、SST, SSE, SSRを計算してみます。

ここからは、houseを説明変数にした回帰分析モデルだけで計算します。

SSTはtotal sum of squaresというもので、それぞれの説明変数と平均値の差を2乗して合計したものです。

計算してみます。

f:id:cross_hyou:20220409163219p:plain

SSTは1155.717と算出されました。

SSRはexplained sum of squaresというもので、OLSで算出される予想値と平均値の差を2乗して合計したものです。予想値は、fitted()関数で算出できます。

f:id:cross_hyou:20220409163658p:plain

SSEは724.9415と算出されました。

SSRはresidual sum of squaresというもので、残差の2乗の合計です。残差はresid()関数で算出できます。

f:id:cross_hyou:20220409163926p:plain

SSRは430.7759と算出されました。

そして、このSST, SSE, SSRの関係は、

SST = SSE + SSR 

という関係が成り立つということです。

SSE + SSRがSSTと同じか、確認します。

f:id:cross_hyou:20220409164205p:plain

SST = 1155.711, SSE + SSR = 1155.7171 と一致しましたね。

そして、回帰分析のR2(決定係数)は

R2 = SSE/SST = 1 - SSR/SST

で計算されます。R2は被説明変数のバらツキ度合いのうち、どのくらいの比率がSSE、つまり回帰分析の式で説明できるのかを表しています。

R2を計算してみます。

f:id:cross_hyou:20220409164646p:plain

0R2は、0.6272654と算出されました。summary()関数で出力されていたR2と同じです。

f:id:cross_hyou:20220409164935p:plain

続いて、SER(standard error of regression)を計算しましょう。

SERは回帰モデルの誤差項の標準偏差の推定値です。

sqrt(SSR/(n - 2))

で求めることができます。nは観測数です。

f:id:cross_hyou:20220409165845p:plain

SERは0.9656165と算出されました。summary()関数の出力結果と同じですね。

f:id:cross_hyou:20220409170055p:plain

SERが算出できたので、houseのstandard error(標準誤差)を求めることができます。

SER/sqrt(SST_house)

で計算できます。
SST_houseはそれぞれのhouseの値と平均値の値を2乗したものの合計です。

f:id:cross_hyou:20220409171713p:plain

houseの標準誤差は、0.004922938と算出できました。summary()関数の出力結果と同じです。

f:id:cross_hyou:20220409171938p:plain

houseの係数を標準誤差で割り算すればt値が算出されます。

coef()関数でOLSの係数は取り出せます。

f:id:cross_hyou:20220409172400p:plain

houseのt値は27.88347となりました。summary()関数の出力結果と一致します。

f:id:cross_hyou:20220409172625p:plain

t値が27.88と2よりも遥かに大きい値なので、houseの係数のp値は0です。

f:id:cross_hyou:20220409174117p:plain

(Intercept)のStandard Errorも算出します。

f:id:cross_hyou:20220410170023p:plain

f:id:cross_hyou:20220410170044p:plain

今回は以上です。

次回は、

 

www.crosshyou.info

です。

初めから読むには、

 

www.crosshyou.info

です。

 









 

小売物価統計調査のデータ分析3 - R言語のlm()関数で単回帰分析。マニュアル計算でも切片と傾きを求める。

f:id:cross_hyou:20220408212458j:plain

Photo by Simon Berger on Unsplash 

www.crosshyou.info

の続きです。

前回はいろいろなグラフを作成しました。今回は分析的なことをしてみましょう。

まず、sougou: 総合物価指数と他の種類の物価の相関を調べてみます。

f:id:cross_hyou:20220408212946p:plain

一番相関が無いのは、utility: 光熱・水道で、一番相関があるのは、ex_rent: 家賃を除く総合です。家賃を除く総合が総合と相関が高いのは当りまえなので、実質的に一番相関が高いのは、house: 住居です。

そこで、utility: 光熱・水道とhouse: 住居とsougou: 相互の3つについて調べてみましょう。

house: 住居とsougou: 総合の散布図を描いてみます。

f:id:cross_hyou:20220408213503p:plain

f:id:cross_hyou:20220408213514p:plain

house: 住居とsougou: は高い相関があることがわかります。

utility: 光熱・水道とsougou: の散布図を描いてみます。

f:id:cross_hyou:20220408213818p:plain

f:id:cross_hyou:20220408213829p:plain

こちらは相関が無いことがわかります。

次にlm()関数で回帰分析をしてみます。

f:id:cross_hyou:20220408214000p:plain

sougou: 総合を被説明変数に、house: 住居を説明変数にして回帰分析をしました。

house: 住居の係数は0.137でp値はほとんど0なので有意な係数です。

house: 住居の値が1増えると、sougou: 総合の値が0.137増える、ということです。

次は、説明変数にutility: 光熱・水道を使って回帰分析をしてみます。

f:id:cross_hyou:20220408214257p:plain

utility: 光熱・水道の係数は-0.015とかなり小さく、p値は0.299と0.05よりも大きいので有意な係数ではないです。utility: 光熱・水道の水準とは無関係にsougou: 総合の値は決まっているようです。

lm()関数を使えば、回帰モデルの切片や傾きは自動的に計算してくれます。

これをマニュアルで算出してみます。

説明変数が1つの場合の単回帰モデルでは、傾きの式は、

傾き = yとxの共分散 / xの分散

です。

house: 住居を説明変数にしたときの傾きを算出してみます。

f:id:cross_hyou:20220408215548p:plain

0.137となりました。これは、lm()関数で求めた傾きと一致しています。

切片の式は、

切片 = yの平均値 - 傾き*xの平均値

です。

傾きは0.1372686とわかりましたので、切片を算出してみます。

f:id:cross_hyou:20220408220048p:plain

86.25となりました。これは、lm()関数で求めた切片と一致しています。

今回は以上です。

次回は

 

www.crosshyou.info

です。

初めから読むには、

 

www.crosshyou.info

です。