www.crosshyou.info

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

都道府県別の被服及び履物費のデータの分析3 - R言語で回帰分析。まずは単回帰分析。15~64歳の人口割合を説明変数にする。

UnsplashBoris Smokrovicが撮影した写真 

www.crosshyou.info

の続きです。

今回はR言語で回帰分析をしてみます。

被説明変数は、wear_shoe: 被服及び履物費です。回帰分析をはじめる前にwear_shoeとその他の説明変数の候補との相関関係を確認しておきましょう。

一番相関係数が大きいのは、percapita23: 一人当たりの県民所得(平成23年基準)ですね。

これらの散布図も描いてみます。

まずは、上のようにそれぞれの散布図をオブジェクトとして保存しておきます。

そして、gridExtraパッケージのgrid.arrange()関数をつかって9つの散布図を同時に表示します。

こうやって散布図にすると、wear_shoeとの関係性がよくわかりますね。

wariai: 15~64歳の人口割合、を説明変数にして回帰分析してみます。説明変数が一つなので、単回帰分析です。

lm()関数で回帰分析のオブジェクトを作ります。

回帰分析の結果を見るのに、通常はsummary()関数を使いますが、今回はmoderndiveパッケージのget_regression_table()関数を使います。

式で表すと、

wear_shoe = -3283 + 264 * wariai + u

となります。wariaiが1増えると264だけwear_shoeが増えるという式です。

wear_shoeとwariaiの散布図に回帰直線と重ねてみます。

geom_abline()関数で直線を引きました。

lm_wariai$coefficient[1] で切片、lm_wariai$coefficient[2]で傾きを取り出せます。

今回は以上です。

次回は、

www.crosshyou.info

です。

 

はじめから読むには、

www.crosshyou.info

です。

 

読書記録 - 「生物学探偵セオ・クレイ 街の狩人」 アンドリュー・メイン著 ハヤカワ・ミステリ文庫

生物学探偵セオ・クレイのシリーズ2冊目です。探偵というよりは、シリアルキラー・ハンターです。街の狩人というタイトルのとおり今回はロスアンゼルスやアトランタという街が舞台です。

セオが科学的な推論と知識を使って今回の標的を追い詰めていきます。

 

都道府県別の被服及び履物費のデータの分析2 - R言語でグラフを描く。The Five Named Graphs でデータを視覚化する。

UnsplashPierre Van Crombruggheが撮影した写真 

www.crosshyou.info

の続きです。

今回は、Chapter 2 Data Visualization | Statistical Inference via Data Science (moderndive.com) を参考にしてR言語のggplot2パッケージを使っていくつかグラフを描きます。こちらのサイトには、The Five Named Graphsと称して5種類のグラフが紹介されています。ひとつずつ、実行してみます。

1つ目は散布図です。これはggplot() + geom_point()関数で作成できます。

wear_shoeをY軸に、mitsudoをX軸にした散布図です。両者はあまり関連なさそうです。

geom_point()の中に、alpha = 0.2 として点を薄くして、theme_bw()を追加して背景を白にしました。

Y軸をwear_shoe, X軸をpercapita17にした散布図です。filter()関数を使ってpercapita17がNAの行を除外してからggplot()関数とgeom_point()関数を使いました。percapita17の値が大きいほうがwear_shoeの値も大きい傾向があるようです。

2番目のグラフはライングラフです。ggplot()関数とgeom_line()関数を使います。

ggplot()の中で、group = pref としているので都道府県ごとのライングラフになります。最近のwear_shoeの値は低下傾向ですね。

3番目のグラフは、ヒストグラムです。ggplot()関数とgeom_histogram()関数を使います。

wear_shoeのヒストグラムです。左右対称な感じの分布ですね。bins = 25 としているので25個の棒があります。

4番目のグラフは箱ひげ図です。ggplot()関数とgeom_boxplot()関数を使います。

group = year として年ごとに箱ひげ図を描きました。1991年、1992年頃が一番金額が多いですね。

そして、5番目のグラフは棒グラフです。ggplot()関数とgeom_col()関数を使いました。

group_by()関数でyearのグループにしてから、summarize()関数の中でmean()関数を使いyearごとのwear_shoeの平均値を算出してグラフにしました。

以上の5種類のグラフから、

wear_shoeは1975年から1992年くらいまでは上昇していて、その後は低下傾向になっている、分布は左右対称の正規分布のような分布であることがわかりました。

今回は以上です。

次回は、

www.crosshyou.info

です。

 

初めから読むには、

www.crosshyou.info

です。

都道府県別の被服及び履物費のデータの分析1 - R言語にCSVファイルのデータを読み込む。

UnsplashNico Knaackが撮影した写真 

今回は、都道府県別の被服及び履物費のデータを分析してみようと思います。

データは、政府統計の総合窓口(www.e-stat.go.jp)から取得しました。

被服及び履物費を被説明変数として、人口密度や15~64歳人口割合、1人当たり県民所得を説明変数にしようと思いますので、一緒にダウンロードします。

こんな感じのCSVファイルをダウンロードしました。このファイルをR言語に取り込んで分析します。

まずは、library()関数を使ってtidyverseライブラリを読み込んでおきます。

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

glimpse()関数でデータが読み込まれているかどうかみてみます。

yearとprefが文字化けしています。何故か日本語のCSVファイルを読み込むと文字化けしてしまうんですよね。。とりあえず、1. year_codeを4桁の西暦に直す。2. yearとprefを削除する。この2つをmutate()関数とselect()関数を使って実行します。

うまくいったかglimpse()関数で見てみます。

うまくいきました。

次に、あらかじめ用意してある下図のようなCSVファイルを読み込みます。

これはcodeと英語の都道府県名があるCSVファイルです。東日本なら1のダミー変数、east, 東京都、大阪府、愛知県、千葉県、埼玉県、神奈川県なら1をとるダミー変数のbig6, 海が無い県なら1をとるダミー変数のnoseaという変数もついています。

そうしたら、dfのpref_codeとpref_codeのcodeを鍵にして結合します。inner_join()関数を使います。

glimpse()関数で結果を確認します。

うまくいきました。

つぎは、wear_shoeのデータが無い行をfilter()関数を使って削除します。

glimpse()関数を使って結果を確認します。

pref_codeはもう必要ないので削除して、year_codeという変数名をyearに変更して、変数の表示順を並び替えます。select()関数とrename()関数を使います。

summary()関数で各変数の基本統計値をみてみます。

year: 調査年は1975が最小で2007が最大です。およそ30年間のデータがありますね。

pref: 都道府県名はlength:1551と表示されていますので、このデータフレームは1551の行があることがわかります。

wear_shoe: 被服及び履物費(円)の最小値は6518円、最大値は30969円、平均値は18452円です。NAはありません。

mitsudo: 可住面積1km2当たり人口密度(人)の最小値は251.2人、最大値は9200.9人、平均値は1334.2人です。NAはありません。

wariai: 15~64歳の人口の割合(%)の最小値は58.90%、最大値は69.40%、平均値は63.65%です。NAが1410あります。

percapita17: 平成17年基準の1人当たり県民所得(千円)の最小値は203万9千円、最大値は526万6千円、平均値は282万4千円です。NAは1222個あります。

percapita23: 平成23年基準の1人当たり県民所得(千円)の最小値は200万3千円、最大値は596万6千円、平均値は283万5千円です。NAは1457個あります。

east: 東日本なら1を取るダミー変数で平均値は0.5106です。NAはありません。

big6: Tokyo, Osaka, Aichi, Chiba, Saiatama, Kanagwaなら1を取るダミー変数で平均値は0.1277です。NAはありません。

nosea: 海無し県なら1を取るダミー変数で平均値は0.1702です。NAはありません。

以上で、分析する準備が整いました。

今回は以上です。

次回は、

www.crosshyou.info

です。

J. Leagueのデータの分析 - R言語で「攻撃は最大の防御なり」か「防御は最大の攻撃なり」かを調べる。

UnsplashWesley Tingeyが撮影した写真 

今回は、J. Leagueのデータを分析してみます。勝ち点と得点、失点の関係を調べます。

まず、データをJ. Leagueの公式サイトから取得しました。

J. League Data Site (j-league.or.jp)

Webスクレイピングできればいいのですが、そこまでの技術が無いのでコピペしてExcelに貼り付けました。

こんな感じでまとめました。

これを、R言語で分析します。

まず、必要なライブラリを読み込みます。

tidyverseパッケージはRを使うときはとりあえず読み込んでおいたほうがいいですよね。modendiveパッケージは、現在、読んでいるデータ分析の本、Statistical Inference via Data Science (moderndive.com)のパッケージです。今回はこの本の内容を参考にして分析します。readxlパッケージはExcelのデータを読み込むときに使います。

では、read_xlsx()関数でデータ読み込みます。

glimpse()関数でデータが読み込まれたかどうか見てみます。

無事に読み込みできたようです。日本語のチーム名が何故か文字化けしてしまうんですよね。。今回はチーム名は分析には使わないので大丈夫です。

変数を簡単に説明します。year: 西暦、rank: 順位、name: 英語のチーム名、team: 日本語のチーム名、point: 勝点、game: 試合数、win: 勝数、lose: 敗数、tile: 引分数、goal_get: 得点、goal_give: 失点、get_give: 得失点差、です。

今回の分析では、pointが被説明変数で、goal_getとgoal_giveが説明変数です。

point, goal_get, goal_giveの基本データをみてみます。summary()関数を使います。

pointの平均値は46.65, goal_getとgoal_giveの平均値は46.16と同じです。

最小値は、pointは14, goal_getは16, goal_giveは24。

最大値は、pointは92, goal_getは88, goal_giveは88となっています。

3つの値の範囲はだいたい同じですね。

この3つの変数の相関係数を見てみます。cor()関数を使います。

pointとgoal_getは0.77の相関係数なので、比較的強い正の相関です。

pointとgoal_giveは-0.70の相関なので、比較的強い負の相関です。

goal_getとgoal_giveは-0.25の相関なので、比較的弱い負の相関です。

この結果は直感と矛盾しませんね、得点が多ければ勝点は多く、失点が多ければ勝点は少なく、得点が多ければ失点は少ない、勝点が多ければ得点は多く、勝点が少なければ失点は多く、失点が少なければ得点が多い、ということです。

plot()関数で散布図を描いてみます。

pointとgoal_get, pointとgoal_giveは相関が強く、goal_getとgoal_giveは相関が弱いことがわかります。

それでは、lm()関数で回帰分析をしてみます。

続いて、moderndiveパッケージのget_regression_table()関数で結果を表示します。

goal_getの係数は0.688でgoal_giveの係数は-0.656です。

つまり、goal_giveが変わらないとすると、1点得点すると勝点が0.688増えるということ、goal_getが変わらないとすると、1点失点すると勝点が0.656減るということです。

ということは、1得点(=勝点0.688) > 1失点(=勝点0.656)ということなので、「攻撃は最大の防御なり」ですね。

以上の結果はJ. League全部のチームの結果です。もしかしたら、上位チームと下位チームでは違うかもしれません。勝数のほうが敗数よりも多いチームは、1, そうでないときは0というダミー変数を作って、これを回帰モデルに加えてみましょう。

まずは、ダミー変数を作成します。

こうして作成した、winnerを加えた回帰モデルを作成します。

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

すこしややこしいですが winner = 1 のときとwinner = 0 のときでわけて考えます。

winner = 1, つまり勝ち越しチームは

point = 40.1 + 9.23 + (0.594 + 0.022) * goal_get - (0.507 + 0.13) * goal_give

= 49.33 + 0.616 * goal_get - 0.520 * goal_give となります。

winner = 0, つまり負け越しチームは、

point = 40.1 + 0.594 * goal_get - 0.507 + goal_give となります。

どちらの場合でも、goal_getの係数の絶対値ほうがgoal_giveの係数の絶対値よりも大きいですから、1得点の価値のほうが1失点の価値よりも大きいことがわかりました。

今度は、J. Leagueの開催された年を前半、後半に区切ったダミー変数を作成して、分析してみましょう。

まずは、前半の年なら1のダミー変数を作成します。

回帰モデルを作成します。

結果をget_regression_table()関数で表示します。

old_year = 1 のときと、old_year = 0 のときでわけて推計式を考えます。

old_year = 1 のとき、

point = 46.6 - 4.03 + (0.707 - 0.029) * goal_get - (0.687 - 0.079) * goal_give

= 43.57 + 0.736 * goal_get - 0.608 * goal_give

old_year = 0 のとき、

point = 46.6 + 0.707 * goal_get - 0.687 * goal_give

です。どちらのケースでもgoal_getの係数の絶対値のほうがgola_giveの係数の絶対値よりも大きいですね。1得点の価値のほうが1失点の価値よりも大きいです。

以上の分析結果から、J. Leagueでは「攻撃は最大の防御なり」ということが言えるようです。

 

読書記録 - 「人類と病 - 国際政治から見る感染症と健康格差」 詫磨佳代 著 中公新書

人類が感染症や生活習慣病とどのように対処してきたのかをまとめた本です。

中世のペスト(黒死病)から現在の新型コロナウイルス、糖尿病などの生活習慣病などに人類がどのように向き合ってきたのか、この本を読んで感じたことは、病気と闘うには医学の力だけでは足りなくて、政治や産業界の力も必要だということです。

 

OECD Researchers data analysis 6 - Multiple Linear Regression in R

UnsplashEvi T.が撮影した写真 

www.crosshyou.info

This post is following of the above post.
In the previous post, I did sinple linear regression, it menas there is only one explanatory vatiable. In this post I will do multiple linear regression, it means there are multiple explanatory variables.

Firstly, I make time variable which is factor type variable.

Let's see table of time.

Then, I add time for explanatory vatiable.

Let's see regrettion table with moderndive package's get_regression_table() function.

p_value of time: 2017 or l_tot_1000employed:time is greater than 0.05, so I know those variables are not statistically significant.

Let's visualize this regression.

Next, I use l_women_pc_researcher instead of time.

Let's see the regression table.

I see l_tot_1000employed, l_women_pc_researcher and l_tot_1000employed:l_women_pc_researcher have small p-value less than 0.05.
So, those three are statistically significant.

Let's visialize actual l_usd_cap and estimated l_usd_cap.

Firstly, I show regression data point with moderndive's get_regression_package

I used ggplot() + geom_point() function to draw a scatter plot.

If actual l_usd_cap and estimated l_usd_cap are exactly same, all scatter pooints shold be on the line.

Let's see the histogram of the residuals.

Now, I have three linear regression models. Let's compare them with moderndive's get_regression_summaries() function.

model3, which has two explanatory variables, l_tot_1000employed and l_women_pc_researcher, has the smallest emse and the largest adj_r_squared.

That's it. Thank you!

To read from the fist post,

www.crosshyou.info